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_ :: Cplx -> Cplx -> Cplx -> Double -> ((Double, Double, Double), Cplx)
rf_ Cplx
x Cplx
y Cplx
z Double
err =
  let a :: Cplx
a = (Cplx
xforall a. Num a => a -> a -> a
+Cplx
yforall a. Num a => a -> a -> a
+Cplx
z)forall a. Fractional a => a -> a -> a
/Cplx
3 in
  let delta :: [Double]
delta = forall a b. (a -> b) -> [a] -> [b]
map (\Cplx
u -> forall a. RealFloat a => Complex a -> a
magnitude(Cplx
1forall a. Num a => a -> a -> a
-Cplx
uforall a. Fractional a => a -> a -> a
/Cplx
a)) [Cplx
x,Cplx
y,Cplx
z] in
  if forall (t :: * -> *) a. (Foldable t, Ord a) => t a -> a
maximum [Double]
delta forall a. Ord a => a -> a -> Bool
< Double
err
    then (([Double]
delta forall a. [a] -> Int -> a
!! Int
0, [Double]
delta forall a. [a] -> Int -> a
!! Int
1, [Double]
delta forall a. [a] -> Int -> a
!! Int
2), Cplx
a)
    else
      let (Cplx
sqrtx, Cplx
sqrty, Cplx
sqrtz) = (forall a. Floating a => a -> a
sqrt Cplx
x, forall a. Floating a => a -> a
sqrt Cplx
y, forall a. Floating a => a -> a
sqrt Cplx
z) in
      let lambda :: Cplx
lambda = Cplx
sqrtxforall a. Num a => a -> a -> a
*Cplx
sqrty forall a. Num a => a -> a -> a
+ Cplx
sqrtyforall a. Num a => a -> a -> a
*Cplx
sqrtz forall a. Num a => a -> a -> a
+ Cplx
sqrtzforall a. Num a => a -> a -> a
*Cplx
sqrtx in
      Cplx -> Cplx -> Cplx -> Double -> ((Double, Double, Double), Cplx)
rf_ ((Cplx
xforall a. Num a => a -> a -> a
+Cplx
lambda)forall a. Fractional a => a -> a -> a
/Cplx
4) ((Cplx
yforall a. Num a => a -> a -> a
+Cplx
lambda)forall a. Fractional a => a -> a -> a
/Cplx
4) ((Cplx
zforall a. Num a => a -> a -> a
+Cplx
lambda)forall a. Fractional a => a -> a -> a
/Cplx
4) Double
err

-- | Carlson integral RF.
carlsonRF' :: 
     Double -- ^ bound on the relative error
  -> Cplx   -- ^ first variable
  -> Cplx   -- ^ second variable 
  -> Cplx   -- ^ third variable 
  -> Cplx
carlsonRF' :: Double -> Cplx -> Cplx -> Cplx -> Cplx
carlsonRF' Double
err Cplx
x Cplx
y Cplx
z =
  if Int
zeros forall a. Ord a => a -> a -> Bool
> Int
1
    then forall a. HasCallStack => [Char] -> a
error [Char]
"At most one of x, y, z can be 0"
    else
      let ((Double
dx,Double
dy,Double
dz), Cplx
a) = Cplx -> Cplx -> Cplx -> Double -> ((Double, Double, Double), Cplx)
rf_ Cplx
x Cplx
y Cplx
z Double
err in
      let (Double
e2,Double
e3) = (Double
dxforall a. Num a => a -> a -> a
*Double
dy forall a. Num a => a -> a -> a
+ Double
dyforall a. Num a => a -> a -> a
*Double
dz forall a. Num a => a -> a -> a
+ Double
dzforall a. Num a => a -> a -> a
*Double
dx, Double
dxforall a. Num a => a -> a -> a
*Double
dyforall a. Num a => a -> a -> a
*Double
dz) in
      Double -> Cplx
toCplx(Double
1 forall a. Num a => a -> a -> a
- Double
e2forall a. Fractional a => a -> a -> a
/Double
10 forall a. Num a => a -> a -> a
+ Double
e3forall a. Fractional a => a -> a -> a
/Double
14 forall a. Num a => a -> a -> a
+ Double
e2forall a. Num a => a -> a -> a
*Double
e2forall a. Fractional a => a -> a -> a
/Double
24 forall a. Num a => a -> a -> a
- Double
3forall a. Num a => a -> a -> a
*Double
e2forall a. Num a => a -> a -> a
*Double
e3forall a. Fractional a => a -> a -> a
/Double
44 forall a. Num a => a -> a -> a
- Double
5forall a. Num a => a -> a -> a
*Double
e2forall a. Num a => a -> a -> a
*Double
e2forall a. Num a => a -> a -> a
*Double
e2forall a. Fractional a => a -> a -> a
/Double
208 forall a. Num a => a -> a -> a
+
          Double
3forall a. Num a => a -> a -> a
*Double
e3forall a. Num a => a -> a -> a
*Double
e3forall a. Fractional a => a -> a -> a
/Double
104 forall a. Num a => a -> a -> a
+ Double
e2forall a. Num a => a -> a -> a
*Double
e2forall a. Num a => a -> a -> a
*Double
e3forall a. Fractional a => a -> a -> a
/Double
16) forall a. Fractional a => a -> a -> a
/ forall a. Floating a => a -> a
sqrt Cplx
a
  where
    zeros :: Int
