{-# LANGUAGE ForeignFunctionInterface #-}

-- SPDX-License-Identifier: MPL-2.0

{- | Real special functions with checked domains. 'SFResult' reports an
/estimated/ absolute truncation/rounding error, never an interval certificate.
Series and quadrature have finite work budgets and reject nonconvergence.
Complex hypergeometric functions use principal powers and their power-series
disk; no unspecified continuation across a branch cut is performed.
Equations and normalization: https://dlmf.nist.gov/ (chapters 5--33).
The C99 libm gamma/error/Bessel entry points are used as platform elementary
kernels. They do not change the process error handler or floating-point mode.
-}
module HMath.Special (
    SFResult (..),
    gamma,
    logGamma,
    gammaP,
    gammaQ,
    incompleteGamma,
    beta,
    logBeta,
    betaI,
    incompleteBeta,
    digamma,
    trigamma,
    polygamma,
    erf,
    erfc,
    inverseErf,
    erfcx,
    normalCDF,
    normalQuantile,
    dawson,
    exponentialIntegral,
    logarithmicIntegral,
    sineIntegral,
    cosineIntegral,
    fresnel,
    ellipticF,
    ellipticE,
    ellipticPi,
    jacobiElliptic,
    besselJ,
    besselY,
    besselI,
    besselK,
    sphericalBesselJ,
    sphericalBesselY,
    hankel,
    airy,
    struve,
    hypergeometric,
    hypergeometric1F1,
    hypergeometric2F1,
    complex1F1,
    legendre,
    associatedLegendre,
    sphericalHarmonic,
    riemannZeta,
    hurwitzZeta,
    polylogarithm,
    dilogarithm,
    lambertW,
    debye,
    clausen,
    fermiDirac,
    boseEinstein,
    coulombWave,
    gegenbauer,
    parabolicCylinder,
    mathieu,
    bernoulliPolynomial,
    eulerPolynomial,
    wigner3j,
    wigner6j,
    wigner9j,
    clebschGordan,
    scaledBesselI,
    scaledBesselK,
) where

import Data.Complex hiding (phase)
import qualified Data.Foldable as Foldable (foldl')
import Foreign.C.Types
import HMath.Error
import HMath.Matrix (fromListChecked, toList)
import HMath.Matrix.LAPACK (eigenvalues, eigenvectors, symmetricEigen)

foreign import ccall unsafe "tgamma" cGamma :: CDouble -> CDouble
foreign import ccall unsafe "lgamma" cLogGamma :: CDouble -> CDouble
foreign import ccall unsafe "erf" cErf :: CDouble -> CDouble
foreign import ccall unsafe "erfc" cErfc :: CDouble -> CDouble
foreign import ccall unsafe "jn" cJn :: CInt -> CDouble -> CDouble
foreign import ccall unsafe "yn" cYn :: CInt -> CDouble -> CDouble

-- | Finite value and estimated absolute numerical error. The estimate includes truncation/discretization and a small rounding allowance; it is not a certified enclosure.
data SFResult = SFResult
    { SFResult -> Double
sfValue :: Double
    -- ^ Returned special-function numerical estimate.
    , SFResult -> Double
sfError :: Double
    -- ^ Absolute numerical error estimate accompanying sfValue.
    }
    deriving (SFResult -> SFResult -> Bool
(SFResult -> SFResult -> Bool)
-> (SFResult -> SFResult -> Bool) -> Eq SFResult
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
$c== :: SFResult -> SFResult -> Bool
== :: SFResult -> SFResult -> Bool
$c/= :: SFResult -> SFResult -> Bool
/= :: SFResult -> SFResult -> Bool
Eq, Int -> SFResult -> ShowS
[SFResult] -> ShowS
SFResult -> String
(Int -> SFResult -> ShowS)
-> (SFResult -> String) -> ([SFResult] -> ShowS) -> Show SFResult
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
$cshowsPrec :: Int -> SFResult -> ShowS
showsPrec :: Int -> SFResult -> ShowS
$cshow :: SFResult -> String
show :: SFResult -> String
$cshowList :: [SFResult] -> ShowS
showList :: [SFResult] -> ShowS
Show)

finite :: Double -> Bool
finite :: Double -> Bool
finite Double
x = Bool -> Bool
not (Double -> Bool
forall a. RealFloat a => a -> Bool
isNaN Double
x Bool -> Bool -> Bool
|| Double -> Bool
forall a. RealFloat a => a -> Bool
isInfinite Double
x)
checked :: Double -> Either MathError Double
checked :: Double -> Either MathError Double
checked Double
x = if Double -> Bool
finite Double
x then Double -> Either MathError Double
forall a b. b -> Either a b
Right Double
x else MathError -> Either MathError Double
forall a b. a -> Either a b
Left MathError
NonFiniteResult
estimate :: Double -> Double -> Either MathError SFResult
estimate :: Double -> Double -> Either MathError SFResult
estimate Double
x Double
e = if Double -> Bool
finite Double
x Bool -> Bool -> Bool
&& Double -> Bool
finite Double
e then SFResult -> Either MathError SFResult
forall a b. b -> Either a b
Right (Double -> Double -> SFResult
SFResult Double
x (Double -> Double
forall a. Num a => a -> a
abs Double
e Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
8e-16 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double
forall a. Num a => a -> a
abs Double
x)) else MathError -> Either MathError SFResult
forall a b. a -> Either a b
Left MathError
NonFiniteResult
require :: Bool -> String -> Either MathError ()
require :: Bool -> String -> Either MathError ()
require Bool
p String
msg = if Bool
p then () -> Either MathError ()
forall a b. b -> Either a b
Right () else MathError -> Either MathError ()
forall a b. a -> Either a b
Left (String -> MathError
DomainError String
msg)

-- | Real gamma using C99 libm; rejects nonfinite arguments and nonpositive integer poles. Overflow is reported as NonFiniteResult.
gamma :: Double -> Either MathError Double
gamma :: Double -> Either MathError Double
gamma Double
x = do
    Bool -> String -> Either MathError ()
require (Double -> Bool
finite Double
x Bool -> Bool -> Bool
&& (Double
x Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double
0 Bool -> Bool -> Bool
|| Double
x Double -> Double -> Bool
forall a. Eq a => a -> a -> Bool
/= Integer -> Double
forall a. Num a => Integer -> a
fromInteger (Double -> Integer
forall b. Integral b => Double -> b
forall a b. (RealFrac a, Integral b) => a -> b
round Double
x))) String
"gamma pole or non-finite argument"
    Double -> Either MathError Double
checked (CDouble -> Double
forall a b. (Real a, Fractional b) => a -> b
realToFrac (CDouble -> CDouble
cGamma (Double -> CDouble
forall a b. (Real a, Fractional b) => a -> b
realToFrac Double
x)))

-- | Natural logarithm of gamma for strictly positive finite arguments; avoids overflow of gamma itself.
logGamma :: Double -> Either MathError Double
logGamma :: Double -> Either MathError Double
logGamma Double
x = do
    Bool -> String -> Either MathError ()
require (Double -> Bool
finite Double
x Bool -> Bool -> Bool
&& Double
x Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double
0) String
"logGamma requires x>0"; Double -> Either MathError Double
checked (CDouble -> Double
forall a b. (Real a, Fractional b) => a -> b
realToFrac (CDouble -> CDouble
cLogGamma (Double -> CDouble
forall a b. (Real a, Fractional b) => a -> b
realToFrac Double
x)))

-- | Log beta for strictly positive finite parameters.
logBeta :: Double -> Double -> Either MathError Double
logBeta :: Double -> Double -> Either MathError Double
logBeta Double
a Double
b = do
    Bool -> String -> Either MathError ()
require (Double
a Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double
0 Bool -> Bool -> Bool
&& Double
b Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double
0 Bool -> Bool -> Bool
&& Double -> Bool
finite (Double
a Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
b)) String
"beta positive parameters required"
    x <- Double -> Either MathError Double
logGamma Double
a
    y <- logGamma b
    z <- logGamma (a + b)
    checked (x + y - z)

-- | Beta function for strictly positive parameters, evaluated through log gamma.
beta :: Double -> Double -> Either MathError Double
beta :: Double -> Double -> Either MathError Double
beta Double
a Double
b = Double -> Double
forall a. Floating a => a -> a
exp (Double -> Double)
-> Either MathError Double -> Either MathError Double
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
<$> Double -> Double -> Either MathError Double
logBeta Double
a Double
b

{- | Regularized lower incomplete gamma. Series below a+1, Lentz continued
fraction above; upper tail is evaluated directly instead of subtracting 1.
-}
gammaP :: Double -> Double -> Either MathError Double
gammaP :: Double -> Double -> Either MathError Double
gammaP Double
a Double
x = do (p, _) <- Double -> Double -> Either MathError (Double, Double)
gammaPQ Double
a Double
x; pure p

-- | Regularized upper gamma Q(a,x); the upper-tail continued fraction is evaluated directly when x>=a+1.
gammaQ :: Double -> Double -> Either MathError Double
gammaQ :: Double -> Double -> Either MathError Double
gammaQ Double
a Double
x = do (_, q) <- Double -> Double -> Either MathError (Double, Double)
gammaPQ Double
a Double
x; pure q

gammaPQ :: Double -> Double -> Either MathError (Double, Double)
gammaPQ :: Double -> Double -> Either MathError (Double, Double)
gammaPQ Double
a Double
x = do
    Bool -> String -> Either MathError ()
require (Double -> Bool
finite Double
a Bool -> Bool -> Bool
&& Double -> Bool
finite Double
x Bool -> Bool -> Bool
&& Double
a Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double
0 Bool -> Bool -> Bool
&& Double
x Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
>= Double
0) String
"incomplete gamma needs a>0, x>=0"
    lg <- Double -> Either MathError Double
logGamma Double
a
    if x == 0
        then Right (0, 1)
        else
            if x < a + 1
                then do
                    s <- series 10000 (1 / a) (\Int
n Double
t -> Double
t Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
x Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
a Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
n))
                    p <- checked (exp (a * log x - x - lg) * sfValue s)
                    pure (max 0 (min 1 p), max 0 (1 - p))
                else do
                    let b = Double
x Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
a
                        d = Double
1 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
b
                        go :: Int -> Double -> Double -> Double -> Double -> Either MathError Double
                        go Int
n Double
bb Double
cc Double
dd Double
h
                            | Int
n Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
> Int
10000 = MathError -> Either MathError Double
forall a b. a -> Either a b
Left (String -> Int -> MathError
NonConvergence String
"incomplete gamma fraction" Int
10000)
                            | Bool
otherwise =
                                let an :: Double
an = Double -> Double
forall a. Num a => a -> a
negate (Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
n Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
n Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
a))
                                    bn :: Double
bn = Double
bb Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
2
                                    dn :: Double
dn = Double -> Double
forall {a}. (Ord a, Fractional a) => a -> a
nonzero (Double
an Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
dd Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
bn)
                                    cn :: Double
cn = Double -> Double
forall {a}. (Ord a, Fractional a) => a -> a
nonzero (Double
bn Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
an Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
cc)
                                    di :: Double
di = Double
1 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
dn
                                    change :: Double
change = Double
di Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
cn
                                    hn :: Double
hn = Double
h Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
change
                                 in if Double -> Double
forall a. Num a => a -> a
abs (Double
change Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
1) Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
2e-15 then Double -> Either MathError Double
checked Double
hn else Int
-> Double -> Double -> Double -> Double -> Either MathError Double
go (Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1) Double
bn Double
cn Double
di Double
hn
                    h <- go 1 b 1e300 d d
                    q <- checked (exp (a * log x - x - lg) * h)
                    pure (max 0 (1 - q), max 0 (min 1 q))
  where
    nonzero :: a -> a
nonzero a
z = if a -> a
forall a. Num a => a -> a
abs a
z a -> a -> Bool
forall a. Ord a => a -> a -> Bool
< a
1e-300 then a
1e-300 else a
z

-- | Unregularized lower gamma gamma(a,x)=Gamma(a)*P(a,x). Overflow is an error.
incompleteGamma :: Double -> Double -> Either MathError Double
incompleteGamma :: Double -> Double -> Either MathError Double
incompleteGamma Double
a Double
x = do p <- Double -> Double -> Either MathError Double
gammaP Double
a Double
x; g <- gamma a; checked (p * g)

-- | Regularized lower incomplete beta I_x(a,b), positive parameters and x in [0,1], via reflected Lentz continued fraction.
betaI :: Double -> Double -> Double -> Either MathError Double
betaI :: Double -> Double -> Double -> Either MathError Double
betaI Double
a Double
b Double
x = do
    Bool -> String -> Either MathError ()
require
        (Double -> Bool
finite (Double
a Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
b Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
x) Bool -> Bool -> Bool
&& Double
a Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double
0 Bool -> Bool -> Bool
&& Double
b Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double
0 Bool -> Bool -> Bool
&& Double
x Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
>= Double
0 Bool -> Bool -> Bool
&& Double
x Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
<= Double
1)
        String
"regularized beta requires positive parameters and x in [0,1]"
    if Double
x Double -> Double -> Bool
forall a. Eq a => a -> a -> Bool
== Double
0 Bool -> Bool -> Bool
|| Double
x Double -> Double -> Bool
forall a. Eq a => a -> a -> Bool
== Double
1
        then Double -> Either MathError Double
forall a b. b -> Either a b
Right Double
x
        else
            if Double
x Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> (Double
a Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
1) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
a Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
b Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
2)
                then (Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
-) (Double -> Double)
-> Either MathError Double -> Either MathError Double
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
<$> Double -> Double -> Double -> Either MathError Double
betaI Double
b Double
a (Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
x)
                else do
                    lb <- Double -> Double -> Either MathError Double
logBeta Double
a Double
b
                    let qab = Double
a Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
b
                        qap = Double
a Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
1
                        qam = Double
a Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
1
                        nz a
v = if a -> a
forall a. Num a => a -> a
abs a
v a -> a -> Bool
forall a. Ord a => a -> a -> Bool
< a
1e-300 then a
1e-300 else a
v
                        go :: Int -> Double -> Double -> Double -> Either MathError Double
                        go Int
m Double
c Double
d Double
h
                            | Int
m Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
> Int
10000 = MathError -> Either MathError Double
forall a b. a -> Either a b
Left (String -> Int -> MathError
NonConvergence String
"beta continued fraction" Int
10000)
                            | Bool
otherwise =
                                let fm :: Double
fm = Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
m
                                    m2 :: Double
m2 = Double
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
fm
                                    aa :: Double
aa = Double
fm Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
b Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
fm) Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
x Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ ((Double
qam Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
m2) Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
a Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
m2))
                                    d1 :: Double
d1 = Double
1 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double -> Double
forall {a}. (Ord a, Fractional a) => a -> a
nz (Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
aa Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
d)
                                    c1 :: Double
c1 = Double -> Double
forall {a}. (Ord a, Fractional a) => a -> a
nz (Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
aa Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
c)
                                    h1 :: Double
h1 = Double
h Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
d1 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
c1
                                    ab :: Double
ab = -(Double
a Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
fm) Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
qab Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
fm) Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
x Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ ((Double
a Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
m2) Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
qap Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
m2))
                                    d2 :: Double
d2 = Double
1 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double -> Double
forall {a}. (Ord a, Fractional a) => a -> a
nz (Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
ab Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
d1)
                                    c2 :: Double
c2 = Double -> Double
forall {a}. (Ord a, Fractional a) => a -> a
nz (Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
ab Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
c1)
                                    change :: Double
change = Double
d2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
c2
                                    h2 :: Double
h2 = Double
h1 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
change
                                 in if Double -> Double
forall a. Num a => a -> a
abs (Double
change Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
1) Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
2e-15 then Double -> Either MathError Double
checked Double
h2 else Int -> Double -> Double -> Double -> Either MathError Double
go (Int
m Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1) Double
c2 Double
d2 Double
h2
                    cf <- go 1 1 (1 / (1 - qab * x / qap)) (1 / (1 - qab * x / qap))
                    checked (max 0 (min 1 (exp (a * log x + b * log (1 - x) - lb) * cf / a)))

-- | Unregularized lower incomplete beta B_x(a,b).
incompleteBeta :: Double -> Double -> Double -> Either MathError Double
incompleteBeta :: Double -> Double -> Double -> Either MathError Double
incompleteBeta Double
a Double
b Double
x = do p <- Double -> Double -> Double -> Either MathError Double
betaI Double
a Double
b Double
x; v <- beta a b; checked (p * v)

series :: Int -> Double -> (Int -> Double -> Double) -> Either MathError SFResult
series :: Int
-> Double -> (Int -> Double -> Double) -> Either MathError SFResult
series Int
budget Double
first Int -> Double -> Double
step = Int -> Double -> Double -> Double -> Either MathError SFResult
go Int
1 Double
first Double
first (Double -> Double
forall a. Num a => a -> a
abs Double
first)
  where
    go :: Int -> Double -> Double -> Double -> Either MathError SFResult
go Int
n Double
t Double
s Double
absoluteSum
        | Int
n Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
> Int
budget = MathError -> Either MathError SFResult
forall a b. a -> Either a b
Left (String -> Int -> MathError
NonConvergence String
"special-function series" Int
budget)
        | Bool
otherwise =
            let u :: Double
u = Int -> Double -> Double
step Int
n Double
t
                q :: Double
q = Double
s Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
u
                absoluteNext :: Double
absoluteNext = Double
absoluteSum Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double -> Double
forall a. Num a => a -> a
abs Double
u
                ratio :: Double
ratio = if Double
t Double -> Double -> Bool
forall a. Eq a => a -> a -> Bool
== Double
0 then Double
0 else Double -> Double
forall a. Num a => a -> a
abs (Double
u Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
t)
                tailEstimate :: Double
tailEstimate = if Double
ratio Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
1 then Double -> Double
forall a. Num a => a -> a
abs Double
u Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double -> Double -> Double
forall a. Ord a => a -> a -> a
max Double
1e-16 (Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
ratio) else Double -> Double
forall a. Num a => a -> a
abs Double
u
             in if Bool -> Bool
