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

>>>>> "van" == van Nek <van.nek at arcor.de> writes:

    van> 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.

    van> Ray,

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

I was too brief.  Let a=18231....  Then compute isqrt(a*2^(2*n)).
That gives the numerator of the root.  The denominator is 2^n.  Then
do what you did at the end of fprt18231:

(fpquotient (intofp <isqrt>) (intofp (ash 1 n)))

This could probably be optimized since dividing by a power of 2 just
subtracts from the exponent.

    van> It took me some minutes to remember what the idea of my implementation of fppi1 was. 
    van> Perhaps I should have written a better documentation.

You at least gave the formula you used!  That's a lot better than
many, many parts of maxima.  Some of the magic constants in comppi
could have been documented though.  But with the formula, it doesn't
take long to figure out what the constants are.

    van> 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.

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

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

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

Well, you said you wanted to do some timing.  It just seemed to me the
worst case timing would happen if fpprec was just one bit over 32*2^n.
Then at iteration n, prec would be 32*2^n, which is just one bit too
small, so an additional iteration would be needed.  The extra
iteration would now operate on numbers with twice as many bits, and
hence take much longer in comparison to the case where fpprec equal to
32*2^n.  And then you would throw away almost all of the extra bits
produced.

Perhaps using isqrt or iroot doesn't have this issue of possibly
computing almost twice the number of bits that are needed.

Just a guess.  I didn't actually try any of this.

Ray