zeros = forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum (forall a b. (a -> b) -> [a] -> [b]
map (\Cplx
u -> forall a. Enum a => a -> Int
fromEnum (Cplx
u forall a. Eq a => a -> a -> Bool
== Cplx
0)) [Cplx
x,Cplx
y,Cplx
z])

-- | Carlson integral RF.
carlsonRF :: 
     Cplx   -- ^ first variable
  -> Cplx   -- ^ second variable 
  -> Cplx   -- ^ third variable 
  -> Cplx
carlsonRF :: Cplx -> Cplx -> Cplx -> Cplx
carlsonRF = Double -> Cplx -> Cplx -> Cplx -> Cplx
carlsonRF' Double
1e-15

rd_ :: Cplx -> Cplx -> Cplx -> Cplx -> Cplx -> Double ->
       ((Double,Double,Double), Cplx, Cplx, Cplx)
rd_ :: Cplx
-> Cplx
-> Cplx
-> Cplx
-> Cplx
-> Double
-> ((Double, Double, Double), Cplx, Cplx, Cplx)
rd_ Cplx
x Cplx
y Cplx
z Cplx
s Cplx
fac Double
err =
  let a :: Cplx
a = (Cplx
xforall a. Num a => a -> a -> a
+Cplx
yforall a. Num a => a -> a -> a
+Cplx
zforall a. Num a => a -> a -> a
+Cplx
zforall a. Num a => a -> a -> a
+Cplx
z)forall a. Fractional a => a -> a -> a
/Cplx
5 in
  let delta :: [Double]
delta = forall a b. (a -> b) -> [a] -> [b]
map (\Cplx
u -> forall a. RealFloat a => Complex a -> a
magnitude(Cplx
1forall a. Num a => a -> a -> a
-Cplx
uforall a. Fractional a => a -> a -> a
/Cplx
a)) [Cplx
x,Cplx
y,Cplx
z] in
  if forall (t :: * -> *) a. (Foldable t, Ord a) => t a -> a
maximum [Double]
delta forall a. Ord a => a -> a -> Bool
< Double
err
    then (([Double]
delta forall a. [a] -> Int -> a
!! Int
0, [Double]
delta forall a. [a] -> Int -> a
!! Int
1, [Double]
delta forall a. [a] -> Int -> a
!! Int
2), Cplx
a, Cplx
s, Cplx
fac)
    else
      let (Cplx
sqrtx, Cplx
sqrty, Cplx
sqrtz) = (forall a. Floating a => a -> a
sqrt Cplx
x, forall a. Floating a => a -> a
sqrt Cplx
y, forall a. Floating a => a -> a
sqrt Cplx
z) in
      let lambda :: Cplx
lambda = Cplx
sqrtxforall a. Num a => a -> a -> a
*Cplx
sqrty forall a. Num a => a -> a -> a
+ Cplx
sqrtyforall a. Num a => a -> a -> a
*Cplx
sqrtz forall a. Num a => a -> a -> a
+ Cplx
sqrtzforall a. Num a => a -> a -> a
*Cplx
sqrtx in
      let s' :: Cplx
s' = Cplx
s forall a. Num a => a -> a -> a
+ Cplx
fac forall a. Fractional a => a -> a -> a
/ (forall a. Floating a => a -> a
sqrt Cplx
z forall a. Num a => a -> a -> a
* (Cplx
z forall a. Num a => a -> a -> a
+ Cplx
lambda)) in
      Cplx
-> Cplx
-> Cplx
-> Cplx
-> Cplx
-> Double
-> ((Double, Double, Double), Cplx, Cplx, Cplx)
rd_ ((Cplx
xforall a. Num a => a -> a -> a
+Cplx
lambda)forall a. Fractional a => a -> a -> a
/Cplx
4) ((Cplx
yforall a. Num a => a -> a -> a
+Cplx
lambda)forall a. Fractional a => a -> a -> a
/Cplx
4) ((Cplx
zforall a. Num a => a -> a -> a
+Cplx
lambda)forall a. Fractional a => a -> a -> a
/Cplx
4) Cplx
s' (Cplx
facforall a. Fractional a => a -> a -> a
/Cplx
4) Double
err

-- | Carlson integral RD.
carlsonRD' ::
     Double -- ^ bound on the relative error
  -> Cplx   -- ^ first variable
  -> Cplx   -- ^ second variable 
  -> Cplx   -- ^ third variable 
  -> Cplx