not (Double -> Bool
finite Double
q Bool -> Bool -> Bool
&& Double -> Bool
finite Double
absoluteNext)
                    then MathError -> Either MathError SFResult
forall a b. a -> Either a b
Left MathError
NonFiniteResult
                    else
                        if Double
u Double -> Double -> Bool
forall a. Eq a => a -> a -> Bool
== Double
0 Bool -> Bool -> Bool
|| Double -> Double
forall a. Num a => a -> a
abs Double
u Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
2e-15 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double
forall a. Num a => a -> a
abs Double
q
                            then Double -> Double -> Either MathError SFResult
estimate Double
q (Double
tailEstimate Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
8e-16 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
absoluteNext)
                            else Int -> Double -> Double -> Double -> Either MathError SFResult
go (Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1) Double
u Double
q Double
absoluteNext

-- Adaptive Simpson with endpoint-finite integrands. Error is Richardson's
-- difference, not a bound. Budget counts binary subdivision depth (20).
integrate :: (Double -> Double) -> Double -> Double -> Either MathError SFResult
integrate :: (Double -> Double) -> Double -> Double -> Either MathError SFResult
integrate Double -> Double
f Double
a Double
b
    | Double
a Double -> Double -> Bool
forall a. Eq a => a -> a -> Bool
== Double
b = SFResult -> Either MathError SFResult
forall a b. b -> Either a b
Right (Double -> Double -> SFResult
SFResult Double
0 Double
0)
    | Bool
otherwise = Int
-> Double
-> Double
-> Double
-> Double
-> Double
-> Either MathError SFResult
go Int
20 Double
a Double
b (Double -> Double
f Double
a) (Double -> Double
f ((Double
a Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
b) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
2)) (Double -> Double
f Double
b)
  where
    sim :: a -> a -> a -> a -> a -> a
sim a
x a
y a
u a
v a
w = (a
y a -> a -> a
forall a. Num a => a -> a -> a
- a
x) a -> a -> a
forall a. Num a => a -> a -> a
* (a
u a -> a -> a
forall a. Num a => a -> a -> a
+ a
4 a -> a -> a
forall a. Num a => a -> a -> a
* a
v a -> a -> a
forall a. Num a => a -> a -> a
+ a
w) a -> a -> a
forall a. Fractional a => a -> a -> a
/ a
6
    go :: Int -> Double -> Double -> Double -> Double -> Double -> Either MathError SFResult
    go :: Int
-> Double
-> Double
-> Double
-> Double
-> Double
-> Either MathError SFResult
go Int
depth Double
x Double
y Double
u Double
v Double
w
        | (Double -> Bool) -> [Double] -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
any (Bool -> Bool
not (Bool -> Bool) -> (Double -> Bool) -> Double -> Bool
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Double -> Bool
finite) [Double
u, Double
v, Double
w] = MathError -> Either MathError SFResult
forall a b. a -> Either a b
Left MathError
NonFiniteResult
        | Bool
otherwise =
            let m :: Double
m = (Double
x Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
y) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
2
                l :: Double
l = Double -> Double
f ((Double
x Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
m) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
2)
                r :: Double
r = Double -> Double
f ((Double
m Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
y) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
2)
                old :: Double
old = Double -> Double -> Double -> Double -> Double -> Double
forall {a}. Fractional a => a -> a -> a -> a -> a -> a
sim Double
x Double
y Double
u Double
v Double
w
                new :: Double
new = Double -> Double -> Double -> Double -> Double -> Double
forall {a}. Fractional a => a -> a -> a -> a -> a -> a
sim Double
x Double
m Double
u Double
l Double
v Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double -> Double -> Double -> Double -> Double -> Double
forall {a}. Fractional a => a -> a -> a -> a -> a -> a
sim Double
m Double
y Double
v Double
r Double
w
                e :: Double
e = Double -> Double
forall a. Num a => a -> a
abs (Double
new Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
old) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
15
             in if Bool -> Bool
not (Double -> Bool
finite Double
new)
                    then MathError -> Either MathError SFResult
forall a b. a -> Either a b
Left MathError
NonFiniteResult
                    else
                        if Double
e Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
1e-11 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double -> Double
forall a. Ord a => a -> a -> a
max Double
1 (Double -> Double
forall a. Num a => a -> a
abs Double
new)
                            then Double -> Double -> Either MathError SFResult
estimate (Double
new Double -> Double -> Double
forall a. Num a => a -> a -> a
+ (Double
new Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
old) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
15) Double
e
                            else
                                if Int
depth Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
<= Int
0
                                    then MathError -> Either MathError SFResult
forall a b. a -> Either a b
Left (String -> Int -> MathError
NonConvergence String
"special quadrature" Int
20)
                                    else do
                                        p <- Int
-> Double
-> Double
-> Double
-> Double
-> Double
-> Either MathError SFResult
go (Int
depth Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1) Double
x Double
m Double
u Double
l Double
v
                                        q <- go (depth - 1) m y v r w
                                        estimate (sfValue p + sfValue q) (sfError p + sfError q)

-- | Log-gamma derivative for x>0, using recurrence to x>=12 and an inverse-power asymptotic expansion.
digamma :: Double -> Either MathError Double
digamma :: Double -> Either MathError Double
digamma Double
x = do
    Bool -> String -> Either MathError ()
require (Double -> Bool
finite Double
x Bool -> Bool -> Bool
&& Double
x Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double
0) String
"digamma x>0 required"
    let (Double
y, Double
s) = Double -> Double -> (Double, Double)
forall {b}. (Ord b, Fractional b) => b -> b -> (b, b)
shift Double
x Double
0
        z :: Double
z = Double
1 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
y Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
y)
        v :: Double
v = Double -> Double
forall a. Floating a => a -> a
log Double
y Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
1 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
y) Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
z Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
1 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
12 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
z Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
1 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
120 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
z Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
1 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
252 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
z Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
240))) Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
s
     in Double -> Either MathError Double
checked Double
v
  where
    shift :: b -> b -> (b, b)
shift b
y b
s | b
y b -> b -> Bool
forall a. Ord a => a -> a -> Bool
< b
12 = b -> b -> (b, b)
shift (b
y b -> b -> b
forall a. Num a => a -> a -> a
+ b
1) (b
s b -> b -> b
forall a. Num a => a -> a -> a
- b
1 b -> b -> b
forall a. Fractional a => a -> a -> a
/ b
y) | Bool
otherwise = (b
y, b
s)

-- | First polygamma derivative for x>0, using the Hurwitz-zeta identity.
trigamma :: Double -> Either MathError Double
trigamma :: Double -> Either MathError Double
trigamma = Int -> Double -> Either MathError Double
polygamma Int
1

-- | Derivative order 0..32 of digamma for x>0; higher orders use signed factorial times Hurwitz zeta.
polygamma :: Int -> Double -> Either MathError Double
polygamma :: Int -> Double -> Either MathError Double
polygamma Int
n Double
x
    | Int
n Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
0 = Double -> Either MathError Double
digamma Double
x
    | Bool
otherwise = do
        Bool -> String -> Either MathError ()
require (Int
n Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
>= Int
1 Bool -> Bool -> Bool
&& Int
n Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
<= Int
32 Bool -> Bool -> Bool
&& Double -> Bool
finite Double
x Bool -> Bool -> Bool
&& Double
x Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double
0) String
"polygamma order 0..32 and x>0 required"
        z <- Double -> Double -> Either MathError SFResult
hurwitzZeta (Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
n Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
1) Double
x
        checked ((if odd n then 1 else -1) * fact n * sfValue z)

-- | Real error function from C99 libm; finite arguments required.
erf :: Double -> Either MathError Double
erf :: Double -> Either MathError Double
erf Double
x = do Bool -> String -> Either MathError ()
require (Double -> Bool
finite Double
x) String
"erf finite input required"; Double -> Either MathError Double
checked (CDouble -> Double
forall a b. (Real a, Fractional b) => a -> b
realToFrac (CDouble -> CDouble
cErf (Double -> CDouble
forall a b. (Real a, Fractional b) => a -> b
realToFrac Double
x)))

-- | Real complementary error function from C99 libm; upper tails avoid subtracting erf from one.
erfc :: Double -> Either MathError Double
erfc :: Double -> Either MathError Double
erfc Double
x = do Bool -> String -> Either MathError ()
require (Double -> Bool
finite Double
x) String
"erfc finite input required"; Double -> Either MathError Double
checked (CDouble -> Double
forall a b. (Real a, Fractional b) => a -> b
realToFrac (CDouble -> CDouble
cErfc (Double -> CDouble
forall a b. (Real a, Fractional b) => a -> b
realToFrac Double
x)))

-- | exp(x^2)*erfc(x), with a large-positive-x inverse-power asymptotic expansion to avoid overflow.
erfcx :: Double -> Either MathError SFResult
erfcx :: Double -> Either MathError SFResult
erfcx Double
x
    | Bool -> Bool
not (Double -> Bool
finite Double
x) = MathError -> Either MathError SFResult
forall a b. a -> Either a b
Left (String -> MathError
DomainError String
"erfcx finite input required")
    | Double
x Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double
26 = do
        s <- Int
-> Double -> (Int -> Double -> Double) -> Either MathError SFResult
series Int
100 Double
1 (\Int
n Double
t -> Double
t Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double
forall a. Num a => a -> a
negate (Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Int
2 Int -> Int -> Int
forall a. Num a => a -> a -> a
* Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1)) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
x Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
x))
        estimate (sfValue s / (sqrt pi * x)) (sfError s / (sqrt pi * x))
    | Bool
otherwise = do y <- Double -> Either MathError Double
erfc Double
x; estimate (exp (x * x) * y) (abs y * 2e-15)

-- | Standard-normal CDF through erfc, with stable small lower tails.
normalCDF :: Double -> Either MathError Double
normalCDF :: Double -> Either MathError Double
normalCDF Double
x = (Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
2) (Double -> Double)
-> Either MathError Double -> Either MathError Double
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
<$> Double -> Either MathError Double
erfc (-Double
x Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double -> Double
forall a. Floating a => a -> a
sqrt Double
2)

-- | Standard-normal inverse CDF for 0\<p\<1. Symmetry preserves tiny upper-tail probabilities before bounded bisection.
normalQuantile :: Double -> Either MathError Double
normalQuantile :: Double -> Either MathError Double
normalQuantile Double
p = do
    Bool -> String -> Either MathError ()
require (Double -> Bool
finite Double
p Bool -> Bool -> Bool
&& Double
p Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double
0 Bool -> Bool -> Bool
&& Double
p Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
1) String
"normal quantile requires 0<p<1"
    if Double
p Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double
0.5 then Double -> Double
forall a. Num a => a -> a
negate (Double -> Double)
-> Either MathError Double -> Either MathError Double
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
<$> Double -> Either MathError Double
normalQuantile (Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
p) else Int
-> Double
-> Double
-> Double
-> (Double -> Either MathError Double)
-> Either MathError Double
bisect Int
100 (-Double
40) Double
0 Double
p Double -> Either MathError Double
normalCDF

-- | Real inverse error function for -1\<x\<1, using the normal inverse identity.
inverseErf :: Double -> Either MathError Double
inverseErf :: Double -> Either MathError Double
inverseErf Double
x = do
    Bool -> String -> Either MathError ()
require (Double -> Bool
finite Double
x Bool -> Bool -> Bool
&& Double -> Double
forall a. Num a => a -> a
abs Double
x Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
1) String
"inverseErf requires -1<x<1"
    (Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double -> Double
forall a. Floating a => a -> a
sqrt Double
2) (Double -> Double)
-> Either MathError Double -> Either MathError Double
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
<$> Double -> Either MathError Double
normalQuantile ((Double
x Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
1) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
2)

bisect ::
    Int -> Double -> Double -> Double -> (Double -> Either MathError Double) -> Either MathError Double
bisect :: Int
-> Double
-> Double
-> Double
-> (Double -> Either MathError Double)
-> Either MathError Double
bisect Int
n Double
lo Double
hi Double
target Double -> Either MathError Double
f
    | Int
n Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
<= Int
0 = Double -> Either MathError Double
forall a b. b -> Either a b
Right ((Double
lo Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
hi) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
2)
    | Bool
otherwise = do
        let m :: Double
m = (Double
lo Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
hi) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
2
        v <- Double -> Either MathError Double
f Double
m
        if v < target then bisect (n - 1) m hi target f else bisect (n - 1) lo m target f
fact :: Int -> Double
fact :: Int -> Double
fact Int
n = [Double] -> Double
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
product [Double
1 .. Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
n]

-- | Dawson integral for |x|\<=100; scaled quadrature below |x|=10, inverse-power asymptotics above.
dawson :: Double -> Either MathError SFResult
dawson :: Double -> Either MathError SFResult
dawson Double
x = do
    Bool -> String -> Either MathError ()
require (Double -> Bool
finite Double
x Bool -> Bool -> Bool
&& Double -> Double
forall a. Num a => a -> a
abs Double
x Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
<= Double
100) String
"Dawson supports |x|<=100"
    if Double -> Double
forall a. Num a => a -> a
abs Double
x Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double
10
        then Int
-> Double -> (Int -> Double -> Double) -> Either MathError SFResult
series Int
100 (Double
1 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
x)) (\Int
n Double
t -> Double
t Double -> Double -> Double
forall a. Num a => a -> a -> a
* Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Int
2 Int -> Int -> Int
forall a. Num a => a -> a -> a
* Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
x Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
x))
        else (Double -> Double) -> Double -> Double -> Either MathError SFResult
integrate (\Double
t -> Double -> Double
forall a. Floating a => a -> a
exp (Double
t Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
t Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
x Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
x)) Double
0 Double
x

-- | E_n(x), integer n=0..100000,x>0. Uses a continued fraction for x>=1, the logarithmic E1 series and stable upward recurrence below one.
exponentialIntegral :: Int -> Double -> Either MathError SFResult
exponentialIntegral :: Int -> Double -> Either MathError SFResult
exponentialIntegral Int
n Double
x = do
    Bool -> String -> Either MathError ()
require (Int
n Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
>= Int
0 Bool -> Bool -> Bool
&& Int
n Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
<= Int
100000 Bool -> Bool -> Bool
&& Double -> Bool
finite Double
x Bool -> Bool -> Bool
&& Double
x Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double
0) String
"E_n needs integer n in0..100000 and x>0"
    if Int
n Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
0
        then Double -> Double -> Either MathError SFResult
estimate (Double -> Double
forall a. Floating a => a -> a
exp (-Double
x) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
x) Double
0
        else
            if Double
x Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
>= Double
1
                then do
                    let nonzero :: a -> a
nonzero a
y = if a -> a
forall a. Num a => a -> a
abs a
y a -> a -> Bool
forall a. Ord a => a -> a -> Bool
< a
1e-300 then a
1e-300 else a
y
                        go :: Int -> Double -> Double -> Double -> Double -> Either MathError SFResult
                        go :: Int
-> Double
-> Double
-> Double
-> Double
-> Either MathError SFResult
go Int
k Double
b Double
c Double
d Double
h
                            | Int
k Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
> Int
10000 = MathError -> Either MathError SFResult
forall a b. a -> Either a b
Left (String -> Int -> MathError
NonConvergence String
"exponential integral fraction" Int
10000)
                            | Bool
otherwise =
                                let an :: Double
an = -Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
k Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
k Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
n Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
1)
                                    bn :: Double
bn = Double
b Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
2
                                    dn :: Double
dn = Double
1 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double -> Double
forall {a}. (Ord a, Fractional a) => a -> a
nonzero (Double
bn Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
an Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
d)
                                    cn :: Double
cn = Double -> Double
forall {a}. (Ord a, Fractional a) => a -> a
nonzero (Double
bn Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
an Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
c)
                                    change :: Double
change = Double
dn Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
cn
                                    hn :: Double
hn = Double
h Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
change
                                 in if Double -> Double
forall a. Num a => a -> a
abs (Double
change Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
1) Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
2e-15
                                        then Double -> Double -> Either MathError SFResult
estimate (Double -> Double
forall a. Floating a => a -> a
exp (-Double
x) Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
hn) (Double -> Double
forall a. Num a => a -> a
abs (Double -> Double
forall a. Floating a => a -> a
exp (-Double
x) Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
hn Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
change Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
1)))
                                        else Int
-> Double
-> Double
-> Double
-> Double
-> Either MathError SFResult
go (Int
k Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1) Double
bn Double
cn Double
dn Double
hn
                    Int
-> Double
-> Double
-> Double
-> Double
-> Either MathError SFResult
go Int
1 (Double
x Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
n) Double
1e300 (Double
1 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
x Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
n)) (Double
1 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
x Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
n))
                else do
                    value <- Int
-> Double -> (Int -> Double -> Double) -> Either MathError SFResult
series Int
10000 (-Double
x) (\Int
k Double
t -> -Double
t Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
x Double -> Double -> Double
forall a. Num a => a -> a -> a
* Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
k Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral ((Int
k Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1) Int -> Int -> Int
forall a. Num a => a -> a -> a
* (Int
k Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1)))
                    let first = -Double
0.5772156649015329 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double -> Double
forall a. Floating a => a -> a
log Double
x Double -> Double -> Double
forall a. Num a => a -> a -> a
- SFResult -> Double
sfValue SFResult
value
                        advance (Double
v, Double
e) p
k =
                            let fk :: Double
fk = p -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral p
k; term :: Double
term = Double -> Double
forall a. Floating a => a -> a
exp (-Double
x) Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
x Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
v
                             in (Double
term Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
fk, (Double
x Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
e Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
4e-16 Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double -> Double
forall a. Num a => a -> a
abs (Double -> Double
forall a. Floating a => a -> a
exp (-Double
x)) Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double -> Double
forall a. Num a => a -> a
abs (Double
x Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
v))) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
fk)
                        (answer, errorEstimate) = Foldable.foldl' advance (first, sfError value) [1 .. n - 1]
                    estimate answer errorEstimate

{- | Real principal-value Ei for x>0; logarithmic integral li(x)=Ei(log x),
with the real PV continuation on 0\<x\<1.
-}
logarithmicIntegral :: Double -> Either MathError SFResult
logarithmicIntegral :: Double -> Either MathError SFResult
logarithmicIntegral Double
x = do
    Bool -> String -> Either MathError ()
