## [Maxima-commits] CVS: maxima/src ellipt.lisp,1.34,1.35

 [Maxima-commits] CVS: maxima/src ellipt.lisp,1.34,1.35 From: Raymond Toy - 2009-01-28 04:24:41 ```Update of /cvsroot/maxima/maxima/src In directory 23jxhf1.ch3.sourceforge.com:/tmp/cvs-serv9882 Modified Files: ellipt.lisp Log Message: First cut at adding bigfloat evaluation for elliptic integrals. The key algorithms are translations of Jim FitzSimons' bigfloat implemention of elliptic integrals. Preliminary support added for inverse_jacobi_sn, inverse_jacobi_cn, inverse_jacobi_dn and elliptic_f functions. Not well tested yet. Index: ellipt.lisp =================================================================== RCS file: /cvsroot/maxima/maxima/src/ellipt.lisp,v retrieving revision 1.34 retrieving revision 1.35 diff -u -d -r1.34 -r1.35 --- ellipt.lisp 26 Jan 2009 00:53:39 -0000 1.34 +++ ellipt.lisp 28 Jan 2009 04:24:33 -0000 1.35 @@ -891,6 +891,10 @@ (complex-number-p m)) (complexify (elliptic-f (cl:asin (complex (\$realpart u) (\$imagpart u))) (complex (\$realpart m) (\$imagpart m))))) + ((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)))) ((zerop1 u) ;; asn(0,m) = 0 0) @@ -925,6 +929,10 @@ (complex-number-p m)) (complexify (elliptic-f (cl:acos (complex (\$realpart u) (\$imagpart u))) (complex (\$realpart m) (\$imagpart m))))) + ((or (bigfloat-numerical-eval-p u m) + (complex-bigfloat-numerical-eval-p u m)) + (to (bigfloat::bf-elliptic-f (bigfloat:acos (bigfloat:to u)) + (bigfloat:to m)))) ((zerop1 m) ;; asn(x,0) = F(acos(x),0) = acos(x) `((%elliptic_f) ((%acos) ,u) 0)) @@ -967,6 +975,14 @@ (sqrt m)))) (complexify (elliptic-f (cl:asin phi) (complex (\$realpart m) (\$imagpart m)))))) + ((or (bigfloat-numerical-eval-p u m) + (complex-bigfloat-numerical-eval-p u m)) + (let* ((u (bigfloat:to (\$bfloat u))) + (phi (bigfloat:/ (bigfloat:* (bigfloat:sqrt (bigfloat:- 1 u)) + (bigfloat:sqrt (bigfloat:+ 1 u))) + (bigfloat:sqrt (bigfloat:to (\$bfloat m)))))) + (to (bigfloat::bf-elliptic-f (bigfloat:asin phi) + (bigfloat:to m))))) ((onep1 m) ;; x = dn(u,1) = sech(u). so u = asech(x) `((%asech) ,u)) @@ -1656,6 +1672,10 @@ (complex-number-p m)) (complexify (elliptic-f (complex (\$realpart phi) (\$imagpart phi)) (complex (\$realpart m) (\$imagpart m))))) + ((or (bigfloat-numerical-eval-p phi m) + (complex-bigfloat-numerical-eval-p phi m)) + (to (bigfloat::bf-elliptic-f (bigfloat:to (\$bfloat phi)) + (bigfloat:to (\$bfloat m))))) ((zerop1 phi) 0) ((zerop1 m) @@ -2155,6 +2175,286 @@ (values (+ drj-1 drj-2 (* (- a b) drj-3) (/ 3 (sqrt a))) drj-4))) +(in-package #-gcl #:bigfloat #+gcl "BIGFLOAT") +;; Translation of Jim FitzSimons' bigfloat implementation of elliptic +;; integrals from http://www.getnet.com/~cherry/elliptbf3.mac. +;; +;; The algorithms are based on B.C. Carlson's "Numerical Computation +;; of Real or Complex Elliptic Integrals". These are updated to the +;; algorithms in Journal of Computational and Applied Mathematics 118 +;; (2000) 71-85 "Reduction Theorems for Elliptic Integrands with the +;; Square Root of two quadritic factors" +;; + + +(defvar *bferrtol* (bigfloat 1d-6)) + +;; rc(x,y) = integrate(1/2*(t+x)^(-1/2)/(t+y), t, 0, inf) +;; +;; log(x) = (x-1)*rc(((1+x)/2)^2, x), x > 0 +;; asin(x) = x * rc(1-x^2, 1), |x|<= 1 +;; acos(x) = sqrt(1-x^2)*rc(x^2,1), 0 <= x <=1 +;; atan(x) = x * rc(1,1+x^2) +;; asinh(x) = x * rc(1+x^2,1) +;; acosh(x) = sqrt(x^2-1) * rc(x^2,1), x >= 1 +;; atanh(x) = x * rc(1,1-x^2), |x|<=1 + +(defun bf-rc (x y) + (let ((yn y) + xn z w a an pwr4 n epslon lambda sn s) + (cond ((and (zerop (imagpart yn)) + (minusp (realpart yn))) + (setf xn (- x y)) + (setf yn (- yn)) + (setf z yn) + (setf w (sqrt (/ x xn)))) + (t + (setf xn x) + (setf z yn) + (setf w (bigfloat 1)))) + (setf a (/ (+ xn yn yn) 3)) + (setf epslon (/ (abs (- a xn)) *bferrtol*)) + (setf an a) + (setf pwr4 (bigfloat 1)) + (setf n 0) + (loop while (> (* epslon pwr4) (abs an)) + do + (setf pwr4 (/ pwr4 4)) + (setf lambda (+ (* 2 (sqrt xn) (sqrt yn)) yn)) + (setf an (/ (+ an lambda) 4)) + (setf xn (/ (+ xn lambda) 4)) + (setf yn (/ (+ yn lambda) 4)) + (incf n)) + ;; c2=3/10,c3=1/7,c4=3/8,c5=9/22,c6=159/208,c7=9/8 + (setf sn (/ (* pwr4 (- z a)) an)) + (setf s (* sn sn (+ (bigfloat 3/10) + (* sn (+ (bigfloat 1/7) + (* sn (+ (bigfloat 3/8) + (* sn (+ (bigfloat 9/22) + (* sn (+ (bigfloat 159/208) + (* sn (bigfloat 9/8))))))))))))) + (/ (* w (+ 1 s)) + (sqrt an)))) + + + +;; rd(x,y,z) = integrate(3/2/sqrt(t+x)/sqrt(t+y)/sqrt(t+z), t, 0, inf) +;; +;; E(K) = rf(0, 1-K^2, 1) - (K^2/3)*rd(0,1-K^2,1) +;; +;; B = integrate(s^2/sqrt(1-s^4), s, 0 ,1) +;; = beta(3/4,1/2)/4 +;; = sqrt(%pi)*gamma(3/4)/gamma(1/4) +;; = 1/3*rd(0,2,1) +(defun bf-rd (x y z) + (let* ((xn x) + (yn y) + (zn z) + (a (/ (+ xn yn (* 3 zn)) 5)) + (epslon (/ (max (abs (- a xn)) + (abs (- a yn)) + (abs (- a zn))) + *bferrtol*)) + (an a) + (sigma 0) + (power4 1) + (n 0) + xnroot ynroot znroot lam) + (loop while (> (* power4 epslon) (abs an)) + do + (setf xnroot (sqrt xn)) + (setf ynroot (sqrt yn)) + (setf znroot (sqrt zn)) + (setf lam (+ (* xnroot ynroot) + (* xnroot znroot) + (* ynroot znroot))) + (setf sigma (+ sigma (/ power4 + (* znroot (+ zn lam))))) + (setf power4 (* power4 1/4)) + (setf xn (* (+ xn lam) 1/4)) + (setf yn (* (+ yn lam) 1/4)) + (setf zn (* (+ zn lam) 1/4)) + (setf an (* (+ an lam) 1/4)) + (incf n)) + ;; c1=-3/14,c2=1/6,c3=9/88,c4=9/22,c5=-3/22,c6=-9/52,c7=3/26 + (let* ((xndev (/ (* (- a x) power4) an)) + (yndev (/ (* (- a y) power4) an)) + (zndev (- (* (+ xndev yndev) 1/3))) + (ee2 (- (* xndev yndev) (* 6 zndev zndev))) + (ee3 (* (- (* 3 xndev yndev) + (* 8 zndev zndev)) + zndev)) + (ee4 (* 3 (- (* xndev yndev) (* zndev zndev)) zndev zndev)) + (ee5 (* xndev yndev zndev zndev zndev)) + (s (+ 1 + (* -3/14 ee2) + (* 1/6 ee3) + (* 9/88 ee2 ee2) + (* -3/22 ee4) + (* -9/52 ee2 ee3) + (* 3/26 ee5) + (* -1/16 ee2 ee2 ee2) + (* 3/10 ee3 ee3) + (* 3/20 ee2 ee4) + (* 45/272 ee2 ee2 ee3) + (* -9/68 (+ (* ee2 ee5) (* ee3 ee4)))))) + (+ (* 3 sigma) + (/ (* power4 s) + (expt an 3/2)))))) + +(defun bf-rf (x y z) + (let* ((xn x) + (yn y) + (zn z) + (a (/ (+ xn yn zn) 3)) + (epslon (/ (max (abs (- a xn)) + (abs (- a yn)) + (abs (- a zn))) + *bferrtol*)) + (an a) + (power4 1) + (n 0) + xnroot ynroot znroot lam) + (loop while (> (* power4 epslon) (abs an)) + do + (setf xnroot (sqrt xn)) + (setf ynroot (sqrt yn)) + (setf znroot (sqrt zn)) + (setf lam (+ (* xnroot ynroot) + (* xnroot znroot) + (* ynroot znroot))) + (setf power4 (* power4 1/4)) + (setf xn (* (+ xn lam) 1/4)) + (setf yn (* (+ yn lam) 1/4)) + (setf zn (* (+ zn lam) 1/4)) + (setf an (* (+ an lam) 1/4)) + (incf n)) + ;; c1=-3/14,c2=1/6,c3=9/88,c4=9/22,c5=-3/22,c6=-9/52,c7=3/26 + (let* ((xndev (/ (* (- a x) power4) an)) + (yndev (/ (* (- a y) power4) an)) + (zndev (- (+ xndev yndev))) + (ee2 (- (* xndev yndev) (* 6 zndev zndev))) + (ee3 (* xndev yndev zndev)) + (s (+ 1 + (* -1/10 ee2) + (* 1/14 ee3) + (* 1/24 ee2 ee2) + (* -3/44 ee2 ee3)))) + (/ s (sqrt an))))) + +(defun bf-rj1 (x y z p) + (let* ((xn x) + (yn y) + (zn z) + (pn p) + (en (* (- pn xn) + (- pn yn) + (- pn zn))) + (sigma (bigfloat 0)) + (power4 (bigfloat 1)) + (k 0) + (a (/ (+ xn yn zn pn pn) 5)) + (epslon (/ (max (abs (- a xn)) + (abs (- a yn)) + (abs (- a zn)) + (abs (- a pn))) + *bferrtol*)) + (an a) + xnroot ynroot znroot pnroot lam dn) + (loop while (> (* power4 epslon) (abs an)) + do + (setf xnroot (sqrt xn)) + (setf ynroot (sqrt yn)) + (setf znroot (sqrt zn)) + (setf pnroot (sqrt pn)) + (setf lam (+ (* xnroot ynroot) + (* xnroot znroot) + (* ynroot znroot))) + (setf dn (* (+ pnroot xnroot) + (+ pnroot ynroot) + (+ pnroot znroot))) + (setf sigma (+ sigma + (/ (* power4 + (bf-rc 1 (+ 1 (/ en (* dn dn))))) + dn))) + (setf power4 (* power4 1/4)) + (setf en (/ en 64)) + (setf xn (* (+ xn lam) 1/4)) + (setf yn (* (+ yn lam) 1/4)) + (setf zn (* (+ zn lam) 1/4)) + (setf pn (* (+ pn lam) 1/4)) + (setf an (* (+ an lam) 1/4)) + (incf k)) + (let* ((xndev (/ (* (- a x) power4) an)) + (yndev (/ (* (- a y) power4) an)) + (zndev (/ (* (- a z) power4) an)) + (pndev (* -0.5 (+ xndev yndev zndev))) + (ee2 (+ (* xndev yndev) + (* xndev zndev) + (* yndev zndev) + (* -3 pndev pndev))) + (ee3 (+ (* xndev yndev zndev) + (* 2 ee2 pndev) + (* 4 pndev pndev pndev))) + (ee4 (* (+ (* 2 xndev yndev zndev) + (* ee2 pndev) + (* 3 pndev pndev pndev)) + pndev)) + (ee5 (* xndev yndev zndev pndev pndev)) + (s (+ 1 + (* -3/14 ee2) + (* 1/6 ee3) + (* 9/88 ee2 ee2) + (* -3/22 ee4) + (* -9/52 ee2 ee3) + (* 3/26 ee5) + (* -1/16 ee2 ee2 ee2) + (* 3/10 ee3 ee3) + (* 3/20 ee2 ee4) + (* 45/272 ee2 ee2 ee3) + (* -9/68 (+ (* ee2 ee5) (* ee3 ee4)))))) + (+ (* 6 sigma) + (/ (* power4 s) + (sqrt (* an an an))))))) + +(defun bf-rj (x y z p) + (let* ((xn x) + (yn y) + (zn z) + (qn (- p))) + (cond ((and (and (zerop (imagpart xn)) (>= (realpart xn) 0)) + (and (zerop (imagpart yn)) (>= (realpart yn) 0)) + (and (zerop (imagpart zn)) (>= (realpart zn) 0)) + (and (zerop (imagpart qn)) (> (realpart qn) 0))) + (destructuring-bind (xn yn zn) + (sort (list xn yn zn) #'<) + (let* ((pn (+ yn (* (- zn yn) (/ (- yn xn) (+ yn qn))))) + (s (- (* (- pn yn) (bf-rj1 xn yn zn pn)) + (* 3 (bf-rf xn yn zn))))) + (setf s (+ s (* 3 (sqrt (/ (* xn yn zn) + (+ (* xn zn) (* pn qn)))) + (bf-rc (+ (* xn zn) (* pn qn)) (* pn qn))))) + (/ s (+ yn qn))))) + (t + (bf-rj1 x y z p))))) + +(defun bf-rg (x y z) + (* 0.5 + (+ (* z (bf-rf x y z)) + (* (- z x) + (- y z) + (bf-rd x y z) + 1/3) + (sqrt (/ (* x y) z))))) + +;; elliptic_f(phi,m) = sin(phi)*rf(cos(phi)^2, 1-m*sin(phi)^2,1) +(defun bf-elliptic-f (phi m) + (let ((s (sin phi)) + (c (cos phi))) + (* s (bf-rf (* c c) (- 1 (* m s s)) 1)))) + +(in-package :maxima) + ;;; Other Jacobian elliptic functions ;; jacobi_ns(u,m) = 1/jacobi_sn(u,m) @@ -3236,6 +3536,10 @@ (complex-number-p m)) (complexify (elliptic-f (cl:asin (/ (complex (\$realpart u) (\$imagpart u)))) (complex (\$realpart m) (\$imagpart m))))) + ((or (bigfloat-numerical-eval-p u m) + (complex-bigfloat-numerical-eval-p u m)) + (to (bigfloat::bf-elliptic-f (bigfloat:asin (bigfloat:/ (bigfloat:to (\$bfloat u)))) + (bigfloat:to (\$bfloat m))))) ((zerop1 m) ;; ans(x,0) = F(asin(1/x),0) = asin(1/x) `((%elliptic_f) ((%asin) ((mexpt) ,u -1)) 0)) ```

 [Maxima-commits] CVS: maxima/src ellipt.lisp,1.34,1.35 From: Raymond Toy - 2009-01-28 04:24:41 ```Update of /cvsroot/maxima/maxima/src In directory 23jxhf1.ch3.sourceforge.com:/tmp/cvs-serv9882 Modified Files: ellipt.lisp Log Message: First cut at adding bigfloat evaluation for elliptic integrals. The key algorithms are translations of Jim FitzSimons' bigfloat implemention of elliptic integrals. Preliminary support added for inverse_jacobi_sn, inverse_jacobi_cn, inverse_jacobi_dn and elliptic_f functions. Not well tested yet. Index: ellipt.lisp =================================================================== RCS file: /cvsroot/maxima/maxima/src/ellipt.lisp,v retrieving revision 1.34 retrieving revision 1.35 diff -u -d -r1.34 -r1.35 --- ellipt.lisp 26 Jan 2009 00:53:39 -0000 1.34 +++ ellipt.lisp 28 Jan 2009 04:24:33 -0000 1.35 @@ -891,6 +891,10 @@ (complex-number-p m)) (complexify (elliptic-f (cl:asin (complex (\$realpart u) (\$imagpart u))) (complex (\$realpart m) (\$imagpart m))))) + ((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)))) ((zerop1 u) ;; asn(0,m) = 0 0) @@ -925,6 +929,10 @@ (complex-number-p m)) (complexify (elliptic-f (cl:acos (complex (\$realpart u) (\$imagpart u))) (complex (\$realpart m) (\$imagpart m))))) + ((or (bigfloat-numerical-eval-p u m) + (complex-bigfloat-numerical-eval-p u m)) + (to (bigfloat::bf-elliptic-f (bigfloat:acos (bigfloat:to u)) + (bigfloat:to m)))) ((zerop1 m) ;; asn(x,0) = F(acos(x),0) = acos(x) `((%elliptic_f) ((%acos) ,u) 0)) @@ -967,6 +975,14 @@ (sqrt m)))) (complexify (elliptic-f (cl:asin phi) (complex (\$realpart m) (\$imagpart m)))))) + ((or (bigfloat-numerical-eval-p u m) + (complex-bigfloat-numerical-eval-p u m)) + (let* ((u (bigfloat:to (\$bfloat u))) + (phi (bigfloat:/ (bigfloat:* (bigfloat:sqrt (bigfloat:- 1 u)) + (bigfloat:sqrt (bigfloat:+ 1 u))) + (bigfloat:sqrt (bigfloat:to (\$bfloat m)))))) + (to (bigfloat::bf-elliptic-f (bigfloat:asin phi) + (bigfloat:to m))))) ((onep1 m) ;; x = dn(u,1) = sech(u). so u = asech(x) `((%asech) ,u)) @@ -1656,6 +1672,10 @@ (complex-number-p m)) (complexify (elliptic-f (complex (\$realpart phi) (\$imagpart phi)) (complex (\$realpart m) (\$imagpart m))))) + ((or (bigfloat-numerical-eval-p phi m) + (complex-bigfloat-numerical-eval-p phi m)) + (to (bigfloat::bf-elliptic-f (bigfloat:to (\$bfloat phi)) + (bigfloat:to (\$bfloat m))))) ((zerop1 phi) 0) ((zerop1 m) @@ -2155,6 +2175,286 @@ (values (+ drj-1 drj-2 (* (- a b) drj-3) (/ 3 (sqrt a))) drj-4))) +(in-package #-gcl #:bigfloat #+gcl "BIGFLOAT") +;; Translation of Jim FitzSimons' bigfloat implementation of elliptic +;; integrals from http://www.getnet.com/~cherry/elliptbf3.mac. +;; +;; The algorithms are based on B.C. Carlson's "Numerical Computation +;; of Real or Complex Elliptic Integrals". These are updated to the +;; algorithms in Journal of Computational and Applied Mathematics 118 +;; (2000) 71-85 "Reduction Theorems for Elliptic Integrands with the +;; Square Root of two quadritic factors" +;; + + +(defvar *bferrtol* (bigfloat 1d-6)) + +;; rc(x,y) = integrate(1/2*(t+x)^(-1/2)/(t+y), t, 0, inf) +;; +;; log(x) = (x-1)*rc(((1+x)/2)^2, x), x > 0 +;; asin(x) = x * rc(1-x^2, 1), |x|<= 1 +;; acos(x) = sqrt(1-x^2)*rc(x^2,1), 0 <= x <=1 +;; atan(x) = x * rc(1,1+x^2) +;; asinh(x) = x * rc(1+x^2,1) +;; acosh(x) = sqrt(x^2-1) * rc(x^2,1), x >= 1 +;; atanh(x) = x * rc(1,1-x^2), |x|<=1 + +(defun bf-rc (x y) + (let ((yn y) + xn z w a an pwr4 n epslon lambda sn s) + (cond ((and (zerop (imagpart yn)) + (minusp (realpart yn))) + (setf xn (- x y)) + (setf yn (- yn)) + (setf z yn) + (setf w (sqrt (/ x xn)))) + (t + (setf xn x) + (setf z yn) + (setf w (bigfloat 1)))) + (setf a (/ (+ xn yn yn) 3)) + (setf epslon (/ (abs (- a xn)) *bferrtol*)) + (setf an a) + (setf pwr4 (bigfloat 1)) + (setf n 0) + (loop while (> (* epslon pwr4) (abs an)) + do + (setf pwr4 (/ pwr4 4)) + (setf lambda (+ (* 2 (sqrt xn) (sqrt yn)) yn)) + (setf an (/ (+ an lambda) 4)) + (setf xn (/ (+ xn lambda) 4)) + (setf yn (/ (+ yn lambda) 4)) + (incf n)) + ;; c2=3/10,c3=1/7,c4=3/8,c5=9/22,c6=159/208,c7=9/8 + (setf sn (/ (* pwr4 (- z a)) an)) + (setf s (* sn sn (+ (bigfloat 3/10) + (* sn (+ (bigfloat 1/7) + (* sn (+ (bigfloat 3/8) + (* sn (+ (bigfloat 9/22) + (* sn (+ (bigfloat 159/208) + (* sn (bigfloat 9/8))))))))))))) + (/ (* w (+ 1 s)) + (sqrt an)))) + + + +;; rd(x,y,z) = integrate(3/2/sqrt(t+x)/sqrt(t+y)/sqrt(t+z), t, 0, inf) +;; +;; E(K) = rf(0, 1-K^2, 1) - (K^2/3)*rd(0,1-K^2,1) +;; +;; B = integrate(s^2/sqrt(1-s^4), s, 0 ,1) +;; = beta(3/4,1/2)/4 +;; = sqrt(%pi)*gamma(3/4)/gamma(1/4) +;; = 1/3*rd(0,2,1) +(defun bf-rd (x y z) + (let* ((xn x) + (yn y) + (zn z) + (a (/ (+ xn yn (* 3 zn)) 5)) + (epslon (/ (max (abs (- a xn)) + (abs (- a yn)) + (abs (- a zn))) + *bferrtol*)) + (an a) + (sigma 0) + (power4 1) + (n 0) + xnroot ynroot znroot lam) + (loop while (> (* power4 epslon) (abs an)) + do + (setf xnroot (sqrt xn)) + (setf ynroot (sqrt yn)) + (setf znroot (sqrt zn)) + (setf lam (+ (* xnroot ynroot) + (* xnroot znroot) + (* ynroot znroot))) + (setf sigma (+ sigma (/ power4 + (* znroot (+ zn lam))))) + (setf power4 (* power4 1/4)) + (setf xn (* (+ xn lam) 1/4)) + (setf yn (* (+ yn lam) 1/4)) + (setf zn (* (+ zn lam) 1/4)) + (setf an (* (+ an lam) 1/4)) + (incf n)) + ;; c1=-3/14,c2=1/6,c3=9/88,c4=9/22,c5=-3/22,c6=-9/52,c7=3/26 + (let* ((xndev (/ (* (- a x) power4) an)) + (yndev (/ (* (- a y) power4) an)) + (zndev (- (* (+ xndev yndev) 1/3))) + (ee2 (- (* xndev yndev) (* 6 zndev zndev))) + (ee3 (* (- (* 3 xndev yndev) + (* 8 zndev zndev)) + zndev)) + (ee4 (* 3 (- (* xndev yndev) (* zndev zndev)) zndev zndev)) + (ee5 (* xndev yndev zndev zndev zndev)) + (s (+ 1 + (* -3/14 ee2) + (* 1/6 ee3) + (* 9/88 ee2 ee2) + (* -3/22 ee4) + (* -9/52 ee2 ee3) + (* 3/26 ee5) + (* -1/16 ee2 ee2 ee2) + (* 3/10 ee3 ee3) + (* 3/20 ee2 ee4) + (* 45/272 ee2 ee2 ee3) + (* -9/68 (+ (* ee2 ee5) (* ee3 ee4)))))) + (+ (* 3 sigma) + (/ (* power4 s) + (expt an 3/2)))))) + +(defun bf-rf (x y z) + (let* ((xn x) + (yn y) + (zn z) + (a (/ (+ xn yn zn) 3)) + (epslon (/ (max (abs (- a xn)) + (abs (- a yn)) + (abs (- a zn))) + *bferrtol*)) + (an a) + (power4 1) + (n 0) + xnroot ynroot znroot lam) + (loop while (> (* power4 epslon) (abs an)) + do + (setf xnroot (sqrt xn)) + (setf ynroot (sqrt yn)) + (setf znroot (sqrt zn)) + (setf lam (+ (* xnroot ynroot) + (* xnroot znroot) + (* ynroot znroot))) + (setf power4 (* power4 1/4)) + (setf xn (* (+ xn lam) 1/4)) + (setf yn (* (+ yn lam) 1/4)) + (setf zn (* (+ zn lam) 1/4)) + (setf an (* (+ an lam) 1/4)) + (incf n)) + ;; c1=-3/14,c2=1/6,c3=9/88,c4=9/22,c5=-3/22,c6=-9/52,c7=3/26 + (let* ((xndev (/ (* (- a x) power4) an)) + (yndev (/ (* (- a y) power4) an)) + (zndev (- (+ xndev yndev))) + (ee2 (- (* xndev yndev) (* 6 zndev zndev))) + (ee3 (* xndev yndev zndev)) + (s (+ 1 + (* -1/10 ee2) + (* 1/14 ee3) + (* 1/24 ee2 ee2) + (* -3/44 ee2 ee3)))) + (/ s (sqrt an))))) + +(defun bf-rj1 (x y z p) + (let* ((xn x) + (yn y) + (zn z) + (pn p) + (en (* (- pn xn) + (- pn yn) + (- pn zn))) + (sigma (bigfloat 0)) + (power4 (bigfloat 1)) + (k 0) + (a (/ (+ xn yn zn pn pn) 5)) + (epslon (/ (max (abs (- a xn)) + (abs (- a yn)) + (abs (- a zn)) + (abs (- a pn))) + *bferrtol*)) + (an a) + xnroot ynroot znroot pnroot lam dn) + (loop while (> (* power4 epslon) (abs an)) + do + (setf xnroot (sqrt xn)) + (setf ynroot (sqrt yn)) + (setf znroot (sqrt zn)) + (setf pnroot (sqrt pn)) + (setf lam (+ (* xnroot ynroot) + (* xnroot znroot) + (* ynroot znroot))) + (setf dn (* (+ pnroot xnroot) + (+ pnroot ynroot) + (+ pnroot znroot))) + (setf sigma (+ sigma + (/ (* power4 + (bf-rc 1 (+ 1 (/ en (* dn dn))))) + dn))) + (setf power4 (* power4 1/4)) + (setf en (/ en 64)) + (setf xn (* (+ xn lam) 1/4)) + (setf yn (* (+ yn lam) 1/4)) + (setf zn (* (+ zn lam) 1/4)) + (setf pn (* (+ pn lam) 1/4)) + (setf an (* (+ an lam) 1/4)) + (incf k)) + (let* ((xndev (/ (* (- a x) power4) an)) + (yndev (/ (* (- a y) power4) an)) + (zndev (/ (* (- a z) power4) an)) + (pndev (* -0.5 (+ xndev yndev zndev))) + (ee2 (+ (* xndev yndev) + (* xndev zndev) + (* yndev zndev) + (* -3 pndev pndev))) + (ee3 (+ (* xndev yndev zndev) + (* 2 ee2 pndev) + (* 4 pndev pndev pndev))) + (ee4 (* (+ (* 2 xndev yndev zndev) + (* ee2 pndev) + (* 3 pndev pndev pndev)) + pndev)) + (ee5 (* xndev yndev zndev pndev pndev)) + (s (+ 1 + (* -3/14 ee2) + (* 1/6 ee3) + (* 9/88 ee2 ee2) + (* -3/22 ee4) + (* -9/52 ee2 ee3) + (* 3/26 ee5) + (* -1/16 ee2 ee2 ee2) + (* 3/10 ee3 ee3) + (* 3/20 ee2 ee4) + (* 45/272 ee2 ee2 ee3) + (* -9/68 (+ (* ee2 ee5) (* ee3 ee4)))))) + (+ (* 6 sigma) + (/ (* power4 s) + (sqrt (* an an an))))))) + +(defun bf-rj (x y z p) + (let* ((xn x) + (yn y) + (zn z) + (qn (- p))) + (cond ((and (and (zerop (imagpart xn)) (>= (realpart xn) 0)) + (and (zerop (imagpart yn)) (>= (realpart yn) 0)) + (and (zerop (imagpart zn)) (>= (realpart zn) 0)) + (and (zerop (imagpart qn)) (> (realpart qn) 0))) + (destructuring-bind (xn yn zn) + (sort (list xn yn zn) #'<) + (let* ((pn (+ yn (* (- zn yn) (/ (- yn xn) (+ yn qn))))) + (s (- (* (- pn yn) (bf-rj1 xn yn zn pn)) + (* 3 (bf-rf xn yn zn))))) + (setf s (+ s (* 3 (sqrt (/ (* xn yn zn) + (+ (* xn zn) (* pn qn)))) + (bf-rc (+ (* xn zn) (* pn qn)) (* pn qn))))) + (/ s (+ yn qn))))) + (t + (bf-rj1 x y z p))))) + +(defun bf-rg (x y z) + (* 0.5 + (+ (* z (bf-rf x y z)) + (* (- z x) + (- y z) + (bf-rd x y z) + 1/3) + (sqrt (/ (* x y) z))))) + +;; elliptic_f(phi,m) = sin(phi)*rf(cos(phi)^2, 1-m*sin(phi)^2,1) +(defun bf-elliptic-f (phi m) + (let ((s (sin phi)) + (c (cos phi))) + (* s (bf-rf (* c c) (- 1 (* m s s)) 1)))) + +(in-package :maxima) + ;;; Other Jacobian elliptic functions ;; jacobi_ns(u,m) = 1/jacobi_sn(u,m) @@ -3236,6 +3536,10 @@ (complex-number-p m)) (complexify (elliptic-f (cl:asin (/ (complex (\$realpart u) (\$imagpart u)))) (complex (\$realpart m) (\$imagpart m))))) + ((or (bigfloat-numerical-eval-p u m) + (complex-bigfloat-numerical-eval-p u m)) + (to (bigfloat::bf-elliptic-f (bigfloat:asin (bigfloat:/ (bigfloat:to (\$bfloat u)))) + (bigfloat:to (\$bfloat m))))) ((zerop1 m) ;; ans(x,0) = F(asin(1/x),0) = asin(1/x) `((%elliptic_f) ((%asin) ((mexpt) ,u -1)) 0)) ```

No, thanks