carlsonRD' :: Double -> Cplx -> Cplx -> Cplx -> Cplx
carlsonRD' Double
err Cplx
x Cplx
y Cplx
z =
  if Int
zeros forall a. Ord a => a -> a -> Bool
> Int
1
    then forall a. HasCallStack => [Char] -> a
error [Char]
"At most one of x, y, z can be 0"
    else
      let ((Double
dx,Double
dy,Double
dz), Cplx
a, Cplx
s, Cplx
fac) = Cplx
-> Cplx
-> Cplx
-> Cplx
-> Cplx
-> Double
-> ((Double, Double, Double), Cplx, Cplx, Cplx)
rd_ Cplx
x Cplx
y Cplx
z Cplx
0 Cplx
1 Double
err in
      let
        (Double
e2,Double
e3,Double
e4,Double
e5) = (Double
dxforall a. Num a => a -> a -> a
*Double
dy forall a. Num a => a -> a -> a
+ Double
dyforall a. Num a => a -> a -> a
*Double
dz forall a. Num a => a -> a -> a
+ Double
3forall a. Num a => a -> a -> a
*Double
dzforall a. Num a => a -> a -> a
*Double
dz forall a. Num a => a -> a -> a
+ Double
2forall a. Num a => a -> a -> a
*Double
dzforall a. Num a => a -> a -> a
*Double
dx forall a. Num a => a -> a -> a
+ Double
dxforall a. Num a => a -> a -> a
*Double
dz forall a. Num a => a -> a -> a
+ Double
2forall a. Num a => a -> a -> a
*Double
dyforall a. Num a => a -> a -> a
*Double
dz,
          Double
dzforall a. Num a => a -> a -> a
*Double
dzforall a. Num a => a -> a -> a
*Double
dz forall a. Num a => a -> a -> a
+ Double
dxforall a. Num a => a -> a -> a
*Double
dzforall a. Num a => a -> a -> a
*Double
dz forall a. Num a => a -> a -> a
+ Double
3forall a. Num a => a -> a -> a
*Double
dxforall a. Num a => a -> a -> a
*Double
dyforall a. Num a => a -> a -> a
*Double
dz forall a. Num a => a -> a -> a
+ Double
2forall a. Num a => a -> a -> a
*Double
dyforall a. Num a => a -> a -> a
*Double
dzforall a. Num a => a -> a -> a
*Double
dz forall a. Num a => a -> a -> a
+ Double
dyforall a. Num a => a -> a -> a
*Double
dzforall a. Num a => a -> a -> a
*Double
dz forall a. Num a => a -> a -> a
+ Double
2forall a. Num a => a -> a -> a
*Double
dxforall a. Num a => a -> a -> a
*Double
dzforall a. Num a => a -> a -> a
*Double
dz,
          Double
dyforall a. Num a => a -> a -> a
*Double
dzforall a. Num a => a -> a -> a
*Double
dzforall a. Num a => a -> a -> a
*Double
dz forall a. Num a => a -> a -> a
+ Double
dxforall a. Num a => a -> a -> a
*Double
dzforall a. Num a => a -> a -> a
*Double
dzforall a. Num a => a -> a -> a
*Double
dz forall a. Num a => a -> a -> a
+ Double
dxforall a. Num a => a -> a -> a
*Double
dyforall a. Num a => a -> a -> a
*Double
dzforall a. Num a => a -> a -> a
*Double
dz forall a. Num a => a -> a -> a
+ Double
2forall a. Num a => a -> a -> a
*Double
dxforall a. Num a => a -> a -> a
*Double
dyforall a. Num a => a -> a -> a
*Double
dzforall a. Num a => a -> a -> a
*Double
dz,
          Double
dxforall a. Num a => a -> a -> a
*Double
dyforall a. Num a => a -> a -> a
*Double
dzforall a. Num a => a -> a -> a
*Double
dzforall a. Num a => a -> a -> a
*Double
dz) in
      Cplx