require (Double -> Bool
finite Double
x Bool -> Bool -> Bool
&& Double
x Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double
0 Bool -> Bool -> Bool
&& Double
x Double -> Double -> Bool
forall a. Eq a => a -> a -> Bool
/= Double
1 Bool -> Bool -> Bool
&& Double -> Double
forall a. Num a => a -> a
abs (Double -> Double
forall a. Floating a => a -> a
log Double
x) Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
500) String
"li requires x>0, x/=1, |log x|<500"
    let y :: Double
y = Double -> Double
forall a. Floating a => a -> a
log Double
x
    if Double
y Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
0
        then do
            value <- Int -> Double -> Either MathError SFResult
exponentialIntegral Int
1 (-Double
y)
            estimate (-sfValue value) (sfError value)
        else do
            value <- Int
-> Double -> (Int -> Double -> Double) -> Either MathError SFResult
series Int
10000 Double
y (\Int
n Double
t -> Double
t Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
y Double -> Double -> Double
forall a. Num a => a -> a -> a
* Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
n Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral ((Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1) Int -> Int -> Int
forall a. Num a => a -> a -> a
* (Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1)))
            estimate (0.5772156649015329 + log y + sfValue value) (sfError value)

-- | Real Si(x) on |x|\<=1000, with its removable singularity explicitly filled at zero.
sineIntegral :: Double -> Either MathError SFResult
sineIntegral :: Double -> Either MathError SFResult
sineIntegral Double
x = do
    Bool -> String -> Either MathError ()
require (Double -> Bool
finite Double
x Bool -> Bool -> Bool
&& Double -> Double
forall a. Num a => a -> a
abs Double
x Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
<= Double
1000) String
"Si supports |x|<=1000"
    if Double -> Double
forall a. Num a => a -> a
abs Double
x Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double
50
        then do
            (f, g) <- Double -> Either MathError (SFResult, SFResult)
oscillatoryFG (Double -> Double
forall a. Num a => a -> a
abs Double
x)
            let t = Double -> Double
forall a. Num a => a -> a
abs Double
x; value = Double
forall a. Floating a => a
pi Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double -> Double
forall a. Floating a => a -> a
cos Double
t Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
t Double -> Double -> Double
forall a. Num a => a -> a -> a
* SFResult -> Double
sfValue SFResult
f Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double -> Double
forall a. Floating a => a -> a
sin Double
t Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
t Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
t) Double -> Double -> Double
forall a. Num a => a -> a -> a
* SFResult -> Double
sfValue SFResult
g
            estimate ((if x < 0 then -1 else 1) * value) (sfError f / t + sfError g / (t * t))
        else
            (Double -> Double) -> [Double] -> Either MathError SFResult
integrateSegments
                (\Double
t -> if Double
t Double -> Double -> Bool
forall a. Eq a => a -> a -> Bool
== Double
0 then Double
1 else Double -> Double
forall a. Floating a => a -> a
sin Double
t Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
t)
                (Double -> Double -> Int -> [Double]
linearSegments Double
0 Double
x (Int -> Int -> Int
forall a. Ord a => a -> a -> a
max Int
1 (Double -> Int
forall b. Integral b => Double -> b
forall a b. (RealFrac a, Integral b) => a -> b
ceiling (Double -> Double
forall a. Num a => a -> a
abs Double
x))))

{- | Real Ci(x) on 0\<x\<=1000, using phase-resolved quadrature near zero and
inverse-power oscillatory asymptotics above x=50.
-}
cosineIntegral :: Double -> Either MathError SFResult
cosineIntegral :: Double -> Either MathError SFResult
cosineIntegral Double
x = do
    Bool -> String -> Either MathError ()
require (Double -> Bool
finite Double
x Bool -> Bool -> Bool
&& Double
x Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double
0 Bool -> Bool -> Bool
&& Double
x Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
<= Double
1000) String
"Ci supports 0<x<=1000"
    if Double
x Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double
50
        then do
            (f, g) <- Double -> Either MathError (SFResult, SFResult)
oscillatoryFG Double
x
            estimate (sin x / x * sfValue f - cos x / (x * x) * sfValue g) (sfError f / x + sfError g / (x * x))
        else do
            value <-
                (Double -> Double) -> [Double] -> Either MathError SFResult
integrateSegments
                    (\Double
t -> if Double
t Double -> Double -> Bool
forall a. Eq a => a -> a -> Bool
== Double
0 then Double
0 else (Double -> Double
forall a. Floating a => a -> a
cos Double
t Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
1) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
t)
                    (Double -> Double -> Int -> [Double]
linearSegments Double
0 Double
x (Int -> Int -> Int
forall a. Ord a => a -> a -> a
max Int
1 (Double -> Int
forall b. Integral b => Double -> b
forall a b. (RealFrac a, Integral b) => a -> b
ceiling Double
x)))
            estimate (0.5772156649015329 + log x + sfValue value) (sfError value + 1e-14 * abs (log x))

oscillatoryFG :: Double -> Either MathError (SFResult, SFResult)
oscillatoryFG :: Double -> Either MathError (SFResult, SFResult)
oscillatoryFG Double
x = do
    f <- Int
-> Double -> (Int -> Double -> Double) -> Either MathError SFResult
series Int
100 Double
1 (\Int
k Double
t -> -Double
t Double -> Double -> Double
forall a. Num a => a -> a -> a
* Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral ((Int
2 Int -> Int -> Int
forall a. Num a => a -> a -> a
* Int
k) Int -> Int -> Int
forall a. Num a => a -> a -> a
* (Int
2 Int -> Int -> Int
forall a. Num a => a -> a -> a
* Int
k Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1)) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
x Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
x))
    g <- series 100 1 (\Int
k Double
t -> -Double
t Double -> Double -> Double
forall a. Num a => a -> a -> a
* Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral ((Int
2 Int -> Int -> Int
forall a. Num a => a -> a -> a
* Int
k Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1) Int -> Int -> Int
forall a. Num a => a -> a -> a
* (Int
2 Int -> Int -> Int
forall a. Num a => a -> a -> a
* Int
k)) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
x Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
x))
    pure (f, g)
linearSegments :: Double -> Double -> Int -> [Double]
linearSegments :: Double -> Double -> Int -> [Double]
linearSegments Double
a Double
b Int
n = [Double
a Double -> Double -> Double
forall a. Num a => a -> a -> a
+ (Double
b Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
a) Double -> Double -> Double
forall a. Num a => a -> a -> a
* Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
i Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
n | Int
i <- [Int
0 .. Int
n]]
integrateSegments :: (Double -> Double) -> [Double] -> Either MathError SFResult
integrateSegments :: (Double -> Double) -> [Double] -> Either MathError SFResult
integrateSegments Double -> Double
f [Double]
knots = do
    parts <- ((Double, Double) -> Either MathError SFResult)
-> [(Double, Double)] -> Either MathError [SFResult]
forall (t :: * -> *) (m :: * -> *) a b.
(Traversable t, Monad m) =>
(a -> m b) -> t a -> m (t b)
forall (m :: * -> *) a b. Monad m => (a -> m b) -> [a] -> m [b]
mapM ((Double -> Double -> Either MathError SFResult)
-> (Double, Double) -> Either MathError SFResult
forall a b c. (a -> b -> c) -> (a, b) -> c
uncurry ((Double -> Double) -> Double -> Double -> Either MathError SFResult
integrate Double -> Double
f)) ([Double] -> [Double] -> [(Double, Double)]
forall a b. [a] -> [b] -> [(a, b)]
zip [Double]
knots (Int -> [Double] -> [Double]
forall a. Int -> [a] -> [a]
drop Int
1 [Double]
knots))
    estimate (sum (map sfValue parts)) (sum (map sfError parts))

{- | Fresnel (C,S), phase pi*t^2/2, |x|\<=100. Phase-sized subintervals prevent
oscillatory quadrature aliasing; above |x|=8 inverse-power asymptotics avoid it.
-}
fresnel :: Double -> Either MathError (SFResult, SFResult)
fresnel :: Double -> Either MathError (SFResult, SFResult)
fresnel Double
x = do
    Bool -> String -> Either MathError ()
require (Double -> Bool
finite Double
x Bool -> Bool -> Bool
&& Double -> Double
forall a. Num a => a -> a
abs Double
x Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
<= Double
100) String
"Fresnel supports |x|<=100"
    let t :: Double
t = Double -> Double
forall a. Num a => a -> a
abs Double
x; sign :: Double
sign = if Double
x Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
0 then -Double
1 else Double
1
    if Double
t Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
>= Double
8
        then do
            let scale :: Double
scale = Double
forall a. Floating a => a
pi Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
forall a. Floating a => a
pi Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
t Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
t Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
t Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
t; theta :: Double
theta = Double
forall a. Floating a => a
pi Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
t Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
t Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
2
            f <- Int
-> Double -> (Int -> Double -> Double) -> Either MathError SFResult
series Int
100 Double
1 (\Int
k Double
v -> -Double
v Double -> Double -> Double
forall a. Num a => a -> a -> a
* Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral ((Int
4 Int -> Int -> Int
forall a. Num a => a -> a -> a
* Int
k Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1) Int -> Int -> Int
forall a. Num a => a -> a -> a
* (Int
4 Int -> Int -> Int
forall a. Num a => a -> a -> a
* Int
k Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
3)) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
scale)
            g <- series 100 1 (\Int
k Double
v -> -Double
v Double -> Double -> Double
forall a. Num a => a -> a -> a
* Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral ((Int
4 Int -> Int -> Int
forall a. Num a => a -> a -> a
* Int
k Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1) Int -> Int -> Int
forall a. Num a => a -> a -> a
* (Int
4 Int -> Int -> Int
forall a. Num a => a -> a -> a
* Int
k Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1)) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
scale)
            let fv = SFResult -> Double
sfValue SFResult
f Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
forall a. Floating a => a
pi Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
t)
                gv = SFResult -> Double
sfValue SFResult
g Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
forall a. Floating a => a
pi Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
forall a. Floating a => a
pi Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
t Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
t Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
t)
                err = SFResult -> Double
sfError SFResult
f Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
forall a. Floating a => a
pi Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
t) Double -> Double -> Double
forall a. Num a => a -> a -> a
+ SFResult -> Double
sfError SFResult
g Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
forall a. Floating a => a
pi Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
forall a. Floating a => a
pi Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
t Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
t Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
t) Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
1e-14
            c <- estimate (sign * (0.5 + fv * sin theta - gv * cos theta)) err
            ss <- estimate (sign * (0.5 - fv * cos theta - gv * sin theta)) err
            pure (c, ss)
        else do
            let n :: Int
n = Int -> Int -> Int
forall a. Ord a => a -> a -> a
max Int
1 (Double -> Int
forall b. Integral b => Double -> b
forall a b. (RealFrac a, Integral b) => a -> b
ceiling (Double
t Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
t)) :: Int
                knots :: [Double]
knots = [Double
sign Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double
forall a. Floating a => a -> a
sqrt (Double
t Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
t Double -> Double -> Double
forall a. Num a => a -> a -> a
* Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
i Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
n) | Int
i <- [Int
0 .. Int
n]]
            c <- (Double -> Double) -> [Double] -> Either MathError SFResult
integrateSegments (\Double
u -> Double -> Double
forall a. Floating a => a -> a
cos (Double
forall a. Floating a => a
pi Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
u Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
u Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
2)) [Double]
knots
            ss <- integrateSegments (\Double
u -> Double -> Double
forall a. Floating a => a -> a
sin (Double
forall a. Floating a => a
pi Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
u Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
u Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
2)) knots
            pure (c, ss)

{- | Legendre incomplete elliptic integrals, parameter m (not modulus k),
real branch |phi|\<=pi/2, m\<1; Pi additionally n\<1.
-}
ellipticF :: Double -> Double -> Either MathError SFResult
ellipticF :: Double -> Double -> Either MathError SFResult
ellipticF Double
phi Double
m = Double -> Double -> Double -> Int -> Either MathError SFResult
elliptic Double
phi Double
m Double
0 Int
0

-- | Legendre incomplete second-kind elliptic integral with the same amplitude/parameter convention.
ellipticE :: Double -> Double -> Either MathError SFResult
ellipticE :: Double -> Double -> Either MathError SFResult
ellipticE Double
phi Double
m = Double -> Double -> Double -> Int -> Either MathError SFResult
elliptic Double
phi Double
m Double
0 Int
1

-- | Legendre third-kind integral with arguments characteristic n, amplitude phi, parameter m; n\<1, m\<1, |phi|\<=pi/2.
ellipticPi :: Double -> Double -> Double -> Either MathError SFResult
ellipticPi :: Double -> Double -> Double -> Either MathError SFResult
ellipticPi Double
n Double
phi Double
m = Double -> Double -> Double -> Int -> Either MathError SFResult
elliptic Double
phi Double
m Double
n Int
2

elliptic :: Double -> Double -> Double -> Int -> Either MathError SFResult
elliptic :: Double -> Double -> Double -> Int -> Either MathError SFResult
elliptic Double
phi Double
m Double
n Int
kind = do
    Bool -> String -> Either MathError ()
require
        ((Double -> Bool) -> [Double] -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
all Double -> Bool
finite [Double
phi, Double
m, Double
n] Bool -> Bool -> Bool
&& Double -> Double
forall a. Num a => a -> a
abs Double
phi Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
<= Double
forall a. Floating a => a
pi Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
2 Bool -> Bool -> Bool
&& Double
m Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
1 Bool -> Bool -> Bool
&& Double
n Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
1)
        String
"real elliptic branch requires |phi|<=pi/2,m<1,n<1"
    (Double -> Double) -> Double -> Double -> Either MathError SFResult
integrate
        ( \Double
t ->
            let q :: Double
q = Double -> Double
forall a. Floating a => a -> a
sqrt (Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
m Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double
forall a. Floating a => a -> a
sin Double
t Double -> Int -> Double
forall a b. (Num a, Integral b) => a -> b -> a
^ (Int
2 :: Int))
             in if Int
kind Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
1 then Double
q else if Int
kind Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
2 then Double
1 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ ((Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
n Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double
forall a. Floating a => a -> a
sin Double
t Double -> Int -> Double
forall a b. (Num a, Integral b) => a -> b -> a
^ (Int
2 :: Int)) Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
q) else Double
1 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
q
        )
        Double
0
        Double
phi

{- | Real (sn,cn,dn), |u|\<=10^6 and 0\<=m\<1, using periodic amplitude
inversion. The complete period uses AGM; the argument cap prevents accumulated
period uncertainty from being presented as a resolved large-argument phase.
-}
jacobiElliptic :: Double -> Double -> Either MathError (Double, Double, Double)
jacobiElliptic :: Double -> Double -> Either MathError (Double, Double, Double)
jacobiElliptic Double
u Double
m = do
    Bool -> String -> Either MathError ()
require
        (Double -> Bool
finite Double
u Bool -> Bool -> Bool
&& Double -> Double
forall a. Num a => a -> a
abs Double
u Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
<= Double
1e6 Bool -> Bool -> Bool
&& Double -> Bool
finite Double
m Bool -> Bool -> Bool
&& Double
m Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
>= Double
0 Bool -> Bool -> Bool
&& Double
m Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
1)
        String
"Jacobi requires |u|<=1e6 and 0<=m<1"
    let agm :: Int -> Double -> Double -> Either MathError Double
        agm :: Int -> Double -> Double -> Either MathError Double
agm Int
0 Double
_ Double
_ = MathError -> Either MathError Double
forall a b. a -> Either a b
Left (String -> Int -> MathError
NonConvergence String
"Jacobi complete elliptic AGM" Int
100)
        agm Int
count Double
a Double
b =
            let arithmetic :: Double
arithmetic = (Double
a Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
b) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
2
                geometric :: Double
geometric = Double -> Double
forall a. Floating a => a -> a
sqrt (Double
a Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
b)
                -- Rounded AGM iterates can settle one ulp apart. Requiring a
                -- smaller difference would exhaust the budget at a fixed point.
                stagnant :: Bool
stagnant = Double
arithmetic Double -> Double -> Bool
forall a. Eq a => a -> a -> Bool
== Double
a Bool -> Bool -> Bool
&& Double
geometric Double -> Double -> Bool
forall a. Eq a => a -> a -> Bool
== Double
b
             in if Double -> Double
forall a. Num a => a -> a
abs (Double
a Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
b) Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
<= Double
2e-16 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
a Bool -> Bool -> Bool
|| Bool
stagnant
                    then Double -> Either MathError Double
checked (Double
forall a. Floating a => a
pi Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
a Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
b))
                    else Int -> Double -> Double -> Either MathError Double
agm (Int
count Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1) Double
arithmetic Double
geometric
    k <- Int -> Double -> Double -> Either MathError Double
agm Int
100 Double
1 (Double -> Double
forall a. Floating a => a -> a
sqrt (Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
m))
    let period = Double
4 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
k
        v = Double
u Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
period Double -> Double -> Double
forall a. Num a => a -> a -> a
* Integer -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Double -> Integer
forall b. Integral b => Double -> b
forall a b. (RealFrac a, Integral b) => a -> b
floor ((Double
u Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
k) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
period) :: Integer)
        sgn = if Double
v Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
0 then -Double
1 else Double
1
        w = Double -> Double
forall a. Num a => a -> a
abs Double
v
        reflected = Double
w Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double
k
        target = if Bool
reflected then Double
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
k Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
w else Double
w
    phi <- bisect 60 0 (pi / 2) target (\Double
t -> SFResult -> Double
sfValue (SFResult -> Double)
-> Either MathError SFResult -> Either MathError Double
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
<$> Double -> Double -> Either MathError SFResult
ellipticF Double
t Double
m)
    let sn = Double
sgn Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double
forall a. Floating a => a -> a
sin Double
phi; cn = (if Bool
reflected then -Double
1 else Double
1) Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double
forall a. Floating a => a -> a
cos Double
phi
    pure (sn, cn, sqrt (1 - m * sn * sn))

-- | Integer-order J_n(x) from libm, |n|\<=100000 and finite real x.
besselJ :: Int -> Double -> Either MathError Double
besselJ :: Int -> Double -> Either MathError Double
besselJ Int
n Double
x = do
    Bool -> String -> Either MathError ()
