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.

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.

"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.

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.

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.

> 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

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_*).

  
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.

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 ... )

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.

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?

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.

(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.

(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.

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

(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".

(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.

(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.

(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.

(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.

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.

(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.

@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. )

(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.

(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