Improve precision for incomplete gamma function

Properly evaluate constant factor before Taylor series for P. This markedly improves precision. We're still quite far from 1-ulp precision but at very least we don't lose 5(!) digits of precision in worst case Fixes #42

Author: Shimuuar

Numeric/SpecFunctions/Internal.hs

@@ -409,24 +409,18 @@ incompleteGamma a x
  || (abs mu < 0.4)
  -- Gautschi's algorithm.
  --
- -- Evaluate series for P(a,x). See [Temme1994] Eq. 5.5
- --
- -- Evaluation of constant multiplier is somewhat tricky. On one
- -- hand working in linear domain gives markedly better precision
- -- than log domain (in log domain we can easily lose 2-3
- -- digits). On other it's prone to overflow: we're working with
- -- exponents and gamma function here!. Solution is to try to
- -- compute value normally and revert to log domain it anything
- -- overflowed
+ -- Evaluate series for P(a,x). See [Temme1994] Eq. 5.5 and [NOTE:
+ -- incompleteGamma.taylorP]
  factorP
- | isInfinite num = logDom
- | isInfinite denom = logDom
- | otherwise = num / denom
- where
- num = x ** a
- denom = exp x * exp (logGamma (a + 1))
- logDom = exp (log x * a - x - logGamma (a+1))
-
+ | a < 10 = x ** a
+ / (exp x * exp (logGamma (a + 1)))
+ | a < 1182.5 = (x * exp 1 / a) ** a
+ / exp x
+ / sqrt (2*pi*a)
+ / exp (logGammaCorrection a)
+ | otherwise = (x * exp 1 / a * exp (-x/a)) ** a
+ / sqrt (2*pi*a)
+ / exp (logGammaCorrection a)
  taylorSeriesP
  = sumPowerSeries x (scanSequence (/) 1 $ enumSequenceFrom (a+1))
  * factorP
@@ -1337,3 +1331,53 @@ factorialTable = U.fromListN 171
  , 4.269068009004705e304
  , 7.257415615307998e306
  ]
+
+
+-- [NOTE: incompleteGamma.taylorP]
+--
+-- Incompltete gamma uses several algorithms for different parts of
+-- parameter space. Most troublesome is P(a,x) Taylor series
+-- [Temme1994,Eq.5.5] which requires to evaluate rather nasty
+-- expression:
+--
+-- x^a x^a
+-- ------------- = -------------
+-- exp(x)·Γ(a+1) exp(x)·a·Γ(a)
+--
+-- Conditions:
+-- | 0.5<x<1.1 = x < 4/3*a
+-- | otherwise = x < a
+--
+-- For small `a` computation could be performed directly. However for
+-- largish values of `a` it's possible some of factor in the
+-- expression overflow. Values below take into account ranges for
+-- Taylor P approximation:
+--
+-- · a > 155 - x^a could overflow
+-- · a > 1182.5 - exp(x) could overflow
+--
+-- Usual way to avoid overflow problem is to perform calculations in
+-- the log domain. It however doesn't work very well in this case
+-- since we encounter catastrophic cancellations and could easily lose
+-- up to 6(!) digits for large `a`.
+--
+-- So we take another approach and use Stirling approximation with
+-- correction (logGammaCorrection).
+--
+-- x^a / x·e \^a 1
+-- ≈ ------------------------- = | --- | · ----------------
+-- exp(x)·sqrt(2πa)·(a/e)^a) \ a / exp(x)·sqrt(2πa)
+--
+-- We're using this approach as soon as logGammaCorrection starts
+-- working (a>10) because we don't have implementation for gamma
+-- function and exp(logGamma z) results in errors for large a.
+--
+-- Once we get into region when exp(x) could overflow we rewrite
+-- expression above once more:
+--
+-- / x·e \^a 1
+-- | --- · e^(-x/a) | · ---------
+-- \ a / sqrt(2πa)
+--
+-- This approach doesn't work very well but it's still big improvement
+-- over calculations in the log domain.

tests/Tests/SpecFunctions.hs

