[Maxima] SF [2159499] Full bigfloat precision for Gamma after the second call
van Nek
van.nek at arcor.de
Sat Oct 18 06:39:31 CDT 2008
Am 16 Oct 2008 um 10:35 hat Raymond Toy geschrieben:
> 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)))
Now I got it. What a clever idea. Again I have learned something.
Well, I did some timing experiments:
I checked four versions of fprt18231_
(defun fprt18231_1 nil
(let ((a 1823176476672000)
(n 42698670666)
(d 1000)
h )
(do ((prec 32 (* 2 prec)))
((> prec fpprec))
(setq h n)
(setq n (+ (* n n) (* a d d)))
(setq d (* 2 h d)) )
(fpquotient (intofp n) (intofp d)) ))
(defun fprt18231_2 nil
(let ((a 1823176476672000))
(setq a (ash a (ash fpprec 1)))
(fpquotient (intofp (car (iroot a 2))) (intofp (ash 1 fpprec))) ))
(defun fprt18231_3 nil
(let ((a 1823176476672000))
(setq a (ash a (ash fpprec 1)))
(fpquotient (intofp (isqrt a)) (intofp (ash 1 fpprec))) ))
(defun fprt18231_4 nil
(let ((a 1823176476672000))
(setq a (ash a (- (ash fpprec 1) 52))) ;; we already have at least 52 bits precision
(fpquotient (intofp (isqrt a)) (intofp (ash 1 (- fpprec 26)))) ))
result: version 3 and 4 are about the same concerning time, version 2 is significantly slower,
so I tested version 1 und 3 more in detail.
I always used precisions $fpprec, where fpprec is just 32*2^n + some bits, which is the worst
case for version 1.
On Windows / GCL:
precs: [616, 1233, 2466, 4932, 9864, 19728, 39456]$
time1: [0.17, 0.49, 1.22, 3.18, 8.6, 23.92, 70.59]$
time3: [0.48, 1.17, 3.52, 11.54, 37.43, 116.63, 368.76]$
(in all cases 1000 runs, so time for a single run in millisec)
So version 1 is of type O(x^1.45) and version 3 of type O(x^1.6).
With GCL there is a significant advantage for version 1.
On Ubuntu / clisp (same computer):
precs: [616, 1233, 2466, 4932, 9864, 19728, 39456]$
time1: [0.356022, 0.748047, 2.016126, 5.820365, 18.08513, 58.195637, 195.188199]$
time3: [0.192013, 0.308018, 0.816051, 2.632164, 9.720607, 37.12232, 148.445278]$
Here version 1 is of type O(x^1.53) and version 3 of type O(x^1.66).
If I assume that one might not always hit the worst case 32*2^n+1 for version 1, on clisp
both versions are about the same. 1 is better in high precs and 3 in small precs.
Volker
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://www.math.utexas.edu/pipermail/maxima/attachments/20081018/f7b9c148/attachment.htm
More information about the Maxima
mailing list