require (Integer -> Integer
forall a. Num a => a -> a
abs (Int -> Integer
forall a. Integral a => a -> Integer
toInteger Int
n) Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
<= Integer
100000 Bool -> Bool -> Bool
&& Double -> Bool
finite Double
x) String
"Bessel J integer order budget exceeded"
    Double -> Either MathError Double
checked (CDouble -> Double
forall a b. (Real a, Fractional b) => a -> b
realToFrac (CInt -> CDouble -> CDouble
cJn (Int -> CInt
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
n) (Double -> CDouble
forall a b. (Real a, Fractional b) => a -> b
realToFrac Double
x)))

-- | Integer-order Y_n(x) from libm, |n|\<=100000 and x>0.
besselY :: Int -> Double -> Either MathError Double
besselY :: Int -> Double -> Either MathError Double
besselY Int
n Double
x = do
    Bool -> String -> Either MathError ()
require (Integer -> Integer
forall a. Num a => a -> a
abs (Int -> Integer
forall a. Integral a => a -> Integer
toInteger Int
n) Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
<= Integer
100000 Bool -> Bool -> Bool
&& Double -> Bool
finite Double
x Bool -> Bool -> Bool
&& Double
x Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double
0) String
"Bessel Y integer order, x>0 required"
    Double -> Either MathError Double
checked (CDouble -> Double
forall a b. (Real a, Fractional b) => a -> b
realToFrac (CInt -> CDouble -> CDouble
cYn (Int -> CInt
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
n) (Double -> CDouble
forall a b. (Real a, Fractional b) => a -> b
realToFrac Double
x)))

-- | Modified Bessel I_n(x) via its convergent series, integer order 0..100 and |x|\<=700. Overflow is an error.
besselI :: Int -> Double -> Either MathError SFResult
besselI :: Int -> Double -> Either MathError SFResult
besselI Int
n Double
x = do
    Bool -> String -> Either MathError ()
require (Int
n Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
>= Int
0 Bool -> Bool -> Bool
&& Int
n Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
<= Int
100 Bool -> Bool -> Bool
&& Double -> Bool
finite Double
x Bool -> Bool -> Bool
&& Double -> Double
forall a. Num a => a -> a
abs Double
x Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
<= Double
700) String
"Bessel I order 0..100, |x|<=700"
    g <- Double -> Either MathError Double
gamma (Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
n Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
1)
    series
        10000
        ((x / 2) ^ n / g)
        (\Int
k Double
t -> Double
t Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
x Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
x Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
4 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
k Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
k Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
n)))

-- | exp(-|x|)*I_n(x). Uses the power series for |x|\<=50 and asymptotics when n\<=20 and |x|>n^2.
scaledBesselI :: Int -> Double -> Either MathError SFResult
scaledBesselI :: Int -> Double -> Either MathError SFResult
scaledBesselI Int
n Double
x
    | Double -> Double
forall a. Num a => a -> a
abs Double
x Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
<= Double
50 = do s <- Int -> Double -> Either MathError SFResult
besselI Int
n Double
x; estimate (exp (-abs x) * sfValue s) (exp (-abs x) * sfError s)
    | Bool
otherwise = do
        Bool -> String -> Either MathError ()
require
            (Int
n Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
>= Int
0 Bool -> Bool -> Bool
&& Int
n Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
<= Int
20 Bool -> Bool -> Bool
&& Double -> Bool
finite Double
x Bool -> Bool -> Bool
&& Double -> Double
forall a. Num a => a -> a
abs Double
x Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
* Int
n))
            String
"scaled I asymptotic requires |x|>n^2, n<=20"
        s <-
            Int
-> Double -> (Int -> Double -> Double) -> Either MathError SFResult
series
                Int
100
                Double
1
                ( \Int
k Double
t -> Double
t Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double
forall a. Num a => a -> a
negate (Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Int
4 Int -> Int -> Int
forall a. Num a => a -> a -> a
* Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
* Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- (Int
2 Int -> Int -> Int
forall a. Num a => a -> a -> a
* Int
k Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1) Int -> Int -> Int
forall a b. (Num a, Integral b) => a -> b -> a
^ (Int
2 :: Int))) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
8 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double
forall a. Num a => a -> a
abs Double
x Double -> Double -> Double
forall a. Num a => a -> a -> a
* Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
k)
                )
        let sign = if Double
x Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
0 Bool -> Bool -> Bool
&& Int -> Bool
forall a. Integral a => a -> Bool
odd Int
n then -Double
1 else Double
1
         in estimate
                (sign * sfValue s / (sqrt (2 * pi) * sqrt (abs x)))
                (sfError s / (sqrt (2 * pi) * sqrt (abs x)))

-- | Modified Bessel K_n(x)=exp(-x)*scaledBesselK(n,x), integer order 0..100,x>0.
besselK :: Int -> Double -> Either MathError SFResult
besselK :: Int -> Double -> Either MathError SFResult
besselK Int
n Double
x = do s <- Int -> Double -> Either MathError SFResult
scaledBesselK Int
n Double
x; estimate (exp (-x) * sfValue s) (exp (-x) * sfError s)

-- | exp(x)*K_n(x) from its cosh integral, integer order 0..100,x>0; scaling is inside the integrand.
scaledBesselK :: Int -> Double -> Either MathError SFResult
scaledBesselK :: Int -> Double -> Either MathError SFResult
scaledBesselK Int
n Double
x = do
    Bool -> String -> Either MathError ()
require (Int
n Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
>= Int
0 Bool -> Bool -> Bool
&& Int
n Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
<= Int
100 Bool -> Bool -> Bool
&& Double -> Bool
finite Double
x Bool -> Bool -> Bool
&& Double
x Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double
0) String
"Bessel K needs order 0..100,x>0"
    if Double
x Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double -> Double -> Double
forall a. Ord a => a -> a -> a
max Double
50 (Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
* Int
n))
        then do
            value <-
                Int
-> Double -> (Int -> Double -> Double) -> Either MathError SFResult
series
                    Int
100
                    Double
1
                    (\Int
k Double
t -> Double
t Double -> Double -> Double
forall a. Num a => a -> a -> a
* Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Int
4 Int -> Int -> Int
forall a. Num a => a -> a -> a
* Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
* Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- (Int
2 Int -> Int -> Int
forall a. Num a => a -> a -> a
* Int
k Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1) Int -> Int -> Int
forall a b. (Num a, Integral b) => a -> b -> a
^ (Int
2 :: Int)) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
8 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
x Double -> Double -> Double
forall a. Num a => a -> a -> a
* Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
k))
            let scale = Double -> Double
forall a. Floating a => a -> a
sqrt (Double
forall a. Floating a => a
pi Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
2) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double -> Double
forall a. Floating a => a -> a
sqrt Double
x
            estimate (scale * sfValue value) (scale * sfError value)
        else
            if Double
x Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
1e-4
                then do
                    -- The logarithmic small-x expansion avoids a truncated cosh integral whose
                    -- plateau extends beyond any practical finite integration limit.
                    let z :: Double
z = Double
x Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
x Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
4
                        terms :: [Double]
terms = Int -> [Double] -> [Double]
forall a. Int -> [a] -> [a]
take Int
12 ((Double -> Int -> Double) -> Double -> [Int] -> [Double]
forall b a. (b -> a -> b) -> b -> [a] -> [b]
scanl (\Double
t Int
k -> Double
t Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
z Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Int
k Int -> Int -> Int
forall a. Num a => a -> a -> a
* Int
k)) Double
1 [Int
1 :: Int ..])
                        hs :: [Double]
hs = (Double -> Double -> Double) -> Double -> [Double] -> [Double]
forall b a. (b -> a -> b) -> b -> [a] -> [b]
scanl Double -> Double -> Double
forall a. Num a => a -> a -> a
(+) Double
0 [Double
1 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
k | Int
k <- [Int
1 :: Int ..]]
                        i0 :: Double
i0 = [Double] -> Double
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum [Double]
terms
                        logTerm :: Double
logTerm = Double -> Double
forall a. Floating a => a -> a
log Double
x Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double -> Double
forall a. Floating a => a -> a
log Double
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
0.5772156649015329
                        i1 :: Double
i1 = Double
x Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* [Double] -> Double
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum (Int -> [Double] -> [Double]
forall a. Int -> [a] -> [a]
take Int
12 ((Double -> Int -> Double) -> Double -> [Int] -> [Double]
forall b a. (b -> a -> b) -> b -> [a] -> [b]
scanl (\Double
t Int
k -> Double
t Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
z Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Int
k Int -> Int -> Int
forall a. Num a => a -> a -> a
* (Int
k Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1))) Double
1 [Int
1 :: Int ..]))
                        k0 :: Double
k0 = -Double
logTerm Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
i0 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ [Double] -> Double
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum ((Double -> Double -> Double) -> [Double] -> [Double] -> [Double]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith Double -> Double -> Double
forall a. Num a => a -> a -> a
(*) [Double]
terms [Double]
hs)
                        k1 :: Double
k1 =
                            Double
i0 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
x
                                Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
logTerm Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
i1
                                Double -> Double -> Double
forall a. Num a => a -> a -> a
- [Double] -> Double
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum [Double
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
k Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
x Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
h Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
t | (Int
k, (Double
h, Double
t)) <- [Int] -> [(Double, Double)] -> [(Int, (Double, Double))]
forall a b. [a] -> [b] -> [(a, b)]
zip [Int
1 :: Int ..] (Int -> [(Double, Double)] -> [(Double, Double)]
forall a. Int -> [a] -> [a]
drop Int
1 ([Double] -> [Double] -> [(Double, Double)]
forall a b. [a] -> [b] -> [(a, b)]
zip [Double]
hs [Double]
terms))]
                        go :: Int -> Double -> Double -> Double
go Int
j Double
a Double
b | Int
j Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
n = Double
a | Bool
otherwise = Int -> Double -> Double -> Double
go (Int
j Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1) Double
b (Double
a Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Int
j Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
x Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
b)
                        answer :: Double
answer = if Int
n Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
0 then Double
k0 else Int -> Double -> Double -> Double
go Int
0 Double
k0 Double
k1
                    Double -> Double -> Either MathError SFResult
estimate (Double -> Double
forall a. Floating a => a -> a
exp Double
x Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
answer) (Double
1e-13 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double
forall a. Num a => a -> a
abs Double
answer)
                else
                    (Double -> Double) -> Double -> Double -> Either MathError SFResult
integrate
                        ( \Double
t ->
                            let y :: Double
y = Double -> Double
forall a. Num a => a -> a
abs (Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
n Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
t); logCosh :: Double
logCosh = Double
y Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double -> Double
forall a. Floating a => a -> a
log (Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double -> Double
forall a. Floating a => a -> a
exp (-Double
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
y)) Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double -> Double
forall a. Floating a => a -> a
log Double
2
                             in Double -> Double
forall a. Floating a => a -> a
exp (-Double
x Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double -> Double
forall a. Floating a => a -> a
cosh Double
t Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
1) Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
logCosh)
                        )
                        Double
0
                        (Double -> Double -> Double
forall a. Ord a => a -> a -> a
min Double
40 (Double -> Double
forall a. Floating a => a -> a
acosh (Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double -> Double -> Double
forall a. Ord a => a -> a -> a
max Double
50 (Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
* Int
n) Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
50) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
x)))

-- | Spherical j_n(x), order 0..100; uses a series for small x relative to order and upward recurrence otherwise.
sphericalBesselJ :: Int -> Double -> Either MathError Double
sphericalBesselJ :: Int -> Double -> Either MathError Double
sphericalBesselJ Int
n Double
x = do
    Bool -> String -> Either MathError ()
require (Int
n Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
>= Int
0 Bool -> Bool -> Bool
&& Int
n Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
<= Int
100 Bool -> Bool -> Bool
&& Double -> Bool
finite Double
x) String
"spherical Bessel J order 0..100"
    if Double
x Double -> Double -> Bool
forall a. Eq a => a -> a -> Bool
== Double
0
        then Double -> Either MathError Double
forall a b. b -> Either a b
Right (if Int
n Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
0 then Double
1 else Double
0)
        else
            if Double -> Double
forall a. Num a => a -> a
abs Double
x Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
n Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
1
                then do
                    s <-
                        Int
-> Double -> (Int -> Double -> Double) -> Either MathError SFResult
series
                            Int
10000
                            (Double
x Double -> Int -> Double
forall a b. (Num a, Integral b) => a -> b -> a
^ Int
n Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ [Double] -> Double
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
product [Double
1, Double
3 .. Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Int
2 Int -> Int -> Int
forall a. Num a => a -> a -> a
* Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1)])
                            (\Int
k Double
t -> -Double
t Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
x Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
x Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
k Double -> Double -> Double
forall a. Num a => a -> a -> a
* Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Int
2 Int -> Int -> Int
forall a. Num a => a -> a -> a
* Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
2 Int -> Int -> Int
forall a. Num a => a -> a -> a
* Int
k Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1)))
                    pure (sfValue s)
                else Double -> Either MathError Double
forall a. a -> Either MathError a
forall (f :: * -> *) a. Applicative f => a -> f a
pure (Int -> Double -> Double -> Double
forall {a}. Integral a => a -> Double -> Double -> Double
recur Int
n (Double -> Double
forall a. Floating a => a -> a
sin Double
x Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
x) (Double -> Double
forall a. Floating a => a -> a
sin Double
x Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
x Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
x) Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double -> Double
forall a. Floating a => a -> a
cos Double
x Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
x))
  where
    recur :: a -> Double -> Double -> Double
recur a
0 Double
a Double
_ = Double
a
    recur a
1 Double
_ Double
b = Double
b
    recur a
k Double
a Double
b = a -> Double -> Double -> Double
go a
1 Double
a Double
b
      where
        go :: a -> Double -> Double -> Double
go a
j Double
p Double
q | a
j a -> a -> Bool
forall a. Ord a => a -> a -> Bool
>= a
k = Double
q | Bool
otherwise = a -> Double -> Double -> Double
go (a
j a -> a -> a
forall a. Num a => a -> a -> a
+ a
1) Double
q (a -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral (a
2 a -> a -> a
forall a. Num a => a -> a -> a
* a
j a -> a -> a
forall a. Num a => a -> a -> a
+ a
1) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
x Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
q Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
p)

-- | Spherical y_n(x), integer order 0..100,x>0, by upward recurrence.
sphericalBesselY :: Int -> Double -> Either MathError Double
sphericalBesselY :: Int -> Double -> Either MathError Double
sphericalBesselY Int
n Double
x = do
    Bool -> String -> Either MathError ()
require (Int
n Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
>= Int
0 Bool -> Bool -> Bool
&& Int
n Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
<= Int
100 Bool -> Bool -> Bool
&& Double -> Bool
finite Double
x Bool -> Bool -> Bool
&& Double
x Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double
0) String
"spherical Y order 0..100,x>0"
    Double -> Either MathError Double
checked (Int -> Double -> Double -> Double
go Int
0 (-Double -> Double
forall a. Floating a => a -> a
cos Double
x Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
x) (-Double -> Double
forall a. Floating a => a -> a
cos Double
x Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
x Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
x) Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double -> Double
forall a. Floating a => a -> a
sin Double
x Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
x))
  where
    go :: Int -> Double -> Double -> Double
go Int
j Double
a Double
b | Int
j Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
n = Double
a | Bool
otherwise = Int -> Double -> Double -> Double
go (Int
j Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1) Double
b (Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Int
2 Int -> Int -> Int
forall a. Num a => a -> a -> a
* Int
j Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
3) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
x Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
b Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
a)

-- | Complex H_n^(kind)(x)=J_n(x)+/-iY_n(x), kind 1 or 2, integer order and x>0.
hankel :: Int -> Int -> Double -> Either MathError (Complex Double)
hankel :: Int -> Int -> Double -> Either MathError (Complex Double)
hankel Int
kind Int
n Double
x = do
    Bool -> String -> Either MathError ()
require (Int
kind Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
1 Bool -> Bool -> Bool
|| Int
kind Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
2) String
"Hankel kind 1 or 2"
    j <- Int -> Double -> Either MathError Double
besselJ Int
n Double
x
    y <- besselY n x
    pure (j :+ (if kind == 1 then y else -y))

{- | Airy Ai,Bi and derivatives solve y''=x y by convergent Maclaurin series.
The bounded real reference domain |x|\<=8 avoids severe recessive cancellation.
-}
airy :: Double -> Either MathError (SFResult, SFResult, SFResult, SFResult)
airy :: Double -> Either MathError (SFResult, SFResult, SFResult, SFResult)
airy Double
x = do
    Bool -> String -> Either MathError ()
require (Double -> Bool
finite Double
x Bool -> Bool -> Bool
&& Double -> Double
forall a. Num a => a -> a
abs Double
x Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
<= Double
8) String
"Airy reference domain |x|<=8"
    g1 <- Double -> Either MathError Double
gamma (Double
1 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
3)
    g2 <- gamma (2 / 3)
    let a0 = Double
1 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
3 Double -> Double -> Double
forall a. Floating a => a -> a -> a
** (Double
2 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
3) Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
g2)
        a1 = -Double
1 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
3 Double -> Double -> Double
forall a. Floating a => a -> a -> a
** (Double
1 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
3) Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
g1)
        eval Double
y0 Double
y1 =
            let cs :: [Double]
cs = Int -> [Double] -> [Double]
forall a. Int -> [a] -> [a]
take Int
150 (Double
y0 Double -> [Double] -> [Double]
forall a. a -> [a] -> [a]
: Double
y1 Double -> [Double] -> [Double]
forall a. a -> [a] -> [a]
: Double
0 Double -> [Double] -> [Double]
forall a. a -> [a] -> [a]
: (Double -> Int -> Double) -> [Double] -> [Int] -> [Double]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith (\Double
c Int
n -> Double
c Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
* (Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1))) [Double]
cs [Int
3 :: Int ..])
                v :: Double
v = [Double] -> Double
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum ((Double -> Double -> Double) -> [Double] -> [Double] -> [Double]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith Double -> Double -> Double
forall a. Num a => a -> a -> a
(*) [Double]
cs ((Double -> Double) -> Double -> [Double]
forall a. (a -> a) -> a -> [a]
iterate (Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
x) Double
1))
                d :: Double
