JuliaMath / SpecialFunctions.jl / 814

Committed 25 Sep 2020 - 15:06 coverage: 85.701% (+0.05%) from 85.647%

Build # 814

Build Type

Pull #242

travis-ci

Committed by

web-flow

Commit Message

Merge 225521a95 into 1bb0913fc

Pull Request Pull Request #242: Adds Owen's T function

Run Details

32 of 33 new or added lines in 1 file covered. (96.97%)

62 existing lines in 10 files now uncovered.

3644 of 4252 relevant lines covered (85.7%)

2014.15 hits per line

Source File
Press 'n' to go to next uncovered line, 'b' for previous

86.34

/src/gamma_inc.jl

1	using Base.Math: @horner
2	using Base.MPFR: ROUNDING_MODE
3	#useful constants
4	const acc0 = [5.0e-15 , 5.0e-7 , 5.0e-4] #accuracy options
5	const big1 = [25.0 , 14.0, 10.0]
6	const e0 = [0.25e-3 , 0.25e-1 , 0.14]
7	const x0=[31.0 , 17.0 , 9.7]
8	const alog10 = log(10)
9	const rt2pin = 1.0/sqrt(2*pi)
10	const rtpi = sqrt(pi)
11	const exparg = -745.1
12	const stirling_coef = [1.996379051590076518221, -0.17971032528832887213e-2, 0.131292857963846713e-4, -0.2340875228178749e-6, 0.72291210671127e-8, -0.3280997607821e-9, 0.198750709010e-10, -0.15092141830e-11, 0.1375340084e-12, -0.145728923e-13, 0.17532367e-14, -0.2351465e-15, 0.346551e-16, -0.55471e-17, 0.9548e-18, -0.1748e-18, 0.332e-19, -0.58e-20]
13	const auxgam_coef = [-1.013609258009865776949, 0.784903531024782283535e-1, 0.67588668743258315530e-2, -0.12790434869623468120e-2, 0.462939838642739585e-4, 0.43381681744740352e-5, -0.5326872422618006e-6, 0.172233457410539e-7, 0.8300542107118e-9, -0.10553994239968e-9, 0.39415842851e-11, 0.362068537e-13, -0.107440229e-13, 0.5000413e-15, -0.62452e-17, -0.5185e-18, 0.347e-19, -0.9e-21]
14
15	#----------------COEFFICIENTS FOR TEMME EXPANSION------------------
16
17	const d00 = -.333333333333333E+00
18	const d0=[.833333333333333E-01 , -.148148148148148E-01 , .115740740740741E-02 , .352733686067019E-03 , -.178755144032922E-03 , .391926317852244E-04]
19	const d10 = -.185185185185185E-02
20	const d1=[-.347222222222222E-02 , .264550264550265E-02 , -.990226337448560E-03 , .205761316872428E-03]
21	const d20 = .413359788359788E-02
22	const d2=[-.268132716049383E-02 , .771604938271605E-03]
23	const d30 = .649434156378601E-03
24	const d3=[.229472093621399E-03 , -.469189494395256E-03]
25	const d40 = -.861888290916712E-03
26	const d4=[.784039221720067E-03]
27	const d50 = -.336798553366358E-03
28	const d5=[-.697281375836586E-04]
29	const d60 = .531307936463992E-03
30	const d6=[-.592166437353694E-03]
31	const d70 = .344367606892378E-03
32	const d80 = -.652623918595309E-03
33	#Source of logmxp1(x): https://github.com/JuliaStats/StatsFuns.jl/blob/master/src/basicfuns.jl
34	# The kernel of log1pmx
35	# Accuracy within ~2ulps for -0.227 < x < 0.315
36	function _log1pmx_ker(x::Float64)
37	r = x/(x+2.0)	2,568×
38	t = r*r	2,568×
39	w = @horner(t,	2,568×
40	6.66666666666666667e-1, # 2/3
41	4.00000000000000000e-1, # 2/5
42	2.85714285714285714e-1, # 2/7
43	2.22222222222222222e-1, # 2/9
44	1.81818181818181818e-1, # 2/11
45	1.53846153846153846e-1, # 2/13
46	1.33333333333333333e-1, # 2/15
47	1.17647058823529412e-1) # 2/17
48	hxsq = 0.5xx	2,568×
49	r(hxsq+wt)-hxsq	2,568×
50	end
51	"""
52	log1pmx(x::Float64)
53	Return `log(1 + x) - x`
54	Use naive calculation or range reduction outside kernel range. Accurate ~2ulps for all `x`.
55	"""
56	function log1pmx(x::Float64)
57	if !(-0.7 < x < 0.9)	2,568×
58	return log1p(x) - x	!
59	elseif x > 0.315	2,568×
60	u = (x-0.5)/1.5	128×
61	return _log1pmx_ker(u) - 9.45348918918356180e-2 - 0.5*u	128×
62	elseif x > -0.227	2,440×
63	return _log1pmx_ker(x)	2,216×
64	elseif x > -0.4	224×
65	u = (x+0.25)/0.75	120×
66	return _log1pmx_ker(u) - 3.76820724517809274e-2 + 0.25*u	120×
67	elseif x > -0.6	104×
68	u = (x+0.5)*2.0	100×
69	return _log1pmx_ker(u) - 1.93147180559945309e-1 + 0.5*u	100×
70	else
71	u = (x+0.625)/0.375	4×
72	return _log1pmx_ker(u) - 3.55829253011726237e-1 + 0.625*u	4×
73	end
74	end
75	"""
76	logmxp1(x::Float64)
77	Return `log(x) - x + 1` carefully evaluated.
78	"""
79	function logmxp1(x::Float64)
80	if x <= 0.3	1,824×
81	return (log(x) + 1.0) - x	!
82	elseif x <= 0.4	1,824×
83	u = (x-0.375)/0.375	!
84	return _log1pmx_ker(u) - 3.55829253011726237e-1 + 0.625*u	!
85	elseif x <= 0.6	1,824×
86	u = 2.0*(x-0.5)	!
87	return _log1pmx_ker(u) - 1.93147180559945309e-1 + 0.5*u	!
88	else
89	return log1pmx(x - 1.0)	1,824×
90	end
91	end
92
93	"""
94	rgamma1pm1(a)
95
96	Computation of ``1/Gamma(a+1) - 1`` for `-0.5<=a<=1.5` : ``1/\\Gamma (a+1) - 1``
97	Uses the relation `gamma(a+1) = a*gamma(a)`.
98	"""
99	function rgamma1pm1(a::Float64)
100	t=a	8,444×
101	rangereduce = a > 0.5	8,444×
102	t = rangereduce ? a-1 : a #-0.5<= t <= 0.5	8,444×
103	if t == 0.0	8,444×
104	return 0.0	640×
105	elseif t < 0.0	7,804×
106	top = @horner(t , -.422784335098468E+00 , -.771330383816272E+00 , -.244757765222226E+00 , .118378989872749E+00 , .930357293360349E-03 , -.118290993445146E-01 , .223047661158249E-02 , .266505979058923E-03 , -.132674909766242E-03)	3,780×
107	bot = @horner(t , 1.0 , .273076135303957E+00 , .559398236957378E-01)	3,780×
108	w = top/bot	3,780×
109	return rangereduce ? tw/a : a(w+1)	3,780×
110	else
111	top = @horner(t , .577215664901533E+00 , -.409078193005776E+00 , -.230975380857675E+00 , .597275330452234E-01 , .766968181649490E-02 , -.514889771323592E-02 , .589597428611429E-03)	4,024×
112	bot = @horner(t , 1.0 , .427569613095214E+00 , .158451672430138E+00 , .261132021441447E-01 , .423244297896961E-02)	4,024×
113	w = top/bot	4,024×
114	return rangereduce ? (t/a)(w - 1.0) : aw	4,024×
115	end
116	end
117	"""
118	rgammax(a,x)
119
120	Evaluation of ``1/\\Gamma(a) e^{-x} x^{a}``.
121	"""
122	function rgammax(a::Float64,x::Float64)
123	if x == 0.0	3,112×
124	return 0.0	!
125	elseif a >= 20.0	3,112×
126	u =x/a	8×
127	if u == 0.0	8×
128	return 0.0	!
129	end
130	t = (1.0/a)^2	8×
131	t1 = (((0.75t - 1.0)t + 3.5)t - 105.0)/(a1260.0) #Using stirling Series : https://dlmf.nist.gov/5.11.1	8×
132	t1 = t1 + a*logmxp1(u)	8×
133	if t1 >= exparg	8×
134	return rt2pinsqrt(a)exp(t1)	8×
135	end
136	else
137	t = a*log(x) - x	3,104×
138	if t < exparg	3,104×
139	return 0.0	!
140	elseif a >= 1.0	3,104×
141	return exp(t)/gamma(a)	1,760×
142	else
143	return (aexp(t))(1.0 + rgamma1pm1(a))	1,344×
144	end
145	end
146	end
147	"""
148	auxgam(x)
149
150	Compute function `g` in ``1/\\Gamma(x+1) = 1+x(x-1)g(x)``, `-1 <= x <= 1`.
151	"""
152	function auxgam(x::Float64)
153	if x < 0	1,160×
154	return -(1.0+(1+x)(1+x)auxgam(1+x))/(1.0-x)	4×
155	else
156	t = 2*x - 1.0	1,156×
157	return chepolsum(t, auxgam_coef)	1,156×
158	end
159	end
160
161	"""
162	loggamma1p(a)
163
164	Compute ``log(\\Gamma(1+a))`` for `-1 < a <= 1`.
165	"""
166	function loggamma1p(x::Float64)
167	return -log1p(x(x-1.0)auxgam(x))	1,156×
168	end
169
170	"""
171	chepolsum(n,x,a)
172
173	Computes a series of Chebyshev Polynomials given by: `a[1]/2 + a[2]T1(x) + .... + a[n]T{n-1}(X)`.
174	"""
175	function chepolsum(x::Float64, a::Array{Float64,1})
176	n = length(a)	1,168×
177	if n == 1	1,168×
178	return a[1]/2.0	!
179	elseif n == 2	1,168×
180	return a[1]/2.0 + a[2]*x	!
181	else
182	tx = 2*x	1,168×
183	r = a[n]	1,168×
184	h = a[n-1] + r*tx	1,168×
185	for k = n-2:-1:2	2,336×
186	s=r	17,520×
187	r=h	17,520×
188	h=a[k]+r*tx-s	33,872×
189	end
190	return a[1]/2.0 - r + h*x	1,168×
191	end
192	end
193
194	"""
195	stirling(x)
196
197	Compute ``\\ln{\\Gamma(x)} - (x-0.5)*\\ln{x} + x - \\ln{(2\\pi)}/2`` using stirling series for asymptotic part
198	in https://dlmf.nist.gov/5.11.1
199	"""
200	function stirling(x::Float64)
201	if x < floatmin(Float64)*1000.0	484×
202	return floatmax(Float64)/1000.0	!
203	elseif x < 1	484×
204	return loggamma1p(x) - (x+0.5)log(x) + x - log((2pi))/2.0	4×
205	elseif x < 2	480×
206	return loggamma1p(x-1) - (x-0.5)log(x) + x - log((2pi))/2.0	!
207	elseif x < 3	480×
208	return loggamma1p(x-2) - (x-0.5)log(x) + x - log((2pi))/2.0 + log(x-1)	!
209	elseif x < 12	480×
210	z=18.0/(x*x)-1.0	12×
211	return chepolsum(z, stirling_coef)/(12.0*x)	12×
212	else
213	z = 1.0/(x*x)	468×
214	if x < 1000	468×
215	return @horner(z, 0.25721014990011306473e-1, 0.82475966166999631057e-1, -0.25328157302663562668e-2, 0.60992926669463371e-3, -0.33543297638406e-3, 0.250505279903e-3)/((0.30865217988013567769 + z)*x)	468×
216	else
217	return (((-z/1680.0+1.0/1260.0)z-1.0/360.0)z+1.0/12.0)/x	!
218	end
219	end
220	end
221	@doc raw"""
222	gammax(x)
223
224	```math
225	\operatorname{gammax}(x) = \begin{cases}e^{\operatorname{stirling}(x)}\quad\quad\quad \text{if} \quad x>0,\\
226	\frac{\Gamma(x)}{\sqrt{2 \pi}e^{-x + (x-0.5)\operatorname{log}(x)}},\quad \text{if} \quad x\leq 0.
227	\end{cases}
228	```
229	"""
230	function gammax(x::Float64)
231	if x >= 3	708×
232	return exp(stirling(x))	392×
233	elseif x > 0	316×
234	return gamma(x)/(exp(-x+(x-0.5)log(x))sqrt(2*pi))	316×
235	else
236	return floatmax(Float64)/1000.0	!
237	end
238	end
239	"""
240	lambdaeta(eta)
241
242	Compute the value of eta satisfying ``eta^{2}/2 = \\lambda-1-\\ln{\\lambda}``.
243	"""
244	function lambdaeta(eta::Float64)
245	s = etaeta0.5	736×
246	if eta == 0.0	736×
247	la = 1	!
248	elseif eta < -1.0	736×
249	r = exp(-1-s)	!
250	la = @horner(r, 0.0, 1.0, 1.0, 1.5, 8.0/3.0, 125.0/24.0, 15.0/5.0)	!
251	elseif eta < 1.0	736×
252	r = eta	680×
253	la = @horner(r, 1.0, 1.0, 1.0/3.0, 1.0/36.0, -1.0/270.0, 1.0/4320.0, 1.0/17010.0)	680×
254	else
255	r = 11 + s	56×
256	l = log(r)	56×
257	la = r+l	56×
258	r = 1.0/r	56×
259	ak1 = 1.0	56×
260	ak2 = (2 - l)*0.5	56×
261	ak3 = (@horner(l,6,-9,2))/6.0	56×
262	ak4 = -(@horner(l,-12,36,-22,3))/12.0	56×
263	ak5 = (@horner(l,60,-300,350,-125,12))/60.0	56×
264	ak6 = -(@horner(l,-120,900,-1700,1125,-274,20))/120.0	56×
265	la = la + l*@horner(r,0.0, ak1, ak2, ak3, ak4, ak5, ak6)	56×
266	end
267	r = 1	736×
268	if (eta>-3.5 && eta<-0.03) \|\| (eta>0.03 && eta<40)	1,168×
269	r=1	640×
270	q=la	640×
271	while r>1.0e-8	2,308×
272	la = q*(s+log(q))/(q-1.0)	1,028×
273	r = abs(q/la-1)	1,028×
274	q=la	1,028×
275	end
276	end
277	return la	736×
278	end
279	"""
280	Computing the first coefficient for the expansion :
281	```math
282	\\epsilon (\\eta_{0},a) = \\epsilon_{1} (\\eta_{0},a)/a + \\epsilon_{2} (\\eta_{0},a)/a^{2} + \\epsilon_{3} (\\eta_{0},a)/a^{3}
283	```
284	Refer Eqn (3.12) in the paper
285	"""
286	function coeff1(eta::Float64)
287	if abs(eta) < 1.0	696×
288	coeff1 = @horner(eta, -3.333333333438e-1, -2.070740359969e-1, -5.041806657154e-2, -4.923635739372e-3, -4.293658292782e-5) / @horner(eta, 1.000000000000e+0, 7.045554412463e-1, 2.118190062224e-1, 3.048648397436e-2, 1.605037988091e-3)	668×
289	else
290	la = lambdaeta(eta)	28×
291	coeff1 = log(eta/(la-1.0))/eta	28×
292	end
293	return coeff1	696×
294	end
295	"""
296	Computing the second coefficient for the expansion :
297	```math
298	\\epsilon (\\eta_{0},a) = \\epsilon_{1} (\\eta_{0},a)/a + \\epsilon_{2} (\\eta_{0},a)/a^{2} + \\epsilon_{3} (\\eta_{0},a)/a^{3}
299	```
300	Refer Eqn (3.12) in the paper
301	"""
302	function coeff2(eta::Float64)
303
304	if eta < -5.0	696×
305	x=eta*eta	!
306	lnmeta = log(-eta)	!
307	coeff2 = (12.0 - x - 6.0lnmetalnmeta)/(12.0xeta)	!
308	elseif eta < -2.0	696×
309	coeff2 = @horner(eta, -1.72847633523e-2, -1.59372646475e-2, -4.64910887221e-3, -6.06834887760e-4, -6.14830384279e-6) / @horner(eta, 1.00000000000e+0, 7.64050615669e-1, 2.97143406325e-1, 5.79490176079e-2, 5.74558524851e-3)	!
310	elseif eta < 2.0	696×
311	coeff2 = @horner(eta, -1.72839517431e-2, -1.46362417966e-2, -3.57406772616e-3, -3.91032032692e-4, 2.49634036069e-6) / @horner(eta, 1.00000000000e+0, 6.90560400696e-1, 2.49962384741e-1, 4.43843438769e-2, 4.24073217211e-3)	696×
312	elseif eta < 1000.0	!
313	x = 1.0/eta	!
314	coeff2 = @horner(x, 9.99944669480e-1, 1.04649839762e+2, 8.57204033806e+2, 7.31901559577e+2, 4.55174411671e+1) / @horner(x, 1.00000000000e+0, 1.04526456943e+2, 8.23313447808e+2, 3.11993802124e+3, 3.97003311219e+3)	!
315	else
316	coeff2 = -1.0/(12.0*eta)	!
317	end
318	return coeff2	696×
319	end
320	"""
321	Computing the third and last coefficient for the expansion :
322	```math
323	\\epsilon (\\eta_{0},a) = \\epsilon_{1} (\\eta_{0},a)/a + \\epsilon_{2} (\\eta_{0},a)/a^{2} + \\epsilon_{3} (\\eta_{0},a)/a^{3}
324	```
325	Refer Eqn (3.12) in the paper
326	"""
327	function coeff3(eta::Float64)
328	if eta < -8.0	696×
329	x=eta*eta	!
330	y = log(-eta)/eta	!
331	coeff3=(-30.0+etay(6.0xyy-12.0+x))/(12.0etaxx)	!
332	elseif eta < -4.0	696×
333	coeff3 = (@horner(eta, 4.95346498136e-2, 2.99521337141e-2, 6.88296911516e-3, 5.12634846317e-4, -2.01411722031e-5) / @horner(eta, 1.00000000000e+0, 7.59803615283e-1, 2.61547111595e-1, 4.64854522477e-2, 4.03751193496e-3))/(eta*eta)	!
334	elseif eta < -2.0	696×
335	coeff3 = @horner(eta, 4.52313583942e-3, 1.20744920113e-3, -7.89724156582e-5, -5.04476066942e-5, -5.35770949796e-6) / @horner(eta, 1.00000000000e+0, 9.12203410349e-1, 4.05368773071e-1, 9.01638932349e-2, 9.48935714996e-3)	!
336	elseif eta < 2.0	696×
337	coeff3 = @horner(eta, 4.39937562904e-3, 4.87225670639e-4, -1.28470657374e-4, 5.29110969589e-6, 1.57166771750e-7) / @horner(eta, 1.00000000000e+0, 7.94435257415e-1, 3.33094721709e-1, 7.03527806143e-2, 8.06110846078e-3)	696×
338	elseif eta < 10.0	!
339	x=1.0/eta	!
340	coeff3 = (@horner(x, -1.14811912320e-3, -1.12850923276e-1, 1.51623048511e+0, -2.18472031183e-1, 7.30002451555e-2) / @horner(x, 1.00000000000e+0, 1.42482206905e+1, 6.97360396285e+1, 2.18938950816e+2, 2.77067027185e+2))/(eta*eta)	!
341	elseif eta < 100.0	!
342	x=1.0/eta	!
343	coeff3 = (@horner(x, -1.45727889667e-4, -2.90806748131e-1, -1.33085045450e+1, 1.99722374056e+2, -1.14311378756e+1) / @horner(x, 1.00000000000e+0, 1.39612587808e+2, 2.18901116348e+3, 7.11524019009e+3, 4.55746081453e+4))/(eta*eta)	!
344	else
345	eta3 = etaetaeta	!
346	coeff3 = -log(eta)/(12.0*eta3)	!
347	end
348	return coeff3	696×
349	end
350
351
352	"""
353	gamma_inc_cf(a, x, ind)
354
355	Computes ``P(a,x)`` by continued fraction expansion given by :
356	```math
357	R(a,x) * \\frac{1}{1-\\frac{z}{a+1+\\frac{z}{a+2-\\frac{(a+1)z}{a+3+\\frac{2z}{a+4-\\frac{(a+2)z}{a+5+\\frac{3z}{a+6-\\dots}}}}}}}
358	```
359	Used when `1 <= a <= BIG` and `x < x0`.
360
361	External links: [DLMF](https://dlmf.nist.gov/8.9.2)
362
363	See also: [`gamma_inc(a,x,ind)`](@ref SpecialFunctions.gamma_inc)
364	"""
365	function gamma_inc_cf(a::Float64, x::Float64, ind::Integer)
366	acc = acc0[ind + 1]	892×
367	tol = 4.0*acc	892×
368	a2nm1 = 1.0	892×
369	a2n = 1.0	892×
370	b2nm1 = x	892×
371	b2n = x + (1.0 - a)	892×
372	c = 1.0	892×
373	while true	28,564×
374	a2nm1 = xa2n + ca2nm1	28,564×
375	b2nm1 = xb2n + cb2nm1	28,564×
376	c = c + 1.0	28,564×
377	t = c - a	28,564×
378	a2n = a2nm1 + t*a2n	28,564×
379	b2n = b2nm1 + t*b2n	28,564×
380	a2nm1 = a2nm1/b2n	28,564×
381	b2nm1 = b2nm1/b2n	28,564×
382	a2n = a2n/b2n	28,564×
383	b2n = 1.0	28,564×
384	if abs(a2n - a2nm1/b2nm1) < tol*a2n	28,564×
385	break	28,564×
386	end
387	end
388	q = rgammax(a,x)*a2n	892×
389	return (1.0 - q, q)	892×
390	end
391	"""
392	gamma_inc_taylor(a, x, ind)
393
394	Compute ``P(a,x)`` using Taylor Series for P/R given by :
395	```math
396	R(a,x)/a * (1 + \\sum_{n=1}^{\\infty} x^{n}/((a+1)(a+2)...(a+n)))
397	```
398	Used when `1 <= a <= BIG` and `x <= max{a, ln 10}`.
399
400	External links: [DLMF](https://dlmf.nist.gov/8.11.2)
401
402	See also: [`gamma_inc(a,x,ind)`](@ref SpecialFunctions.gamma_inc)
403	"""
404	function gamma_inc_taylor(a::Float64, x::Float64, ind::Integer)
405	acc = acc0[ind + 1]	660×
406	wk = zeros(30)	19,800×
407	flag = false	660×
408	apn = a + 1.0	660×
409	t = x/apn	660×
410	wk[1] = t	660×
411	loop=2	660×
412	for indx = 2:20	660×
413	apn = apn + 1.0	2,896×
414	t = t*(x/apn)	2,896×
415	if t <= 1.0e-3	2,896×
416	loop = indx	660×
417	flag = true	660×
418	break	660×
419	end
420	wk[indx] = t	2,236×
421	end
422	if !flag	660×
UNCOV 423	loop = 20	!
424	end
425	sm = t	660×
426	tol = 0.5*acc #tolerance	660×
427	while true	7,096×
428	apn = apn+1.0	7,096×
429	t = t*(x/apn)	7,096×
430	sm = sm + t	7,096×
431	if t <= tol	7,096×
432	break	7,096×
433	end
434	end
435	for j = loop-1:-1:1	1,320×
436	sm += wk[j]	5,132×
437	end
438	p = (rgammax(a,x)/a)*(1.0 + sm)	660×
439	return (p, 1.0-p)	660×
440	end
441	"""
442	gamma_inc_asym(a, x, ind)
443
444	Compute ``P(a,x)`` using asymptotic expansion given by:
445	```math
446	R(a,x)/a * (1 + \\sum_{n=1}^{N-1}(a_{n}/x^{n} + \\Theta _{n}a_{n}/x^{n}))
447	```
448	where `R(a,x) = rgammax(a,x)`. Used when `1 <= a <= BIG` and `x >= x0`.
449
450	External links: [DLMF](https://dlmf.nist.gov/8.11.2)
451
452	See also: [`gamma_inc(a,x,ind)`](@ref SpecialFunctions.gamma_inc)
453	"""
454	function gamma_inc_asym(a::Float64, x::Float64, ind::Integer)
455	wk = zeros(30)	120×
456	flag = false	4×
457	acc = acc0[ind + 1]	4×
458	amn = a-1.0	4×
459	t=amn/x	4×
460	wk[1]=t	4×
461	loop=2	4×
462	for indx = 2 : 20	4×
463	amn = amn-1.0	12×
464	t=t*(amn/x)	12×
465	if abs(t) <= 1.0e-3	12×
466	loop = indx	4×
467	flag = true	4×
468	break	4×
469	end
470	wk[indx]=t	8×
471	end
472	if !flag	4×
UNCOV 473	loop = 20	!
474	end
475	sm = t	4×
476	while true	20×
477	if abs(t) < acc	20×
478	break	4×
479	end
480	amn=amn-1.0	16×
481	t=t*(amn/x)	16×
482	sm=sm+t	16×
483	end
484	for j = loop-1:-1:1	8×
485	sm += wk[j]	20×
486	end
487	q = (rgammax(a,x)/x)*(1.0 + sm)	4×
488	return (1.0 - q, q)	4×
489	end
490	"""
491	gamma_inc_taylor_x(a,x,ind)
492
493	Computes ``P(a,x)`` based on Taylor expansion of ``P(a,x)/x**a`` given by:
494	```math
495	J = -a * \\sum_{1}^{\\infty} (-x)^{n}/((a+n)n!)
496	```
497	and ``P(a,x)/x**a`` is given by :
498	```math
499	(1 - J)/ \\Gamma(a+1)
500	```
501	resulting from term-by-term integration of `gamma_inc(a,x,ind)`.
502	This is used when `a < 1` and `x < 1.1` - Refer Eqn (9) in the paper.
503
504	See also: [`gamma_inc(a,x,ind)`](@ref SpecialFunctions.gamma_inc)
505	"""
506	function gamma_inc_taylor_x(a::Float64, x::Float64, ind::Integer)
507	acc = acc0[ind + 1]	2,664×
508	l=3.0	2,664×
509	c=x	2,664×
510	sm= x/(a + 3.0)	2,664×
511	tol = 3.0*acc/(a + 1.0)	2,664×
512	while true	18,584×
513	l=l+1.0	18,584×
514	c=-c*(x/l)	18,584×
515	t=c/(a+l)	18,584×
516	sm=sm+t	18,584×
517	if abs(t) <= tol	18,584×
518	break	18,584×
519	end
520	end
521	temp = ax((sm/6.0 - 0.5/(a + 2.0))*x + 1.0/(a + 1.0))	2,664×
522	z = a*log(x)	2,664×
523	h = rgamma1pm1(a)	2,664×
524	g = 1.0 + h	2,664×
525	if (x < 0.25 && z > -.13394) \|\| a < x/2.59	4,976×
526	l = expm1(z)	368×
527	w = 1.0+l	368×
528	q = max((wtemp - l)g - h, 0.0)	368×
529	return (1.0 - q, q)	368×
530	else
531	w = exp(z)	2,296×
532	p = wg(1.0 - temp)	2,296×
533	return (p, 1.0 - p)	2,296×
534	end
535	end
536	"""
537	gamma_inc_minimax(a,x,z)
538
539	Compute ``P(a,x)`` using minimax approximations given by :
540	```math
541	1/2 * erfc(\\sqrt{y}) - e^{-y}/\\sqrt{2\\pia} T(a,\\lambda)
542	``` where
543	```math
544	T(a,\\lambda) = \\sum_{0}^{N} c_{k}(z)a^{-k}
545	```
546
547	This is a higher accuracy approximation of Temme expansion, which deals with the region near `a ≈ x` with `a` large.
548	Refer Appendix F in the paper for the extensive set of coefficients calculated using Brent's multiple precision arithmetic(set at 50 digits) in BRENT, R. P. A FORTRAN multiple-precision arithmetic package, ACM Trans. Math. Softw. 4(1978), 57-70 .
549
550	External links: [DLMF](https://dlmf.nist.gov/8.12.8)
551
552	See also: [`gamma_inc(a,x,ind)`](@ref SpecialFunctions.gamma_inc)
553	"""
554	function gamma_inc_minimax(a::Float64, x::Float64, z::Float64)
555	l = x/a	816×
556	s = 1.0 - l	816×
557	y = -a*logmxp1(l)	816×
558	c = exp(-y)	816×
559	w = 0.5*erfcx(sqrt(y))	816×
560
561	if abs(s) <= 1.0e-3	816×
562	c0 = @horner(z , d00 , d0[1] , d0[2] , d0[3])	8×
563	c1 = @horner(z , d10 , d1[1] , d1[2] , d1[3])	8×
564	c2 = @horner(z , d20 , d2[1] , d2[2])	8×
565	c3 = @horner(z , d30 , d3[1] , d3[2])	8×
566	c4 = @horner(z , d40 , d4[1])	8×
567	c5 = @horner(z , d50 , d5[1])	8×
568	c6 = @horner(z , d60 , d6[1])	8×
569
570	t = @horner(z , c0 , c1 , c2 , c3 , c4 , c5 , c6 , d70 , d80)	8×
571	if l < 1.0	8×
572	p = c(w - rt2pint/sqrt(a))	!
573	return (p , 1.0 - p)	!
574	else
575	q = c(w + rt2pint/sqrt(a))	8×
576	return (1.0 - q , q)	8×
577	end
578	end
579	#---USING THE MINIMAX APPROXIMATIONS---
580	c0 = @horner(z , -.333333333333333E+00 , -.159840143443990E+00 , -.335378520024220E-01 , -.231272501940775E-02)/(@horner(z , 1.0 , .729520430331981E+00 , .238549219145773E+00, .376245718289389E-01 , .239521354917408E-02 , -.939001940478355E-05 , .633763414209504E-06) )	808×
581	c1 = @horner(z , -.185185185184291E-02 , -.491687131726920E-02 , -.587926036018402E-03 , -.398783924370770E-05)/(@horner(z , 1.0 , .780110511677243E+00 , .283344278023803E+00 , .506042559238939E-01 , .386325038602125E-02))	808×
582	c2 = @horner(z , .413359788442192E-02 , .669564126155663E-03)/(@horner(z , 1.0 , .810647620703045E+00 , .339173452092224E+00 , .682034997401259E-01 , .650837693041777E-02 , -.421924263980656E-03))	808×
583	c3 = @horner(z , .649434157619770E-03 , .810586158563431E-03)/(@horner(z , 1.0 , .894800593794972E+00, .406288930253881E+00 , .906610359762969E-01 , .905375887385478E-02 , -.632276587352120E-03))	808×
584	c4 = @horner(z , -.861888301199388E-03 , -.105014537920131E-03)/(@horner(z , 1.0 , .103151890792185E+01 , .591353097931237E+00 , .178295773562970E+00 , .322609381345173E-01))	808×
585	c5 = @horner(z , -.336806989710598E-03, -.435211415445014E-03)/(@horner(z , 1.0 , .108515217314415E+01 , .600380376956324E+00 , .178716720452422E+00))	808×
586	c6 = @horner(z , .531279816209452E-03 , -.182503596367782E-03)/(@horner(z , 1.0 , .770341682526774E+00 , .345608222411837E+00))	808×
587	c7 = @horner(z , .344430064306926E-03, .443219646726422E-03)/(@horner(z , 1.0 , .115029088777769E+01 , .821824741357866E+00))	808×
588	c8 = @horner(z , -.686013280418038E-03 , .878371203603888E-03)	808×
589
590	t = @horner(1.0/a , c0 , c1 , c2 , c3 , c4 , c5 , c6 , c7 , c8)	808×
591	if l < 1.0	808×
592	p = c(w - rt2pint/sqrt(a))	420×
593	return (p , 1.0 - p)	420×
594	else
595	q = c(w + rt2pint/sqrt(a))	388×
596	return (1.0 - q , q)	388×
597	end
598	end
599	"""
600	gamma_inc_temme(a, x, z)
601
602	Compute ``P(a,x)`` using Temme's expansion given by :
603	```math
604	1/2 * erfc(\\sqrt{y}) - e^{-y}/\\sqrt{2\\pia} T(a,\\lambda)
605	``` where
606	```math
607	T(a,\\lambda) = \\sum_{0}^{N} c_{k}(z)a^{-k}
608	```
609	This mainly solves the problem near the region when `a ≈ x` with a large, and is a lower accuracy version of the minimax approximation.
610
611	External links: [DLMF](https://dlmf.nist.gov/8.12.8)
612
613	See also: [`gamma_inc(a,x,ind)`](@ref SpecialFunctions.gamma_inc)
614	"""
615	function gamma_inc_temme(a::Float64, x::Float64, z::Float64)
616	l = x/a	4×
617	y = -a*logmxp1(x/a)	4×
618	c = exp(-y)	4×
619	w = 0.5*erfcx(sqrt(y))	4×
620	c0 = @horner(z , d00 , d0[1] , d0[2] , d0[3] , d0[4] , d0[5] , d0[6])	4×
621	c1 = @horner(z , d10 , d1[1] , d1[2] , d1[3] , d1[4])	4×
622	c2 = @horner(z , d20 , d2[1])	4×
623	t = @horner(1.0/a , c0 , c1 , c2)	4×
624	if l < 1.0	4×
625	p = c(w - rt2pint/sqrt(a))	4×
626	return (p , 1.0 - p)	4×
627	else
628	q = c(w + rt2pint/sqrt(a))	!
629	return (1.0 - q , q)	!
630	end
631	end
632	"""
633	gamma_inc_temme_1(a, x, z, ind)
634
635	Computes ``P(a,x)`` using simplified Temme expansion near ``y=0`` by :
636	```math
637	E(y) - (1 - y)/\\sqrt{2\\pia} T(a,\\lambda)
638	```
639	where
640	```math
641	E(y) = 1/2 - (1 - y/3)*(\\sqrt(y/\\pi))
642	```
643	Used instead of it's previous function when ``\\sigma <= e_{0}/\\sqrt{a}``.
644
645	External links: [DLMF](https://dlmf.nist.gov/8.12.8)
646	"""
647	function gamma_inc_temme_1(a::Float64, x::Float64, z::Float64, ind::Integer)
648	iop = ind + 1	88×
649	l = x/a	88×
650	y = -a * logmxp1(l)	88×
651	if aeps()eps() > 3.28e-3	88×
652	throw(DomainError((a,x,ind,"P(a,x) or Q(a,x) is computationally indeterminant in this case.")))	!
653	end
654	c = exp(-y)	88×
655	w = 0.5*erfcx(sqrt(y))	88×
656	u = 1.0/a	88×
657	if iop == 1	88×
658	c0 = @horner(z , d00 , d0[1] , d0[2] , d0[3])	84×
659	c1 = @horner(z , d10 , d1[1] , d1[2] , d1[3])	84×
660	c2 = @horner(z , d20 , d2[1] , d2[2])	84×
661	c3 = @horner(z , d30 , d3[1] , d3[2])	84×
662	c4 = @horner(z , d40 , d4[1])	84×
663	c5 = @horner(z , d50 , d5[1])	84×
664	c6 = @horner(z , d60 , d6[1])	84×
665
666	t = @horner(u , c0 , c1 , c2 , c3 , c4 , c5 , c6 , d70 , d80)	84×
667
668	elseif iop == 2	4×
669	c0 = @horner(d00 , d0[1] , d0[2])	!
670	c1 = @horner(d10 , d1[1])	!
671	t = @horner(u , c0 , c1 , d20)	!
672
673	else
674	t = @horner(z , d00 , d0[1])	4×
675
676	end
677	if l < 1.0	88×
678	p = c(w - rt2pint/sqrt(a))	8×
679	return (p , 1.0 - p)	8×
680	else
681	q = c(w + rt2pint/sqrt(a))	80×
682	return (1.0 - q , q)	80×
683	end
684	end
685	"""
686	gamma_inc_fsum(a,x)
687
688	Compute using Finite Sums for ``Q(a,x)`` when `a >= 1 && 2a` is integer.
689	Used when `a <= x <= x0` and `a = n/2`.
690	Refer Eqn (14) in the paper.
691
692	See also: [`gamma_inc(a,x,ind)`](@ref SpecialFunctions.gamma_inc)
693	"""
694	function gamma_inc_fsum(a::Float64, x::Float64)
695	if isinteger(a)	140×
696	sm = exp(-x)	132×
697	t = sm	132×
698	N = 1	132×
699	c=0.0	132×
700	i = trunc(Int,a )	132×
701	else
702	rtx = sqrt(x)	8×
703	sm = erfc(rtx)	8×
704	t = exp(-x)/(rtpi*rtx)	8×
705	N=0	8×
706	c=-0.5	8×
707	i = trunc(Int,a )	8×
708	end
709	for lp = N : i-1	160×
710	if i == 0	252×
711	break	!
712	end
713	c += 1.0	252×
714	t = (x*t)/c	252×
715	sm += t	252×
716	end
717	q = sm	140×
718	return (1.0 - q, q)	140×
719
720	end
721	"""
722	gamma_inc_inv_psmall(a,p)
723
724	Compute `x0` - initial approximation when `p` is small.
725	Here we invert the series in Eqn (2.20) in the paper and write the inversion problem as:
726	```math
727	x = r\\left[1 + a\\sum_{k=1}^{\\infty}\\frac{(-x)^{n}}{(a+n)n!}\\right]^{-1/a},
728	```
729	where ``r = (p\\Gamma(1+a))^{1/a}``
730	Inverting this relation we obtain ``x = r + \\sum_{k=2}^{\\infty}c_{k}r^{k}``.
731	"""
732	function gamma_inc_inv_psmall(a::Float64, p::Float64)
733	logr = (1.0/a)*(log(p) + logabsgamma(a + 1.0)[1])	936×
734	r = exp(logr)	936×
735	ap1 = a + 1.0	936×
736	ap1² = (ap1)*ap1	936×
737	ap1³ = (ap1)*ap1²	936×
738	ap1⁴ = ap1²*ap1²	936×
739	ap2 = a+2.0	936×
740	ap2² = ap2*ap2	936×
741	ck1 = 1.0	936×
742	ck2 = 1.0/(1.0+a)	936×
743	ck3 = 0.5(3a+5)/(ap1²*(a+2))	936×
744	ck4 = (1.0/3.0)(@horner(a,31,33,8))/(ap1³ap2*(a+3))	936×
745	ck5 = (1.0/24.0)(@horner(a, 2888, 5661, 3971, 1179, 125))/(ap1⁴ap2²(a+3)(a+4))	936×
746	x0 = @horner(r, 0.0, ck1, ck2, ck3, ck4, ck5)	936×
747	return x0	936×
748	end
749	"""
750	gamma_inc_inv_qsmall(a,q)
751
752	Compute `x0` - initial approximation when `q` is small from ``e^{-x_{0}} x_{0}^{a} = q \\Gamma(a)``.
753	Asymptotic expansions Eqn (2.29) in the paper is used here and higher approximations are obtained using
754	```math
755	x \\sim x_{0} - L + b \\sum_{k=1}^{\\infty} d_{k}/x_{0}^{k}
756	```
757	where ``b = 1-a``, ``L = \\ln{x_0}``.
758	"""
759	function gamma_inc_inv_qsmall(a::Float64, q::Float64)
760	b=1.0-a	12×
761	eta=sqrt(-2/alog(qgammax(a)sqrt(2pi)/sqrt(a)))	12×
762	x0 = a*lambdaeta(eta)	12×
763	l = log(x0)	12×
764
765	if a > 0.12 \|\| x0 > 5	12×
766	r = 1.0/x0	12×
767	ck1 = l - 1.0	12×
768	ck2 = (@horner(l, @horner(b, 2, 3), @horner(b, -2, -2), 1))/2.0	12×
769	ck3 = (@horner(l, @horner(b, -12, -24, -11), @horner(b, 12, 24, 6), @horner(b, -6, -9), 2))/6.0	12×
770	ck4 = (@horner(l, @horner(b, 72, 162, 120, 25), @horner(b, -72, -168, -114, -12), @horner(b, 36, 84, 36), @horner(b, -12, -22), 3))/12.0	12×
771	x0 = x0 - l + b * r * @horner(r, ck1, ck2, ck3, ck4)	12×
772	else
773	r = 1.0/x0	!
UNCOV 774	l² = l*l	!
775	ck1 = l - 1.0	!
776	x0 = x0 - l + b * r * ck1	!
777	end
778	return x0	12×
779	end
780	"""
781	gamma_inc_inv_asmall(a,p,q,pcase)
782
783	Compute `x0` - initial approximation when `a` is small.
784	Here the solution `x` of ``P(a,x)=p`` satisfies ``x_{l} < x < x_{u}``
785	where ``x_{l} = (p\\Gamma(a+1))^{1/a}`` and ``x_{u} = -\\log{(1 - p\\Gamma(a+1))}``, and is used as starting value for Newton iteration.
786	"""
787	function gamma_inc_inv_asmall(a::Float64, p::Float64, q::Float64, pcase::Bool)
788	logp = (pcase) ? log(p) : log1p(-q)	640×
789	return exp((1.0/a)*(logp +loggamma1p(a)))	324×
790	end
791	"""
792	gamma_inc_inv_alarge(a,porq,s)
793
794	Compute `x0` - initial approximation when `a` is large.
795	The inversion problem is rewritten as :
796	```math
797	0.5 \\operatorname{erfc}(\\eta \\sqrt{a/2}) + R_{a}(\\eta) = q
798	```
799	For large values of `a` we can write: ``\\eta(q,a) = \\eta_{0}(q,a) + \\epsilon(\\eta_{0},a)``
800	and it is possible to expand:
801	```math
802	\\epsilon(\\eta_{0},a) = \\epsilon_{1}(\\eta_{0},a)/a + \\epsilon_{2}(\\eta_{0},a)/a^{2} + \\epsilon_{3}(\\eta_{0},a)/a^{3} + ...
803	```
804	which is calculated by coeff1, coeff2 and coeff3 functions below.
805	This returns a tuple `(x0,fp)`, where `fp` is computed since it's an approximation for the coefficient after inverting the original power series.
806	"""
807	function gamma_inc_inv_alarge(a::Float64, porq::Float64, s::Integer)
808	r = erfcinv(2*porq)	696×
809	eta = sr/sqrt(a0.5)	696×
810	eta += (coeff1(eta) + (coeff2(eta) + coeff3(eta)/a)/a)/a	696×
811	x0 = a*lambdaeta(eta)	696×
812	fp = -sqrt(a/(2pi))exp(-0.5aeta*eta)/gammax(a)	696×
813	return (x0,fp)	696×
814	end
815	# Reference : 'Computation of the incomplete gamma function ratios and their inverse' by Armido R DiDonato , Alfred H Morris.
816	# Published in Journal: ACM Transactions on Mathematical Software (TOMS)
817	# Volume 12 Issue 4, Dec. 1986 Pages 377-393
818	# doi>10.1145/22721.23109
819
820	"""
821	gamma_inc(a,x,IND)
822
823	Returns a tuple ``(p, q)`` where ``p + q = 1``, and
824	``p=P(a,x)`` is the Incomplete gamma function ratio given by:
825	```math
826	P(a,x)=\\frac{1}{\\Gamma (a)} \\int_{0}^{x} e^{-t}t^{a-1} dt.
827	```
828	and ``q=Q(a,x)`` is the Incomplete gamma function ratio given by:
829	```math
830	Q(x,a)=\\frac{1}{\\Gamma (a)} \\int_{x}^{\\infty} e^{-t}t^{a-1} dt.
831	```
832
833	`IND ∈ [0,1,2]` sets accuracy: `IND=0` means 14 significant digits accuracy, `IND=1` means 6 significant digit, and `IND=2` means only 3 digit accuracy.
834
835	External links: [DLMF](https://dlmf.nist.gov/8.2.4), [Wikipedia](https://en.wikipedia.org/wiki/Incomplete_gamma_function)
836
837	See also [`gamma(z)`](@ref SpecialFunctions.gamma), [`gamma_inc_inv(a,p,q)`](@ref SpecialFunctions.gamma_inc_inv)
838	"""
839	function gamma_inc(a::Float64,x::Float64,ind::Integer)
840	iop = ind + 1	5,524×
841	acc = acc0[iop]	5,524×
842	if a<0.0 \|\| x<0.0	11,044×
843	throw(DomainError((a,x,ind,"`a` and `x` must be greater than 0 ---- Domain : (0,inf)")))	4×
844	elseif a==0.0 && x==0.0	5,520×
845	throw(DomainError((a,x,ind,"`a` and `x` must be greater than 0 ---- Domain : (0,inf)")))	4×
846	elseif a*x==0.0	5,516×
847	if x<=a	188×
848	return (0.0,1.0)	188×
849	else
850	return (1.0,0.0)	!
851	end
852	end
853
854	if a >= 1.0	5,328×
855	if a >= big1[iop]	1,932×
856	l = x/a	908×
857	if l == 0.0	908×
858	return (0.0,1.0)	!
859	end
860	s = 1.0 - l	908×
861	z = -logmxp1(l)	908×
862	if z >= 700.0/a	908×
863	return (0.0,1.0)	!
864	end
865	y = a*z	908×
866	rta = sqrt(a)	908×
867	if abs(s) <= e0[iop]/rta	908×
868	z = sqrt(z + z)	84×
869	if l < 1.0	84×
870	z=-z	4×
871	end
872	return gamma_inc_temme_1(a, x, z, ind)	84×
873	end
874
875	if abs(s) <= 0.4	824×
876	if abs(s) <= 2.0eps() && aeps()*eps() > 3.28e-3	824×
877	throw(DomainError((a,x,ind,"P(a,x) or Q(a,x) is computationally indeterminant in this case.")))	!
878	end
879	c = exp(-y)	824×
880	w = 0.5*erfcx(sqrt(y))	824×
881	u = 1.0/a	824×
882	z = sqrt(z + z)	824×
883	if l < 1.0	824×
884	z=-z	428×
885	end
886	if iop == 1	824×
887	return gamma_inc_minimax(a,x,z)	816×
888	elseif iop == 2	8×
889	return gamma_inc_temme(a,x,z)	4×
890	else
891	t = @horner(z , d00 , d0[1] , d0[2] , d0[3])	4×
892	return gamma_inc_temme_1(a, x, z, ind)	4×
893	end
894	end
895	elseif a > x \|\| x >= x0[iop] \|\| !isinteger(2*a)	1,492×
896	r = rgammax(a,x)	884×
897	if r == 0.0	1,105×
898	if x <= a	!
899	return (0.0,1.0)	!
900	else
901	return (1.0,0.0)	!
902	end
903	end
904	if x <= max(a,alog10)	884×
905	return gamma_inc_taylor(a, x, ind)	660×
906	elseif x < x0[iop]	224×
907	return gamma_inc_cf(a, x, ind)	220×
908	else
909	return gamma_inc_asym(a, x, ind)	4×
910	end
911	else
912	return gamma_inc_fsum(a,x)	140×
913
914	end
915	elseif a == 0.5	3,396×
916	if x >= 0.25	60×
917	q = erfc(sqrt(x))	56×
918	return ( 1.0 - q , q )	56×
919	end
920	p = erf(sqrt(x))	4×
921	return ( p , 1.0 - p )	4×
922	elseif x < 1.1	3,336×
923	return gamma_inc_taylor_x(a, x, ind)	2,664×
924	end
925	r = rgammax(a,x)	672×
926	if r == 0.0	840×
927	return (1.0, 0.0)	!
928	else
929	return gamma_inc_cf(a, x, ind)	672×
930	end
931
932
933
934	end
935
936	function gamma_inc(a::BigFloat,x::BigFloat,ind::Integer) #BigFloat version from GNU MPFR wrapped via ccall
937	z = BigFloat()	4×
938	ccall((:mpfr_gamma_inc, :libmpfr), Int32 , (Ref{BigFloat} , Ref{BigFloat} , Ref{BigFloat} , Int32) , z , a , x , ROUNDING_MODE[])	4×
939	q = z/gamma(a)	4×
940	return (1.0 - q, q)	4×
941	end
942	gamma_inc(a::Float32,x::Float32,ind::Integer) = ( Float32(gamma_inc(Float64(a),Float64(x),ind)[1]) , Float32(gamma_inc(Float64(a),Float64(x),ind)[2]) )	4×
943	gamma_inc(a::Float16,x::Float16,ind::Integer) = ( Float16(gamma_inc(Float64(a),Float64(x),ind)[1]) , Float16(gamma_inc(Float64(a),Float64(x),ind)[2]) )	4×
944	gamma_inc(a::Real,x::Real,ind::Integer) = (gamma_inc(float(a),float(x),ind))	24×
945	gamma_inc(a::Integer,x::Integer,ind::Integer) = gamma_inc(Float64(a),Float64(x),ind)	44×
946	gamma_inc(a::AbstractFloat,x::AbstractFloat,ind::Integer) = throw(MethodError(gamma_inc,(a,x,ind,"")))	!
947
948	#EFFICIENT AND ACCURATE ALGORITHMS FOR THECOMPUTATION AND INVERSION OF THE INCOMPLETEGAMMA FUNCTION RATIOS by Amparo Gil, Javier Segura, Nico M. Temme
949	#SIAM Journal on Scientific Computing 34(6) (2012), A2965-A2981
950	# arXiv:1306.1754
951
952	"""
953	gamma_inc_inv(a,p,q)
954
955	Inverts the `gamma_inc(a,x)` function, by computing `x` given `a`,`p`,`q` in ``P(a,x)=p`` and ``Q(a,x)=q``.
956
957	External links: [DLMF](https://dlmf.nist.gov/8.2.4), [Wikipedia](https://en.wikipedia.org/wiki/Incomplete_gamma_function)
958
959	See also: [`gamma_inc(a,x,ind)`](@ref SpecialFunctions.gamma_inc).
960	"""
961	function gamma_inc_inv(a::Float64, p::Float64, q::Float64)
962	if p < 0.5	1,984×
963	pcase = true	980×
964	porq = p	980×
965	s = -1	980×
966	else
967	pcase = false	1,004×
968	porq = q	1,004×
969	s = 1	1,004×
970	end
971	haseta = false	1,984×
972
973	logr = (1.0/a)*( log(p) + logabsgamma(a + 1.0)[1] )	1,984×
974	if logr < log(0.2*(1+a)) #small value of p	1,984×
975	x0 = gamma_inc_inv_psmall(a,p)	936×
976	elseif ((q < min(0.02, exp(-1.5*a)/gamma(a))) && (a < 10)) #small q	1,048×
977	x0 = gamma_inc_inv_qsmall(a, q)	12×
978	elseif abs(porq - 0.5) < 1.0e-05	1,036×
979	x0 = a - 1.0/3.0 + (8.0/405.0 + 184.0/25515.0/a) / a	16×
980	elseif abs(a-1.0) < 1.0e-4	1,020×
981	x0 = pcase ? -log1p(-p) : -log(q)	!
982	elseif a < 1.0 # small value of a	1,020×
983	x0 = gamma_inc_inv_asmall(a, p, q, pcase)	324×
984	else #large a
985	haseta = true	696×
986	(x0,fp) = gamma_inc_inv_alarge(a,porq,s)	696×
987	end
988
989	t = 1	1,984×
990	x = x0	1,984×
991	n = 1	1,984×
992	logabsgam = logabsgamma(a)[1]	1,984×
993	#Newton-like higher order iteration with upper limit as 15.
994	while t > 1.0e-15 && n < 15	8,908×
995	x = x0	4,940×
996	x² = x*x	4,940×
997	if !haseta	4,940×
998	dlnr = (1.0-a)*log(x) + x + logabsgam	3,548×
999	if dlnr > log(floatmax(Float64)/1000.0)	3,548×
1000	n=20	!
1001	else
1002	r = exp(dlnr)	7,096×
1003	end
1004	else
1005	r = -(1/fp)*x	1,392×
1006	end
1007
1008	(px,qx) = gamma_inc(a,x,0)	4,940×
1009	ck1 = (pcase) ? -r(px-p) : r(qx-q)	8,040×
1010	ck2 = (x-a+1.0)/(2.0*x)	4,940×
1011	ck3 = (@horner(x, @horner(a, 1, -3, 2), @horner(a, 4, -4), 2))/(6*x²)	4,940×
1012	r = ck1	4,940×
1013	if a > 0.1	4,940×
1014	x0 = @horner(r, x, 1.0, ck2, ck3)	3,972×
1015	else
1016	if a > 0.05	968×
1017	x0 = @horner(r, x, 1.0, ck2)	!
1018	else
1019	x0 = x + r	968×
1020	end
1021	end
1022
1023	t=abs(x/x0 - 1.0)	4,940×
1024	n+=1	4,940×
1025	x=x0	4,940×
1026
1027	end
1028	return x	1,984×
1029	end
1030
1031	gamma_inc_inv(a::Float32, p::Float32, q::Float32) = Float32( gamma_inc_inv(Float64(a), Float64(p), Float64(q)) )	!
1032	gamma_inc_inv(a::Float16, p::Float16, q::Float16) = Float16( gamma_inc_inv(Float64(a), Float64(p), Float64(q)) )	!
1033	gamma_inc_inv(a::Real, p::Real, q::Real) = ( gamma_inc_inv(float(a), float(p), float(q)) )	!
1034	gamma_inc_inv(a::Integer, p::Integer, q::Integer) = ( gamma_inc_inv(Float64(a), Float64(p), Float64(q)) )	!
1035	gamma_inc_inv(a::AbstractFloat, p::AbstractFloat, q::AbstractFloat) = throw(MethodError(gamma_inc_inv,(a,p,q,"")))	!

JuliaMath / SpecialFunctions.jl / 814

Source File Press 'n' to go to next uncovered line, 'b' for previous

Source File
Press 'n' to go to next uncovered line, 'b' for previous