[Maxima-commits] CVS: maxima/src ellipt.lisp,1.48,1.49

[Raymond T. <rt...@us...>]

Update of /cvsroot/maxima/maxima/src
In directory 23jxhf1.ch3.sourceforge.com:/tmp/cvs-serv16616/src

Modified Files:
	ellipt.lisp 
Log Message:
src/ellipt.lisp:
o Fix typo.  It's A&S 16.14.3, not 16.14.13
o Remove *bferrtol*
o Replace *bferrtol* with the function bferrtol that computes the
  tolerance as a function of fpprec.  sqrt(2^(-fpprec)) seems to work.
o Use (bferrtol) instead of *bferrtol*.

tests/rtest_elliptic.mac:
o Fix bug in test_table: Don't set numer to true.
o For some reason jacobi_cn tests needs 6b-14 instead of 3b-14.
o Add test for elliptic_kc for known analytic values using 100 digits
  of precision.


Index: ellipt.lisp
===================================================================
RCS file: /cvsroot/maxima/maxima/src/ellipt.lisp,v
retrieving revision 1.48
retrieving revision 1.49
diff -u -d -r1.48 -r1.49
--- ellipt.lisp	20 Feb 2009 04:10:47 -0000	1.48
+++ ellipt.lisp	21 Feb 2009 21:51:44 -0000	1.49
@@ -293,7 +293,7 @@
 	 ;; cn(u,m) = sech(u) - 1/4*(1-m)*(sinh(u)*cosh(u)-u)*tanh(u)*sech(u)
 	 (/ (cosh u)))
 	(t
-	 ;; Use the ascending Landen transformation, A&S 16.14.13.
+	 ;; Use the ascending Landen transformation, A&S 16.14.3.
 	 (multiple-value-bind (v mu root-mu1)
 	     (ascending-transform u m)
 	   (let ((d (dn v mu)))
@@ -2268,7 +2268,11 @@
 ;;
 
 
-(defvar *bferrtol* (bigfloat 1d-6))
+(defun bferrtol ()
+  ;; Compute error tolerance as sqrt(2^(-fpprec)).  Not sure this is
+  ;; quite right, but it makes the routines more accurate as fpprec
+  ;; increases.
+  (bigfloat (sqrt (scale-float (bigfloat 1) (- maxima::fpprec)))))
 
 ;; rc(x,y) = integrate(1/2*(t+x)^(-1/2)/(t+y), t, 0, inf)
 ;;
@@ -2294,7 +2298,7 @@
 	   (setf z yn)
 	   (setf w (bigfloat 1))))
     (setf a (/ (+ xn yn yn) 3))
-    (setf epslon (/ (abs (- a xn)) *bferrtol*))
+    (setf epslon (/ (abs (- a xn)) (bferrtol)))
     (setf an a)
     (setf pwr4 (bigfloat 1))
     (setf n 0)
@@ -2335,7 +2339,7 @@
 	 (epslon (/ (max (abs (- a xn))
 			 (abs (- a yn))
 			 (abs (- a zn)))
-		    *bferrtol*))
+		    (bferrtol)))
 	 (an a)
 	 (sigma 0)
 	 (power4 1)
@@ -2391,7 +2395,7 @@
 	 (epslon (/ (max (abs (- a xn))
 			 (abs (- a yn))
 			 (abs (- a zn)))
-		    *bferrtol*))
+		    (bferrtol)))
 	 (an a)
 	 (power4 1)
 	 (n 0)
@@ -2439,7 +2443,7 @@
 			 (abs (- a yn))
 			 (abs (- a zn))
 			 (abs (- a pn)))
-		    *bferrtol*))
+		    (bferrtol)))
 	 (an a)
 	 xnroot ynroot znroot pnroot lam dn)
     (loop while (> (* power4 epslon) (abs an))