From: Raymond T. <rt...@us...> - 2009-01-06 22:16:36
|
Update of /cvsroot/maxima/maxima/src In directory 23jxhf1.ch3.sourceforge.com:/tmp/cvs-serv15177 Modified Files: ellipt.lisp Log Message: o Fix a compiler warning about unused variable in JACOBI-AGM. o AGM appears to be rather slow and inaccurate in computing Jacobi sn (and cn and dn) when the parameter, m, is 1. But m = 1 simplifies to elementary functions, so use them instead of AGM. Index: ellipt.lisp =================================================================== RCS file: /cvsroot/maxima/maxima/src/ellipt.lisp,v retrieving revision 1.31 retrieving revision 1.32 diff -u -d -r1.31 -r1.32 --- ellipt.lisp 1 Jan 2009 00:52:13 -0000 1.31 +++ ellipt.lisp 6 Jan 2009 21:47:20 -0000 1.32 @@ -182,7 +182,7 @@ (let* ((agm-data (nreverse (rest (agm-scale 1d0 (sqrt (- 1 m)) (sqrt m))))) (phi (destructuring-bind (n a b c) (first agm-data) - (declare (ignore b)) + (declare (ignore b c)) (* a u (ash 1 n)))) (phi1 0d0)) (dolist (agm agm-data) @@ -194,28 +194,49 @@ (values (cl:sin phi) (cl:cos phi) (/ (cl:cos phi) (cl:cos (- phi1 phi)))))) (defun sn (u m) - (multiple-value-bind (s c d) - (jacobi-agm u m) - (declare (ignore c d)) - (if (and (realp u) (realp m)) - (realpart s) - s))) + (cond ((zerop m) + ;; jacobi_sn(u,0) = sin(u) + (sin u)) + ((= m 1) + ;; jacobi_sn(u,1) = tanh(u) + (cl:tanh u)) + (t + (multiple-value-bind (s c d) + (jacobi-agm u m) + (declare (ignore c d)) + (if (and (realp u) (realp m)) + (realpart s) + s))))) (defun cn (u m) - (multiple-value-bind (s c d) - (jacobi-agm u m) - (declare (ignore s d)) - (if (and (realp u) (realp m)) - (realpart c) - c))) + (cond ((zerop m) + ;; jacobi_cn(u,0) = cos(u) + (cos u)) + ((= m 1) + ;; jacobi_cn(u,1) = sech(u) + (/ (cl:cosh u))) + (t + (multiple-value-bind (s c d) + (jacobi-agm u m) + (declare (ignore s d)) + (if (and (realp u) (realp m)) + (realpart c) + c))))) (defun dn (u m) - (multiple-value-bind (s c d) - (jacobi-agm u m) - (declare (ignore s c)) - (if (and (realp u) (realp m)) - (realpart d) - d))) + (cond ((zerop m) + ;; jacobi_dn(u,0) = 1 + 1) + ((= m 1) + ;; jacobi_dn(u,1) = sech(u) + (/ (cl:cosh u))) + (t + (multiple-value-bind (s c d) + (jacobi-agm u m) + (declare (ignore s c)) + (if (and (realp u) (realp m)) + (realpart d) + d))))) ;; ;; How this works, I think. ;; |