* Bill Allombert [2016-11-02 14:43]:
> On Tue, Jul 05, 2016 at 12:15:35PM -0400, Charles Greathouse wrote:
> > Package: pari
> > Version: 2.8.0
> > Severity: wishlist
> >
> > 2.6.2 introduced the GP function expm1. It would be nice to have a
> > companion function log1p which computes log(1 + x) with better precision
> > when |x| is small.
>
> I want to comment that both expm1 and log1p are available in math.h, and
> there is a few places in PARI we might want to use log1p(x) instead of
> log(1+x).  This is assuming log1p is portable enough.
>
> All the more reason to provide log1p in GP too.

Sorry about the delay...

I started implementing this today and noticed that the implementation
would be as simple as log(1 + x) because we already have better
precision when |x| is small: indeed, 1 + x increases the accuracy by
max(0, -expo(x)) bits, exactly the amount needed to make up for the
cancellation.

 default(realprecision, 10000); x = 1e-100;
 a = log1p(x);
 b = log(1.0+x);
 default(realprecision, 10040); x = 1e-100;
 c = log1p(x);  \\ reference point

 ? abs(a-c)/c  \\ relative error in log1p(x)
 %6 = 1.4027986153764837216 E-10019

So the C implementation

GEN
glog1p(GEN x, long prec)
{
  pari_sp av = avma;
  return gerepileupto(av, glog(gaddgs(x,1), prec));
}

would actually be enough. There is a very minor speed gain to be obtained
for tiny real x and low accuracy (else we use AGM anyway) because we could
replace

  y := ((1 + x) - 1) / ((1 + x) + 1)  followed by summation on y

by

  y := x / (2 + x)  followed by summation on y

i.e. one addition saved (the other, 1 + x, being essentially a copy),
this is a wee bit stabler as well (1ulp...).


To conclude, the only reason I see we might want to provide log1p is for
symmetry with expm1. And possibly prevent misguided users from using
log(1.0 + x)  [ which *does* suffer from cancellation ]

Any comments ?

Cheers,

    K.B.
--
Karim Belabas, IMB (UMR 5251)  Tel: (+33) (0)5 40 00 26 17
Universite de Bordeaux         Fax: (+33) (0)5 40 00 21 23
351, cours de la Liberation    http://www.math.u-bordeaux.fr/~kbelabas/
F-33405 Talence (France)       http://pari.math.u-bordeaux.fr/  [PARI/GP]
`