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

van Nek wrote:
> Am 14 Oct 2008 um 18:18 hat Raymond Toy geschrieben:
> 
> 
>> *** At this point maxima wants 438 bits of pi.  But fprt18231_ sees
>>     $fpprec = 64 and therefore only computes 64 digits or about 215
>>     bits.  We really want 438 bits.  I think it was a misunderstanding
>>     about fpprec and $fpprec.  fpprec changes throughout the code and
>>     is dynamically changed to produce the desired precision.  I think
>>     sin is especially likely to increase fpprec so that it can do
>>     range reduction accurately.
> 
> Ray, 
> 
> I think you are right here. When I wrote the code I wasn't aware of possible interferences 
> with other functions concerning fpprec. More or less I only thought of a standalone 
> computation of pi. Sorry for this stupidness! (Just for the record: The standalone 

No need to apologize.  This stuff isn't very well documented and it
would be easy to make this mistake.
> 
> I believe the problem is fixed when $fpprec is replaced by something like fpprec/3.3
> So I propose the following fix. Please check it.
> 
> (defun fprt18231_ nil  
>   (let ((a 1823176476672000)
>         (n 42698670666) 
>         (d 1000)     
>         (p (ceiling fpprec 3.3)))    
>      (do ((prec 10 (* 2 prec)))
>          ((> prec p))
>        (setq h n)
>        (setq n (+ (* n n) (* a d d)))
>        (setq d (* 2 h d)) )
>      (fpquotient (intofp n) (intofp d)) ))

Why divide fpprec by 3.3?  To convert to back to digits?

In that case, why not start with prec = 36 (10 digits) and make the
condition (> prec fpprec)?

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