@@ -80,7 +80,10 @@ tests = testGroup "Special functions"
  , testGroup "incomplete gamma"
  [ testCase "incompleteGamma table" $
  forTable "tests/tables/igamma.dat" $ \[a,x,exact] -> do
- checkTabularPure 16 (show (a,x)) exact (incompleteGamma a x)
+ let err | a < 10 = 16
+ | a == 201 = 200
+ | otherwise = 32
+ checkTabularPure err (show (a,x)) exact (incompleteGamma a x)
  , testProperty "incomplete gamma - increases" $
  \(abs -> s) (abs -> x) (abs -> y) -> s > 0 ==> monotonicallyIncreases (incompleteGamma s) x y
  , testProperty "0 <= gamma <= 1" incompleteGammaInRange

tests/tables/generate.py

@@ -1,4 +1,4 @@
-#!/usr/bin/python
+#!/usr/bin/env python3
 """
 """
 import itertools

tests/tables/igamma.dat

@@ -7,6 +7,7 @@
 0.000100000000000000005	10	0.999999999584179852012789127241
 0.000100000000000000005	100	1.0
 0.000100000000000000005	1000	1.0
+0.000100000000000000005	3301	1.0
 0.00100000000000000002	9.99999999999999955e-07	0.986848133694076665339510700432
 0.00100000000000000002	1.00000000000000008e-05	0.989123044695782668850940465364
 0.00100000000000000002	0.00100000000000000002	0.993687646708860290096219637037
@@ -16,6 +17,7 @@
 0.00100000000000000002	10	0.999999995830692182809738396874
 0.00100000000000000002	100	1.0
 0.00100000000000000002	1000	1.0
+0.00100000000000000002	3301	1.0
 0.0100000000000000002	9.99999999999999955e-07	0.875933759832353305070970748339
 0.0100000000000000002	1.00000000000000008e-05	0.896336798267197200971193210546
 0.0100000000000000002	0.00100000000000000002	0.938570652526128985382985766694
@@ -25,6 +27,7 @@
 0.0100000000000000002	10	0.99999995718295022590061122615
 0.0100000000000000002	100	1.0
 0.0100000000000000002	1000	1.0
+0.0100000000000000002	3301	1.0
 0.100000000000000006	9.99999999999999955e-07	0.264033654327922324857187514752
 0.100000000000000006	1.00000000000000008e-05	0.332398405040503295401903974218
 0.100000000000000006	0.00100000000000000002	0.526768568392445111817869842601
@@ -34,6 +37,7 @@
 0.100000000000000006	10	0.999999445201428209809392042838
 0.100000000000000006	100	1.0
 0.100000000000000006	1000	1.0
+0.100000000000000006	3301	1.0
 0.200000000000000011	9.99999999999999955e-07	0.0687190937987684780583487585964
 0.200000000000000011	1.00000000000000008e-05	0.108912260585591831357806184556
 0.200000000000000011	0.00100000000000000002	0.273530102033034019134013607911
@@ -43,6 +47,7 @@
 0.200000000000000011	10	0.999998540143010797930830697127
 0.200000000000000011	100	1.0
 0.200000000000000011	1000	1.0
+0.200000000000000011	3301	1.0
 0.299999999999999989	9.99999999999999955e-07	0.0176595495901936665671778397068
 0.299999999999999989	1.00000000000000008e-05	0.0352353606155625769158079986437
 0.299999999999999989	0.00100000000000000002	0.140242458924867370902963425856
@@ -52,6 +57,7 @@
 0.299999999999999989	10	0.99999715515533278951644372428
 0.299999999999999989	100	1.0
 0.299999999999999989	1000	1.0
+0.299999999999999989	3301	1.0
 0.400000000000000022	9.99999999999999955e-07	0.00448690737698480048018769908146
 0.400000000000000022	1.00000000000000008e-05	0.0112705727782256817005130809278
 0.400000000000000022	0.00100000000000000002	0.0710923978953330760218610475958
@@ -61,6 +67,7 @@
 0.400000000000000022	10	0.999995127544578487259167603836
 0.400000000000000022	100	1.0
 0.400000000000000022	1000	1.0
+0.400000000000000022	3301	1.0
 0.5	9.99999999999999955e-07	0.00112837879096923635441785924383
 0.5	1.00000000000000008e-05	0.00356823633818045042058077366094
 0.5	0.00100000000000000002	0.0356705917296798854171108747554
@@ -70,6 +77,7 @@
 0.5	10	0.999992255783568955916362323619
 0.5	100	1.0
 0.5	1000	1.0
+0.5	3301	1.0
 0.599999999999999978	9.99999999999999955e-07	0.000281123932739927969642852394278
 0.599999999999999978	1.00000000000000008e-05	0.00111917075717695837092315556803
 0.599999999999999978	0.00100000000000000002	0.0177310780570833845378940622156
@@ -79,6 +87,7 @@
 0.599999999999999978	10	0.999988293084421631090774556099
 0.599999999999999978	100	1.0
 0.599999999999999978	1000	1.0
+0.599999999999999978	3301	1.0
 2	9.99999999999999955e-07	4.9999966666679162138149041803e-13
 2	1.00000000000000008e-05	4.99996666679166715135638665241e-11
 2	0.00100000000000000002	0.000000499666791633340297383350611252
@@ -88,6 +97,7 @@
 2	10	0.999500600772612666633108493329
 2	100	1.0
 2	1000	1.0
+2	3301	1.0
 3	9.99999999999999955e-07	1.6666654166671664402685631963e-19
 3	1.00000000000000008e-05	1.66665416671666693678925483422e-16
 3	0.00100000000000000002	1.66541716652780763845435502936e-10
@@ -97,6 +107,7 @@
 3	10	0.997230604284488424056328917551
 3	100	1.0
 3	1000	1.0
+3	3301	1.0
 6	9.99999999999999955e-07	1.38888769841321886877873406008e-39
 6	1.00000000000000008e-05	1.38887698417906798768211894792e-33
 6	0.00100000000000000002	1.38769893337745993330072800382e-21
@@ -106,6 +117,7 @@
 6	10	0.932914037120968217714240937173
 6	100	1.0
 6	1000	1.0
+6	3301	1.0
 12	9.99999999999999955e-07	2.08767377170244306355144979982e-81
 12	1.00000000000000008e-05	2.08765642802367929860815380041e-69
 12	0.00100000000000000002	2.08574950796625691581696824468e-45
@@ -115,6 +127,7 @@
 12	10	0.30322385369689331180582927404
 12	100	0.99999999999999999999999999999
 12	1000	1.0
+12	3301	1.0
 18	9.99999999999999955e-07	1.56191921714497984569652172566e-124
 18	1.00000000000000008e-05	1.56190589978546723049068336087e-106
 18	0.00100000000000000002	1.56044168515546614690893456571e-70
@@ -124,6 +137,7 @@
 18	10	0.0142776135970496129089830965836
 18	100	0.999999999999999999999998742916
 18	1000	1.0
+18	3301	1.0
 101	9.99999999999999955e-07	1.06090022487198225461869933883e-766
 101	1.00000000000000008e-05	1.06089077042094832542309156458e-665
 101	0.00100000000000000002	1.05985129506822713235416044578e-463
@@ -133,6 +147,7 @@
 101	10	5.33940546071971052337450543217e-64
 101	100	0.473437801470001529623393607105
 101	1000	1.0
+101	3301	1.0
 201	9.99999999999999955e-07	6.30833677503352627325382205745e-1584
 201	1.00000000000000008e-05	6.30828028132019299365575041933e-1383
 201	0.00100000000000000002	6.30206906057329746426353105998e-981
@@ -142,6 +157,7 @@
 201	10	3.01310889066541622340100430602e-181
 201	100	4.62617947019577288096442044516e-19
 201	1000	1.0
+201	3301	1.0
 1000	9.99999999999999955e-07	2.4851656605824546988202041697e-8568
 1000	1.00000000000000008e-05	2.48514331653642003877116895294e-7568
 1000	0.00100000000000000002	2.48268669750003876617472467308e-5568
@@ -151,3 +167,14 @@
 1000	10	1.13964958763725134069338684487e-1572
 1000	100	1.0270971815582277877665615585e-611
 1000	1000	0.504205244180215508503777843602
+1000	3301	1.0
+3003	9.99999999999999955e-07	8.90812855529069811196158732723e-27160
+3003	1.00000000000000008e-05	8.90804840918643793907789460332e-24157
+3003	0.00100000000000000002	8.89923673804629728415411068459e-18151
+3003	0.0100000000000000002	8.81952937102869844156580901528e-15148
+3003	0.100000000000000006	8.0606844309353722733298755906e-12145
+3003	1	3.27821191297260784144201433956e-9142
+3003	10	4.05779611158376897690039719814e-6143
+3003	100	3.42800829838332693574242322961e-3179
+3003	1000	6.77752709164469153719156665409e-567
+3003	3301	0.999999933219976432069280011406

tests/tables/inputs/igamma.dat

@@ -16,6 +16,7 @@ a =
 101
 201
 1000
+3003
 
 x =
 1e-6
@@ -26,4 +27,5 @@ x =
 1
 10
 100
-1000
+1000
+3301