d = [Double] -> Double
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum ((Int -> Double -> Double) -> [Int] -> [Double] -> [Double]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith (\Int
n Double
c -> Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
n Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
c Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
x Double -> Int -> Double
forall a b. (Num a, Integral b) => a -> b -> a
^ (Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1)) [Int
1 :: Int ..] (Int -> [Double] -> [Double]
forall a. Int -> [a] -> [a]
drop Int
1 [Double]
cs))
                e :: Double
e = Double
1e-14 Double -> Double -> Double
forall a. Num a => a -> a -> a
* [Double] -> Double
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum ((Double -> Double) -> [Double] -> [Double]
forall a b. (a -> b) -> [a] -> [b]
map Double -> Double
forall a. Num a => a -> a
abs ((Double -> Double -> Double) -> [Double] -> [Double] -> [Double]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith Double -> Double -> Double
forall a. Num a => a -> a -> a
(*) [Double]
cs ((Double -> Double) -> Double -> [Double]
forall a. (a -> a) -> a -> [a]
iterate (Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
x) Double
1)))
             in (Double -> Double -> SFResult
SFResult Double
v Double
e, Double -> Double -> SFResult
SFResult Double
d (Double
e Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
150))
        (a, ad) = eval a0 a1
        (b, bd) = eval (sqrt 3 * a0) (-sqrt 3 * a1)
    pure (a, ad, b, bd)

-- | Real Struve H_nu(x) by a convergent series; nu>-3/2 and 0\<=x\<=30.
struve :: Double -> Double -> Either MathError SFResult
struve :: Double -> Double -> Either MathError SFResult
struve Double
nu Double
x = do
    Bool -> String -> Either MathError ()
require (Double -> Bool
finite (Double
nu Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
x) Bool -> Bool -> Bool
&& Double
nu Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> -Double
1.5 Bool -> Bool -> Bool
&& Double
x Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
>= Double
0 Bool -> Bool -> Bool
&& Double
x Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
<= Double
30) String
"Struve H requires nu>-1.5,0<=x<=30"
    g <- Double -> Either MathError Double
gamma (Double
nu Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
1.5)
    series
        10000
        ((x / 2) ** (nu + 1) / (sqrt pi / 2 * g))
        (\Int
k Double
t -> -Double
t Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
x Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
x Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
4 Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
k Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
0.5) Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
k Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
nu Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
0.5)))

-- | Generalized pFq power series. Requires p\<=q or p=q+1 with |z|\<1 and no nonpositive-integer denominator parameters; finite budget 10000 terms.
hypergeometric :: [Double] -> [Double] -> Double -> Either MathError SFResult
hypergeometric :: [Double] -> [Double] -> Double -> Either MathError SFResult
hypergeometric [Double]
as [Double]
bs Double
z = do
    Bool -> String -> Either MathError ()
require
        ( (Double -> Bool) -> [Double] -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
all Double -> Bool
finite (Double
z Double -> [Double] -> [Double]
forall a. a -> [a] -> [a]
: [Double]
as [Double] -> [Double] -> [Double]
forall a. [a] -> [a] -> [a]
++ [Double]
bs)
            Bool -> Bool -> Bool
&& (Double -> Bool) -> [Double] -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
all (\Double
b -> Double
b Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double
0 Bool -> Bool -> Bool
|| Double
b Double -> Double -> Bool
forall a. Eq a => a -> a -> Bool
/= Integer -> Double
forall a. Num a => Integer -> a
fromInteger (Double -> Integer
forall b. Integral b => Double -> b
forall a b. (RealFrac a, Integral b) => a -> b
round Double
b)) [Double]
bs
            Bool -> Bool -> Bool
&& [Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
as Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
<= [Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
bs Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1
            Bool -> Bool -> Bool
&& ([Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
as Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
<= [Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
bs Bool -> Bool -> Bool
|| Double -> Double
forall a. Num a => a -> a
abs Double
z Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
1)
        )
        String
"hypergeometric power-series domain or denominator pole"
    Int
-> Double -> (Int -> Double -> Double) -> Either MathError SFResult
series
        Int
10000
        Double
1
        ( \Int
k Double
t ->
            Double
t
                Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
z
                Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
k
                Double -> Double -> Double
forall a. Num a => a -> a -> a
* [Double] -> Double
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
product [Double
a Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Int
k Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1) | Double
a <- [Double]
as]
                Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ [Double] -> Double
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
product [Double
b Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Int
k Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1) | Double
b <- [Double]
bs]
        )

-- | Kummer confluent hypergeometric 1F1(a;b;z), real power-series evaluation with denominator-pole rejection.
hypergeometric1F1 :: Double -> Double -> Double -> Either MathError SFResult
hypergeometric1F1 :: Double -> Double -> Double -> Either MathError SFResult
hypergeometric1F1 Double
a Double
b Double
z
    | Double
z Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
0 = do
        -- Kummer reflection turns a catastrophically alternating negative-axis
        -- series into positive-axis evaluation, with scaling applied afterwards.
        value <- [Double] -> [Double] -> Double -> Either MathError SFResult
hypergeometric [Double
b Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
a] [Double
b] (-Double
z)
        estimate (exp z * sfValue value) (exp z * sfError value)
    | Bool
otherwise = [Double] -> [Double] -> Double -> Either MathError SFResult
hypergeometric [Double
a] [Double
b] Double
z

-- | Gauss 2F1(a,b;c;z) in the real power-series disk |z|\<1.
hypergeometric2F1 :: Double -> Double -> Double -> Double -> Either MathError SFResult
hypergeometric2F1 :: Double -> Double -> Double -> Double -> Either MathError SFResult
hypergeometric2F1 Double
a Double
b Double
c Double
z = [Double] -> [Double] -> Double -> Either MathError SFResult
hypergeometric [Double
a, Double
b] [Double
c] Double
z

-- | Complex Kummer 1F1 in its entire z plane, rejecting denominator poles and nonfinite arithmetic; cancellation can limit accuracy for large imaginary z.
complex1F1 ::
    Complex Double -> Complex Double -> Complex Double -> Either MathError (Complex Double)
complex1F1 :: Complex Double
-> Complex Double
-> Complex Double
-> Either MathError (Complex Double)
complex1F1 Complex Double
a Complex Double
b Complex Double
z = Int
-> Complex Double
-> Complex Double
-> Either MathError (Complex Double)
go Int
1 Complex Double
1 Complex Double
1
  where
    go :: Int -> Complex Double -> Complex Double -> Either MathError (Complex Double)
    go :: Int
-> Complex Double
-> Complex Double
-> Either MathError (Complex Double)
go Int
k Complex Double
t Complex Double
s
        | Int
k Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
> Int
10000 = MathError -> Either MathError (Complex Double)
forall a b. a -> Either a b
Left (String -> Int -> MathError
NonConvergence String
"complex 1F1" Int
10000)
        | Complex Double -> Double
forall a. RealFloat a => Complex a -> a
magnitude (Complex Double
b Complex Double -> Complex Double -> Complex Double
forall a. Num a => a -> a -> a
+ Int -> Complex Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Int
k Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1)) Double -> Double -> Bool
forall a. Eq a => a -> a -> Bool
== Double
0 = MathError -> Either MathError (Complex Double)
forall a b. a -> Either a b
Left (String -> MathError
DomainError String
"complex 1F1 pole")
        | Bool
otherwise =
            let u :: Complex Double
u = Complex Double
t Complex Double -> Complex Double -> Complex Double
forall a. Num a => a -> a -> a
* (Complex Double
a Complex Double -> Complex Double -> Complex Double
forall a. Num a => a -> a -> a
+ Int -> Complex Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Int
k Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1)) Complex Double -> Complex Double -> Complex Double
forall a. Fractional a => a -> a -> a
/ (Complex Double
b Complex Double -> Complex Double -> Complex Double
forall a. Num a => a -> a -> a
+ Int -> Complex Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Int
k Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1)) Complex Double -> Complex Double -> Complex Double
forall a. Num a => a -> a -> a
* Complex Double
z Complex Double -> Complex Double -> Complex Double
forall a. Fractional a => a -> a -> a
/ Int -> Complex Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
k; q :: Complex Double
q = Complex Double
s Complex Double -> Complex Double -> Complex Double
forall a. Num a => a -> a -> a
+ Complex Double
u
             in if Bool -> Bool
not (Double -> Bool
finite (Complex Double -> Double
forall a. RealFloat a => Complex a -> a
magnitude Complex Double
q))
                    then MathError -> Either MathError (Complex Double)
forall a b. a -> Either a b
Left MathError
NonFiniteResult
                    else if Complex Double -> Double
forall a. RealFloat a => Complex a -> a
magnitude Complex Double
u Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
1e-14 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double -> Double
forall a. Ord a => a -> a -> a
max Double
1 (Complex Double -> Double
forall a. RealFloat a => Complex a -> a
magnitude Complex Double
q) then Complex Double -> Either MathError (Complex Double)
forall a b. b -> Either a b
Right Complex Double
q else Int
-> Complex Double
-> Complex Double
-> Either MathError (Complex Double)
go (Int
k Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1) Complex Double
u Complex Double
q

-- | Legendre P_n(x), n=0..1000 and |x|\<=1.
legendre :: Int -> Double -> Either MathError Double
legendre :: Int -> Double -> Either MathError Double
legendre Int
n Double
x = Int -> Int -> Double -> Either MathError Double
associatedLegendre Int
n Int
0 Double
x

-- | Associated P_l^m(x) with Condon--Shortley phase, 0\<=m\<=l\<=1000 and |x|\<=1.
associatedLegendre :: Int -> Int -> Double -> Either MathError Double
associatedLegendre :: Int -> Int -> Double -> Either MathError Double
associatedLegendre Int
l Int
m Double
x = do
    Bool -> String -> Either MathError ()
require
        (Int
l Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
>= Int
0 Bool -> Bool -> Bool
&& Int
l Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
<= Int
1000 Bool -> Bool -> Bool
&& Int
m Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
>= Int
0 Bool -> Bool -> Bool
&& Int
m Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
<= Int
l Bool -> Bool -> Bool
&& Double -> Bool
finite Double
x Bool -> Bool -> Bool
&& Double -> Double
forall a. Num a => a -> a
abs Double
x Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
<= Double
1)
        String
"associated Legendre needs 0<=m<=l<=1000, |x|<=1"
    let p0 :: Double
p0 = [Double] -> Double
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
product [Double
1, Double
3 .. Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Int
2 Int -> Int -> Int
forall a. Num a => a -> a -> a
* Int
m Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1)] Double -> Double -> Double
forall a. Num a => a -> a -> a
* (-Double -> Double
forall a. Floating a => a -> a
sqrt (Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
x Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
x)) Double -> Int -> Double
forall a b. (Num a, Integral b) => a -> b -> a
^ Int
m
        p1 :: Double
p1 = Double
x Double -> Double -> Double
forall a. Num a => a -> a -> a
* Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Int
2 Int -> Int -> Int
forall a. Num a => a -> a -> a
* Int
m Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1) Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
p0
        go :: Int -> Double -> Double -> Double
go Int
j Double
a Double
b
            | Int
j Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
l = Double
a
            | Bool
otherwise =
                Int -> Double -> Double -> Double
go
                    (Int
j Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1)
                    Double
b
                    ((Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Int
2 Int -> Int -> Int
forall a. Num a => a -> a -> a
* Int
j Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
3) Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
x Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
b Double -> Double -> Double
forall a. Num a => a -> a -> a
- Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Int
j Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
m Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1) Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
a) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Int
j Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
m Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
2))
     in Double -> Either MathError Double
checked (if Int
l Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
m then Double
p0 else Int -> Double -> Double -> Double
go Int
m Double
p0 Double
p1)

-- | Orthonormal complex spherical harmonic Y_l^m(theta,phi), theta in [0,pi], including negative-m conjugation symmetry.
sphericalHarmonic :: Int -> Int -> Double -> Double -> Either MathError (Complex Double)
sphericalHarmonic :: Int -> Int -> Double -> Double -> Either MathError (Complex Double)
sphericalHarmonic Int
l Int
m Double
theta Double
phi = do
    Bool -> String -> Either MathError ()
require
        (Int
l Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
>= Int
0 Bool -> Bool -> Bool
&& Int -> Int
forall a. Num a => a -> a
abs Int
m Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
<= Int
l Bool -> Bool -> Bool
&& Double -> Bool
finite (Double
theta Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
phi) Bool -> Bool -> Bool
&& Double
theta Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
>= Double
0 Bool -> Bool -> Bool
&& Double
theta Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
<= Double
forall a. Floating a => a
pi)
        String
"spherical harmonic invalid degree/angle"
    p <- Int -> Int -> Double -> Either MathError Double
associatedLegendre Int
l (Int -> Int
forall a. Num a => a -> a
abs Int
m) (Double -> Double
forall a. Floating a => a -> a
cos Double
theta)
    a <- logGamma (fromIntegral (l - abs m + 1))
    b <- logGamma (fromIntegral (l + abs m + 1))
    let y =
            (Double -> Double
forall a. Floating a => a -> a
sqrt (Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Int
2 Int -> Int -> Int
forall a. Num a => a -> a -> a
* Int
l Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
4 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
forall a. Floating a => a
pi) Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double
forall a. Floating a => a -> a
exp (Double
a Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
b)) Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
p Double -> Double -> Complex Double
forall a. a -> a -> Complex a
:+ Double
0)
                Complex Double -> Complex Double -> Complex Double
forall a. Num a => a -> a -> a
* Double -> Complex Double
forall a. Floating a => a -> Complex a
cis (Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Int -> Int
forall a. Num a => a -> a
abs Int
m) Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
phi)
    pure (if m < 0 then (if odd m then -1 else 1) * conjugate y else y)

-- | Hurwitz zeta for s>1,q>0 using 64 explicit terms plus seven Euler--Maclaurin corrections.
hurwitzZeta :: Double -> Double -> Either MathError SFResult
hurwitzZeta :: Double -> Double -> Either MathError SFResult
hurwitzZeta Double
s Double
q = do
    Bool -> String -> Either MathError ()
require (Double -> Bool
finite (Double
s Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
q) Bool -> Bool -> Bool
&& Double
s Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double
1 Bool -> Bool -> Bool
&& Double
q Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double
0) String
"Hurwitz zeta requires s>1,q>0"
    let n :: Int
n = Int
64 :: Int
        a :: Double
a = Double
q Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
n
        terms :: [Double]
terms = [(Double
q Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
k) Double -> Double -> Double
forall a. Floating a => a -> a -> a
** (-Double
s) | Int
k <- [Int
0 .. Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1]]
        coeff :: [Double]
coeff = [Double
1 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
12, -Double
1 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
720, Double
1 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
30240, -Double
1 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
1209600, Double
1 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
47900160, -Double
691 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
1307674368000, Double
7 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
523069747200]
        correction :: [Double]
correction =
            [ Double
b Double -> Double -> Double
forall a. Num a => a -> a -> a
* [Double] -> Double
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
product [Double
s Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
j | Int
j <- [Int
0 .. Int
2 Int -> Int -> Int
forall a. Num a => a -> a -> a
* Int
k Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
2]] Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
a Double -> Double -> Double
forall a. Floating a => a -> a -> a
** (-Double
s Double -> Double -> Double
forall a. Num a => a -> a -> a
- Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Int
2 Int -> Int -> Int
forall a. Num a => a -> a -> a
* Int
k Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1))
            | (Int
k, Double
b) <- [Int] -> [Double] -> [(Int, Double)]
forall a b. [a] -> [b] -> [(a, b)]
zip [Int
1 :: Int ..] [Double]
coeff
            ]
    Double -> Double -> Either MathError SFResult
estimate
        ([Double] -> Double
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum [Double]
terms Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
a Double -> Double -> Double
forall a. Floating a => a -> a -> a
** (Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
s) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
s Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
1) Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
0.5 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
a Double -> Double -> Double
forall a. Floating a => a -> a -> a
** (-Double
s) Double -> Double -> Double
forall a. Num a => a -> a -> a
+ [Double] -> Double
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum [Double]
correction)
        (Double -> Double
forall a. Num a => a -> a
abs ([Double] -> Double
forall a. HasCallStack => [a] -> a
last [Double]
correction))

-- | Real zeta: Euler--Maclaurin for s>1, Euler-transformed eta for 0\<s\<1 and reflection for -170\<s\<0. The pole s=1 is rejected.
riemannZeta :: Double -> Either MathError SFResult
riemannZeta :: Double -> Either MathError SFResult
riemannZeta Double
s
    | Double
s Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double
1 = Double -> Double -> Either MathError SFResult
hurwitzZeta Double
s Double
1
    | Double
s Double -> Double -> Bool
forall a. Eq a => a -> a -> Bool
== Double
0 = SFResult -> Either MathError SFResult
forall a b. b -> Either a b
Right (Double -> Double -> SFResult
SFResult (-Double
0.5) Double
0)
    | Double
s Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
0 = do
        Bool -> String -> Either MathError ()
require (Double -> Bool
finite Double
s Bool -> Bool -> Bool
&& Double
s Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> -Double
170) String
"zeta reflection domain -170<s<0"
        if Double
s Double -> Double -> Bool
forall a. Eq a => a -> a -> Bool
== Integer -> Double
forall a. Num a => Integer -> a
fromInteger (Double -> Integer
forall b. Integral b => Double -> b
forall a b. (RealFrac a, Integral b) => a -> b
round Double
s) Bool -> Bool -> Bool
&& Integer -> Bool
forall a. Integral a => a -> Bool
even (Double -> Integer
forall b. Integral b => Double -> b
forall a b. (RealFrac a, Integral b) => a -> b
round Double
s :: Integer)
            then SFResult -> Either MathError SFResult
forall a b. b -> Either a b
Right (Double -> Double -> SFResult
SFResult Double
0 Double
0)
            else do
                z <- Double -> Double -> Either MathError SFResult
hurwitzZeta (Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
s) Double
1
                g <- gamma (1 - s)
                estimate
                    (2 ** s * pi ** (s - 1) * sin (pi * s / 2) * g * sfValue z)
                    (abs (2 ** s * pi ** (s - 1) * sin (pi * s / 2) * g) * sfError z)
    | Bool
