packages feed

elliptic-integrals (empty) → 0.1.0.0

raw patch · 11 files changed

+587/−0 lines, 11 filesdep +basedep +elliptic-integralsdep +tastysetup-changed

Dependencies added: base, elliptic-integrals, tasty, tasty-hunit

Files

+ CHANGELOG.md view
@@ -0,0 +1,4 @@+0.1.0.0+-------+* initial release+
+ LICENSE view
@@ -0,0 +1,30 @@+Copyright Stéphane Laurent (c) 2023++All rights reserved.++Redistribution and use in source and binary forms, with or without+modification, are permitted provided that the following conditions are met:++    * Redistributions of source code must retain the above copyright+      notice, this list of conditions and the following disclaimer.++    * Redistributions in binary form must reproduce the above+      copyright notice, this list of conditions and the following+      disclaimer in the documentation and/or other materials provided+      with the distribution.++    * Neither the name of Stéphane Laurent nor the names of other+      contributors may be used to endorse or promote products derived+      from this software without specific prior written permission.++THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS+"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT+LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR+A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT+OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,+SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT+LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,+DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY+THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT+(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE+OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+ README.md view
@@ -0,0 +1,3 @@+# elliptic-integrals++Evaluation of the Carlson elliptic integrals and the incomplete elliptic integrals with complex arguments.
+ Setup.hs view
@@ -0,0 +1,2 @@+import Distribution.Simple+main = defaultMain
+ elliptic-integrals.cabal view
@@ -0,0 +1,40 @@+name:                elliptic-integrals+version:             0.1.0.0+synopsis:            Carlson Elliptic Integrals and Incomplete Elliptic Integrals+description:         Evaluation of the Carlson elliptic integrals and the incomplete elliptic integrals with complex arguments.+homepage:            https://github.com/stla/elliptic-integrals#readme+license:             BSD3+license-file:        LICENSE+author:              Stéphane Laurent+maintainer:          laurent_step@outlook.fr+copyright:           2023 Stéphane Laurent+category:            Math, Numeric+build-type:          Simple+extra-source-files:  README.md+                     CHANGELOG.md+cabal-version:       >=1.10++library+  hs-source-dirs:      src+  exposed-modules:     Math.EllipticIntegrals+  other-modules:       Math.EllipticIntegrals.Internal+                     , Math.EllipticIntegrals.Carlson+                     , Math.EllipticIntegrals.Elliptic+  build-depends:       base >= 4.7 && < 5+  default-language:    Haskell2010+  ghc-options:         -Wall++test-suite unit-tests+  type:                 exitcode-stdio-1.0+  main-is:              Main.hs+  hs-source-dirs:       tests/+  other-modules:        Approx+  Build-Depends:        base >= 4.7 && < 5+                      , tasty+                      , tasty-hunit+                      , elliptic-integrals+  Default-Language:     Haskell2010++source-repository head+  type:     git+  location: https://github.com/stla/elliptic-integrals
+ src/Math/EllipticIntegrals.hs view
@@ -0,0 +1,3 @@+module Math.EllipticIntegrals (module X) where+import Math.EllipticIntegrals.Carlson       as X+import Math.EllipticIntegrals.Elliptic      as X
+ src/Math/EllipticIntegrals/Carlson.hs view
@@ -0,0 +1,224 @@+module Math.EllipticIntegrals.Carlson+  (carlsonRF, carlsonRF',+  carlsonRD, carlsonRD',+  carlsonRJ, carlsonRJ',+  carlsonRC, carlsonRC',+  carlsonRG, carlsonRG')+  where+import           Data.Complex+import           Math.EllipticIntegrals.Internal++rf_ :: Cplx -> Cplx -> Cplx -> Double -> ((Double,Double,Double), Cplx)+rf_ x y z err =+  let a = (x+y+z)/3 in+  let delta = map (\u -> magnitude(1-u/a)) [x,y,z] in+  if maximum delta < err+    then ((delta !! 0, delta !! 1, delta !! 2), a)+    else+      let (sqrtx, sqrty, sqrtz) = (sqrt x, sqrt y, sqrt z) in+      let lambda = sqrtx*sqrty + sqrty*sqrtz + sqrtz*sqrtx in+      rf_ ((x+lambda)/4) ((y+lambda)/4) ((z+lambda)/4) err++-- | Carlson integral RF.+carlsonRF' :: +     Double -- ^ bound on the relative error+  -> Cplx   -- ^ first variable+  -> Cplx   -- ^ second variable +  -> Cplx   -- ^ third variable +  -> Cplx+carlsonRF' err x y z =+  if zeros > 1+    then error "At most one of x, y, z can be 0"+    else+      let ((dx,dy,dz), a) = rf_ x y z err in+      let (e2,e3) = (dx*dy + dy*dz + dz*dx, dx*dy*dz) in+      toCplx(1 - e2/10 + e3/14 + e2*e2/24 - 3*e2*e3/44 - 5*e2*e2*e2/208 ++          3*e3*e3/104 + e2*e2*e3/16) / sqrt a+  where+    zeros = sum (map (\u -> fromEnum (u == 0)) [x,y,z])++-- | Carlson integral RF.+carlsonRF :: +     Cplx   -- ^ first variable+  -> Cplx   -- ^ second variable +  -> Cplx   -- ^ third variable +  -> Cplx+carlsonRF = carlsonRF' 1e-15++rd_ :: Cplx -> Cplx -> Cplx -> Cplx -> Cplx -> Double ->+       ((Double,Double,Double), Cplx, Cplx, Cplx)+rd_ x y z s fac err =+  let a = (x+y+z+z+z)/5 in+  let delta = map (\u -> magnitude(1-u/a)) [x,y,z] in+  if maximum delta < err+    then ((delta !! 0, delta !! 1, delta !! 2), a, s, fac)+    else+      let (sqrtx, sqrty, sqrtz) = (sqrt x, sqrt y, sqrt z) in+      let lambda = sqrtx*sqrty + sqrty*sqrtz + sqrtz*sqrtx in+      let s' = s + fac / (sqrt z * (z + lambda)) in+      rd_ ((x+lambda)/4) ((y+lambda)/4) ((z+lambda)/4) s' (fac/4) err++-- | Carlson integral RD.+carlsonRD' ::+     Double -- ^ bound on the relative error+  -> Cplx   -- ^ first variable+  -> Cplx   -- ^ second variable +  -> Cplx   -- ^ third variable +  -> Cplx+carlsonRD' err x y z =+  if zeros > 1+    then error "At most one of x, y, z can be 0"+    else+      let ((dx,dy,dz), a, s, fac) = rd_ x y z 0 1 err in+      let+        (e2,e3,e4,e5) = (dx*dy + dy*dz + 3*dz*dz + 2*dz*dx + dx*dz + 2*dy*dz,+          dz*dz*dz + dx*dz*dz + 3*dx*dy*dz + 2*dy*dz*dz + dy*dz*dz + 2*dx*dz*dz,+          dy*dz*dz*dz + dx*dz*dz*dz + dx*dy*dz*dz + 2*dx*dy*dz*dz,+          dx*dy*dz*dz*dz) in+      3*s + fac * toCplx(1 - 3*e2/14 + e3/6 + 9*e2*e2/88 - 3*e4/22 - 9*e2*e3/52 ++        3*e5/26 - e2*e2*e2/16 + 3*e3*e3/40 + 3*e2*e4/20 + 45*e2*e2*e3/272 -+        9*(e3*e4 + e2*e5)/68) / a / sqrt a+  where+    zeros = sum (map (\u -> fromEnum (u == 0)) [x,y,z])++-- | Carlson integral RD.+carlsonRD ::+     Cplx   -- ^ first variable+  -> Cplx   -- ^ second variable +  -> Cplx   -- ^ third variable +  -> Cplx+carlsonRD = carlsonRD' 1e-15++rj_ :: Cplx -> Cplx -> Cplx -> Cplx -> Cplx -> Double -> Cplx -> Int ->+       Double -> [Cplx] -> [Cplx] -> Double -> (Cplx, Int, [Cplx], [Cplx])+rj_ x y z p a maxmagns delta f fac d e err =+  let q = (4/err)**(1/6) * maxmagns / fromIntegral f in+  if magnitude a > q+    then (a, f, d, e)+    else+      let dnew = (sqrt p + sqrt x)*(sqrt p + sqrt y)*(sqrt p + sqrt z)+          d' = (fromIntegral f * dnew) : d+          e' = (toCplx fac * delta / dnew / dnew) : e+          lambda = sqrt x * sqrt y + sqrt y * sqrt z + sqrt z * sqrt x+          x' = (x + lambda) / 4+          y' = (y + lambda) / 4+          z' = (z + lambda) / 4+          p' = (p + lambda) / 4+          a' = (a + lambda) / 4+      in+      rj_ x' y' z' p' a' maxmagns delta (4*f) (fac/64) d' e' err++-- | Carlson integral RJ.+carlsonRJ' ::+     Double -- ^ bound on the relative error+  -> Cplx   -- ^ first variable+  -> Cplx   -- ^ second variable +  -> Cplx   -- ^ third variable +  -> Cplx   -- ^ fourth variable +  -> Cplx+carlsonRJ' err x y z p =+  if zeros > 1+    then error "At most one of x, y, z, p can be 0"+    else+      let a0 = (x + y + z + p + p) / 5+          maxmagns = maximum $ map (\u -> magnitude(a0-u)) [x, y, z, p]+          delta = (p-x)*(p-y)*(p-z)+      in+      let (a, f, d, e) = rj_ x y z p a0 maxmagns delta 1 1 [] [] err+          f' = fromIntegral f+      in+      let x' = (a0 - x) / f' / a+          y' = (a0 - y) / f' / a+          z' = (a0 - z) / f' / a+          p' = -(x'+y'+z') / 2+          e2 = x'*y' + x'*z' + y'*z' - 3*p'*p'+          e3 = x'*y'*z' + 2*e2*p' + 4*p'*p'*p'+          e4 = p'*(2*x'*y'*z' + e2*p' + 3*p'*p'*p')+          e5 = x'*y'*z'*p'*p'+          h = zipWith (\u v -> atanx_over_x(sqrt u) / v) e d+      in+      (1 - 3*e2/14 + e3/6 + 9*e2*e2/88 - 3*e4/22 - 9*e2*e3/52 + 3*e5/26) /+        f' / a / sqrt a + 6 * sum h+  where+    zeros = sum (map (\u -> fromEnum (u == 0)) [x,y,z,p])+    atanx_over_x w = if w == 0 then 1 else atanC w / w++-- | Carlson integral RJ.+carlsonRJ ::+     Cplx   -- ^ first variable+  -> Cplx   -- ^ second variable +  -> Cplx   -- ^ third variable +  -> Cplx   -- ^ fourth variable +  -> Cplx+carlsonRJ = carlsonRJ' 1e-15+++rc_ :: Cplx -> Cplx -> Cplx -> Double -> Int -> Double -> (Cplx, Int)+rc_ x y a magn f err =+  let q = (1/3/err)**(1/8) * magn / fromIntegral f in+  if magnitude a > q+    then (a, f)+    else+      let lambda = 2 * sqrt x * sqrt y + y+          a' = (a + lambda) / 4+          x' = (x + lambda) / 4+          y' = (y + lambda) / 4+      in+      rc_ x' y' a' magn (4*f) err++-- | Carlson integral RC.+carlsonRC' ::+     Double -- ^ bound on the relative error+  -> Cplx   -- ^ first variable+  -> Cplx   -- ^ second variable +  -> Cplx+carlsonRC' err x y =+  if y == 0+    then error "y cannot be 0"+    else+      let a0 = (x + y + y) / 3+          magn = magnitude(a0-x)+      in+      let (a, f) = rc_ x y a0 magn 1 err+          f' = fromIntegral f+      in+      let s = (y - a0) / f' / a in+      (1 + 3*s*s/10 + s*s*s/7 + 3*s*s*s*s/8 + 9*s*s*s*s*s/22 ++        159*s*s*s*s*s*s/208 + 9*s*s*s*s*s*s*s/8) / sqrt a++-- | Carlson integral RC.+carlsonRC ::+     Cplx   -- ^ first variable+  -> Cplx   -- ^ second variable +  -> Cplx+carlsonRC = carlsonRC' 1e-15+++-- | Carlson integral RG.+carlsonRG' ::+     Double -- ^ bound on the relative error passed to `CarlsonRD'`+  -> Cplx   -- ^ first variable+  -> Cplx   -- ^ second variable +  -> Cplx   -- ^ third variable +  -> Cplx+carlsonRG' err x y z =+  if zeros > 1+    then sqrt(x+y+z) / 2+    else+      if z == 0+        then carlsonRG' err z x y+        else+          (z * carlsonRF' err x y z -+            (x-z) * (y-z) * carlsonRD' err x y z / 3 ++            sqrt x * sqrt y / sqrt z) / 2+  where+    zeros = sum (map (\u -> fromEnum (u == 0)) [x,y,z])++-- | Carlson integral RG.+carlsonRG ::+     Cplx   -- ^ first variable+  -> Cplx   -- ^ second variable +  -> Cplx   -- ^ third variable +  -> Cplx+carlsonRG = carlsonRG' 1e-15+
+ src/Math/EllipticIntegrals/Elliptic.hs view
@@ -0,0 +1,129 @@+module Math.EllipticIntegrals.Elliptic+  where+import Math.EllipticIntegrals.Carlson+import Data.Complex+import Math.EllipticIntegrals.Internal++-- | Elliptic integral of the first kind.+ellipticF' :: +     Double -- ^ bound on the relative error passed to `carlsonRF'`+  -> Cplx   -- ^ amplitude+  -> Cplx   -- ^ parameter+  -> Cplx+ellipticF' err phi m+  | phi == 0 =+    toCplx 0+  | m == 1 && abs(realPart phi) == pi/2 =+    toCplx (0/0)+  | m == 1 && abs(realPart phi) < pi/2 =+    atanh(sin phi)+  | abs(realPart phi) <= pi/2 =+    if m == 0+      then+        phi+      else+        let sine = sin phi in+        let sine2 = sine*sine in+        let (cosine2, oneminusmsine2) = (1 - sine2, 1 - m*sine2) in+        sine * carlsonRF' err cosine2 oneminusmsine2 1+  | otherwise =+    let (phi', k) = getPhiK phi in+    2 * fromIntegral k * ellipticF' err (pi/2) m + ellipticF' err phi' m++-- | Elliptic integral of the first kind.+ellipticF :: +     Cplx -- ^ amplitude+  -> Cplx -- ^ parameter+  -> Cplx+ellipticF = ellipticF' 1e-15++-- | Elliptic integral of the second kind.+ellipticE' :: +     Double -- ^ bound on the relative error passed to `carlsonRF'` and `carlsonRD'` +  -> Cplx   -- ^ amplitude+  -> Cplx   -- ^ parameter+  -> Cplx+ellipticE' err phi m+  | phi == 0 =+    toCplx 0+  | abs(realPart phi) <= pi/2 =+    case m of+      0 -> phi+      1 -> sin phi+      _ ->+        let sine = sin phi in+        let sine2 = sine*sine in+        let (cosine2, oneminusmsine2) = (1 - sine2, 1 - m*sine2) in+        sine * (carlsonRF' err cosine2 oneminusmsine2 1 -+          m * sine2 / 3 * carlsonRD' err cosine2 oneminusmsine2 1)+  | otherwise =+    let (phi', k) = getPhiK phi in+    2 * fromIntegral k * ellipticE' err (pi/2) m + ellipticE' err phi' m++-- | Elliptic integral of the second kind.+ellipticE :: +     Cplx -- ^ amplitude+  -> Cplx -- ^ parameter+  -> Cplx+ellipticE = ellipticE' 1e-15++-- | Elliptic integral of the third kind.+ellipticPI' :: +     Double -- ^ bound on the relative error passed to `carlsonRF'` and `carlsonRJ'` +  -> Cplx   -- ^ amplitude+  -> Cplx   -- ^ characteristic+  -> Cplx   -- ^ parameter+  -> Cplx+ellipticPI' err phi n m+  | phi == 0 =+    toCplx 0+  | phi == pi/2 && n == 1 =+    0/0+  | phi == pi/2 && m == 0 =+    pi/2/sqrt(1-n)+  | phi == pi/2 && m == n =+    ellipticE' err (pi/2) m / (1-m)+  | phi == pi/2 && n == 0 =+    ellipticF' err (pi/2) m+  | abs(realPart phi) <= pi/2 =+    let sine = sin phi in+    let sine2 = sine*sine in+    let (cosine2, oneminusmsine2) = (1 - sine2, 1 - m*sine2) in+    sine * (carlsonRF' err cosine2 oneminusmsine2 1 ++      n * sine2 / 3 * carlsonRJ' err cosine2 oneminusmsine2 1 (1-n*sine2))+  | otherwise =+    let (phi', k) = getPhiK phi in+    2 * fromIntegral k * ellipticPI' err (pi/2) n m + ellipticPI' err phi' n m++-- | Elliptic integral of the third kind.+ellipticPI ::+     Cplx -- ^ amplitude+  -> Cplx -- ^ characteristic+  -> Cplx -- ^ parameter+  -> Cplx+ellipticPI = ellipticPI' 1e-15++-- | Jacobi zeta function.+jacobiZeta' ::+     Double -- ^ bound on the relative error passed to `ellipticF'` and `ellipticE'` +  -> Cplx   -- ^ amplitude+  -> Cplx   -- ^ parameter+  -> Cplx+jacobiZeta' err phi m =+  if m == 1+    then+      if abs(realPart phi) <= pi/2+        then sin phi+        else let (phi',_) = getPhiK phi in sin phi'+    else+      ellipticE' err phi m -+        ellipticE' err (pi/2) m / ellipticF' err (pi/2) m *+        ellipticF' err phi m++-- | Jacobi zeta function.+jacobiZeta ::+     Cplx -- ^ amplitude+  -> Cplx -- ^ parameter+  -> Cplx+jacobiZeta = jacobiZeta' 1e-15+
+ src/Math/EllipticIntegrals/Internal.hs view
@@ -0,0 +1,21 @@+module Math.EllipticIntegrals.Internal+  where+import Data.Complex++type Cplx = Complex Double++toCplx :: Double -> Cplx+toCplx x = x :+ 0.0++getPhiK :: Cplx -> (Cplx, Int)+getPhiK phi+  | realPart phi > pi/2 =+    until (\(x,_) -> realPart x <= pi/2) (\(x,k) -> (x-pi,k+1)) (phi,0)+  | realPart phi < -pi/2 =+    until (\(x,_) -> realPart x >= -pi/2) (\(x,k) -> (x+pi,k-1)) (phi,0)+  | otherwise = (phi,0)++atanC :: Cplx -> Cplx+atanC z = i * (log(1-i*z) - log(1+i*z)) / 2+  where+    i = 0.0 :+ 1.0
+ tests/Approx.hs view
@@ -0,0 +1,8 @@+module Approx where+import Data.Complex++approxDbl :: Int -> Double -> Double+approxDbl n x = fromInteger (round $ x * (10^n)) / (10.0^^n)++approx :: Int -> Complex Double -> Complex Double+approx n z = approxDbl n (realPart z) :+ approxDbl n (imagPart z)
+ tests/Main.hs view
@@ -0,0 +1,123 @@+module Main where+import           Approx+import           Data.Complex+import           Math.EllipticIntegrals+import           Test.Tasty             (defaultMain, testGroup)+import           Test.Tasty.HUnit       (assertEqual, testCase)++i :: Complex Double+i = 0.0 :+ 1.0++main :: IO ()+main = defaultMain $+  testGroup "Tests"+  [ testCase "RF value 1" $+      assertEqual ""+        (approx 12 (carlsonRF 1 2 0))+        (approx 12 1.3110287771461),++    testCase "RF value 2" $+      assertEqual ""+        (approx 12 (carlsonRF i (-i) 0))+        (approx 12 1.8540746773014),++    testCase "RF value 3" $+      assertEqual ""+        (approx 12 (carlsonRF 0.5 1 0))+        (approx 12 1.8540746773014),++    testCase "RF value 4" $+      assertEqual ""+        (approx 13 (carlsonRF (i-1) i 0))+        (approx 13 (0.79612586584234 :+ (-1.2138566698365))),++    testCase "RF value 5" $+      assertEqual ""+        (approx 13 (carlsonRF 2 3 4))+        (approx 13 0.58408284167715),++    testCase "RF value 6" $+      assertEqual ""+        (approx 12 (carlsonRF i (-i) 2))+        (approx 12 1.0441445654064),++    testCase "RF value 7" $+      assertEqual ""+        (approx 13 (carlsonRF (i-1) i (1-i)))+        (approx 13 (0.93912050218619 :+ (-0.53296252018635))),++    testCase "RC value 1" $+      assertEqual ""+        (approx 14 (carlsonRC 0 0.25))+        (approx 14 pi),++    testCase "RC value 2" $+      assertEqual ""+        (approx 14 (carlsonRC 2.25 2))+        (approx 14 (log 2)),++    testCase "RC value 3" $+      assertEqual ""+        (approx 12 (carlsonRC 0 i))+        (approx 12 ((1-i)*1.1107207345396)),++    testCase "RC value 4" $+      assertEqual ""+        (approx 13 (carlsonRC (-i) i))+        (approx 13 (1.2260849569072 :+ (-0.34471136988768))),++    testCase "RC value 5" $+      assertEqual ""+        (approx 14 (carlsonRC 0.25 (-2)))+        (approx 14 ((log 2 :+ (-pi))/ 3)),++    testCase "RC x y = RF x y y" $ do+      let x = 5 :+ 6+          y = 2 :+ (-9)+      assertEqual ""+        (approx 14 (carlsonRC x y))+        (approx 14 (carlsonRF x y y)),++    testCase "RJ x y y p" $ do+      let x = 1 :+ 1+          y = (-2) :+ 3+          p = 0 :+ 4+      assertEqual ""+        (approx 14 (carlsonRJ x y y p))+        (approx 14 (3*(carlsonRC x y - carlsonRC x p) / (p-y))),++    testCase "RJ homogeneity" $ do+      let x = 1 :+ 1+          y = (-2) :+ 3+          z = -3+          p = 0 :+ 4+          kappa = 2 :+ 0+      assertEqual ""+        (approx 14 (carlsonRJ x y z p / kappa / sqrt kappa))+        (approx 14 (carlsonRJ (kappa*x) (kappa*y) (kappa*z) (kappa*p))),++    testCase "Complete elliptic integral K" $ do+      let m = 2 :+ (-3)+      assertEqual ""+        (approx 14 (ellipticF (pi/2) m))+        (approx 14 (carlsonRF 0 (1-m) 1)),++    testCase "Complete ellipticE - RG" $ do+      let m = 2 :+ (-3)+      assertEqual ""+        (approx 14 (ellipticE (pi/2) m))+        (approx 14 (2 * carlsonRG 0 (1-m) 1)),++    testCase "Complete ellipticE - RD" $ do+      let m = 2 :+ (-3)+      assertEqual ""+        (approx 14 (ellipticE (pi/2) m))+        (approx 14 ((1-m) * (carlsonRD 0 (1-m) 1 + carlsonRD 0 1 (1-m)) / 3)),++    testCase "jacobiZeta m=1" $ do+      let z = (-1) :+ 8+      assertEqual ""+        (approx 14 (jacobiZeta z 1))+        (approx 14 (sin z))++  ]