3forall a. Num a => a -> a -> a
*Cplx
s forall a. Num a => a -> a -> a
+ Cplx
fac forall a. Num a => a -> a -> a
* Double -> Cplx
toCplx(Double
1 forall a. Num a => a -> a -> a
- Double
3forall a. Num a => a -> a -> a
*Double
e2forall a. Fractional a => a -> a -> a
/Double
14 forall a. Num a => a -> a -> a
+ Double
e3forall a. Fractional a => a -> a -> a
/Double
6 forall a. Num a => a -> a -> a
+ Double
9forall a. Num a => a -> a -> a
*Double
e2forall a. Num a => a -> a -> a
*Double
e2forall a. Fractional a => a -> a -> a
/Double
88 forall a. Num a => a -> a -> a
- Double
3forall a. Num a => a -> a -> a
*Double
e4forall a. Fractional a => a -> a -> a
/Double
22 forall a. Num a => a -> a -> a
- Double
9forall a. Num a => a -> a -> a
*Double
e2forall a. Num a => a -> a -> a
*Double
e3forall a. Fractional a => a -> a -> a
/Double
52 forall a. Num a => a -> a -> a
+
        Double
3forall a. Num a => a -> a -> a
*Double
e5forall a. Fractional a => a -> a -> a
/Double
26 forall a. Num a => a -> a -> a
- Double
e2forall a. Num a => a -> a -> a
*Double
e2forall a. Num a => a -> a -> a
*Double
e2forall a. Fractional a => a -> a -> a
/Double
16 forall a. Num a => a -> a -> a
+ Double
3forall a. Num a => a -> a -> a
*Double
e3forall a. Num a => a -> a -> a
*Double
e3forall a. Fractional a => a -> a -> a
/Double
40 forall a. Num a => a -> a -> a
+ Double
3forall a. Num a => a -> a -> a
*Double
e2forall a. Num a => a -> a -> a
*Double
e4forall a. Fractional a => a -> a -> a
/Double
20 forall a. Num a => a -> a -> a
+ Double
45forall a. Num a => a -> a -> a
*Double
e2forall a. Num a => a -> a -> a
*Double
e2forall a. Num a => a -> a -> a
*Double
e3forall a. Fractional a => a -> a -> a
/Double
272 forall a. Num a => a -> a -> a
-
        Double
9forall a. Num a => a -> a -> a
*(Double
e3forall a. Num a => a -> a -> a
*Double
e4 forall a. Num a => a -> a -> a
+ Double
e2forall a. Num a => a -> a -> a
*Double
e5)forall a. Fractional a => a -> a -> a
/Double
68) forall a. Fractional a => a -> a -> a
/ Cplx
a forall a. Fractional a => a -> a -> a
/ forall a. Floating a => a -> a
sqrt Cplx
a
  where
    zeros :: Int
zeros = forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum (forall a b. (a -> b) -> [a] -> [b]
map (\Cplx
u -> forall a. Enum a => a -> Int
fromEnum (Cplx
u forall a. Eq a => a -> a -> Bool
== Cplx
0)) [Cplx
x,Cplx
y,Cplx
z])

-- | Carlson integral RD.
carlsonRD ::
     Cplx   -- ^ first variable
  -> Cplx   -- ^ second variable 
  -> Cplx   -- ^ third variable 
  -> Cplx
carlsonRD :: Cplx -> Cplx -> Cplx -> Cplx
carlsonRD = Double -> Cplx -> Cplx -> Cplx -> Cplx
carlsonRD' Double
1e-15

rj_ :: Cplx -> Cplx -> Cplx -> Cplx -> Cplx -> Double -> Cplx -> Int ->
       Double -> [Cplx] -> [Cplx] -> Double -> (Cplx, Int, [Cplx], [Cplx])
rj_ :: Cplx
-> Cplx
-> Cplx
-> Cplx
-> Cplx
-> Double
-> Cplx
-> Int
-> Double
-> [Cplx]
-> [Cplx]
-> Double
-> (Cplx, Int, [Cplx], [Cplx])
rj_ Cplx
x Cplx
y Cplx
z Cplx
p Cplx
a Double
maxmagns Cplx
delta Int
f Double
fac [Cplx]
d [Cplx]
e Double
err =
  let q :: Double
q = (Double
4forall a. Fractional a => a -> a -> a
/Double
err)forall a. Floating a => a -> a -> a
**(Double
1forall a. Fractional a => a -> a -> a
/Double
6) forall a. Num a => a -> a -> a
* Double
maxmagns forall a. Fractional a => a -> a -> a
/ forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
f in
  if forall a. RealFloat a => Complex a -> a
magnitude Cplx
a forall a. Ord a => a -> a -> Bool
> Double
q
    then (Cplx
a, Int
f, [Cplx]
d, [Cplx]
e)
    else
      let dnew :: Cplx
dnew = (forall a. Floating a => a -> a
sqrt Cplx
p forall a. Num a => a -> a -> a
+ forall a. Floating a => a -> a
sqrt Cplx
x)forall a. Num a => a -> a -> a
*(forall a. Floating a => a -> a
sqrt Cplx
p forall a. Num a => a -> a -> a
+ forall a. Floating a => a -> a
sqrt Cplx
y)forall a. Num a => a -> a -> a
*(forall a. Floating a => a -> a
sqrt Cplx
p forall a. Num a => a -> a -> a
+ forall a. Floating a => a -> a
sqrt Cplx
z)
          d' :: [Cplx]
d' = (forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
f forall a. Num a => a -> a -> a
* Cplx
dnew) forall a. a -> [a] -> [a]
: [Cplx]
d
          e' :: [Cplx]
