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

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