Man page issues: logb, significand, cbrt, log2, log10, exp10

Man page issues: logb, significand, cbrt, log2, log10, exp10

[Wilco Dijkstra]

Hi,

> FWIW, it appears that the author of the glibc exp10 implementation
> agrees with me that the implementation is sub-standard:

As Adhemerval pointed out, that was not the implementation used in GLIBC.
The version that was used, wasn't perfect, but still better than MUSL:

> Compare with musl:
>
> https://github.com/rofl0r/musl/blob/master/src/math/exp10.c

The worst case error of MUSL is 4.14 ULP vs 2.01 ULP of the old GLIBC exp10 [1].

And the new exp10 in GLIBC is 0.513 ULP [2].

As I pointed out in the PR, it would be reasonable for compilers to convert
pow (10, x) into exp10 (x) when it is known there is a good implementation
available - in GLIBC, exp10 is both faster and slightly more accurate than pow.

Cheers,
Wilco

[1] https://members.loria.fr/PZimmermann/papers/glibc238-20230921.pdf
[2] https://members.loria.fr/PZimmermann/papers/accuracy.pdf

Re: Man page issues: logb, significand, cbrt, log2, log10, exp10

[Paul Zimmermann]

       Hi Wilco,

> The worst case error of MUSL is 4.14 ULP vs 2.01 ULP of the old GLIBC exp10 [1].
> 
> And the new exp10 in GLIBC is 0.513 ULP [2].

yes the new exp10 is more accurate, nice work!

> As I pointed out in the PR, it would be reasonable for compilers to convert
> pow (10, x) into exp10 (x) when it is known there is a good implementation
> available - in GLIBC, exp10 is both faster and slightly more accurate than pow.

ok at compile time, but there should be a way to disable this optimization,
for example if you want to compare pow(10,x) and exp10(x). (At run time
such a conversion would slow down the critical path of pow.)

Paul

Re: Man page issues: logb, significand, cbrt, log2, log10, exp10

[Adhemerval Zanella Netto]

On 02/03/24 23:02, Morten Welinder wrote:
> Thanks.
> 
> There is (was?) already crlibm out there.
> https://core-math.gitlabpages.inria.fr/  No particular need for wheel
> reinvention here.
> 
> FWIW, it appears that the author of the glibc exp10 implementation
> agrees with me that the implementation is sub-standard:
> 
> https://codebrowser.dev/glibc/glibc/math/e_exp10.c.html

This code was not used by any port and we recently removed to avoid this
very issue [1]. The exp10 implementation used by all ports
is at sysdeps/ieee754/dbl-64/e_exp10.c (i386/m68k are exceptions and my
plan to eventually remove this implementation in favor the generic
one [2]).

And the exp10 implementation was recently improved [3], with the
author suggesting the worst-case error in round-to- should be
nearest to 0.513 ULP (I am not sure if he did some empirical testing
to verify this value, at least with libm-test-ulps the resuts are
bounded to 2 ulp max).

> 
> /* This is a very stupid and inprecise implementation. It'll get
> replaced sometime (soon?). */
> return __ieee754_exp (M_LN10 * arg);
> 
> 
> Compare with musl:
> 
> https://github.com/rofl0r/musl/blob/master/src/math/exp10.c
> 

[1] https://sourceware.org/git/?p=glibc.git;a=commit;h=9c61303ebbdc6e727c89591bff3229c9fbfa438b
[2] https://sourceware.org/pipermail/libc-alpha/2024-January/154107.html
[3] https://sourceware.org/git/?p=glibc.git;a=commit;h=63d0a35d5f223a3f4b68190567b7d4d44545bce5

Re: Man page issues: logb, significand, cbrt, log2, log10, exp10

[Morten Welinder]

Sorry to be a bit of a pain.

Some testing says that the average error from exp10 is 300-500 times
bigger than the average error from pow(10,.). This is consistent
across a large range of arguments.

The maximum error from pow(10,.) in the samples is 1ulp (relative to a
reference rounded value, so the true error is likely less). The
maximum error from exp10 in the samples is 2ulp.

