Karim Belabas on Sun, 08 Oct 2017 10:35:12 +0200

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

* 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]
`