[Bug math/28472] New: pow(10, i) accuracy

[Bug math/28472] New: pow(10, i) accuracy

[mwelinder at gmail dot com]

https://sourceware.org/bugzilla/show_bug.cgi?id=28472

            Bug ID: 28472
           Summary: pow(10, i) accuracy
           Product: glibc
           Version: 2.31
            Status: UNCONFIRMED
          Severity: normal
          Priority: P2
         Component: math
          Assignee: unassigned at sourceware dot org
          Reporter: mwelinder at gmail dot com
  Target Milestone: ---

Created attachment 13725
  --> https://sourceware.org/bugzilla/attachment.cgi?id=13725&action=edit
Test program

When the attached test program is run, I see two cases where pow(10,i)
disagrees with atof(1e$i") for an integer i.

This is a last-bit rounding issue.  The values returned by pow(10,i) are very
close and normally I would cut pow a bit of slack -- it is a lot trickier than
the single-argument functions like sin, cos, or exp.

However, pow(10,i) is a fairly important case that it would make sense to
handle accurately.


$ ./a.out 
i = 23
a = 100000000000000008388608
b = 99999999999999991611392
Midpoint

i = 210
a =
1000000000000000073111882183254852571116159535704205070042237624441112422237792851875363410143857412667610687999697631253349027916052430446705469082528474390439305760542775847335624615778546587814778848485048320
b =
999999999999999927113782419344605574598668153294882673458925392487194643703632279098558059466181044478400725843812838336795121561031396504666917998514458446354143529431921823271795036250068185162804696593727488
b is best after 17 digits

--------------------------------------------------

Note: On another machine I have seen a failure for i=216, but only in 32-bit
mode.  There might be hardware dependencies here.

$ grep Intel /proc/cpuinfo  | head -2
vendor_id       : GenuineIntel
model name      : Intel(R) Core(TM) i7-8550U CPU @ 1.80GHz

$ grep PRETTY /etc/os-release 
PRETTY_NAME="Linux Mint 20.2"

I believe I am using libc version 2.31-0ubuntu9.2cross1


Compilation:

$ gcc -v -g -Wall -O2 eee.c -lm
Using built-in specs.
COLLECT_GCC=gcc
COLLECT_LTO_WRAPPER=/usr/lib/gcc/x86_64-linux-gnu/9/lto-wrapper
OFFLOAD_TARGET_NAMES=nvptx-none:hsa
OFFLOAD_TARGET_DEFAULT=1
Target: x86_64-linux-gnu
Configured with: ../src/configure -v --with-pkgversion='Ubuntu
9.3.0-17ubuntu1~20.04' --with-bugurl=file:///usr/share/doc/gcc-9/README.Bugs
--enable-languages=c,ada,c++,go,brig,d,fortran,objc,obj-c++,gm2 --prefix=/usr
--with-gcc-major-version-only --program-suffix=-9
--program-prefix=x86_64-linux-gnu- --enable-shared --enable-linker-build-id
--libexecdir=/usr/lib --without-included-gettext --enable-threads=posix
--libdir=/usr/lib --enable-nls --enable-clocale=gnu --enable-libstdcxx-debug
--enable-libstdcxx-time=yes --with-default-libstdcxx-abi=new
--enable-gnu-unique-object --disable-vtable-verify --enable-plugin
--enable-default-pie --with-system-zlib --with-target-system-zlib=auto
--enable-objc-gc=auto --enable-multiarch --disable-werror --with-arch-32=i686
--with-abi=m64 --with-multilib-list=m32,m64,mx32 --enable-multilib
--with-tune=generic
--enable-offload-targets=nvptx-none=/build/gcc-9-HskZEa/gcc-9-9.3.0/debian/tmp-nvptx/usr,hsa
--without-cuda-driver --enable-checking=release --build=x86_64-linux-gnu
--host=x86_64-linux-gnu --target=x86_64-linux-gnu
Thread model: posix
gcc version 9.3.0 (Ubuntu 9.3.0-17ubuntu1~20.04) 
COLLECT_GCC_OPTIONS='-v' '-g' '-Wall' '-O2' '-mtune=generic' '-march=x86-64'
 /usr/lib/gcc/x86_64-linux-gnu/9/cc1 -quiet -v -imultiarch x86_64-linux-gnu
eee.c -quiet -dumpbase eee.c -mtune=generic -march=x86-64 -auxbase eee -g -O2
-Wall -version -fasynchronous-unwind-tables -fstack-protector-strong
-Wformat-security -fstack-clash-protection -fcf-protection -o /tmp/ccbP42NQ.s
GNU C17 (Ubuntu 9.3.0-17ubuntu1~20.04) version 9.3.0 (x86_64-linux-gnu)
        compiled by GNU C version 9.3.0, GMP version 6.2.0, MPFR version 4.0.2,
MPC version 1.1.0, isl version isl-0.22.1-GMP

GGC heuristics: --param ggc-min-expand=100 --param ggc-min-heapsize=131072
ignoring nonexistent directory "/usr/local/include/x86_64-linux-gnu"
ignoring nonexistent directory "/usr/lib/gcc/x86_64-linux-gnu/9/include-fixed"
ignoring nonexistent directory
"/usr/lib/gcc/x86_64-linux-gnu/9/../../../../x86_64-linux-gnu/include"
#include "..." search starts here:
#include <...> search starts here:
 /usr/lib/gcc/x86_64-linux-gnu/9/include
 /usr/local/include
 /usr/include/x86_64-linux-gnu
 /usr/include
End of search list.
GNU C17 (Ubuntu 9.3.0-17ubuntu1~20.04) version 9.3.0 (x86_64-linux-gnu)
        compiled by GNU C version 9.3.0, GMP version 6.2.0, MPFR version 4.0.2,
MPC version 1.1.0, isl version isl-0.22.1-GMP

GGC heuristics: --param ggc-min-expand=100 --param ggc-min-heapsize=131072
Compiler executable checksum: bbf13931d8de1abe14040c9909cb6969
COLLECT_GCC_OPTIONS='-v' '-g' '-Wall' '-O2' '-mtune=generic' '-march=x86-64'
 as -v --64 -o /tmp/ccSWUKiU.o /tmp/ccbP42NQ.s
GNU assembler version 2.34 (x86_64-linux-gnu) using BFD version (GNU Binutils
for Ubuntu) 2.34
COMPILER_PATH=/usr/lib/gcc/x86_64-linux-gnu/9/:/usr/lib/gcc/x86_64-linux-gnu/9/:/usr/lib/gcc/x86_64-linux-gnu/:/usr/lib/gcc/x86_64-linux-gnu/9/:/usr/lib/gcc/x86_64-linux-gnu/
LIBRARY_PATH=/usr/lib/gcc/x86_64-linux-gnu/9/:/usr/lib/gcc/x86_64-linux-gnu/9/../../../x86_64-linux-gnu/:/usr/lib/gcc/x86_64-linux-gnu/9/../../../../lib/:/lib/x86_64-linux-gnu/:/lib/../lib/:/usr/lib/x86_64-linux-gnu/:/usr/lib/../lib/:/usr/lib/gcc/x86_64-linux-gnu/9/../../../:/lib/:/usr/lib/
COLLECT_GCC_OPTIONS='-v' '-g' '-Wall' '-O2' '-mtune=generic' '-march=x86-64'
 /usr/lib/gcc/x86_64-linux-gnu/9/collect2 -plugin
/usr/lib/gcc/x86_64-linux-gnu/9/liblto_plugin.so
-plugin-opt=/usr/lib/gcc/x86_64-linux-gnu/9/lto-wrapper
-plugin-opt=-fresolution=/tmp/ccMHqq8T.res -plugin-opt=-pass-through=-lgcc
-plugin-opt=-pass-through=-lgcc_s -plugin-opt=-pass-through=-lc
-plugin-opt=-pass-through=-lgcc -plugin-opt=-pass-through=-lgcc_s --build-id
--eh-frame-hdr -m elf_x86_64 --hash-style=gnu --as-needed -dynamic-linker
/lib64/ld-linux-x86-64.so.2 -pie -z now -z relro
/usr/lib/gcc/x86_64-linux-gnu/9/../../../x86_64-linux-gnu/Scrt1.o
/usr/lib/gcc/x86_64-linux-gnu/9/../../../x86_64-linux-gnu/crti.o
/usr/lib/gcc/x86_64-linux-gnu/9/crtbeginS.o -L/usr/lib/gcc/x86_64-linux-gnu/9
-L/usr/lib/gcc/x86_64-linux-gnu/9/../../../x86_64-linux-gnu
-L/usr/lib/gcc/x86_64-linux-gnu/9/../../../../lib -L/lib/x86_64-linux-gnu
-L/lib/../lib -L/usr/lib/x86_64-linux-gnu -L/usr/lib/../lib
-L/usr/lib/gcc/x86_64-linux-gnu/9/../../.. /tmp/ccSWUKiU.o -lm -lgcc
--push-state --as-needed -lgcc_s --pop-state -lc -lgcc --push-state --as-needed
-lgcc_s --pop-state /usr/lib/gcc/x86_64-linux-gnu/9/crtendS.o
/usr/lib/gcc/x86_64-linux-gnu/9/../../../x86_64-linux-gnu/crtn.o
COLLECT_GCC_OPTIONS='-v' '-g' '-Wall' '-O2' '-mtune=generic' '-march=x86-64'

-- 
You are receiving this mail because:
You are on the CC list for the bug.

[Bug math/28472] pow(10, i) accuracy

[joseph at codesourcery dot com]

https://sourceware.org/bugzilla/show_bug.cgi?id=28472

--- Comment #1 from joseph at codesourcery dot com <joseph at codesourcery dot com> ---
Last-bit rounding issues aren't bugs except for functions such as sqrt or 
fma that are fully defined by bindings to IEEE 754 operations.

I'd encourage using exp10 (in glibc for a very long time, standardized in 
C23) instead of pow with first argument 10; that probably won't improve 
last-bit accuracy, but is likely to be faster.