e' = (Double -> Cplx
toCplx Double
fac forall a. Num a => a -> a -> a
* Cplx
delta forall a. Fractional a => a -> a -> a
/ Cplx
dnew forall a. Fractional a => a -> a -> a
/ Cplx
dnew) forall a. a -> [a] -> [a]
: [Cplx]
e
          lambda :: Cplx
lambda = forall a. Floating a => a -> a
sqrt Cplx
x forall a. Num a => a -> a -> a
* forall a. Floating a => a -> a
sqrt Cplx
y forall a. Num a => a -> a -> a
+ forall a. Floating a => a -> a
sqrt Cplx
y forall a. Num a => a -> a -> a
* forall a. Floating a => a -> a
sqrt Cplx
z forall a. Num a => a -> a -> a
+ forall a. Floating a => a -> a
sqrt Cplx
z forall a. Num a => a -> a -> a
* forall a. Floating a => a -> a
sqrt Cplx
x
          x' :: Cplx
x' = (Cplx
x forall a. Num a => a -> a -> a
+ Cplx
lambda) forall a. Fractional a => a -> a -> a
/ Cplx
4
          y' :: Cplx
y' = (Cplx
y forall a. Num a => a -> a -> a
+ Cplx
lambda) forall a. Fractional a => a -> a -> a
/ Cplx
4
          z' :: Cplx
z' = (Cplx
z forall a. Num a => a -> a -> a
+ Cplx
lambda) forall a. Fractional a => a -> a -> a
/ Cplx
4
          p' :: Cplx
p' = (Cplx
p forall a. Num a => a -> a -> a
+ Cplx
lambda) forall a. Fractional a => a -> a -> a
/ Cplx
4
          a' :: Cplx
a' = (Cplx
a forall a. Num a => a -> a -> a
+ Cplx
lambda) forall a. Fractional a => a -> a -> a
/ Cplx
4
      in
      Cplx
-> Cplx
-> Cplx
-> Cplx
-> Cplx
-> Double
-> Cplx
-> Int
-> Double
-> [Cplx]
-> [Cplx]
-> Double
-> (Cplx, Int, [Cplx], [Cplx])
rj_ Cplx
x' Cplx
y' Cplx
z' Cplx
p' Cplx
a' Double
maxmagns Cplx
delta (Int
4forall a. Num a => a -> a -> a
*Int
f) (Double
facforall a. Fractional a => a -> a -> a
/Double
64) [Cplx]
d' [Cplx]
e' Double
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' :: Double -> Cplx -> Cplx -> Cplx -> Cplx -> Cplx
carlsonRJ' Double
err Cplx
x Cplx
y Cplx
z Cplx
p =
  if Int
zeros forall a. Ord a => a -> a -> Bool
> Int
1
    then forall a. HasCallStack => [Char] -> a
error [Char]
"At most one of x, y, z, p can be 0"
    else
      let a0 :: Cplx
a0 = (Cplx
x forall a. Num a => a -> a -> a
+ Cplx
y forall a. Num a => a -> a -> a
+ Cplx
z forall a. Num a => a -> a -> a
+ Cplx
p forall a. Num a => a -> a -> a
+ Cplx
p) forall a. Fractional a => a -> a -> a
/ Cplx
5
          maxmagns :: Double
maxmagns = forall (t :: * -> *) a. (Foldable t, Ord a) => t a -> a
maximum forall a b. (a -> b) -> a -> b
$ forall a b. (a -> b) -> [a] -> [b]
map (\Cplx
u -> forall a. RealFloat a => Complex a -> a
magnitude(Cplx
a0forall a. Num a => a -> a -> a
-Cplx
u)) [Cplx
x, Cplx
y, Cplx
z, Cplx
p]
          delta :: Cplx
delta = (Cplx
pforall a. Num a => a -> a -> a
-Cplx
x)forall a. Num a => a -> a -> a
*(Cplx
pforall a. Num a => a -> a -> a
-Cplx
y)forall a. Num a => a -> a -> a
*(Cplx
pforall a. Num a => a -> a -> a
-Cplx
z)
      in
      let (Cplx
a, Int
f, [Cplx]
d, [Cplx]
e) = Cplx
-> Cplx
-> Cplx
-> Cplx
-> Cplx
-> Double
-> Cplx
-> Int
-> Double
-> [Cplx]
-> [Cplx]
-> Double
-> (Cplx, Int, [Cplx], [Cplx])
rj_ Cplx
x Cplx
y Cplx
z Cplx
p Cplx
a0 Double
maxmagns Cplx
delta Int
1 Double
1 [] [] Double
err
          f' :: Cplx
f' = forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
f
      in
      let x' :: Cplx
x' = (Cplx
a0 forall a. Num a => a -> a -> a
- Cplx
x) forall a. Fractional a => a -> a -> a
/ Cplx
f' forall a. Fractional a => a -> a -> a
/ Cplx
a
          y' :: Cplx
