module Math.Gamma.Incomplete
( lowerGammaCF, pCF
, lowerGammaHypGeom, lnLowerGammaHypGeom, pHypGeom
, upperGammaCF, lnUpperGammaConvergents, qCF
) where
import Math.Gamma
import Math.ContinuedFraction
import Math.Sequence.Converge
lowerGammaCF :: (Floating a, Ord a) => a -> a -> Math.ContinuedFraction.CF a
lowerGammaCF s z = gcf 0
[ (p,q)
| p <- pow_x_s_div_exp_x s z
: interleave
[negate spn * z | spn <- iterate (1+) s]
[n * z | n <- iterate (1+) 1]
| q <- iterate (1+) s
]
lowerGammaHypGeom :: Floating b => b -> b -> b
lowerGammaHypGeom 0 0 = 0/0
lowerGammaHypGeom s x = x ** s * exp (negate x) / s * m_1_sp1 s x
lnLowerGammaHypGeom :: Floating a => a -> a -> a
lnLowerGammaHypGeom 0 0 = 0/0
lnLowerGammaHypGeom s x
= log ((signum x)**s * sign_m / signum s)
+ s*log (abs x) x log (abs s) + log_m
where
(sign_m, log_m) = log_m_1_sp1 s x
pCF :: (Gamma a, Ord a, Enum a) => a -> a -> CF a
pCF s x = gcf 0
[ (p,q)
| p <- pow_x_s_div_gamma_s_div_exp_x s x
: interleave
[negate spn * x | spn <- [s..]]
[n * x | n <- [1..]]
| q <- [s..]
]
pHypGeom :: (Gamma a, Ord a) => a -> a -> a
pHypGeom 0 0 = 0/0
pHypGeom s x
| s < 0
= signum x ** s * sin (pi*s) / (pi)
* exp (s * log (abs x) x + lnGamma (s)) * m_1_sp1 s x
| s == 0 || x == 0
= 0
| otherwise
= signum x ** s * exp (s * log (abs x) x lnGamma (s+1)) * m_1_sp1 s x
qCF :: (Gamma a, Ord a, Enum a) => a -> a -> CF a
qCF s x = gcf 0
[ (p,q)
| p <- pow_x_s_div_gamma_s_div_exp_x s x
: zipWith (*) [1..] (iterate (subtract 1) (s1))
| q <- [n + x s | n <- [1,3..]]
]
upperGammaCF :: (Floating a, Ord a) => a -> a -> CF a
upperGammaCF s z = gcf 0
[ (p,q)
| p <- pow_x_s_div_exp_x s z
: zipWith (*) (iterate (1+) 1) (iterate (subtract 1) (s1))
| q <- [n + z s | n <- iterate (2+) 1]
]
lnUpperGammaConvergents :: Floating a => a -> a -> [a]
lnUpperGammaConvergents s x = map (a ) (concat (eval theCF))
where
eval = map (map evalSign) . modifiedLentzWith signLog addSignLog negateSignLog 1e-30
a = s * log x x
theCF = gcf (x + 1 s)
[ (p,q)
| p <- zipWith (*) (iterate (1+) 1) (iterate (subtract 1) (s1))
| q <- [n + x s | n <- iterate (2+) 3]
]
evalSign (s,x) = log s + x
signLog x = (signum x, log (abs x))
addSignLog (xS,xL) (yS,yL) = (xS*yS, xL+yL)
negateSignLog (s,l) = (s, negate l)
m_1_sp1 s z = converge . scanl (+) 0 . scanl (*) 1 $
[z / x | x <- iterate (1+) (s+1)]
log_m_1_sp1 s z = converge (concat (log_m_1_sp1_convergents s z))
log_m_1_sp1_convergents s z
= modifiedLentzWith signLog addSignLog negateSignLog 1e-30
$ sumPartialProducts (1:[z / x | x <- iterate (1+) (s+1)])
interleave [] _ = []
interleave _ [] = []
interleave (x:xs) ys = x:interleave ys xs
pow_x_s_div_gamma_s_div_exp_x s x
| x > 0 = exp (log x * s x lnGamma s)
| otherwise = x ** s / (exp x * gamma s)
pow_x_s_div_exp_x s x
| x > 0 = exp (log x * s x)
| otherwise = x ** s / exp x