This is not surprising given that exp10 starts out by introducing a
rounding error due to the multiplication before the call to exp.

I didn't bother looking, but it is almost certainly true that there
are arguments for which exp10 is better than pow(10,.).  However, the
numbers below imply that such arguments are very rare compared to the
other way around.


-------------------------- average ----- max
Binade -7      exp10: 0.3378           1
Binade -7  pow(10,.): 0.0007           1
Binade -6      exp10: 0.3429           2
Binade -6  pow(10,.): 0.0007           1
Binade -5      exp10: 0.3532           2
Binade -5  pow(10,.): 0.0008           1
Binade -4      exp10: 0.3774           2
Binade -4  pow(10,.): 0.0008           1
Binade -3      exp10: 0.4402           2
Binade -3  pow(10,.): 0.0010           1
Binade -2      exp10: 0.4118           2
Binade -2  pow(10,.): 0.0009           1
Binade -1      exp10: 0.4228           2
Binade -1  pow(10,.): 0.0009           1
Binade  0      exp10: 0.4204           2
Binade  0  pow(10,.): 0.0009           1
Binade  1      exp10: 0.4221           2
Binade  1  pow(10,.): 0.0009           1
Binade  2      exp10: 0.4204           2
Binade  2  pow(10,.): 0.0009           1
Binade  3      exp10: 0.4222           2
Binade  3  pow(10,.): 0.0009           1
Binade  4      exp10: 0.4209           2
Binade  4  pow(10,.): 0.0009           1
Binade  5      exp10: 0.4200           2
Binade  5  pow(10,.): 0.0009           1
Binade  6      exp10: 0.4210           2
Binade  6  pow(10,.): 0.0009           1
Binade  7      exp10: 0.4210           2
Binade  7  pow(10,.): 0.0009           1

Notes:
1. Only positive arguments tested.
2. powl (long double version of pow) is used as a reference.  I
double-checked with a double-double version of pow (good to about
100ish bits) that this does not matter.




#define _GNU_SOURCE   1
#include <stdio.h>
#include <math.h>
#include <stdint.h>


static uint64_t
murmur64 (uint64_t h)
{
  h ^= h >> 33;
  h *= 0xff51afd7ed558ccdll;
  h ^= h >> 33;
  h *= 0xc4ceb9fe1a85ec53ll;
  h ^= h >> 33;
  return h;
}

static double
exp10ref (double x)
{
  volatile double y = (double)(powl (10.0l, x));
  return y;
}

static double
exp10viapow (double x)
{
  return pow (10, x);
}


static void
test_binade (int b, double (*f) (double), const char *funcname)
{
  uint64_t h = 0x0123456701234567ll;

  double ulps = 0;
  double mulp = 0;
  int N = 1000000;

  for (int i = 0; i < N; i++) {
    h = murmur64 (h);

    double x = ldexp ((h & 0xfffffffffffffll) | 0x10000000000000ll, b - 52);

    double y = f (x);
    double yref = exp10ref (x);
    double dy = fabs (y - yref);

    double ulp = dy / (nextafter (yref, INFINITY) - yref);
    ulps += ulp;
    if (ulp > mulp) mulp = ulp;
  }

  printf ("Binade %2d %10s: %6.4f  %10.0f\n",
      b, funcname, ulps / N, mulp);
}


int
main ()
{
  for (int b = -7; b <= 7; b++) {
    test_binade (b, exp10, "exp10");
    test_binade (b, exp10viapow, "pow(10,.)");
  }

  return 0;
}

Re: Man page issues: logb, significand, cbrt, log2, log10, exp10

[Alejandro Colomar]