y' = (Cplx
a0 forall a. Num a => a -> a -> a
- Cplx
y) forall a. Fractional a => a -> a -> a
/ Cplx
f' forall a. Fractional a => a -> a -> a
/ Cplx
a
          z' :: Cplx
z' = (Cplx
a0 forall a. Num a => a -> a -> a
- Cplx
z) forall a. Fractional a => a -> a -> a
/ Cplx
f' forall a. Fractional a => a -> a -> a
/ Cplx
a
          p' :: Cplx
p' = -(Cplx
x'forall a. Num a => a -> a -> a
+Cplx
y'forall a. Num a => a -> a -> a
+Cplx
z') forall a. Fractional a => a -> a -> a
/ Cplx
2
          e2 :: Cplx
e2 = Cplx
x'forall a. Num a => a -> a -> a
*Cplx
y' forall a. Num a => a -> a -> a
+ Cplx
x'forall a. Num a => a -> a -> a
*Cplx
z' forall a. Num a => a -> a -> a
+ Cplx
y'forall a. Num a => a -> a -> a
*Cplx
z' forall a. Num a => a -> a -> a
- Cplx
3forall a. Num a => a -> a -> a
*Cplx
p'forall a. Num a => a -> a -> a
*Cplx
p'
          e3 :: Cplx
e3 = Cplx
x'forall a. Num a => a -> a -> a
*Cplx
y'forall a. Num a => a -> a -> a
*Cplx
z' forall a. Num a => a -> a -> a
+ Cplx
2forall a. Num a => a -> a -> a
*Cplx
e2forall a. Num a => a -> a -> a
*Cplx
p' forall a. Num a => a -> a -> a
+ Cplx
4forall a. Num a => a -> a -> a
*Cplx
p'forall a. Num a => a -> a -> a
*Cplx
p'forall a. Num a => a -> a -> a
*Cplx
p'
          e4 :: Cplx
e4 = Cplx
p'forall a. Num a => a -> a -> a
*(Cplx
2forall a. Num a => a -> a -> a
*Cplx
x'forall a. Num a => a -> a -> a
*Cplx
y'forall a. Num a => a -> a -> a
*Cplx
z' forall a. Num a => a -> a -> a
+ Cplx
e2forall a. Num a => a -> a -> a
*Cplx
p' forall a. Num a => a -> a -> a
+ Cplx
3forall a. Num a => a -> a -> a
*Cplx
p'forall a. Num a => a -> a -> a
*Cplx
p'forall a. Num a => a -> a -> a
*Cplx
p')
          e5 :: Cplx
e5 = Cplx
x'forall a. Num a => a -> a -> a
*Cplx
y'forall a. Num a => a -> a -> a
*Cplx
z'forall a. Num a => a -> a -> a
*Cplx
p'forall a. Num a => a -> a -> a
*Cplx
p'
          h :: [Cplx]
h = forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith (\Cplx
u Cplx
v -> Cplx -> Cplx
atanx_over_x(forall a. Floating a => a -> a
sqrt Cplx
u) forall a. Fractional a => a -> a -> a
/ Cplx
v) [Cplx]
e [Cplx]
d
      in
      (Cplx
1 forall a. Num a => a -> a -> a
- Cplx
3forall a. Num a => a -> a -> a
*Cplx
e2forall a. Fractional a => a -> a -> a
/Cplx
14 forall a. Num a => a -> a -> a
+ Cplx
e3forall a. Fractional a => a -> a -> a
/Cplx
6 forall a. Num a => a -> a -> a
+ Cplx
9forall a. Num a => a -> a -> a
*Cplx
e2forall a. Num a => a -> a -> a
*Cplx
e2forall a. Fractional a => a -> a -> a
/Cplx
88 forall a. Num a => a -> a -> a
- Cplx
3forall a. Num a => a -> a -> a
*Cplx
e4forall a. Fractional a => a -> a -> a
/Cplx
22 forall a. Num a => a -> a -> a
- Cplx
9forall a. Num a => a -> a -> a
*Cplx
e2forall a. Num a => a -> a -> a
*Cplx
e3forall a. Fractional a => a -> a -> a
/Cplx
52 forall a. Num a => a -> a -> a
+ Cplx
3forall a. Num a => a -> a -> a
*Cplx
e5forall a. Fractional a => a -> a -> a
/Cplx
26) forall a. Fractional a => a -> a -> a
/
        Cplx
f' forall a. Fractional a => a -> a -> a
/ Cplx
a forall a. Fractional a => a -> a -> a
/ forall a. Floating a => a -> a
sqrt Cplx
a forall a. Num a => a -> a -> a
+ Cplx
6 forall a. Num a => a -> a -> a
* forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum [Cplx]
h
  where
    zeros :: Int
zeros = forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum (forall a b. (a -> b) -> [a] -> [b]
map (\Cplx
u -> forall a. Enum a => a -> Int
fromEnum (Cplx
u forall a. Eq a => a -> a -> Bool
== Cplx
0)) [Cplx
x,Cplx
y,Cplx
z,Cplx
p])
    atanx_over_x :: Cplx -> Cplx
atanx_over_x Cplx
w = if Cplx
w forall a. Eq a => a -> a -> Bool
== Cplx
0 then Cplx
1 else Cplx -> Cplx
atanC Cplx
w forall a. Fractional a => a -> a -> a
/ Cplx
w

-- | Carlson integral RJ.
carlsonRJ ::
     Cplx   -- ^ first variable
  -> Cplx   -- ^ second variable 
  -> Cplx   -- ^ third variable 
  -> Cplx   -- ^ fourth variable 
  -> Cplx
carlsonRJ :: Cplx -> Cplx -> Cplx -> Cplx -> Cplx
carlsonRJ = Double -> Cplx -> Cplx -> Cplx -> Cplx -> Cplx
carlsonRJ' Double
1e-15


rc_ :: Cplx -> Cplx -> Cplx -> Double -> Int -> Double -> (Cplx, Int)
rc_ :: Cplx -> Cplx -> Cplx -> Double -> Int -> Double -> (Cplx, Int)
rc_ Cplx
x Cplx
y Cplx
a Double
magn Int
f Double
err =
  let q :: Double
q = (Double
1forall a. Fractional a => a -> a -> a
/Double
3forall a. Fractional a => a -> a -> a
/Double
err)forall a. Floating a => a -> a -> a
**(Double
1forall a. Fractional a => a -> a -> a
/Double
8) forall a. Num a => a -> a -> a
* Double
magn forall a. Fractional a => a -> a -> a
/ forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
f in
  if forall a. RealFloat a => Complex a -> a
magnitude Cplx
a forall a. Ord a => a -> a -> Bool
> Double
q
    then (Cplx
a, Int
f)
    else
      let lambda :: Cplx
lambda = Cplx
2 forall a. Num a => a -> a -> a
* forall a. Floating a => a -> a
sqrt Cplx
x forall a. Num a => a -> a -> a
* forall a. Floating a => a -> a
sqrt Cplx
y forall a. Num a => a -> a -> a
+ Cplx
y
          a' :: Cplx
a' = (Cplx
a forall a. Num a => a -> a -> a
+ Cplx
lambda) forall a. Fractional a => a -> a -> a
/ Cplx
4
          x' :: Cplx
