The listed bug says bfloat(%gamma) doesn't give a bfloat version of
%gamma, the Euler-Mascheroni constant.
Here is a patch to float.lisp that implements this. Don't know if
this should go into the release or not. Relatively safe since it adds
a missing item without changing anything else. It includes
documentation on the algorithm too. (Gasp!)
On my Ultra 30 300MHz workstation, bfloat(%gamma) with fpprec = 500
takes about 7 sec, so it's reasonably fast.
Ray
Index: float.lisp
===================================================================
RCS file: /cvsroot/maxima/maxima/src/float.lisp,v
retrieving revision 1.3
diff -u -r1.3 float.lisp
--- float.lisp 27 Jun 2002 17:37:20 -0000 1.3
+++ float.lisp 2 Oct 2002 13:49:48 -0000
@@ -88,6 +88,8 @@
"Bigfloat representation of %E")
(DEFMVAR BIGFLOAT%PI '((BIGFLOAT SIMP 56.) 56593902016227522. 2)
"Bigfloat representation of %PI")
+(defmvar bigfloat%gamma '((bigfloat simp 56.) 41592772053807304. 0)
+ "Bigfloat representation of %GAMMA")
;; Internal specials
@@ -111,8 +113,11 @@
(DEFVAR *DECFP NIL)
+;; These hold the most precise values of these contants that have been
+;; computed so far.
(DEFVAR MAX-BFLOAT-%PI BIGFLOAT%PI)
(DEFVAR MAX-BFLOAT-%E BIGFLOAT%E)
+(defvar max-bfloat-%gamma bigfloat%gamma)
(declare-top (SPECIAL *CANCELLED $FLOAT $BFLOAT $RATPRINT $RATEPSILON
$DOMAIN $M1PBRANCH ADJUST)
@@ -410,7 +415,7 @@
(DEFMFUN $BFLOAT (X)
(LET (Y)
(COND ((BIGFLOATP X))
- ((OR (NUMBERP X) (MEMQ X '($%E $%PI)))
+ ((OR (NUMBERP X) (MEMQ X '($%E $%PI $%gamma)))
(BCONS (INTOFP X)))
((OR (ATOM X) (MEMQ 'array (CDAR X)))
(IF (EQ X '$%PHI)
@@ -531,6 +536,7 @@
((EQUAL 0 L) '(0 0))
((EQ L '$%PI) (FPPI))
((EQ L '$%E) (FPE))
+ ((eq l '$%gamma) (fpgamma))
(T (LIST (FPROUND L) (PLUS *M FPPREC)))))
;; It seems to me that this function gets called on an integer
@@ -692,6 +698,15 @@
((< FPPREC (CADDAR MAX-BFLOAT-%PI))
(CDR (SETQ BIGFLOAT%PI (BIGFLOATP MAX-BFLOAT-%PI))))
(T (CDR (SETQ MAX-BFLOAT-%PI (SETQ BIGFLOAT%PI (FPPI1)))))))
+
+(defun fpgamma ()
+ (cond ((= fpprec (caddar bigfloat%gamma)) (cdr bigfloat%gamma))
+ ((< fpprec (caddar bigfloat%gamma))
+ (cdr (setq bigfloat%gamma (bigfloatp bigfloat%gamma))))
+ ((< fpprec (caddar max-bfloat-%gamma))
+ (cdr (setq bigfloat%gamma (bigfloatp max-bfloat-%gamma))))
+ (t (cdr (setq max-bfloat-%gamma (setq bigfloat%gamma (fpgamma1)))))))
+
(DEFUN FPONE NIL
(COND (*DECFP (INTOFP 1)) ((= FPPREC (CADDAR BIGFLOATONE)) (CDR BIGFLOATONE))
@@ -709,12 +724,81 @@
(SETQ B (*QUO (TIMES B (f1- I) (f1- I))
(TIMES 4 I (f1+ I))))
(SETQ C (PLUS C B)))
- (RETURN C)))
+ (RETURN C)))
(DEFUN FPPI1 NIL
(BCONS (LIST (FPROUND (COMPPI (PLUS FPPREC 3))) (PLUS -3 *M))))
-(DEFmfUN FPMAX NA
+;; Compute the main part of the Euler-Mascheroni constant using the
+;; Bessel function approach. See
+;; http://numbers.computation.free.fr/Constants/Gamma/gamma.html for a
+;; description of the algorithm.
+;; Roughly, we have
+;;
+;; %gamma = A(N)/B(N) - log(N) + O(e^(-4*N))
+;;
+;; where
+;;
+;;
+;; b*N
+;; A(N) = sum (N^2/n!)^2*H(n)
+;; n=0
+;;
+;; b*N
+;; B(N) = sum (N^2/n!)^2
+;; n=0
+;;
+;; n
+;; H(n) = sum 1/k
+;; k=1
+;;
+;; with H(0) = 0
+;;
+;; and b = 4.970625759... where b*(log(b)-1) = 3.
+;;
+;; This formula can be easily justified by looking at the value
+;; K0(2*N)/I0(2*N). We see that
+;;
+;; K0(2*N)/I0(2*N) = -log(N) - %gamma + A/B
+;;
+;; where A(N) -> A and B(N) -> B as N -> infinity.
+;;
+;; And assuming N is large, the asymptotic expressions for K0 and I0 show that
+;; K0(2*N)/I0(2*N) ~ exp(-4*N)/N.
+;;
+(defun comp-bf%gamma (prec)
+ ;; Prec is the number of bits we want. We assume the remainder is
+ ;; really e^(-4*N) and not O(e^(-4*N)).
+ ;;
+ ;; We also assume don't need a really precise value of beta because
+ ;; our N's are not so big that we need more.
+ (let* ((big-n (floor (* 1/4 prec (log 2d0))))
+ (big-n-sq (cdr ($bfloat (* big-n big-n))))
+ (beta 4.9706257595442319023d0)
+ (limit (floor (* beta big-n)))
+ (term (cdr bigfloatone))
+ (harmonic (cdr bigfloatzero))
+ (a-sum (cdr bigfloatzero))
+ (b-sum (cdr bigfloatone)))
+ (do ((n 1 (1+ n)))
+ ((> n limit))
+ (let ((bf-n (cdr ($bfloat n))))
+ (setf term (fpquotient (fptimes* term big-n-sq)
+ (fptimes* bf-n bf-n)))
+ (setf harmonic (fpplus harmonic
+ (fpquotient (fpone) bf-n)))
+ (setf a-sum (fpplus a-sum
+ (fptimes* term harmonic)))
+ (setf b-sum (fpplus b-sum term))))
+ (fpplus (fpquotient a-sum b-sum)
+ (fpminus (cdr ($bfloat (list '(%log simp) big-n)))))
+ ))
+
+(defun fpgamma1 ()
+ ;; Use a few extra bits of precision
+ (bcons (list (fpround (first (comp-bf%gamma (plus fpprec 5)))) 0)))
+
+(DEFMfUN FPMAX NA
(PROG (MAX)
(SETQ MAX (ARG 1))
(DO ((I 2 (f1+ I))) ((> I NA))