otherwise = do
        Bool -> String -> Either MathError ()
require (Double -> Bool
finite Double
s Bool -> Bool -> Bool
&& Double
s Double -> Double -> Bool
forall a. Eq a => a -> a -> Bool
/= Double
1) String
"zeta pole at 1"
        let rows :: [[Double]]
rows =
                Int -> [[Double]] -> [[Double]]
forall a. Int -> [a] -> [a]
take
                    Int
64
                    (([Double] -> [Double]) -> [Double] -> [[Double]]
forall a. (a -> a) -> a -> [a]
iterate (\[Double]
xs -> (Double -> Double -> Double) -> [Double] -> [Double] -> [Double]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith (-) [Double]
xs (Int -> [Double] -> [Double]
forall a. Int -> [a] -> [a]
drop Int
1 [Double]
xs)) [Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Int
k Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1) Double -> Double -> Double
forall a. Floating a => a -> a -> a
** (-Double
s) | Int
k <- [Int
0 .. Int
64 :: Int]])
            terms :: [Double]
terms = [([Double]
row [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
0) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
2 Double -> Double -> Double
forall a. Floating a => a -> a -> a
** Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Int
k Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1) | (Int
k, [Double]
row) <- [Int] -> [[Double]] -> [(Int, [Double])]
forall a b. [a] -> [b] -> [(a, b)]
zip [Int
0 :: Int ..] [[Double]]
rows]
            den :: Double
den = Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
2 Double -> Double -> Double
forall a. Floating a => a -> a -> a
** (Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
s)
        Double -> Double -> Either MathError SFResult
estimate ([Double] -> Double
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum [Double]
terms Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
den) (Double -> Double
forall a. Num a => a -> a
abs ([Double] -> Double
forall a. HasCallStack => [a] -> a
last [Double]
terms Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
den) Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
1e-14 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double -> Double
forall a. Num a => a -> a
abs Double
den)

-- | Real polylog Li_s(z) in the convergent power-series disk |z|\<1; no cross-cut continuation is implied.
polylogarithm :: Double -> Double -> Either MathError SFResult
polylogarithm :: Double -> Double -> Either MathError SFResult
polylogarithm Double
s Double
z = do
    Bool -> String -> Either MathError ()
require (Double -> Bool
finite (Double
s Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
z) Bool -> Bool -> Bool
&& Double -> Double
forall a. Num a => a -> a
abs Double
z Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
1) String
"real polylog power series needs |z|<1"
    Int
-> Double -> (Int -> Double -> Double) -> Either MathError SFResult
series Int
100000 Double
z (\Int
k Double
t -> Double
t Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
z Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
k Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Int
k Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1)) Double -> Double -> Double
forall a. Floating a => a -> a -> a
** Double
s)

-- | Real Li_2(x) for x\<=1, with reflection and negative-axis transformations for convergent evaluation.
dilogarithm :: Double -> Either MathError SFResult
dilogarithm :: Double -> Either MathError SFResult
dilogarithm Double
x
    | Double
x Double -> Double -> Bool
forall a. Eq a => a -> a -> Bool
== Double
1 = SFResult -> Either MathError SFResult
forall a b. b -> Either a b
Right (Double -> Double -> SFResult
SFResult (Double
forall a. Floating a => a
pi Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
forall a. Floating a => a
pi Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
6) Double
1e-15)
    | Double
x Double -> Double -> Bool
forall a. Eq a => a -> a -> Bool
== -Double
1 = SFResult -> Either MathError SFResult
forall a b. b -> Either a b
Right (Double -> Double -> SFResult
SFResult (-Double
forall a. Floating a => a
pi Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
forall a. Floating a => a
pi Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
12) Double
1e-15)
    | Double
x Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< -Double
0.5 = do
        Bool -> String -> Either MathError ()
require (Double -> Bool
finite Double
x) String
"dilog finite real x<=1"
        s <- Double -> Either MathError SFResult
dilogarithm (Double
x Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
x Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
1))
        estimate (-sfValue s - log (1 - x) ^ (2 :: Int) / 2) (sfError s)
    | Double
x Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double
0.5 Bool -> Bool -> Bool
&& Double
x Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
1 = do s <- Double -> Either MathError SFResult
dilogarithm (Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
x); estimate (pi * pi / 6 - log x * log (1 - x) - sfValue s) (sfError s)
    | Bool
otherwise = Double -> Double -> Either MathError SFResult
polylogarithm Double
2 Double
x

-- | Real Lambert W branch 0 or -1, x>=-1/e and x\<0 on branch -1, using bounded Halley iteration.
lambertW :: Int -> Double -> Either MathError SFResult
lambertW :: Int -> Double -> Either MathError SFResult
lambertW Int
branch Double
x = do
    Bool -> String -> Either MathError ()
require
        (Double -> Bool
finite Double
x Bool -> Bool -> Bool
&& (Int
branch Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
0 Bool -> Bool -> Bool
|| Int
branch Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== -Int
1) Bool -> Bool -> Bool
&& Double
x Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
>= -Double
1 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double -> Double
forall a. Floating a => a -> a
exp Double
1 Bool -> Bool -> Bool
&& (Int
branch Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
0 Bool -> Bool -> Bool
|| Double
x Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
0))
        String
"real Lambert W branches 0 and -1 require x>=-1/e; -1 needs x<0"
    if Double
x Double -> Double -> Bool
forall a. Eq a => a -> a -> Bool
== -Double
1 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double -> Double
forall a. Floating a => a -> a
exp Double
1
        then SFResult -> Either MathError SFResult
forall a b. b -> Either a b
Right (Double -> Double -> SFResult
SFResult (-Double
1) Double
1e-8)
        else if Double
x Double -> Double -> Bool
forall a. Eq a => a -> a -> Bool
== Double
0 then SFResult -> Either MathError SFResult
forall a b. b -> Either a b
Right (Double -> Double -> SFResult
SFResult Double
0 Double
0) else Int -> Double -> Either MathError SFResult
go Int
100 Double
initial
  where
    initial :: Double
initial | Int
branch Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== -Int
1 = Double -> Double
forall a. Floating a => a -> a
log (-Double
x) Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double -> Double
forall a. Floating a => a -> a
log (-Double -> Double
forall a. Floating a => a -> a
log (-Double
x)) | Double
x Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
1 = Double
x | Bool
otherwise = Double -> Double
forall a. Floating a => a -> a
log Double
x Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double -> Double
forall a. Floating a => a -> a
log (Double -> Double -> Double
forall a. Ord a => a -> a -> a
max Double
1 (Double -> Double
forall a. Floating a => a -> a
log Double
x))
    go :: Int -> Double -> Either MathError SFResult
    go :: Int -> Double -> Either MathError SFResult
go Int
0 Double
_ = MathError -> Either MathError SFResult
forall a b. a -> Either a b
Left (String -> Int -> MathError
NonConvergence String
"Lambert W Halley" Int
100)
    go Int
n Double
w =
        let e :: Double
e = Double -> Double
forall a. Floating a => a -> a
exp Double
w; f :: Double
f = Double
w Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
e Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
x; d :: Double
d = Double
e Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
w Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
1) Double -> Double -> Double
forall a. Num a => a -> a -> a
- (Double
w Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
2) Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
f Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
w Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
2); v :: Double
v = Double
w Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
f Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
d
         in if Double -> Double
forall a. Num a => a -> a
abs (Double
v Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
w) Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
1e-13 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double -> Double
forall a. Ord a => a -> a -> a
max Double
1 (Double -> Double
forall a. Num a => a -> a
abs Double
v) then Double -> Double -> Either MathError SFResult
estimate Double
v (Double -> Double
forall a. Num a => a -> a
abs (Double
v Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
w)) else Int -> Double -> Either MathError SFResult
go (Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1) Double
v

{- | Normalized Debye D_n(x), integer n=1..16,x>=0, with D_n(0)=1.
At large x the complete Bose integral supplies the asymptote and an
exponentially small omitted-tail estimate; logarithmic scaling avoids x^n overflow.
-}
debye :: Int -> Double -> Either MathError SFResult
debye :: Int -> Double -> Either MathError SFResult
debye Int
n Double
x = do
    Bool -> String -> Either MathError ()
require (Int
n Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
>= Int
1 Bool -> Bool -> Bool
&& Int
n Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
<= Int
16 Bool -> Bool -> Bool
&& Double -> Bool
finite Double
x Bool -> Bool -> Bool
&& Double
x Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
>= Double
0) String
"Debye order 1..16,x>=0"
    if Double
x Double -> Double -> Bool
forall a. Eq a => a -> a -> Bool
== Double
0
        then SFResult -> Either MathError SFResult
forall a b. b -> Either a b
Right (Double -> Double -> SFResult
SFResult Double
1 Double
0)
        else
            if Double
x Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double
100 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
n
                then do
                    z <- Double -> Double -> Either MathError SFResult
hurwitzZeta (Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
n Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
1) Double
1
                    lg <- logGamma (fromIntegral n + 1)
                    let value = Double -> Double
forall a. Floating a => a -> a
exp (Double -> Double
forall a. Floating a => a -> a
log (Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
n) Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
lg Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double -> Double
forall a. Floating a => a -> a
log (SFResult -> Double
sfValue SFResult
z) Double -> Double -> Double
forall a. Num a => a -> a -> a
- Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
n Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double
forall a. Floating a => a -> a
log Double
x)
                        tailBound = Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
n Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double
forall a. Floating a => a -> a
exp (-Double
x) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ ((Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double -> Double
forall a. Floating a => a -> a
exp (-Double
x)) Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
n Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
x))
                    estimate value (tailBound + abs value * sfError z / sfValue z)
                else do
                    s <- (Double -> Double) -> Double -> Double -> Either MathError SFResult
integrate (\Double
t -> if Double
t Double -> Double -> Bool
forall a. Eq a => a -> a -> Bool
== Double
0 then if Int
n Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
1 then Double
1 else Double
0 else Double
t Double -> Int -> Double
forall a b. (Num a, Integral b) => a -> b -> a
^ Int
n Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double -> Double
forall a. Floating a => a -> a
exp Double
t Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
1)) Double
0 Double
x
                    estimate (fromIntegral n * sfValue s / x ^ n) (fromIntegral n * sfError s / x ^ n)

-- | Periodic Clausen Cl_2(theta) from an integral with a squared-variable substitution at its logarithmic endpoint.
clausen :: Double -> Either MathError SFResult
clausen :: Double -> Either MathError SFResult
clausen Double
theta = do
    Bool -> String -> Either MathError ()
require (Double -> Bool
finite Double
theta) String
"Clausen finite angle"
    -- libm trigonometric reduction retains phase where theta-2*pi*k would
    -- round away the entire remainder for a huge finite angle.
    let x :: Double
x = Double -> Double -> Double
forall a. RealFloat a => a -> a -> a
atan2 (Double -> Double
forall a. Floating a => a -> a
sin Double
theta) (Double -> Double
forall a. Floating a => a -> a
cos Double
theta)
     in if Double
x Double -> Double -> Bool
forall a. Eq a => a -> a -> Bool
== Double
0
            then SFResult -> Either MathError SFResult
forall a b. b -> Either a b
Right (Double -> Double -> SFResult
SFResult Double
0 Double
0)
            else (Double -> Double) -> Double -> Double -> Either MathError SFResult
integrate (\Double
u -> if Double
u Double -> Double -> Bool
forall a. Eq a => a -> a -> Bool
== Double
0 then Double
0 else -Double
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
x Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
u Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double
forall a. Floating a => a -> a
log (Double -> Double
forall a. Num a => a -> a
abs (Double
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double
forall a. Floating a => a -> a
sin (Double
x Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
u Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
u Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
2)))) Double
0 Double
1

-- | Complete normalized Fermi--Dirac integral F_order(eta), 0\<=order\<=20, |eta|\<=100.
fermiDirac :: Double -> Double -> Either MathError SFResult
fermiDirac :: Double -> Double -> Either MathError SFResult
fermiDirac Double
order Double
eta = Bool -> Double -> Double -> Either MathError SFResult
quantumIntegral Bool
False Double
order Double
eta

-- | Complete normalized Bose--Einstein integral, same order range, -100\<=eta\<=0 and positive order at eta=0.
boseEinstein :: Double -> Double -> Either MathError SFResult
boseEinstein :: Double -> Double -> Either MathError SFResult
boseEinstein Double
order Double
eta = Bool -> Double -> Double -> Either MathError SFResult
quantumIntegral Bool
True Double
order Double
eta

quantumIntegral :: Bool -> Double -> Double -> Either MathError SFResult
quantumIntegral :: Bool -> Double -> Double -> Either MathError SFResult
quantumIntegral Bool
bose Double
order Double
eta = do
    Bool -> String -> Either MathError ()
require
        ( Double -> Bool
finite (Double
order Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
eta)
            Bool -> Bool -> Bool
&& Double
order Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
>= Double
0
            Bool -> Bool -> Bool
&& Double
order Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
<= Double
20
            Bool -> Bool -> Bool
&& Double
eta Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
<= Double
100
            Bool -> Bool -> Bool
&& Double
eta Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
>= -Double
100
            Bool -> Bool -> Bool
&& (Bool -> Bool
not Bool
bose Bool -> Bool -> Bool
|| (Double
eta Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
<= Double
0 Bool -> Bool -> Bool
&& (Double
eta Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
0 Bool -> Bool -> Bool
|| Double
order Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double
0)))
        )
        String
"quantum integral: 0<=order<=20,|eta|<=100; Bose eta<=0 and order>0 at eta=0"
    g <- Double -> Either MathError Double
gamma (Double
order Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
1)
    let top = Double -> Double -> Double
forall a. Ord a => a -> a -> a
max Double
50 (Double
eta Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
50); sgn = if Bool
bose then -Double
1 else Double
1
    s <-
        integrate
            (\Double
u -> if Double
u Double -> Double -> Bool
forall a. Eq a => a -> a -> Bool
== Double
0 then Double
0 else let t :: Double
t = Double
top Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
u Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
u in Double
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
top Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
u Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
t Double -> Double -> Double
forall a. Floating a => a -> a -> a
** Double
order Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double -> Double
forall a. Floating a => a -> a
exp (Double
t Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
eta) Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
sgn))
            0
            1
    estimate (sfValue s / g) (sfError s / g + exp (eta - top) * top ** order / g)

{- | Regular and irregular Coulomb waves F_l,G_l. The regular function follows
the convergent hypergeometric representation. G is integrated backwards from
a first-order asymptotic boundary at R>=300 using RK4; its error estimate
includes the omitted O((l(l+1)+eta^2)/R) term. Real reference domain l\<=8,
|eta|\<=3, 0.5\<=rho\<=30. This is deliberately not a high-accuracy scattering solver.
-}
coulombWave :: Int -> Double -> Double -> Either MathError (SFResult, SFResult)
coulombWave :: Int -> Double -> Double -> Either MathError (SFResult, SFResult)
coulombWave Int
l Double
eta Double
rho = do
    Bool -> String -> Either MathError ()
require
        (Int
l Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
>= Int
0 Bool -> Bool -> Bool
&& Int
l Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
<= Int
8 Bool -> Bool -> Bool
&& Double -> Bool
finite (Double
eta Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
rho) Bool -> Bool -> Bool
&& Double -> Double
forall a. Num a => a -> a
abs Double
eta Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
<= Double
3 Bool -> Bool -> Bool
&& Double
rho Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
>= Double
0.5 Bool -> Bool -> Bool
&& Double
rho Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
<= Double
30)
        String
"Coulomb reference domain l<=8,|eta|<=3,0.5<=rho<=30"
    let z :: Complex Double
z = Complex Double -> Complex Double
complexLogGamma (Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Int
l Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1) Double -> Double -> Complex Double
forall a. a -> a -> Complex a
:+ Double
eta)
        sigma :: Double
sigma = Complex Double -> Double
forall a. Complex a -> a
imagPart Complex Double
z
        lc :: Double
lc = Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
l Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double
forall a. Floating a => a -> a
log Double
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
forall a. Floating a => a
pi Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
eta Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Complex Double -> Double
forall a. Complex a -> a
realPart Complex Double
z Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double -> Double
forall a. Floating a => a -> a
log (Int -> Double
fact (Int
2 Int -> Int -> Int
forall a. Num a => a -> a -> a
* Int
l Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1))
        start :: Double
start = Double
0.5
        aa :: Complex Double
aa = Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Int
l Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1) Double -> Double -> Complex Double
forall a. a -> a -> Complex a
:+ (-Double
eta)
        bb :: Complex Double
bb = Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Int
2 Int -> Int -> Int
forall a. Num a => a -> a -> a
* Int
l Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
2) Double -> Double -> Complex Double
forall a. a -> a -> Complex a
:+ Double
0
    h <- Complex Double
-> Complex Double
-> Complex Double
-> Either MathError (Complex Double)
complex1F1 Complex Double
aa Complex Double
bb (Double
0 Double -> Double -> Complex Double
forall a. a -> a -> Complex a
:+ (Double
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
start))
    hd <- complex1F1 (aa + 1) (bb + 1) (0 :+ (2 * start))
    let pref = Double -> Double
forall a. Floating a => a -> a
exp Double
lc
        f0 = Double
pref Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
start Double -> Int -> Double
forall a b. (Num a, Integral b) => a -> b -> a
^ (Int
l Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1) Double -> Double -> Double
forall a. Num a => a -> a -> a
* Complex Double -> Double
forall a. Complex a -> a
realPart (Double -> Complex Double
forall a. Floating a => a -> Complex a
cis (-Double
start) Complex Double -> Complex Double -> Complex Double
forall a. Num a => a -> a -> a
* Complex Double
h)
        df0 =
            Double
pref
                Double -> Double -> Double