[-- Attachment #1: Type: text/plain, Size: 1977 bytes --]

Hi Vincent,

On Sun, Mar 03, 2024 at 12:46:00PM +0100, Vincent Lefevre wrote:
> On 2024-03-03 03:21:26 +0100, Alejandro Colomar wrote:
> > Maybe just add some headers to core-math, and package it as a
> > standalone library.
> 
> The issue is that it is not portable yet.

Well, one could package it just to the systems to which it is portable,
if that's useful.  That's why a standalone library has more chances of
being available soon than glibc.  You'd need to make it portable (and
other things) to put it in glibc; but if you say "here's libcore-math,
avaiable only in XXX systems", you could get distros to distribute it
already.  And then provide headers that don't clash with glibc, such as
<core-math/*.h> or whatever.

> 
> > > FWIW, it appears that the author of the glibc exp10 implementation
> > > agrees with me that the implementation is sub-standard:
> > > 
> > > https://codebrowser.dev/glibc/glibc/math/e_exp10.c.html
> > > 
> > > /* This is a very stupid and inprecise implementation. It'll get
> > > replaced sometime (soon?). */
> > > return __ieee754_exp (M_LN10 * arg);
> > 
> > Hmmm.  Still, it's simple.  If pow(10, x) is strictly better, maybe one
> > can prove it and send a patch.  Or for something better, it'll take more
> > work.
> 
> If by "strictly better", you mean that for each input, it returns a
> result that is at least as accurate as the one returned by the above
> expression, then, probably no. The reason is that the rounding errors
> in the above expression may partly compensate on a random basis. So,
> for some proportion of inputs, you'll actually get an accurate result.
> And unless pow is designed to be almost correctly rounded, it will
> probably be sometimes worse.

Then glibc's current code is good, I guess.  It's simple, and works for
most programs.

Have a lovely day!
Alex

-- 
<https://www.alejandro-colomar.es/>
Looking for a remote C programming job at the moment.

