Dirk Laurie on Sat, 14 Oct 2017 10:00:22 +0200 |
[Date Prev] [Date Next] [Thread Prev] [Thread Next] [Date Index] [Thread Index]
Re: a(n+1) = log(1+a(n)) |
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,