[Maxima-commits] [git] Maxima CAS branch, master, updated. branch-5_30-base-181-gb9ab71f

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

This is an automated email from the git hooks/post-receive script. It was
generated because a ref change was pushed to the repository containing
the project "Maxima CAS".

The branch, master has been updated
       via  b9ab71f360b9a3e1e575a86bfca81fd4dacb4b50 (commit)
      from  a47e04d91675524ef4c08c05ceb7bf4e6aae0e61 (commit)

Those revisions listed above that are new to this repository have
not appeared on any other notification email; so we list those
revisions in full, below.

- Log -----------------------------------------------------------------
commit b9ab71f360b9a3e1e575a86bfca81fd4dacb4b50
Author: Raymond Toy <toy...@gm...>
Date:   Wed Aug 14 20:04:30 2013 -0700

    Improve accuracy of inverse_jacobi_sn.
    
    Instead of using F(asin(x),m), go directly to Carlson's Rf function.
    This removes the roundoff from computing sin(asin(x)), and related
    identities.

diff --git a/src/ellipt.lisp b/src/ellipt.lisp
index b20a850..d765868 100644
--- a/src/ellipt.lisp
+++ b/src/ellipt.lisp
@@ -1006,18 +1006,39 @@
   (twoargcheck form)
   (let ((u (simpcheck (cadr form) z))
 	(m (simpcheck (caddr form) z)))
+    ;; To numerically evaluate inverse_jacobi_sn (asn), use
+    ;;
+    ;; asn(x,m) = F(asin(x),m)
+    ;;
+    ;; But F(phi,m) = sin(phi)*rf(cos(phi)^2, 1-m*sin(phi)^2,1).  Thus
+    ;;
+    ;; asn(x,m) = F(asin(x),m)
+    ;;          = x*rf(1-x^2,1-m*x^2,1)
+    ;;
+    ;; I (rtoy) am not 100% about the first identity above for all
+    ;; complex values of x and m, but tests seem to indicate that it
+    ;; produces the correct value as verified by verifying
+    ;; jacobi_sn(inverse_jacobi_sn(x,m),m) = x.
     (cond ((float-numerical-eval-p u m)
-	   ;; Numerically evaluate asn
-	   ;;
-	   ;; asn(x,m) = F(asin(x),m)
-	   (complexify (elliptic-f (cl:asin ($float u)) ($float m))))
+	   (complexify (* u (bigfloat::bf-rf (bigfloat:to (float (- 1 (* u u))))
+					     (bigfloat:to (float (- 1 (* m u u))))
+					     1))))
 	  ((complex-float-numerical-eval-p u m)
-	   (complexify (elliptic-f (cl:asin (complex ($realpart u) ($imagpart u)))
-				   (complex ($realpart m) ($imagpart m)))))
+	   (let ((uu (complex ($float ($realpart u))
+			      ($float ($imagpart u))))
+		 (mm (complex ($float ($realpart m))
+			      ($float ($imagpart m)))))
+	     (complexify (* uu (bigfloat::bf-rf (- 1 (* uu uu))
+						(- 1 (* mm uu uu))
+						1)))))
 	  ((or (bigfloat-numerical-eval-p u m)
 	       (complex-bigfloat-numerical-eval-p u m))
-	   (to (bigfloat::bf-elliptic-f (bigfloat:asin (bigfloat:to u))
-					(bigfloat:to m))))
+	   (let ((uu (bigfloat:to u))
+		 (mm (bigfloat:to m)))
+	     (to (bigfloat:* uu
+			     (bigfloat::bf-rf (bigfloat:- 1 (bigfloat:* uu uu))
+					      (bigfloat:- 1 (bigfloat:* mm uu uu))
+					      1)))))
 	  ((zerop1 u)
 	   ;; asn(0,m) = 0
 	   0)

-----------------------------------------------------------------------

Summary of changes:
 src/ellipt.lisp |   37 +++++++++++++++++++++++++++++--------
 1 files changed, 29 insertions(+), 8 deletions(-)


hooks/post-receive
-- 
Maxima CAS