-- 
You are receiving this mail because:
You are on the CC list for the bug.

[Bug math/28472] pow(10, i) accuracy

[mwelinder at gmail dot com]

https://sourceware.org/bugzilla/show_bug.cgi?id=28472

--- Comment #2 from M Welinder <mwelinder at gmail dot com> ---
I am assuming that "Last-bit rounding issues aren't bugs" really means
something like "the standard allow it".  They are still bugs in the "it would
be better to return the right result" sense, even if severity is much lowered
and considerations of speed and man-power enter into the picture.

For what it is worth, IEEE-754 since at least 2008 has recommended
correctly-rounded exp10 and pown, both of which are simpler than the generic
pow.

-- 
You are receiving this mail because:
You are on the CC list for the bug.

[Bug math/28472] pow(10, i) accuracy

[joseph at codesourcery dot com]

https://sourceware.org/bugzilla/show_bug.cgi?id=28472

--- Comment #3 from joseph at codesourcery dot com <joseph at codesourcery dot com> ---
"aren't bugs" means "the accuracy goals for libm functions are documented 
in the glibc manual, and a small error for most functions is within those 
goals" (and even if we strengthened the goals, they'd still allow last-bit 
errors).

The IEEE 754 recommended correctly-rounded operations correspond to 
functions such as cr_exp10 (names reserved in C23), which are probably not 
in scope for glibc considering the difficulty in providing such functions 
across the various floating-point formats supported (determining worst 
cases for correct rounding still requires unfeasible resources for 
binary128) and the maintainability and usability problems with having a 
supported API with arbitrary variations in what functions are supported 
for what types.