[-- Attachment #2: signature.asc --]
[-- Type: application/pgp-signature, Size: 833 bytes --]

Re: Man page issues: logb, significand, cbrt, log2, log10, exp10

[Vincent Lefevre]

On 2024-03-03 03:21:26 +0100, Alejandro Colomar wrote:
> Hi Morten,
> 
> On Sat, Mar 02, 2024 at 09:02:24PM -0500, Morten Welinder wrote:
> > Thanks.
> > 
> > There is (was?) already crlibm out there.
> > https://core-math.gitlabpages.inria.fr/  No particular need for wheel
> > reinvention here.
> 
> crlibm doesn't seem to exist anymore.

The sources are still available at more non-official mirror, but
it is no longer maintained.

> Maybe just add some headers to core-math, and package it as a
> standalone library.

The issue is that it is not portable yet.

> > FWIW, it appears that the author of the glibc exp10 implementation
> > agrees with me that the implementation is sub-standard:
> > 
> > https://codebrowser.dev/glibc/glibc/math/e_exp10.c.html
> > 
> > /* This is a very stupid and inprecise implementation. It'll get
> > replaced sometime (soon?). */
> > return __ieee754_exp (M_LN10 * arg);
> 
> Hmmm.  Still, it's simple.  If pow(10, x) is strictly better, maybe one
> can prove it and send a patch.  Or for something better, it'll take more
> work.

If by "strictly better", you mean that for each input, it returns a
result that is at least as accurate as the one returned by the above
expression, then, probably no. The reason is that the rounding errors
in the above expression may partly compensate on a random basis. So,
for some proportion of inputs, you'll actually get an accurate result.
And unless pow is designed to be almost correctly rounded, it will
probably be sometimes worse.

-- 
Vincent Lefèvre <vincent@vinc17.net> - Web: <https://www.vinc17.net/>
100% accessible validated (X)HTML - Blog: <https://www.vinc17.net/blog/>
Work: CR INRIA - computer arithmetic / AriC project (LIP, ENS-Lyon)

Re: Man page issues: logb, significand, cbrt, log2, log10, exp10

[Alejandro Colomar]

[-- Attachment #1: Type: text/plain, Size: 985 bytes --]

Hi Morten,

On Sat, Mar 02, 2024 at 09:02:24PM -0500, Morten Welinder wrote:
> Thanks.
> 
> There is (was?) already crlibm out there.
> https://core-math.gitlabpages.inria.fr/  No particular need for wheel
> reinvention here.

crlibm doesn't seem to exist anymore.  Maybe just add some headers to
core-math, and package it as a standalone library.

> FWIW, it appears that the author of the glibc exp10 implementation
> agrees with me that the implementation is sub-standard:
> 
> https://codebrowser.dev/glibc/glibc/math/e_exp10.c.html
> 
> /* This is a very stupid and inprecise implementation. It'll get
> replaced sometime (soon?). */
> return __ieee754_exp (M_LN10 * arg);

Hmmm.  Still, it's simple.  If pow(10, x) is strictly better, maybe one
can prove it and send a patch.  Or for something better, it'll take more
work.


Have a lovely night!
Alex

-- 
<https://www.alejandro-colomar.es/>
Looking for a remote C programming job at the moment.

[-- Attachment #2: signature.asc --]
[-- Type: application/pgp-signature, Size: 833 bytes --]

Re: Man page issues: logb, significand, cbrt, log2, log10, exp10

[Morten Welinder]

Thanks.

There is (was?) already crlibm out there.
https://core-math.gitlabpages.inria.fr/  No particular need for wheel
reinvention here.

FWIW, it appears that the author of the glibc exp10 implementation
agrees with me that the implementation is sub-standard:

https://codebrowser.dev/glibc/glibc/math/e_exp10.c.html

/* This is a very stupid and inprecise implementation. It'll get
replaced sometime (soon?). */
return __ieee754_exp (M_LN10 * arg);


Compare with musl:

https://github.com/rofl0r/musl/blob/master/src/math/exp10.c

Re: Man page issues: logb, significand, cbrt, log2, log10, exp10

[Alejandro Colomar]

[-- Attachment #1: Type: text/plain, Size: 2354 bytes --]

[CC += glibc, and those involved in the glibc bug report]

Hi Morten,

On Sat, Mar 02, 2024 at 04:17:36PM -0500, Morten Welinder wrote:
> I think what happens is that the compiler (not glibc) computes that
> exp10 for you and that the compiler happens to be more accurate.
> Here's what I get for the loop:
> 
>   for (int i = 1; i < 20; i++) {
>     printf ("%.20g\n", exp10 (i));
>   }
> 
> welinder@CarbonX1:~$ ./a.out
> 10
> 100
> 1000.0000000000001137
> 10000.000000000001819
> 100000
> 1000000
> 9999999.9999999981374
> 99999999.999999985099
> 999999999.99999988079
> 10000000000
> 100000000000
> 1000000000000
> 10000000000000
> 100000000000000
> 1000000000000000
> 10000000000000000
> 99999999999999984
> 1000000000000000000
> 10000000000000000000
> 
> Here's the bug report to go with this:
> https://sourceware.org/bugzilla/show_bug.cgi?id=28472
> Note comment 6.  It is clearly not a high-priority item for glibc.

Thanks for that link.

I agree with glibc that the standard specifies that these functions need
not be precise.  That lost precission probably results in better
performance.  Most programs won't care that these functions are
inaccurate.

If you need a correctly-rounded version of these functions, which is
perfectly reasonable, the right thing to ask is that libc implements
the cr_ version of these functions.

I also understand that adding functions to glibc isn't straightforward,
so glibc maintainers have reasons to not do it at the moment.  In fact,
lately I've been leaning towards thinking that libc is a huge monster to
which nothing more should be added, at all.

How about writing a new library --maybe call it libm-cr, maybe
libm-cr-pow, maybe libm-cr-exp10, depending on how extensive you want
it-- and add cr_exp10(3) to that library?  You could do that, and just
support the systems you need to support.  The effort would be certainly
smaller than adding the function to glibc.

Regarding the manual pages, I don't remember from the top of my head if
there's any page documenting that libm functions are imprecise.  I would
prefer documenting it in one place, rather than adding caveats to every
libm page.

> 
> M.

Have a lovely night!
Alex

-- 
<https://www.alejandro-colomar.es/>
Looking for a remote C programming job at the moment.

[-- Attachment #2: signature.asc --]
[-- Type: application/pgp-signature, Size: 833 bytes --]