[Maxima] SF  Full bigfloat precision for Gamma after the second call
raymond.toy at ericsson.com
Wed Oct 15 11:34:49 CDT 2008
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.
> 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.
More information about the Maxima