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

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

* Karim Belabas [2017-10-08 10:35]:
> * 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.
[...]
> 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...

N.B. the "simpler" v = log(one + v) is a little faster but less
accurate, due to loss of accuracy in the log [ cf. P.S. above ].
The proposed solution is stabler.

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