Commit 535d2509c2974c113e6c1bfa27ece5d5b8f99bd1

Authored by Raymond Toy
1 parent 258cd829

Add support for Jacobi elliptic functions

qd-elliptic.lisp:
o Implementation for Jacobi sn, cn, and dn functions.

qd-methods.lisp:
o Add EPSILON method to return the floating-point epsilon value for
  the given float.

oct.asd:
o Build qd-elliptic.lisp.
... ... @@ -58,6 +58,8 @@
58 58 :depends-on ("qd-methods"))
59 59 (:file "qd-complex"
60 60 :depends-on ("qd-methods"))
  61 + (:file "qd-elliptic"
  62 + :depends-on ("qd-methods"))
61 63 ))
62 64  
63 65 (defmethod perform ((op test-op) (c (eql (find-system :oct))))
... ...
qd-elliptic.lisp 0 → 100644
  1 +;;;; -*- Mode: lisp -*-
  2 +;;;;
  3 +;;;; Copyright (c) 2011 Raymond Toy
  4 +;;;; Permission is hereby granted, free of charge, to any person
  5 +;;;; obtaining a copy of this software and associated documentation
  6 +;;;; files (the "Software"), to deal in the Software without
  7 +;;;; restriction, including without limitation the rights to use,
  8 +;;;; copy, modify, merge, publish, distribute, sublicense, and/or sell
  9 +;;;; copies of the Software, and to permit persons to whom the
  10 +;;;; Software is furnished to do so, subject to the following
  11 +;;;; conditions:
  12 +;;;;
  13 +;;;; The above copyright notice and this permission notice shall be
  14 +;;;; included in all copies or substantial portions of the Software.
  15 +;;;;
  16 +;;;; THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
  17 +;;;; EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
  18 +;;;; OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
  19 +;;;; NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
  20 +;;;; HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
  21 +;;;; WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
  22 +;;;; FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
  23 +;;;; OTHER DEALINGS IN THE SOFTWARE.
  24 +
  25 +(in-package #:oct)
  26 +
  27 +(declaim (inline descending-transform ascending-transform))
  28 +
  29 +(defun ascending-transform (u m)
  30 + ;; A&S 16.14.1
  31 + ;;
  32 + ;; Take care in computing this transform. For the case where
  33 + ;; m is complex, we should compute sqrt(mu1) first as
  34 + ;; (1-sqrt(m))/(1+sqrt(m)), and then square this to get mu1.
  35 + ;; If not, we may choose the wrong branch when computing
  36 + ;; sqrt(mu1).
  37 + (let* ((root-m (sqrt m))
  38 + (mu (/ (* 4 root-m)
  39 + (expt (1+ root-m) 2)))
  40 + (root-mu1 (/ (- 1 root-m) (+ 1 root-m)))
  41 + (v (/ u (1+ root-mu1))))
  42 + (values v mu root-mu1)))
  43 +
  44 +(defun descending-transform (u m)
  45 + ;; Note: Don't calculate mu first, as given in 16.12.1. We
  46 + ;; should calculate sqrt(mu) = (1-sqrt(m1)/(1+sqrt(m1)), and
  47 + ;; then compute mu = sqrt(mu)^2. If we calculate mu first,
  48 + ;; sqrt(mu) loses information when m or m1 is complex.
  49 + (let* ((root-m1 (sqrt (- 1 m)))
  50 + (root-mu (/ (- 1 root-m1) (+ 1 root-m1)))
  51 + (mu (* root-mu root-mu))
  52 + (v (/ u (1+ root-mu))))
  53 + (values v mu root-mu)))
  54 +
  55 +
  56 +;; Could use the descending transform, but some of my tests show
  57 +;; that it has problems with roundoff errors.
  58 +
  59 +;; WARNING: This doesn't work very well for u > 1000 or so. For
  60 +;; example (elliptic-dn-ascending 1000b0 .5b0) -> 3.228b324, but dn <= 1.
  61 +#+nil
  62 +(defun elliptic-dn-ascending (u m)
  63 + (cond ((zerop m)
  64 + ;; A&S 16.6.3
  65 + 1.0)
  66 + ((< (abs (- 1 m)) (* 4 (epsilon u)))
  67 + ;; A&S 16.6.3
  68 + (/ (cosh u)))
  69 + (t
  70 + (multiple-value-bind (v mu root-mu1)
  71 + (ascending-transform u m)
  72 + ;; A&S 16.14.4
  73 + (let* ((new-dn (elliptic-dn-ascending v mu)))
  74 + (* (/ (- 1 root-mu1) mu)
  75 + (/ (+ root-mu1 (* new-dn new-dn))
  76 + new-dn)))))))
  77 +
  78 +;; Don't use the descending version because it requires cn, dn, and
  79 +;; sn.
  80 +;;
  81 +;; WARNING: This doesn't work very well for large u.
  82 +;; (elliptic-cn-ascending 1000b0 .5b0) -> 4.565b324. But |cn| <= 1.
  83 +#+nil
  84 +(defun elliptic-cn-ascending (u m)
  85 + (cond ((zerop m)
  86 + ;; A&S 16.6.2
  87 + (cos u))
  88 + ((< (abs (- 1 m)) (* 4 (epsilon u)))
  89 + ;; A&S 16.6.2
  90 + (/ (cl:cosh u)))
  91 + (t
  92 + (multiple-value-bind (v mu root-mu1)
  93 + (ascending-transform u m)
  94 + ;; A&S 16.14.3
  95 + (let* ((new-dn (elliptic-dn-ascending v mu)))
  96 + (* (/ (+ 1 root-mu1) mu)
  97 + (/ (- (* new-dn new-dn) root-mu1)
  98 + new-dn)))))))
  99 +
  100 +;;
  101 +;; This appears to work quite well for both real and complex values
  102 +;; of u.
  103 +(defun elliptic-sn-descending (u m)
  104 + (cond ((= m 1)
  105 + ;; A&S 16.6.1
  106 + (tanh u))
  107 + ((< (abs m) (epsilon u))
  108 + ;; A&S 16.6.1
  109 + (sin u))
  110 + (t
  111 + (multiple-value-bind (v mu root-mu)
  112 + (descending-transform u m)
  113 + (let* ((new-sn (elliptic-sn-descending v mu)))
  114 + (/ (* (1+ root-mu) new-sn)
  115 + (1+ (* root-mu new-sn new-sn))))))))
  116 +
  117 +;; We don't use the ascending transform here because it requires
  118 +;; evaluating sn, cn, and dn. The ascending transform only needs
  119 +;; sn.
  120 +#+nil
  121 +(defun elliptic-sn-ascending (u m)
  122 + (if (< (abs (- 1 m)) (* 4 flonum-epsilon))
  123 + ;; A&S 16.6.1
  124 + (tanh u)
  125 + (multiple-value-bind (v mu root-mu1)
  126 + (ascending-transform u m)
  127 + ;; A&S 16.14.2
  128 + (let* ((new-cn (elliptic-cn-ascending v mu))
  129 + (new-dn (elliptic-dn-ascending v mu))
  130 + (new-sn (elliptic-sn-ascending v mu)))
  131 + (/ (* (+ 1 root-mu1) new-sn new-cn)
  132 + new-dn)))))
  133 +
  134 +(defun jacobi-sn (u m)
  135 + (let ((s (elliptic-sn-descending u m)))
  136 + (if (and (realp u) (realp m))
  137 + (realpart s)
  138 + s)))
  139 +
  140 +(defun jacobi-dn (u m)
  141 + ;; Use the Gauss transformation from
  142 + ;; http://functions.wolfram.com/09.29.16.0013.01:
  143 + ;;
  144 + ;;
  145 + ;; dn((1+sqrt(m))*z, 4*sqrt(m)/(1+sqrt(m))^2)
  146 + ;; = (1-sqrt(m)*sn(z, m)^2)/(1+sqrt(m)*sn(z,m)^2)
  147 + ;;
  148 + ;; So
  149 + ;;
  150 + ;; dn(y, mu) = (1-sqrt(m)*sn(z, m)^2)/(1+sqrt(m)*sn(z,m)^2)
  151 + ;;
  152 + ;; where z = y/(1+sqrt(m)) and mu=4*sqrt(m)/(1+sqrt(m))^2.
  153 + ;;
  154 + ;; Solve for m, and we get
  155 + ;;
  156 + ;; sqrt(m) = -(mu+2*sqrt(1-mu)-2)/mu or (-mu+2*sqrt(1-mu)+2)/mu.
  157 + ;;
  158 + ;; I don't think it matters which sqrt we use, so I (rtoy)
  159 + ;; arbitrarily choose the first one above.
  160 + ;;
  161 + ;; Note that (1-sqrt(1-mu))/(1+sqrt(1-mu)) is the same as
  162 + ;; -(mu+2*sqrt(1-mu)-2)/mu. Also, the former is more
  163 + ;; accurate for small mu.
  164 + (let* ((root (let ((root-1-m (sqrt (- 1 m))))
  165 + (/ (- 1 root-1-m)
  166 + (+ 1 root-1-m))))
  167 + (z (/ u (+ 1 root)))
  168 + (s (elliptic-sn-descending z (* root root)))
  169 + (p (* root s s )))
  170 + (/ (- 1 p)
  171 + (+ 1 p))))
  172 +
  173 +(defun jacobi-cn (u m)
  174 + ;; Use the ascending Landen transformation, A&S 16.14.3.
  175 + (multiple-value-bind (v mu root-mu1)
  176 + (ascending-transform u m)
  177 + (let ((d (dn v mu)))
  178 + (* (/ (+ 1 root-mu1) mu)
  179 + (/ (- (* d d) root-mu1)
  180 + d)))))
0 181 \ No newline at end of file
... ...
qd-methods.lisp
... ... @@ -1089,3 +1089,23 @@ underlying floating-point format&quot;
1089 1089 ;; and make a real qd-real float, instead of the hackish
1090 1090 ;; %qd-real.
1091 1091 (set-dispatch-macro-character #\# #\Q #'qd-class-reader)
  1092 +
  1093 +
  1094 +(defmethod epsilon ((m cl:float))
  1095 + (etypecase m
  1096 + (single-float single-float-epsilon)
  1097 + (double-float double-float-epsilon)))
  1098 +
  1099 +(defmethod epsilon ((m cl:complex))
  1100 + (epsilon (realpart m)))
  1101 +
  1102 +(defmethod epsilon ((m qd-real))
  1103 + ;; What is the epsilon value for a quad-double? This is complicated
  1104 + ;; by the fact that things like (+ #q1 #q1q-100) is representable as
  1105 + ;; a quad-double. For most purposes we want epsilon to be close to
  1106 + ;; the 212 bits of precision (4*53 bits) that we normally have with
  1107 + ;; a quad-double.
  1108 + (scale-float #q1 -212))
  1109 +
  1110 +(defmethod epsilon ((m qd-complex))
  1111 + (epsilon (realpart m)))
... ...