SF [2159499] Full bigfloat precision for Gamma after the second call

Am 15 Oct 2008 um 12:34 hat Raymond Toy geschrieben:

> Also, it occurred to me that you are computing sqrt(a).  Could we not
> use the Lisp routine isqrt to compute sqrt(a*2^(2*n))?  (isqrt computes
> an integer such that when squared it is just less than the argument.)
> Then the answer would be isqrt(a*2^(2*n))/2^n.  Of course, for this to
> be fast, the Lisp implementation needs a fast isqrt.  I think cmucl and
> sbcl have fast isqrt.

Ray,

I need sqrt(640320^3/12^2) as a bigfloat number!

It took me some minutes to remember what the idea of my implementation of fppi1 was. 
Perhaps I should have written a better documentation.
I had spent lots of days in checking lots of types of pi computation algorithms. 
There are a great number of algorithms which theoretically should be much faster than 
the one by Chudnovsky I chose. But it turned out that all algorithms which make widely 
use of bigfloat computations (e.g. Borwein iterations) are slow in Maxima. 
Maybe Maxima's bigfloat implementation isn't very fast, maybe bigfloat computations 
are generally slow.
I found that arctan series and series of Ramanujan type which can be based on bignum 
arithmetic fit best to Maxima. The Chudnovsky series I chose is of Ramanujan type and 
I could limit the number of bigfloat computations to 2 calls to fpquotient. 

Your next e-mail:
> Sounds good.  If you're going to do some speed tests, I think it would
> be good to try some with fpprec just over 36*2^n.  This is the case
> where, after n iterations, we're just below the desired precision, and
> we need another iteration, effectively doubling the number of bits in
> the computation.

I just don't get the idea of what you intend here. 

(defun fprt18231_ nil  
  (let ((a 1823176476672000)
        (n 42698670666) 
        (d 1000))   
     (do ((prec 32 (* 2 prec))) ;; 10 digits of precision guarantee 33,... bits
         ((> prec fpprec))
       (setq h n)
       (setq n (+ (* n n) (* a d d)))
       (setq d (* 2 h d)) )
     (fpquotient (intofp n) (intofp d)) ))

42698670666 / 1000 is an approximation of sqrt(1823176476672000) with at least 10 
correct digits. 10 digits guarantee 33,... bits. With a little margin I guarantee 32 bits of 
precision.
Each loop doubles the number of correct digits. The variable prec documents the actual 
precision. E.g. $fpprec=37 and fpprec=125. The loop stops when prec is 128. The precision 
of the result is 128 bit (not 64  - is that what you think? ), which can be proven.

(%i2) fpprec:37$
(%i3) ?fpprec;
(%o3)                                 125
(%i4) :lisp (fprt18231_)
(27063497749562510529621671475639075089 26)
(%i4) 27063497749562510529621671475639075089 * 2.0b0^(26-125);
(%o4)              4.269867066633339581771288916065960827b7
(%i5) sqrt(1823176476672000.0b0);
(%o5)              4.269867066633339581771288916065960827b7

Volker

> Ray
> _______________________________________________
> Maxima mailing list
> Maxima at math.utexas.edu
> http://www.math.utexas.edu/mailman/listinfo/maxima