[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
```