x' = (Cplx
x forall a. Num a => a -> a -> a
+ Cplx
lambda) forall a. Fractional a => a -> a -> a
/ Cplx
4
          y' :: Cplx
y' = (Cplx
y forall a. Num a => a -> a -> a
+ Cplx
lambda) forall a. Fractional a => a -> a -> a
/ Cplx
4
      in
      Cplx -> Cplx -> Cplx -> Double -> Int -> Double -> (Cplx, Int)
rc_ Cplx
x' Cplx
y' Cplx
a' Double
magn (Int
4forall a. Num a => a -> a -> a
*Int
f) Double
err

-- | Carlson integral RC.
carlsonRC' ::
     Double -- ^ bound on the relative error
  -> Cplx   -- ^ first variable
  -> Cplx   -- ^ second variable 
  -> Cplx
carlsonRC' :: Double -> Cplx -> Cplx -> Cplx
carlsonRC' Double
err Cplx
x Cplx
y =
  if Cplx
y forall a. Eq a => a -> a -> Bool
== Cplx
0
    then forall a. HasCallStack => [Char] -> a
error [Char]
"y cannot be 0"
    else
      let a0 :: Cplx
a0 = (Cplx
x forall a. Num a => a -> a -> a
+ Cplx
y forall a. Num a => a -> a -> a
+ Cplx
y) forall a. Fractional a => a -> a -> a
/ Cplx
3
          magn :: Double
magn = forall a. RealFloat a => Complex a -> a
magnitude(Cplx
a0forall a. Num a => a -> a -> a
-Cplx
x)
      in
      let (Cplx
a, Int
f) = Cplx -> Cplx -> Cplx -> Double -> Int -> Double -> (Cplx, Int)
rc_ Cplx
x Cplx
y Cplx
a0 Double
magn Int
1 Double
err
          f' :: Cplx
f' = forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
f
      in
      let s :: Cplx
