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

[Dirk Laurie <dirk.laurie@gmail.com>]

2017-10-08 10:35 GMT+02:00 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.

The functions log(1+x) and (exp(x)-1) are textbook examples of
numerically unstable expressions. Equivalent stable expressions
were designed by W (Velvel) Kahan and built into the instruction
set of the 8087 coprocessor in the 1980s. They are available in C99
as log1p and expm1.

In the case of log(1+x), a simplified version of what you do is:

log1p(x) = my(w=1+x,h=w-1); if(h,log(w)*x/h,x)

If you code that in Python, it will have the correct accuracy.
If you code it in Pari, it will be equivalent to Karim's suggestion
since x/h will equal 1.0 to full accuracy,