diff --git a/CHANGELOG.md b/CHANGELOG.md
new file mode 100644
--- /dev/null
+++ b/CHANGELOG.md
@@ -0,0 +1,4 @@
+0.1.0.0
+-------
+* initial release
+
diff --git a/LICENSE b/LICENSE
new file mode 100644
--- /dev/null
+++ b/LICENSE
@@ -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.
diff --git a/README.md b/README.md
new file mode 100644
--- /dev/null
+++ b/README.md
@@ -0,0 +1,3 @@
+# elliptic-integrals
+
+Evaluation of the Carlson elliptic integrals and the incomplete elliptic integrals with complex arguments.
diff --git a/Setup.hs b/Setup.hs
new file mode 100644
--- /dev/null
+++ b/Setup.hs
@@ -0,0 +1,2 @@
+import Distribution.Simple
+main = defaultMain
diff --git a/elliptic-integrals.cabal b/elliptic-integrals.cabal
new file mode 100644
--- /dev/null
+++ b/elliptic-integrals.cabal
@@ -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
diff --git a/src/Math/EllipticIntegrals.hs b/src/Math/EllipticIntegrals.hs
new file mode 100644
--- /dev/null
+++ b/src/Math/EllipticIntegrals.hs
@@ -0,0 +1,3 @@
+module Math.EllipticIntegrals (module X) where
+import Math.EllipticIntegrals.Carlson       as X
+import Math.EllipticIntegrals.Elliptic      as X
diff --git a/src/Math/EllipticIntegrals/Carlson.hs b/src/Math/EllipticIntegrals/Carlson.hs
new file mode 100644
--- /dev/null
+++ b/src/Math/EllipticIntegrals/Carlson.hs
@@ -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
+
diff --git a/src/Math/EllipticIntegrals/Elliptic.hs b/src/Math/EllipticIntegrals/Elliptic.hs
new file mode 100644
--- /dev/null
+++ b/src/Math/EllipticIntegrals/Elliptic.hs
@@ -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
+
diff --git a/src/Math/EllipticIntegrals/Internal.hs b/src/Math/EllipticIntegrals/Internal.hs
new file mode 100644
--- /dev/null
+++ b/src/Math/EllipticIntegrals/Internal.hs
@@ -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
diff --git a/tests/Approx.hs b/tests/Approx.hs
new file mode 100644
--- /dev/null
+++ b/tests/Approx.hs
@@ -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)
diff --git a/tests/Main.hs b/tests/Main.hs
new file mode 100644
--- /dev/null
+++ b/tests/Main.hs
@@ -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))
+
+  ]