-- 
You are receiving this mail because:
You are on the CC list for the bug.

[Bug math/28472] pow(10, i) accuracy

[newbie-02 at gmx dot de]

https://sourceware.org/bugzilla/show_bug.cgi?id=28472

b. <newbie-02 at gmx dot de> changed:

           What    |Removed                     |Added
----------------------------------------------------------------------------
                 CC|                            |newbie-02 at gmx dot de

--- Comment #4 from b. <newbie-02 at gmx dot de> ---
confirm the error,  
support the desire for improvement,  
powers of ten are important for rounding other values,  
problematic enough that they are exact only from 10^0 ... 10^22 (in doubles),  
they do not need additional deviations,  
IMHO IEEE 754 specifies:  
 - `as calculated with infinite precision and then rounded to representable
precision', 
   and for 'midpoints': 
 - 'to the even' - last bit 0, 
at least with -O0 gcc should fulfill such ...???
would set to new but lack the power or the role.

-- 
You are receiving this mail because:
You are on the CC list for the bug.

[Bug math/28472] pow(10, i) accuracy

[newbie-02 at gmx dot de]

https://sourceware.org/bugzilla/show_bug.cgi?id=28472

--- Comment #5 from b. <newbie-02 at gmx dot de> ---
if ... I'm right in testing ... Intel 'oneAPI' compiler(s) 'icpx' and 'dpcpp'
produce clean powers of 10 ( matching E-string evaluation ) for long double
figures with 'exp10l( x )'.  

since they do not compile themselves but use clang/clang++ for this purpose -
with a huge bunch of options and includes - it should be possible to have a
look at it and gather some inspiration ... ??? hot spots seem to be calling
clang with '-cc1' and '-x c++'.  

( for the bug itself i propose testing and setting to new, besides the
deviation it is a consistency problem which is difficult to explain to users:
that one get's different results in one system for different formulations of
the same  mathematical task. ('pow( 10, x )' != '1Ex' != 'pow10( x )' !=
'exp10( x )' ).  

( and of course gcc shouldn't be second in 'best compilers' ) 

reg. the proposal of using exp10( x ) instead of pow( 10, x ):  
'exp10( x )' compiled with gcc -O2: 230 fails between x = -324 and x = 308, 
'pow( 10, x )': only 2, 
as well with 'exp10( x )' 
( complaining " incompatible implicit declaration of built-in function ‘exp10’
" ) 
as for '__builtin_exp10( x )'.  
as well for gcc, g++, clang++, 
'clang' additional complaints and producing even more fails.

-- 
You are receiving this mail because:
You are on the CC list for the bug.

[Bug math/28472] pow(10, i) accuracy

[mwelinder at gmail dot com]

https://sourceware.org/bugzilla/show_bug.cgi?id=28472

--- Comment #6 from M Welinder <mwelinder at gmail dot com> ---
> I'd encourage using exp10

At present that seems to be a very bad idea.  The accuracy is worse than
pow(10,i).

For example, here are the integer powers of 10 from 0-21.  Problems with 6/22
results.

                                 1.00000000000000000000000000000000000
                                10.00000000000000000000000000000000000
                               100.00000000000000000000000000000000000
                              1000.00000000000011368683772161602973938
                             10000.00000000000181898940354585647583008
                            100000.00000000000000000000000000000000000
                           1000000.00000000000000000000000000000000000
                           9999999.99999999813735485076904296875000000
                          99999999.99999998509883880615234375000000000
                         999999999.99999988079071044921875000000000000
                       10000000000.00000000000000000000000000000000000
                      100000000000.00000000000000000000000000000000000
                     1000000000000.00000000000000000000000000000000000
                    10000000000000.00000000000000000000000000000000000
                   100000000000000.00000000000000000000000000000000000
                  1000000000000000.00000000000000000000000000000000000
                 10000000000000000.00000000000000000000000000000000000
                 99999999999999984.00000000000000000000000000000000000
               1000000000000000000.00000000000000000000000000000000000
              10000000000000000000.00000000000000000000000000000000000
             100000000000000000000.00000000000000000000000000000000000
            1000000000000000000000.00000000000000000000000000000000000

-- 
You are receiving this mail because:
You are on the CC list for the bug.

[Bug math/28472] pow(10, i) accuracy

[joseph at codesourcery dot com]

https://sourceware.org/bugzilla/show_bug.cgi?id=28472

--- Comment #7 from joseph at codesourcery dot com <joseph at codesourcery dot com> ---
Being correctly rounded is out of scope for both pow and exp10 and 
failures to be correctly rounded should not be considered problems, only 
errors of more than a few ulps are bugs in glibc.

We're a lot more likely to get cr_exp10 than cr_pow, but that would depend 
on both CORE-MATH implementing it and someone doing the work to integrate 
it into glibc (note that the present CORE-MATH implementations aren't 
suitably portable for glibc at present - they depend on 64-bit systems, 
for example, via hardcoded assumptions on the size of long and 
requirements in places for 128-bit integer types - are very lacking in 
comments, and may not meet the requirements for correct exceptions from 
cr_*).

-- 
You are receiving this mail because:
You are on the CC list for the bug.

[Bug math/28472] pow(10, i) accuracy

[newbie-02 at gmx dot de]

https://sourceware.org/bugzilla/show_bug.cgi?id=28472

--- Comment #8 from b. <newbie-02 at gmx dot de> ---

are the following differences based on any reasonable idea? is it of help to
dig down this issue?  

printf( "%.25E\n", pow( 10, 23 ) );  -> 9.9999999999999991611392000E+22  

in contrast to:  
int j = 23;  
printf( "%.25E\n", pow( 10, j ) );   -> 1.0000000000000000838860800E+23  

IMHO the results should match!  

same with 210:  

printf( "%.25E\n", pow( 10, 210 ) ); -> 9.9999999999999992711378242E+209  

while:  

j = 210;  
printf( "%.25E\n", pow( 10, j ) );   -> 1.0000000000000000731118822E+210  

seems independent of j defined as integer or double, and constants written as
'x' or 'x.0', but dependent on compiler options:  

no '-Ox' option set or '-O0': above differences,  
'-O1', '-O2' or '-O3' set   : results match, all 9.99999xxx,  

different evaluation at compile- vs. runtime? I seem to remember similar
oddities when a compiler replaces long double constants with doubles at compile
time????  

but why here? me in error? me stupid, me blackout? or is gcc / glibc / ???
kidding with the programmers?  

admit to be confused, esp. about this error being dependent on compiler
options,  while the test program by M Welinder, and gnumeric where the problem
first came up, are not.

-- 
You are receiving this mail because:
You are on the CC list for the bug.

[Bug math/28472] pow(10, i) accuracy

[newbie-02 at gmx dot de]

https://sourceware.org/bugzilla/show_bug.cgi?id=28472

--- Comment #9 from b. <newbie-02 at gmx dot de> ---
gnumeric / goffice have found a solution / workaround to get correct powers of
ten for integer exponents, see
https://gitlab.gnome.org/GNOME/goffice/-/issues/62 .  
It works quite well and performant, for long doubles the table needs some space
( ~260kB ), the effort is justified because the results influence the scaling,
and thus results, in almost all other rounding calculations.  
In this respect, I would recommend implementing similar for glibc until a
better solution is found.  
( having a bug with clear description and test program 'unconfirmed' after more
than 12 month is a questionable performance ... )

-- 
You are receiving this mail because:
You are on the CC list for the bug.

[Bug math/28472] pow(10, i) accuracy

[joseph at codesourcery dot com]

https://sourceware.org/bugzilla/show_bug.cgi?id=28472

--- Comment #10 from joseph at codesourcery dot com <joseph at codesourcery dot com> ---
It's already explained that this is *not a bug*, given the documented 
glibc accuracy goals.  It could be closed as INVALID, but not confirmed as 
a bug.

-- 
You are receiving this mail because:
You are on the CC list for the bug.

[Bug math/28472] pow(10, i) accuracy

[newbie-02 at gmx dot de]

https://sourceware.org/bugzilla/show_bug.cgi?id=28472

--- Comment #11 from b. <newbie-02 at gmx dot de> ---
hello @joseph@codesourcery.com, thanks for the quick answer,  

I consider it problematic when such categorizations block possible progress.  
I see the problem important, and sensibly and reasonably substantiated.  
Is there any option to implement it? enhancement request? patch proposal?

-- 
You are receiving this mail because:
You are on the CC list for the bug.

[Bug math/28472] pow(10, i) accuracy

[joseph at codesourcery dot com]

https://sourceware.org/bugzilla/show_bug.cgi?id=28472

--- Comment #12 from joseph at codesourcery dot com <joseph at codesourcery dot com> ---
I don't consider special-casing 10 in pow to be reasonable.  Integrating 
better function implementations from CORE-MATH would be reasonable, but 
it's important there to take account of portability considerations not 
handled in the CORE-MATH code (glibc supports 32-bit platforms, platforms 
without all the exception macros defined, platforms where there may be 
excess precision for intermediate computations, platforms where various 
__builtin_* functions used in CORE-MATH will be expanded out-of-line and 
slow, etc., and has various ABI and namespace considerations not covered 
in the CORE-MATH code as well, so significant work would be needed to turn 
such an implementation into something suitable for glibc).  Note that for 
any such new function implementatations, having suitable benchmark inputs 
in glibc's benchmarks is important first, to demonstrate that a new 
implementation is at least as fast as well as more accurate.

-- 
You are receiving this mail because:
You are on the CC list for the bug.

[Bug math/28472] pow(10, i) accuracy

[newbie-02 at gmx dot de]

https://sourceware.org/bugzilla/show_bug.cgi?id=28472

--- Comment #13 from b. <newbie-02 at gmx dot de> ---
(In reply to joseph@codesourcery.com from comment #12)
> I don't consider special-casing 10 in pow to be reasonable.  
As M Welinder pointed out in C#2 it's an IEEE recommendation to have these
values correct, and - again - it's important in rounding other values. Rounding
is the only! weapon a user has against FP-math imprecision, this should be
taken important.  
And it fundamentally harms math logic if pow( 10, 23 ) <> 1E+23.  
Can't argue about the other points because I don't know enough about CORE,
_builtin ..., IMHO the workaround is independent of such ...  

The request is to have correct powers of ten.  

How this is achieved is a matter for the experts. It was important to me to
emphasize that it is! important and to show that it is! possible.

-- 
You are receiving this mail because:
You are on the CC list for the bug.

[Bug math/28472] pow(10, i) accuracy

[vincent-srcware at vinc17 dot net]

https://sourceware.org/bugzilla/show_bug.cgi?id=28472

Vincent Lefèvre <vincent-srcware at vinc17 dot net> changed:

           What    |Removed                     |Added
----------------------------------------------------------------------------
                 CC|                            |vincent-srcware at vinc17 dot net

--- Comment #14 from Vincent Lefèvre <vincent-srcware at vinc17 dot net> ---
(In reply to b. from comment #13)
> (In reply to joseph@codesourcery.com from comment #12)
> > I don't consider special-casing 10 in pow to be reasonable.  
> As M Welinder pointed out in C#2 it's an IEEE recommendation to have these
> values correct, and - again - it's important in rounding other values.

What IEEE 754 recommends is that some math functions, like pow and exp10, be
correctly rounded. If a function is not correctly rounded, there are no
additional recommendations at all. In particular, IEEE 754 does not recommend
anything for some particular values, like powers of 10. Having special-casing
may be counter-productive, because this can make the function less smooth at
these points and/or break monotonicity (though monotonicity is not guaranteed).
What could be done is to make such functions more accurate, say ensure faithful
rounding for round-to-nearest (so with an error strictly less than 1 ulp of the
exact result and of the computed result), even though this would make them a
bit slower. Users who do not care about the accuracy could still write faster
(and less accurate) functions than those already provided, anyway.

-- 
You are receiving this mail because:
You are on the CC list for the bug.

[Bug math/28472] pow(10, i) accuracy

[wdijkstr at arm dot com]

https://sourceware.org/bugzilla/show_bug.cgi?id=28472

Wilco <wdijkstr at arm dot com> changed:

           What    |Removed                     |Added
----------------------------------------------------------------------------
                 CC|                            |wdijkstr at arm dot com

--- Comment #15 from Wilco <wdijkstr at arm dot com> ---
GLIBC double precision pow is the most accurate of all libraries tested at
0.523 ULP [1].

The new exp10 is also the most accurate of the 13 tested math libraries.

If you complain about inaccuracies in the most accurate library then maybe your
expectations are a little bit off...

The fact is, binary floating point cannot represent all powers of 10. If you
don't like the rounding behaviour of floating point, don't use floating point.

Note compilers optimize pow (C, x) into exp (x * log (C)) with -Ofast. However
if C is a multiple of 2 or 10, we could use exp2 or exp10 for slightly better
accuracy.

[1] https://members.loria.fr/PZimmermann/papers/accuracy.pdf

-- 
You are receiving this mail because:
You are on the CC list for the bug.

[Bug math/28472] pow(10, i) accuracy

[vincent-srcware at vinc17 dot net]

https://sourceware.org/bugzilla/show_bug.cgi?id=28472

--- Comment #16 from Vincent Lefèvre <vincent-srcware at vinc17 dot net> ---
(In reply to Wilco from comment #15)
> GLIBC double precision pow is the most accurate of all libraries tested at
> 0.523 ULP [1].

What you forget is that this is the accuracy *tested* on arbitrary values. The
actual accuracy may be worse. And this is the case here, with an accuracy
larger than 1 ulp, according to the results in Comment #6!

> The new exp10 is also the most accurate of the 13 tested math libraries.
> 
> If you complain about inaccuracies in the most accurate library then maybe
> your expectations are a little bit off...

In the present case, it may be far worse than the most accurate libraries
(well, it is difficult to say, due to the random tests). The result returned by
glibc is not even faithfully rounded. So the user is right to complain, even
though there is no guarantee from the ISO C standard.

> The fact is, binary floating point cannot represent all powers of 10. If you
> don't like the rounding behaviour of floating point, don't use floating
> point.

Don't blame floating point if this is a poor implementation.

> Note compilers optimize pow (C, x) into exp (x * log (C)) with -Ofast.
> However if C is a multiple of 2 or 10, we could use exp2 or exp10 for
> slightly better accuracy.

I suppose that you mean "if C is a *power* of 2 or 10".

-- 
You are receiving this mail because:
You are on the CC list for the bug.

[Bug math/28472] pow(10, i) accuracy

[wdijkstr at arm dot com]

https://sourceware.org/bugzilla/show_bug.cgi?id=28472

--- Comment #17 from Wilco <wdijkstr at arm dot com> ---
(In reply to Vincent Lefèvre from comment #16)
> (In reply to Wilco from comment #15)
> > GLIBC double precision pow is the most accurate of all libraries tested at
> > 0.523 ULP [1].
> 
> What you forget is that this is the accuracy *tested* on arbitrary values.
> The actual accuracy may be worse. And this is the case here, with an
> accuracy larger than 1 ulp, according to the results in Comment #6!

Please see the implementation - it documents the accuracy across the full input
ranges. The worst-case reported by random testing is slightly lower due to not
being able to test all input values.

And comment #6 discusses exp10, which had a known ULP of 2.01 in previous
GLIBCs.

> > The new exp10 is also the most accurate of the 13 tested math libraries.
> > 
> > If you complain about inaccuracies in the most accurate library then maybe
> > your expectations are a little bit off...
> 
> In the present case, it may be far worse than the most accurate libraries
> (well, it is difficult to say, due to the random tests). The result returned
> by glibc is not even faithfully rounded. So the user is right to complain,
> even though there is no guarantee from the ISO C standard.

No, it's not difficult to say. We computed the accuracy and have *documented*
it in the source code. So it's not only not a "guess", it's actually impossible
to get cases that are worse. Ie. if we have an algorithm that does < 0.55ULP
before rounding, we can't ever get a 2 ULP error.

> > The fact is, binary floating point cannot represent all powers of 10. If you
> > don't like the rounding behaviour of floating point, don't use floating
> > point.
> 
> Don't blame floating point if this is a poor implementation.

Even the old exp10 wasn't disastrously bad like j0/j1/y0/y0/tgamma.

> > Note compilers optimize pow (C, x) into exp (x * log (C)) with -Ofast.
> > However if C is a multiple of 2 or 10, we could use exp2 or exp10 for
> > slightly better accuracy.
> 
> I suppose that you mean "if C is a *power* of 2 or 10".

Correct.

-- 
You are receiving this mail because:
You are on the CC list for the bug.

[Bug math/28472] pow(10, i) accuracy

[vincent-srcware at vinc17 dot net]

https://sourceware.org/bugzilla/show_bug.cgi?id=28472

--- Comment #18 from Vincent Lefèvre <vincent-srcware at vinc17 dot net> ---
(In reply to Wilco from comment #17)
> (In reply to Vincent Lefèvre from comment #16)
> > (In reply to Wilco from comment #15)
> > > GLIBC double precision pow is the most accurate of all libraries tested at
> > > 0.523 ULP [1].
> > 
> > What you forget is that this is the accuracy *tested* on arbitrary values.
> > The actual accuracy may be worse. And this is the case here, with an
> > accuracy larger than 1 ulp, according to the results in Comment #6!
> 
> Please see the implementation - it documents the accuracy across the full
> input ranges. The worst-case reported by random testing is slightly lower
> due to not being able to test all input values.

OK, this is something that is new to me. So perhaps [1] should also document
the proved error bounds.

> And comment #6 discusses exp10, which had a known ULP of 2.01 in previous
> GLIBCs.

Sorry, I got confused by a recent mail (2024-03-02) giving these values and
saying that this was the corresponding bug report, but this bug is really about
pow() according to its title.

> Ie. if we have an algorithm that does < 0.55ULP before rounding, we can't ever > get a 2 ULP error.

OK, but in addition to the source, such proved error bounds (assuming no
mistakes in the analysis) should be given in the glibc manual.

Also, note that an error < 0.55 ULP before rounding does not guarantee that
pow(2,i) will be correct for the integers i. So, perhaps more should be said.

-- 
You are receiving this mail because:
You are on the CC list for the bug.

[Bug math/28472] pow(10, i) accuracy

[wdijkstr at arm dot com]

https://sourceware.org/bugzilla/show_bug.cgi?id=28472

--- Comment #19 from Wilco <wdijkstr at arm dot com> ---
(In reply to Vincent Lefèvre from comment #18)
> (In reply to Wilco from comment #17)
> > (In reply to Vincent Lefèvre from comment #16)
> > > (In reply to Wilco from comment #15)
> > > > GLIBC double precision pow is the most accurate of all libraries tested at
> > > > 0.523 ULP [1].
> > > 
> > > What you forget is that this is the accuracy *tested* on arbitrary values.
> > > The actual accuracy may be worse. And this is the case here, with an
> > > accuracy larger than 1 ulp, according to the results in Comment #6!
> > 
> > Please see the implementation - it documents the accuracy across the full
> > input ranges. The worst-case reported by random testing is slightly lower
> > due to not being able to test all input values.
> 
> OK, this is something that is new to me. So perhaps [1] should also document
> the proved error bounds.

That might be too much work to figure out for every math library (AFAIK only
GLIBC documents these error bounds in the implementations). Note the final
result vary depending on FMA support, evaluation order etc.

> > Ie. if we have an algorithm that does < 0.55ULP before rounding, we can't ever > get a 2 ULP error.
> 
> OK, but in addition to the source, such proved error bounds (assuming no
> mistakes in the analysis) should be given in the glibc manual.
> 
> Also, note that an error < 0.55 ULP before rounding does not guarantee that
> pow(2,i) will be correct for the integers i. So, perhaps more should be said.

As it happens, it does work out perfectly for pow(2, i) and exp2(i). Only exp10
and pow(10, i) have a few cases that are not:

exp10(126): 0x1.7a2ecc414a04p+418 0x1.7a2ecc414a03f7ff6ca1cb527788p+418

pow(10, 23): 0x1.52d02c7e14af7p+76 0x1.52d02c7e14af68p+76 
pow(10, 210): 0x1.8557f31326bbcp+697 0x1.8557f31326bbb7fcd59ff127a01dp+697 

These are very close to the halfway rounding point and would need an extra 16
bits internal precision to get right. pow(10, 23) is exactly on the mid-point,
so is even harder.

Note GLIBC previously had some correctly rounded math functions, but they were
a total nightmare. Whenever you got close to a halfway point, performance
dropped by a factor of a million times, even for cases that were relatively
easy.

That's a lesson learnt that we don't ever want to repeat again.

-- 
You are receiving this mail because:
You are on the CC list for the bug.

[Bug math/28472] pow(10, i) accuracy

[vincent-srcware at vinc17 dot net]

https://sourceware.org/bugzilla/show_bug.cgi?id=28472

--- Comment #20 from Vincent Lefèvre <vincent-srcware at vinc17 dot net> ---
(In reply to Wilco from comment #19)
> (In reply to Vincent Lefèvre from comment #18)
> > OK, this is something that is new to me. So perhaps [1] should also document
> > the proved error bounds.
> 
> That might be too much work to figure out for every math library (AFAIK only
> GLIBC documents these error bounds in the implementations). Note the final
> result vary depending on FMA support, evaluation order etc.

Well, I meant only when the bound is known.

> As it happens, it does work out perfectly for pow(2, i) and exp2(i). Only
> exp10 and pow(10, i) have a few cases that are not:
> 
> exp10(126): 0x1.7a2ecc414a04p+418 0x1.7a2ecc414a03f7ff6ca1cb527788p+418

I meant only for the exact cases (pow(2,i) is always exact for integer i,
except in case of overflow/underflow, of course).

exp10(i) is exact for integers i up to 22 (5^22 < 2^53 < 5^23), but these are
not powers of 2. You need to be careful with powers of 2, because this is where
the value of the ulp changes (I've seen errors in several proofs because of
that); such cases may need a separate check.

> Note GLIBC previously had some correctly rounded math functions, but they
> were a total nightmare. Whenever you got close to a halfway point,
> performance dropped by a factor of a million times, even for cases that were
> relatively easy.

The implementation was much too pessimistic about the probability of failure
for the slow path, which used a precision of 768 bits (much more than what is
assumed for cryptographic hashes, for instance). Moreover, some slow cases (for
pow) were inputs that could easily be found in practice. CORE-MATH should be a
huge improvement.

-- 
You are receiving this mail because:
You are on the CC list for the bug.

[Bug math/28472] pow(10, i) accuracy

[newbie-02 at gmx dot de]

https://sourceware.org/bugzilla/show_bug.cgi?id=28472

--- Comment #21 from b. <newbie-02 at gmx dot de> ---
dear guys,  

in a way you are talking this point to death,  

this issue is not about powers of two, usually they have   
few problems in binary datatypes, but about powers of 10!  

it's important to have them exact ( as good as possible )  
at integral powers of ten,  

and regarding monotonity it's surely better to smoothe the  
nearby values to the exact values than the other way.  

Note, we need 2! results changed for binary64's ( bin80's  
are worse ) and there are no monotonity issues around them.

-- 
You are receiving this mail because:
You are on the CC list for the bug.

[Bug math/28472] pow(10, i) accuracy

[wdijkstr at arm dot com]

https://sourceware.org/bugzilla/show_bug.cgi?id=28472

--- Comment #22 from Wilco <wdijkstr at arm dot com> ---
(In reply to b. from comment #21)
> dear guys,  
>   
> in a way you are talking this point to death,  
>   
> this issue is not about powers of two, usually they have   
> few problems in binary datatypes, but about powers of 10!  
>   
> it's important to have them exact ( as good as possible )  
> at integral powers of ten,  

We do get all of the powers of 10 that are exactly representable correct. And
almost all of the inexact cases are rounded correctly.

> and regarding monotonity it's surely better to smoothe the  
> nearby values to the exact values than the other way.  
>   
> Note, we need 2! results changed for binary64's ( bin80's  
> are worse ) and there are no monotonity issues around them.

There isn't an obvious fix for the current implementation. If you can find a
way without a performance hit, I'd love to hear it. If you want correctly
rounded results (and don't care about performance), just use a correctly
rounded math library.

-- 
You are receiving this mail because:
You are on the CC list for the bug.

[Bug math/28472] pow(10, i) accuracy

[newbie-02 at gmx dot de]

https://sourceware.org/bugzilla/show_bug.cgi?id=28472

--- Comment #23 from b. <newbie-02 at gmx dot de> ---
@Wilco, thanks for your comment, would you mind sharing  
a.) a code pointer,  
b.) a hint which method / tool is used / would be accepted  
to compare performance?  
c.) 'just use a correctly rounded math library',  
a hint which and how to?  
( I'm not! a pro in coding and 'system', just good in  
spotting weak points. )

-- 
You are receiving this mail because:
You are on the CC list for the bug.

[Bug math/28472] pow(10, i) accuracy

[vincent-srcware at vinc17 dot net]

https://sourceware.org/bugzilla/show_bug.cgi?id=28472

--- Comment #24 from Vincent Lefèvre <vincent-srcware at vinc17 dot net> ---
(In reply to b. from comment #23)
> c.) 'just use a correctly rounded math library',  
> a hint which and how to?  

You can look at CORE-MATH, already mentioned several times in this bug:
https://core-math.gitlabpages.inria.fr/

But note that it is currently not portable.

-- 
You are receiving this mail because:
You are on the CC list for the bug.

[Bug math/28472] pow(10, i) accuracy

[wdijkstr at arm dot com]

https://sourceware.org/bugzilla/show_bug.cgi?id=28472

--- Comment #25 from Wilco <wdijkstr at arm dot com> ---
(In reply to b. from comment #23)
> @Wilco, thanks for your comment, would you mind sharing  
> a.) a code pointer,  
> b.) a hint which method / tool is used / would be accepted  
> to compare performance?  
> c.) 'just use a correctly rounded math library',  
> a hint which and how to?  
> ( I'm not! a pro in coding and 'system', just good in  
> spotting weak points. )

The implementation [1] comes from Arm Optimized Routines and is believed to be
the best in the industry. It has been ported to GLIBC, Android and MUSL.

Note pow() doesn't at any point deal with integers. It literally implements pow
(x, y) as exp (log (x) * y) with higher internal precision. You'd have to
increase internal precision further with larger tables and polynomials - but
even that is no guarantee it rounds all cases the way you want. You need
correctly rounded math for that.

Both Optimized Routines and GLIBC have internal tests and benchmarks.

[1] https://github.com/ARM-software/optimized-routines/blob/master/math/pow.c

-- 
You are receiving this mail because:
You are on the CC list for the bug.