forall a. Num a => a -> a -> a
* Complex Double -> Double
forall a. Complex a -> a
realPart
                    ( Double -> Complex Double
forall a. Floating a => a -> Complex a
cis (-Double
start)
                        Complex Double -> Complex Double -> Complex Double
forall a. Num a => a -> a -> a
* ( ((Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Int
l Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1) Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
start Double -> Int -> Double
forall a b. (Num a, Integral b) => a -> b -> a
^ Int
l) Double -> Double -> Complex Double
forall a. a -> a -> Complex a
:+ (-Double
start Double -> Int -> Double
forall a b. (Num a, Integral b) => a -> b -> a
^ (Int
l Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1))) Complex Double -> Complex Double -> Complex Double
forall a. Num a => a -> a -> a
* Complex Double
h
                                Complex Double -> Complex Double -> Complex Double
forall a. Num a => a -> a -> a
+ (Double
start Double -> Int -> Double
forall a b. (Num a, Integral b) => a -> b -> a
^ (Int
l Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1) Double -> Double -> Complex Double
forall a. a -> a -> Complex a
:+ Double
0) Complex Double -> Complex Double -> Complex Double
forall a. Num a => a -> a -> a
* (Double
0 Double -> Double -> Complex Double
forall a. a -> a -> Complex a
:+ Double
2) Complex Double -> Complex Double -> Complex Double
forall a. Num a => a -> a -> a
* Complex Double
aa Complex Double -> Complex Double -> Complex Double
forall a. Fractional a => a -> a -> a
/ Complex Double
bb Complex Double -> Complex Double -> Complex Double
forall a. Num a => a -> a -> a
* Complex Double
hd
                          )
                    )
        accel Double
t Double
y = (Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Int
l Int -> Int -> Int
forall a. Num a => a -> a -> a
* (Int
l Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1)) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
t Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
t) Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
eta Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
t Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
1) Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
y
        advance Double
dt (Double
t, Double
y, Double
v) =
            let a :: Double
a = Double -> Double -> Double
accel Double
t Double
y
                b :: Double
b = Double -> Double -> Double
accel (Double
t Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
dt Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
2) (Double
y Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
dt Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
v Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
2)
                c :: Double
c = Double -> Double -> Double
accel (Double
t Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
dt Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
2) (Double
y Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
dt Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
v Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
dt Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
dt Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
a Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
4)
                d :: Double
d = Double -> Double -> Double
accel (Double
t Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
dt) (Double
y Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
dt Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
v Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
dt Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
dt Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
b Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
2)
             in (Double
t Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
dt, Double
y Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
dt Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
v Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
dt Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
dt Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
a Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
b Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
c) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
6, Double
v Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
dt Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
a Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
b Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
c Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
d) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
6)
        propagate Double
maxStep Double
lo Double
hi Double
y Double
v =
            let direction :: Double
direction = if Double
hi Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
>= Double
lo then Double
1 else -Double
1
                go :: (Double, Double, Double) -> t -> Double
go state :: (Double, Double, Double)
state@(Double
t, Double
_, Double
_) t
steps
                    | Double
direction Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
hi Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
t) Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
<= Double
1e-12 = let (Double
_, Double
answer, Double
_) = (Double, Double, Double)
state in Double
answer
                    | t
steps t -> t -> Bool
forall a. Ord a => a -> a -> Bool
>= t
1000000 = Double
0 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
0
                    | Bool
otherwise =
                        let dt :: Double
dt = Double
direction Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double -> Double
forall a. Ord a => a -> a -> a
min (Double -> Double
forall a. Num a => a -> a
abs (Double
hi Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
t)) (Double -> Double -> Double
forall a. Ord a => a -> a -> a
min Double
maxStep (Double
0.01 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
t))
                         in (Double, Double, Double) -> t -> Double
go (Double -> (Double, Double, Double) -> (Double, Double, Double)
advance Double
dt (Double, Double, Double)
state) (t
steps t -> t -> t
forall a. Num a => a -> a -> a
+ t
1)
             in (Double, Double, Double) -> Int -> Double
forall {t}.
(Ord t, Num t) =>
(Double, Double, Double) -> t -> Double
go (Double
lo, Double
y, Double
v) (Int
0 :: Int)
        -- Integrate from a small regular series value to avoid catastrophic
        -- cancellation of the oscillatory 1F1 representation at large rho.
        f = Double -> Double -> Double -> Double -> Double -> Double
propagate Double
0.0025 Double
start Double
rho Double
f0 Double
df0
        fCoarse = Double -> Double -> Double -> Double -> Double -> Double
propagate Double
0.005 Double
start Double
rho Double
f0 Double
df0
        r = Double
300 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Int
l Int -> Int -> Int
forall a. Num a => a -> a -> a
* Int
l) Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
eta Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
eta
        theta = Double
r Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
eta Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double
forall a. Floating a => a -> a
log (Double
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
r) Double -> Double -> Double
forall a. Num a => a -> a -> a
- Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
l Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
forall a. Floating a => a
pi Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
sigma
        gv = Double -> Double -> Double -> Double -> Double -> Double
propagate Double
0.01 Double
r Double
rho (Double -> Double
forall a. Floating a => a -> a
cos Double
theta) (-Double -> Double
forall a. Floating a => a -> a
sin Double
theta Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
eta Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
r))
    a <- estimate f (abs (f - fCoarse) / 15 + 1e-8 * abs f)
    b <- estimate gv ((1 + abs gv) * (1 + fromIntegral (l * (l + 1)) + eta * eta) / r)
    pure (a, b)

complexLogGamma :: Complex Double -> Complex Double
complexLogGamma :: Complex Double -> Complex Double
complexLogGamma Complex Double
z
    | Complex Double -> Double
forall a. Complex a -> a
realPart Complex Double
z Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
0.5 = (Complex Double -> Complex Double
forall a. Floating a => a -> a
log (Double
forall a. Floating a => a
pi Double -> Double -> Complex Double
forall a. a -> a -> Complex a
:+ Double
0) Complex Double -> Complex Double -> Complex Double
forall a. Num a => a -> a -> a
- Complex Double -> Complex Double
forall a. Floating a => a -> a
log (Complex Double -> Complex Double
forall a. Floating a => a -> a
sin (Complex Double
forall a. Floating a => a
pi Complex Double -> Complex Double -> Complex Double
forall a. Num a => a -> a -> a
* Complex Double
z))) Complex Double -> Complex Double -> Complex Double
forall a. Num a => a -> a -> a
- Complex Double -> Complex Double
complexLogGamma (Complex Double
1 Complex Double -> Complex Double -> Complex Double
forall a. Num a => a -> a -> a
- Complex Double
z)
    | Bool
otherwise =
        let cs :: [Double]
cs =
                [ Double
0.99999999999980993
                , Double
676.5203681218851
                , -Double
1259.1392167224028
                , Double
771.32342877765313
                , -Double
176.61502916214059
                , Double
12.507343278686905
                , -Double
0.13857109526572012
                , Double
9.9843695780195716e-6
                , Double
1.5056327351493116e-7
                ]
            w :: Complex Double
w = Complex Double
z Complex Double -> Complex Double -> Complex Double
forall a. Num a => a -> a -> a
- Complex Double
1
            t :: Complex Double
t = Complex Double
w Complex Double -> Complex Double -> Complex Double
forall a. Num a => a -> a -> a
+ Complex Double
7.5
            a :: Complex Double
a = ([Double]
cs [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
0 Double -> Double -> Complex Double
forall a. a -> a -> Complex a
:+ Double
0) Complex Double -> Complex Double -> Complex Double
forall a. Num a => a -> a -> a
+ [Complex Double] -> Complex Double
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum [(Double
c Double -> Double -> Complex Double
forall a. a -> a -> Complex a
:+ Double
0) Complex Double -> Complex Double -> Complex Double
forall a. Fractional a => a -> a -> a
/ (Complex Double
w Complex Double -> Complex Double -> Complex Double
forall a. Num a => a -> a -> a
+ Int -> Complex Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
k) | (Int
k, Double
c) <- [Int] -> [Double] -> [(Int, Double)]
forall a b. [a] -> [b] -> [(a, b)]
zip [Int
1 :: Int ..] (Int -> [Double] -> [Double]
forall a. Int -> [a] -> [a]
drop Int
1 [Double]
cs)]
         in Complex Double
0.5 Complex Double -> Complex Double -> Complex Double
forall a. Num a => a -> a -> a
* Complex Double -> Complex Double
forall a. Floating a => a -> a
log (Complex Double
2 Complex Double -> Complex Double -> Complex Double
forall a. Num a => a -> a -> a
* Complex Double
forall a. Floating a => a
pi) Complex Double -> Complex Double -> Complex Double
forall a. Num a => a -> a -> a
+ (Complex Double
w Complex Double -> Complex Double -> Complex Double
forall a. Num a => a -> a -> a
+ Complex Double
0.5) Complex Double -> Complex Double -> Complex Double
forall a. Num a => a -> a -> a
* Complex Double -> Complex Double
forall a. Floating a => a -> a
log Complex Double
t Complex Double -> Complex Double -> Complex Double
forall a. Num a => a -> a -> a
- Complex Double
t Complex Double -> Complex Double -> Complex Double
forall a. Num a => a -> a -> a
+ Complex Double -> Complex Double
forall a. Floating a => a -> a
log Complex Double
a

-- | Gegenbauer C_n^alpha(x), degree 0..10000 and alpha>-1/2, alpha/=0, by three-term recurrence.
gegenbauer :: Int -> Double -> Double -> Either MathError Double
gegenbauer :: Int -> Double -> Double -> Either MathError Double
gegenbauer Int
n Double
a Double
x = do
    Bool -> String -> Either MathError ()
require (Int
n Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
>= Int
0 Bool -> Bool -> Bool
&& Int
n Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
<= Int
10000 Bool -> Bool -> Bool
&& Double -> Bool
finite (Double
a Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
x) Bool -> Bool -> Bool
&& Double
a Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> -Double
0.5 Bool -> Bool -> Bool
&& Double
a Double -> Double -> Bool
forall a. Eq a => a -> a -> Bool
/= Double
0) String
"Gegenbauer n>=0,a>-1/2,a/=0"
    Double -> Either MathError Double
checked (Int -> Double -> Double -> Double
go Int
0 Double
1 (Double
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
a Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
x))
  where
    go :: Int -> Double -> Double -> Double
go Int
k Double
p Double
q
        | Int
k Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
n = Double
p
        | Bool
otherwise =
            let j :: Double
j = Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Int
k Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1) in Int -> Double -> Double -> Double
go (Int
k Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1) Double
q ((Double
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
j Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
a) Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
x Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
q Double -> Double -> Double
forall a. Num a => a -> a -> a
- (Double
j Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
a Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
1) Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
p) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
j Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
1))

-- | Real D_nu(x) via Kummer's representation; reference |x|\<=6, |nu|\<=10.
parabolicCylinder :: Double -> Double -> Either MathError SFResult
parabolicCylinder :: Double -> Double -> Either MathError SFResult
parabolicCylinder Double
nu Double
x = do
    Bool -> String -> Either MathError ()
require
        (Double -> Bool
finite (Double
nu Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
x) Bool -> Bool -> Bool
&& Double -> Double
forall a. Num a => a -> a
abs Double
nu Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
<= Double
10 Bool -> Bool -> Bool
&& Double -> Double
forall a. Num a => a -> a
abs Double
x Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
<= Double
6)
        String
"parabolic cylinder reference |nu|<=10,|x|<=6"
    let recipGamma :: Double -> Either MathError Double
recipGamma Double
y = if Double
y Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
<= Double
0 Bool -> Bool -> Bool
&& Double
y Double -> Double -> Bool
forall a. Eq a => a -> a -> Bool
== Integer -> Double
forall a. Num a => Integer -> a
fromInteger (Double -> Integer
forall b. Integral b => Double -> b
forall a b. (RealFrac a, Integral b) => a -> b
round Double
y) then Double -> Either MathError Double
forall a b. b -> Either a b
Right Double
0 else (Double
1 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/) (Double -> Double)
-> Either MathError Double -> Either MathError Double
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
<$> Double -> Either MathError Double
gamma Double
y
    a <- Double -> Either MathError Double
recipGamma ((Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
nu) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
2)
    b <- recipGamma (-nu / 2)
    u <- hypergeometric1F1 (-nu / 2) 0.5 (x * x / 2)
    v <- hypergeometric1F1 ((1 - nu) / 2) 1.5 (x * x / 2)
    let fac = Double
2 Double -> Double -> Double
forall a. Floating a => a -> a -> a
** (Double
nu Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
2) Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double
forall a. Floating a => a -> a
exp (-Double
x Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
x Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
4) Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double
forall a. Floating a => a -> a
sqrt Double
forall a. Floating a => a
pi
    estimate
        (fac * (a * sfValue u - sqrt 2 * x * b * sfValue v))
        (abs fac * (abs a * sfError u + sqrt 2 * abs (x * b) * sfError v))

{- | Mathieu ce_n (even=True) or se_n, from a 32-term Fourier Galerkin
eigenproblem. ce(0)=1 or se'(0)=1. Error estimates the last coefficient;
reference 0\<=n\<=12, |q|\<=10, any real x. Characteristic value is returned.
-}
mathieu :: Bool -> Int -> Double -> Double -> Either MathError (Double, SFResult)
mathieu :: Bool
-> Int -> Double -> Double -> Either MathError (Double, SFResult)
mathieu Bool
evenMode Int
n Double
q Double
x = do
    Bool -> String -> Either MathError ()
require
        (Int
n Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
>= Int
0 Bool -> Bool -> Bool
&& Int
n Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
<= Int
12 Bool -> Bool -> Bool
&& (Bool
evenMode Bool -> Bool -> Bool
|| Int
n Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
> Int
0) Bool -> Bool -> Bool
&& Double -> Bool
finite (Double
q Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
x) Bool -> Bool -> Bool
&& Double -> Double
forall a. Num a => a -> a
abs Double
q Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
<= Double
10)
        String
"Mathieu reference degree 0..12,|q|<=10"
    let parity :: Int
parity = Int
n Int -> Int -> Int
forall a. Integral a => a -> a -> a
`mod` Int
2
        freqs :: [Int]
freqs =
            Int -> [Int] -> [Int]
forall a. Int -> [a] -> [a]
take
                Int
32
                [ if Int
parity Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
0 Bool -> Bool -> Bool
&& Bool -> Bool
not Bool
evenMode then Int
2 else Int
parity
                , if Int
parity Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
0 Bool -> Bool -> Bool
&& Bool -> Bool
not Bool
evenMode then Int
4 else Int
parity Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
2
                ..
                ]
        norm :: Int -> a
norm Int
j = if [Int]
freqs [Int] -> Int -> Int
forall a. HasCallStack => [a] -> Int -> a
!! Int
j Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
0 then a -> a
forall a. Floating a => a -> a
sqrt a
2 else a
1
        entry :: Int -> Int -> Double
entry Int
i Int
j
            | Int
i Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
j =
                Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral (([Int]
freqs [Int] -> Int -> Int
forall a. HasCallStack => [a] -> Int -> a
!! Int
i) Int -> Int -> Int
forall a b. (Num a, Integral b) => a -> b -> a
^ (Int
2 :: Int))
                    Double -> Double -> Double
forall a. Num a => a -> a -> a
+ (if [Int]
freqs [Int] -> Int -> Int
forall a. HasCallStack => [a] -> Int -> a
!! Int
i Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
1 then if Bool
evenMode then Double
q else -Double
q else Double
0)
            | Int -> Int
forall a. Num a => a -> a
abs (Int
i Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
j) Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
1 = Double
q Double -> Double -> Double
forall a. Num a => a -> a -> a
* (if [Int]
freqs [Int] -> Int -> Int
forall a. HasCallStack => [a] -> Int -> a
!! Int -> Int -> Int
forall a. Ord a => a -> a -> a
min Int
i Int
j Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
0 then Double -> Double
forall a. Floating a => a -> a
sqrt Double
2 else Double
1)
            | Bool
otherwise = Double
0
    mat <- Int -> Int -> [Double] -> Either MathError (Matrix Double)
forall a. Int -> Int -> [a] -> Either MathError (Matrix a)
fromListChecked Int
32 Int
32 [Int -> Int -> Double
entry Int
i Int
j | Int
i <- [Int
0 .. Int
31], Int
j <- [Int
0 .. Int
31]]
    eig <- symmetricEigen mat
    let which = if Bool
evenMode then Int
n Int -> Int -> Int
forall a. Integral a => a -> a -> a
`div` Int
2 else (Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1) Int -> Int -> Int
forall a. Integral a => a -> a -> a
`div` Int
2
        coeff = [Complex Double -> Double
forall a. Complex a -> a
realPart (Matrix (Complex Double) -> [Complex Double]
forall a. Matrix a -> [a]
toList (Eigen -> Matrix (Complex Double)
eigenvectors Eigen
eig) [Complex Double] -> Int -> Complex Double
forall a. HasCallStack => [a] -> Int -> a
!! (Int
32 Int -> Int -> Int
forall a. Num a => a -> a -> a
* Int
j Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
which)) | Int
j <- [Int
0 .. Int
31]]
        basis Int
j =
            if Bool
evenMode
                then Double -> Double
forall a. Floating a => a -> a
cos (Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral ([Int]
freqs [Int] -> Int -> Int
forall a. HasCallStack => [a] -> Int -> a
!! Int
j) Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
x) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Int -> Double
forall {a}. Floating a => Int -> a
norm Int
j
                else Double -> Double
forall a. Floating a => a -> a
sin (Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral ([Int]
freqs [Int] -> Int -> Int
forall a. HasCallStack => [a] -> Int -> a
!! Int
j) Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
x)
        denom =
            if Bool