s = (Cplx
y forall a. Num a => a -> a -> a
- Cplx
a0) forall a. Fractional a => a -> a -> a
/ Cplx
f' forall a. Fractional a => a -> a -> a
/ Cplx
a in
      (Cplx
1 forall a. Num a => a -> a -> a
+ Cplx
3forall a. Num a => a -> a -> a
*Cplx
sforall a. Num a => a -> a -> a
*Cplx
sforall a. Fractional a => a -> a -> a
/Cplx
10 forall a. Num a => a -> a -> a
+ Cplx
sforall a. Num a => a -> a -> a
*Cplx
sforall a. Num a => a -> a -> a
*Cplx
sforall a. Fractional a => a -> a -> a
/Cplx
7 forall a. Num a => a -> a -> a
+ Cplx
3forall a. Num a => a -> a -> a
*Cplx
sforall a. Num a => a -> a -> a
*Cplx
sforall a. Num a => a -> a -> a
*Cplx
sforall a. Num a => a -> a -> a
*Cplx
sforall a. Fractional a => a -> a -> a
/Cplx
8 forall a. Num a => a -> a -> a
+ Cplx
9forall a. Num a => a -> a -> a
*Cplx
sforall a. Num a => a -> a -> a
*Cplx
sforall a. Num a => a -> a -> a
*Cplx
sforall a. Num a => a -> a -> a
*Cplx
sforall a. Num a => a -> a -> a
*Cplx
sforall a. Fractional a => a -> a -> a
/Cplx
22 forall a. Num a => a -> a -> a
+
        Cplx
159forall a. Num a => a -> a -> a
*Cplx
sforall a. Num a => a -> a -> a
*Cplx
sforall a. Num a => a -> a -> a
*Cplx
sforall a. Num a => a -> a -> a
*Cplx
sforall a. Num a => a -> a -> a
*Cplx
sforall a. Num a => a -> a -> a
*Cplx
sforall a. Fractional a => a -> a -> a
/Cplx
208 forall a. Num a => a -> a -> a
+ Cplx
9forall a. Num a => a -> a -> a
*Cplx
sforall a. Num a => a -> a -> a
*Cplx
sforall a. Num a => a -> a -> a
*Cplx
sforall a. Num a => a -> a -> a
*Cplx
sforall a. Num a => a -> a -> a
*Cplx
sforall a. Num a => a -> a -> a
*Cplx
sforall a. Num a => a -> a -> a
*Cplx
sforall a. Fractional a => a -> a -> a
/Cplx
8) forall a. Fractional a => a -> a -> a
/ forall a. Floating a => a -> a
sqrt Cplx
a

-- | Carlson integral RC.
carlsonRC ::
     Cplx   -- ^ first variable
  -> Cplx   -- ^ second variable 
  -> Cplx
carlsonRC :: Cplx -> Cplx -> Cplx
carlsonRC = Double -> Cplx -> Cplx -> Cplx
carlsonRC' Double
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' :: Double -> Cplx -> Cplx -> Cplx -> Cplx
carlsonRG' Double
err Cplx
x Cplx
y Cplx
z =
  if Int
zeros forall a. Ord a => a -> a -> Bool
> Int
1
    then forall a. Floating a => a -> a
sqrt(Cplx
xforall a. Num a => a -> a -> a
+Cplx
yforall a. Num a => a -> a -> a
+Cplx
z) forall a. Fractional a => a -> a -> a
/ Cplx
2
    else
      if Cplx
z forall a. Eq a => a -> a -> Bool
== Cplx
0
        then Double -> Cplx -> Cplx -> Cplx -> Cplx
carlsonRG' Double
err Cplx
z Cplx
x Cplx
y
        else
          (Cplx
z forall a. Num a => a -> a -> a
* Double -> Cplx -> Cplx -> Cplx -> Cplx
carlsonRF' Double
err Cplx
x Cplx
y Cplx
z forall a. Num a => a -> a -> a
-
            (Cplx
xforall a. Num a => a -> a -> a
-Cplx
z) forall a. Num a => a -> a -> a
* (Cplx
yforall a. Num a => a -> a -> a
-Cplx
z) forall a. Num a => a -> a -> a
* Double -> Cplx -> Cplx -> Cplx -> Cplx
carlsonRD' Double
err Cplx
x Cplx
y Cplx
z forall a. Fractional a => a -> a -> a
/ Cplx
3 forall a. Num a => a -> a -> a
+
            forall a. Floating a => a -> a
sqrt Cplx
x forall a. Num a => a -> a -> a
* forall a. Floating a => a -> a
sqrt Cplx
y forall a. Fractional a => a -> a -> a
/ forall a. Floating a => a -> a
sqrt Cplx
z) forall a. Fractional a => a -> a -> a
/ Cplx
2
  where
    zeros :: Int
zeros = forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum (forall a b. (a -> b) -> [a] -> [b]
map (\Cplx
u -> forall a. Enum a => a -> Int
fromEnum (Cplx
u forall a. Eq a => a -> a -> Bool
== Cplx
0)) [Cplx
x,Cplx
y,Cplx
z])

-- | Carlson integral RG.
carlsonRG ::
     Cplx   -- ^ first variable
  -> Cplx   -- ^ second variable 
  -> Cplx   -- ^ third variable 
  -> Cplx
carlsonRG :: Cplx -> Cplx -> Cplx -> Cplx
carlsonRG = Double -> Cplx -> Cplx -> Cplx -> Cplx
carlsonRG' Double
1e-15