Re: a(n+1) = log(1+a(n))

[Karim Belabas <Karim.Belabas@math.u-bordeaux.fr>]

* Elim Qiu [2017-10-08 05:44]:
> I'm study the behavior of   n(n a(n) -2) / log(n)
> where a(1) > 0, a(n+1) = log(1+a(n))
> 
> Using Pari:
> 
> f(n) =
> { my(v = 2);
>   for(k=1,n, v = log(1+v));
>   return(n*(n*v -2) / (log(n)));
> }
> 
> It turns out the program runs very slowly. The same logic in python runs
> 100 time faster but not have the accuracy I need.
> 
> Any ideas?

You are hit by PARI's "poor man's interval arithmetic" (always compute
as correctly as the input allows). This works :

  g(n) =
  { my(v = 2, one = 1.0);
    for(k=1,n, v = log(1+v) * one);
    return(n*(n*v -2) / (log(n)));
  }

  (10:25) gp > g(10^6)
  time = 2,376 ms.
  %1 = 0.51488946924922388659082316377748262728

The reason why your original function is very slow is that, since v << 1,
1 + v has more bits of accuracy than v . So that the internal accuracy
increases quickly during your loop, slowing down the computation
immensely. Multiplying by 1.0 (at the original accuracy), I restore the
accuracy we expected.

Cheers,

    K.B.

P.S. It is true that log(1+v) should *lose* some accuracy since 1+v is
close to 1, but we are more conservative when reducing accuracy than
when increasing it. The net "gain" is 1 more word of accuracy per loop
iteration...
--
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]
`