evenMode
                then [Double] -> Double
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum [[Double]
coeff [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
j Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Int -> Double
forall {a}. Floating a => Int -> a
norm Int
j | Int
j <- [Int
0 .. Int
31]]
                else [Double] -> Double
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum [[Double]
coeff [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
j Double -> Double -> Double
forall a. Num a => a -> a -> a
* Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral ([Int]
freqs [Int] -> Int -> Int
forall a. HasCallStack => [a] -> Int -> a
!! Int
j) | Int
j <- [Int
0 .. Int
31]]
        v = [Double] -> Double
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum [[Double]
coeff [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
j Double -> Double -> Double
forall a. Num a => a -> a -> a
* Int -> Double
basis Int
j | Int
j <- [Int
0 .. Int
31]] Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
denom
    value <- estimate v (abs (last coeff / denom))
    pure (realPart (eigenvalues eig !! which), value)

{- | Exact rational Bernoulli/Euler coefficients evaluated as Double. Bernoulli
B1=-1/2; Euler E_n(x) defined by 2 exp(xt)/(exp t+1).
-}
bernoulliPolynomial :: Int -> Double -> Either MathError Double
bernoulliPolynomial :: Int -> Double -> Either MathError Double
bernoulliPolynomial Int
n Double
x = do
    Bool -> String -> Either MathError ()
require (Int
n Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
>= Int
0 Bool -> Bool -> Bool
&& Int
n Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
<= Int
100 Bool -> Bool -> Bool
&& Double -> Bool
finite Double
x) String
"Bernoulli polynomial degree 0..100"
    Double -> Either MathError Double
checked ([Double] -> Double
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum [Integer -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Int -> Int -> Integer
choose Int
n Int
k) Double -> Double -> Double
forall a. Num a => a -> a -> a
* Rational -> Double
forall a. Fractional a => Rational -> a
fromRational ([Rational]
bs [Rational] -> Int -> Rational
forall a. HasCallStack => [a] -> Int -> a
!! Int
k) Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
x Double -> Int -> Double
forall a b. (Num a, Integral b) => a -> b -> a
^ (Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
k) | Int
k <- [Int
0 .. Int
n]])
  where
    bs :: [Rational]
bs =
        Rational
1
            Rational -> [Rational] -> [Rational]
forall a. a -> [a] -> [a]
: [ -[Rational] -> Rational
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum [Integer -> Rational
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Int -> Int -> Integer
choose (Int
m Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1) Int
k) Rational -> Rational -> Rational
forall a. Num a => a -> a -> a
* ([Rational]
bs [Rational] -> Int -> Rational
forall a. HasCallStack => [a] -> Int -> a
!! Int
k) | Int
k <- [Int
0 .. Int
m Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1]] Rational -> Rational -> Rational
forall a. Fractional a => a -> a -> a
/ Int -> Rational
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Int
m Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1)
              | Int
m <- [Int
1 ..]
              ] ::
            [Rational]

-- | Euler E_n(x), n=0..100, from the Bernoulli duplication identity.
eulerPolynomial :: Int -> Double -> Either MathError Double
eulerPolynomial :: Int -> Double -> Either MathError Double
eulerPolynomial Int
n Double
x = do
    Bool -> String -> Either MathError ()
require (Int
n Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
>= Int
0 Bool -> Bool -> Bool
&& Int
n Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
<= Int
100 Bool -> Bool -> Bool
&& Double -> Bool
finite Double
x) String
"Euler polynomial degree 0..100"
    a <- Int -> Double -> Either MathError Double
bernoulliPolynomial (Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1) Double
x
    b <- bernoulliPolynomial (n + 1) (x / 2)
    checked (2 / fromIntegral (n + 1) * (a - 2 ** fromIntegral (n + 1) * b))

choose :: Int -> Int -> Integer
choose :: Int -> Int -> Integer
choose Int
n Int
k
    | Int
k Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< Int
0 Bool -> Bool -> Bool
|| Int
k Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
> Int
n = Integer
0
    | Bool
otherwise =
        [Integer] -> Integer
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
product [Int -> Integer
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
k Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1) .. Int -> Integer
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
n] Integer -> Integer -> Integer
forall a. Integral a => a -> a -> a
`div` [Integer] -> Integer
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
product [Integer
1 .. Int -> Integer
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
k]

{- | Angular momenta and projections are encoded as doubled integers, allowing
half-integers exactly. Factorial Racah sums have j\<=30 and return selection
rule zeros rather than factorials of negative arguments.
-}
wigner3j :: Int -> Int -> Int -> Int -> Int -> Int -> Either MathError Double
wigner3j :: Int -> Int -> Int -> Int -> Int -> Int -> Either MathError Double
wigner3j Int
j1 Int
j2 Int
j3 Int
m1 Int
m2 Int
m3 = do
    Bool -> String -> Either MathError ()
require ((Int -> Bool) -> [Int] -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
all (\Int
j -> Int
j Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
>= Int
0 Bool -> Bool -> Bool
&& Int
j Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
<= Int
60) [Int
j1, Int
j2, Int
j3]) String
"Wigner doubled j must lie in 0..60"
    if Bool -> Bool
not (Int -> Int -> Int -> Bool
triangle Int
j1 Int
j2 Int
j3)
        Bool -> Bool -> Bool
|| Int
m1 Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
m2 Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
m3 Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
/= Int
0
        Bool -> Bool -> Bool
|| [Bool] -> Bool
forall (t :: * -> *). Foldable t => t Bool -> Bool
or [Int -> Int
forall a. Num a => a -> a
abs Int
m Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
> Int
j Bool -> Bool -> Bool
|| Int -> Bool
forall a. Integral a => a -> Bool
odd (Int
j Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
m) | (Int
j, Int
m) <- [Int] -> [Int] -> [(Int, Int)]
forall a b. [a] -> [b] -> [(a, b)]
zip [Int
j1, Int
j2, Int
j3] [Int
m1, Int
m2, Int
m3]]
        then Double -> Either MathError Double
forall a b. b -> Either a b
Right Double
0
        else do
            let a :: Int
a = (Int
j1 Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
j2 Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
j3) Int -> Int -> Int
forall a. Integral a => a -> a -> a
`div` Int
2
                b :: Int
b = (Int
j1 Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
m1) Int -> Int -> Int
forall a. Integral a => a -> a -> a
`div` Int
2
                c :: Int
c = (Int
j2 Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
m2) Int -> Int -> Int
forall a. Integral a => a -> a -> a
`div` Int
2
                d :: Int
d = (Int
j3 Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
j2 Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
m1) Int -> Int -> Int
forall a. Integral a => a -> a -> a
`div` Int
2
                e :: Int
e = (Int
j3 Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
j1 Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
m2) Int -> Int -> Int
forall a. Integral a => a -> a -> a
`div` Int
2
                lo :: Int
lo = [Int] -> Int
forall a. Ord a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Ord a) => t a -> a
maximum [Int
0, -Int
d, -Int
e]
                hi :: Int
hi = [Int] -> Int
forall a. Ord a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Ord a) => t a -> a
minimum [Int
a, Int
b, Int
c]
                sumv :: Rational
sumv =
                    [Rational] -> Rational
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum
                        [ (if Int -> Bool
forall a. Integral a => a -> Bool
even Int
z then Rational
1 else -Rational
1)
                            Rational -> Rational -> Rational
forall a. Fractional a => a -> a -> a
/ Integer -> Rational
forall a. Num a => Integer -> a
fromInteger
                                ( Int -> Integer
integerFact Int
z
                                    Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
* Int -> Integer
integerFact (Int
a Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
z)
                                    Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
* Int -> Integer
integerFact (Int
b Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
z)
                                    Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
* Int -> Integer
integerFact (Int
c Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
z)
                                    Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
* Int -> Integer
integerFact (Int
d Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
z)
                                    Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
* Int -> Integer
integerFact (Int
e Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
z)
                                )
                        | Int
z <- [Int
lo .. Int
hi]
                        ] ::
                        Rational
                pref2 :: Rational
pref2 =
                    Int -> Int -> Int -> Rational
deltaR Int
j1 Int
j2 Int
j3
                        Rational -> Rational -> Rational
forall a. Num a => a -> a -> a
* Integer -> Rational
forall a. Num a => Integer -> a
fromInteger
                            ( [Integer] -> Integer
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
product
                                [ Int -> Integer
integerFact ((Int
j Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
m) Int -> Int -> Int
forall a. Integral a => a -> a -> a
`div` Int
2) Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
* Int -> Integer
integerFact ((Int
j Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
m) Int -> Int -> Int
forall a. Integral a => a -> a -> a
`div` Int
2)
                                | (Int
j, Int
m) <- [Int] -> [Int] -> [(Int, Int)]
forall a b. [a] -> [b] -> [(a, b)]
zip [Int
j1, Int
j2, Int
j3] [Int
m1, Int
m2, Int
m3]
                                ]
                            )
                value :: Double
value =
                    Int -> Double
phase ((Int
j1 Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
j2 Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
m3) Int -> Int -> Int
forall a. Integral a => a -> a -> a
`div` Int
2)
                        Double -> Double -> Double
forall a. Num a => a -> a -> a
* (if Rational
sumv Rational -> Rational -> Bool
forall a. Ord a => a -> a -> Bool
< Rational
0 then -Double
1 else Double
1)
                        Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double
forall a. Floating a => a -> a
sqrt (Rational -> Double
forall a. Fractional a => Rational -> a
fromRational (Rational
pref2 Rational -> Rational -> Rational
forall a. Num a => a -> a -> a
* Rational
sumv Rational -> Rational -> Rational
forall a. Num a => a -> a -> a
* Rational
sumv))
            Double -> Either MathError Double
checked Double
value

phase :: Int -> Double
phase :: Int -> Double
phase Int
n = if Int -> Bool
forall a. Integral a => a -> Bool
even Int
n then Double
1 else -Double
1
triangle :: Int -> Int -> Int -> Bool
triangle :: Int -> Int -> Int -> Bool
triangle Int
a Int
b Int
c = Int
a Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
b Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
>= Int
c Bool -> Bool -> Bool
&& Int
a Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
c Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
>= Int
b Bool -> Bool -> Bool
&& Int
b Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
c Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
>= Int
a Bool -> Bool -> Bool
&& Int -> Bool
forall a. Integral a => a -> Bool
even (Int
a Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
b Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
c)
integerFact :: Int -> Integer
integerFact :: Int -> Integer
integerFact Int
n = [Integer] -> Integer
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
product [Integer
1 .. Int -> Integer
forall a. Integral a => a -> Integer
toInteger Int
n]
deltaR :: Int -> Int -> Int -> Rational
deltaR :: Int -> Int -> Int -> Rational
deltaR Int
a Int
b Int
c =
    Integer -> Rational
forall a. Num a => Integer -> a
fromInteger
        ( Int -> Integer
integerFact ((Int
a Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
b Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
c) Int -> Int -> Int
forall a. Integral a => a -> a -> a
`div` Int
2)
            Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
* Int -> Integer
integerFact ((Int
a Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
b Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
c) Int -> Int -> Int
forall a. Integral a => a -> a -> a
`div` Int
2)
            Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
* Int -> Integer
integerFact ((-Int
a Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
b Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
c) Int -> Int -> Int
forall a. Integral a => a -> a -> a
`div` Int
2)
        )
        Rational -> Rational -> Rational
forall a. Fractional a => a -> a -> a
/ Integer -> Rational
forall a. Num a => Integer -> a
fromInteger (Int -> Integer
integerFact ((Int
a Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
b Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
c) Int -> Int -> Int
forall a. Integral a => a -> a -> a
`div` Int
2 Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1))

-- | Racah factorial sum for a Wigner 6j symbol; each angular momentum is a doubled integer in 0..60.
wigner6j :: Int -> Int -> Int -> Int -> Int -> Int -> Either MathError Double
wigner6j :: Int -> Int -> Int -> Int -> Int -> Int -> Either MathError Double
wigner6j Int
a Int
b Int
c Int
d Int
e Int
f = do
    Bool -> String -> Either MathError ()
require ((Int -> Bool) -> [Int] -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
all (\Int
j -> Int
j Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
>= Int
0 Bool -> Bool -> Bool
&& Int
j Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
<= Int
60) [Int
a, Int
b, Int
c, Int
d, Int
e, Int
f]) String
"Wigner doubled j 0..60"
    if Bool -> Bool
not ([Bool] -> Bool
forall (t :: * -> *). Foldable t => t Bool -> Bool
and [Int -> Int -> Int -> Bool
triangle Int
a Int
b Int
c, Int -> Int -> Int -> Bool
triangle Int
a Int
e Int
f, Int -> Int -> Int -> Bool
triangle Int
d Int
b Int
f, Int -> Int -> Int -> Bool
triangle Int
d Int
e Int
c])
        then Double -> Either MathError Double
forall a b. b -> Either a b
Right Double
0
        else do
            let low :: [Int]
low = (Int -> Int) -> [Int] -> [Int]
forall a b. (a -> b) -> [a] -> [b]
map (Int -> Int -> Int
forall a. Integral a => a -> a -> a
`div` Int
2) [Int
a Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
b Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
c, Int
a Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
e Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
f, Int
d Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
b Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
f, Int
d Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
e Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
c]
                high :: [Int]
high = (Int -> Int) -> [Int] -> [Int]
forall a b. (a -> b) -> [a] -> [b]
map (Int -> Int -> Int
forall a. Integral a => a -> a -> a
`div` Int
2) [Int
a Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
b Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
d Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
e, Int
a Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
c Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
d Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
f, Int
b Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
c Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
e Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
f]
                s :: Rational
s =
                    [Rational] -> Rational
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum
                        [ (if Int -> Bool
forall a. Integral a => a -> Bool
even Int
z then Rational
1 else -Rational
1)
                            Rational -> Rational -> Rational
forall a. Num a => a -> a -> a
* Integer -> Rational
forall a. Num a => Integer -> a
fromInteger (Int -> Integer
integerFact (Int
z Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1))
                            Rational -> Rational -> Rational
forall a. Fractional a => a -> a -> a
/ Integer -> Rational
forall a. Num a => Integer -> a
fromInteger ([Integer] -> Integer
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
product [Int -> Integer
integerFact (Int
z Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
t) | Int
t <- [Int]
low] Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
* [Integer] -> Integer
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
product [Int -> Integer
integerFact (Int
t Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
z) | Int
t <- [Int]
high])
                        | Int
z <- [[Int] -> Int
forall a. Ord a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Ord a) => t a -> a
maximum [Int]
low .. [Int] -> Int
forall a. Ord a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Ord a) => t a -> a
minimum [Int]
high]
                        ] ::
                        Rational
                squared :: Rational
squared = Int -> Int -> Int -> Rational
deltaR Int
a Int
b Int
c Rational -> Rational -> Rational
forall a. Num a => a -> a -> a
* Int -> Int -> Int -> Rational
deltaR Int
a Int
e Int
f Rational -> Rational -> Rational
forall a. Num a => a -> a -> a
* Int -> Int -> Int -> Rational
deltaR Int
d Int
b Int
f Rational -> Rational -> Rational
forall a. Num a => a -> a -> a
* Int -> Int -> Int -> Rational
deltaR Int
d Int
e Int
c Rational -> Rational -> Rational
forall a. Num a => a -> a -> a
* Rational
s Rational -> Rational -> Rational
forall a. Num a => a -> a -> a
* Rational
s
            Double -> Either MathError Double
checked ((if Rational
s Rational -> Rational -> Bool
forall a. Ord a => a -> a -> Bool
< Rational
0 then -Double
1 else Double
1) Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double
forall a. Floating a => a -> a
sqrt (Rational -> Double
forall a. Fractional a => Rational -> a
fromRational Rational
squared))

-- | Wigner 9j as a finite sum of three 6j symbols; accepts nine row-major doubled angular momenta in 0..30.
wigner9j :: [Int] -> Either MathError Double
wigner9j :: [Int] -> Either MathError Double
wigner9j [Int]
js = do
    Bool -> String -> Either MathError ()
require
        ([Int] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Int]
js Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
9 Bool -> Bool -> Bool
&& (Int -> Bool) -> [Int] -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
all (\Int
j -> Int
j Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
>= Int
0 Bool -> Bool -> Bool
&& Int
j Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
<= Int
30) [Int]
js)
        String
"9j needs nine doubled angular momenta 0..30"
    let a :: Int
a = [Int]
js [Int] -> Int -> Int
forall a. HasCallStack => [a] -> Int -> a
!! Int
0
        b :: Int
b = [Int]
js [Int] -> Int -> Int
forall a. HasCallStack => [a] -> Int -> a
!! Int
1
        c :: Int
c = [Int]
js [Int] -> Int -> Int
forall a. HasCallStack => [a] -> Int -> a
!! Int
2
        d :: Int
d = [Int]
js [Int] -> Int -> Int
forall a. HasCallStack => [a] -> Int -> a
!! Int
3
        e :: Int
e = [Int]
js [Int] -> Int -> Int
forall a. HasCallStack => [a] -> Int -> a
!! Int
4
        f :: Int
f = [Int]
js [Int] -> Int -> Int
forall a. HasCallStack => [a] -> Int -> a
!! Int
5
        g :: Int
g = [Int]
js [Int] -> Int -> Int
forall a. HasCallStack => [a] -> Int -> a
!! Int
6
        h :: Int
h = [Int]
js [Int] -> Int -> Int
forall a. HasCallStack => [a] -> Int -> a
!! Int
7
        i :: Int
i = [Int]
js [Int] -> Int -> Int
forall a. HasCallStack => [a] -> Int -> a
!! Int
8
    terms <-
        (Int -> Either MathError Double)
-> [Int] -> Either MathError [Double]
forall (t :: * -> *) (m :: * -> *) a b.
(Traversable t, Monad m) =>
(a -> m b) -> t a -> m (t b)
forall (m :: * -> *) a b. Monad m => (a -> m b) -> [a] -> m [b]
mapM
            ( \Int
x -> do
                u <- Int -> Int -> Int -> Int -> Int -> Int -> Either MathError Double
wigner6j Int
a Int
d Int
g Int
h Int
i Int
x
                v <- wigner6j b e h d x f
                w <- wigner6j c f i x a b
                pure (phase x * fromIntegral (x + 1) * u * v * w)
            )
            [Int
0 .. Int
60]
    checked (sum terms)

-- | Clebsch--Gordan coefficient \<j1 m1,j2 m2|j m> through the Condon--Shortley 3j convention; all arguments doubled integers.
clebschGordan :: Int -> Int -> Int -> Int -> Int -> Int -> Either MathError Double
clebschGordan :: Int -> Int -> Int -> Int -> Int -> Int -> Either MathError Double
clebschGordan Int
j1 Int
m1 Int
j2 Int
m2 Int
j Int
m = do
    v <- Int -> Int -> Int -> Int -> Int -> Int -> Either MathError Double
wigner3j Int
j1 Int
j2 Int
j Int
m1 Int
m2 (-Int
m)
    checked (phase ((j1 - j2 + m) `div` 2) * sqrt (fromIntegral (j + 1)) * v)