Line data Source code
1 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
3 : !!!! !!!!
4 : !!!! ParaMonte: Parallel Monte Carlo and Machine Learning Library. !!!!
5 : !!!! !!!!
6 : !!!! Copyright (C) 2012-present, The Computational Data Science Lab !!!!
7 : !!!! !!!!
8 : !!!! This file is part of the ParaMonte library. !!!!
9 : !!!! !!!!
10 : !!!! LICENSE !!!!
11 : !!!! !!!!
12 : !!!! https://github.com/cdslaborg/paramonte/blob/main/LICENSE.md !!!!
13 : !!!! !!!!
14 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
15 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
16 :
17 : !> \brief
18 : !> This include file contains the implementation of procedures in [pm_mathErf](@ref pm_mathErf).
19 : !>
20 : !> \finmain
21 : !>
22 : !> \author
23 : !> \AmirShahmoradi, Oct 16, 2009, 11:14 AM, Michigan
24 :
25 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
26 :
27 : #if getErfInv_ENABLED
28 8013 : if (present(abserr)) then
29 0 : call setErfInv(erfinv, x, abserr)
30 : else
31 8013 : call setErfInv(erfinv, x, 100 * epsilon(0._RKC))
32 : end if
33 : #elif setErfInv_ENABLED
34 : real(RKC) :: temp, absx
35 : real(RKC), parameter :: TWO_OVER_SQRTPI = 2._RKC / sqrt(acos(-1._RKC))
36 : ! Chebyshev declarations.
37 : integer(IK) , parameter :: MINDEL = 0, NDELTA = 37, NLAMDA = 26, NMU = 25, NXI = 38
38 : integer(IK) , parameter :: MAXDEL = MINDEL + NDELTA ! 37
39 : integer(IK) , parameter :: MINLAM = MAXDEL + 1 ! 38
40 : integer(IK) , parameter :: MAXLAM = MINLAM + NLAMDA ! 64
41 : integer(IK) , parameter :: MINMU = MAXLAM + 1 ! 65
42 : integer(IK) , parameter :: MAXMU = MINMU + NMU ! 90
43 : integer(IK) , parameter :: MINXI = MAXMU + 1 ! 91
44 : integer(IK) , parameter :: MAXXI = MINXI + NXI ! 129
45 : ! The vector `chebyshScaler` is used to scale the argument of the Chebyshev polynomial expansion in the range 1.D-300 < 1 - X < 0.2.
46 : real(RKC), parameter :: chebyshScaler(6) = [ -1.548813042373261659512742_RKC &
47 : , +2.565490123147816151928163_RKC &
48 : , -.5594576313298323225436913_RKC &
49 : , +2.287915716263357638965891_RKC &
50 : , -9.199992358830151031278420_RKC &
51 : , +2.794990820124599493768426_RKC ]
52 : ! The coefficients of polynomial expansions, stored in the order DELTA(0..37), LAMDA(0..26), MU(0..25), XI(0..38), where
53 : ! DELTA are coefficients of the Chebyshev polynomial expansion of erfinv(X) for 0.9975 < X <= 1-5.0D-16.
54 : ! LAMDA are coefficients of the Chebyshev polynomial expansion of erfinv(X) for 0.8 < X <= 0.9975.
55 : ! MU are coefficients of the Chebyshev polynomial expansion of erfinv(X) for 1 - 5.0D-16 < X <= 1 - 1.D-300.
56 : ! XI are coefficients of the Chebyshev polynomial expansion of erfinv(X) in the range 0.0 <= X <= 0.8.
57 : real(RKC), parameter :: coef(0:MAXXI) = [ +0.9566797090204925274526373_RKC &
58 : , -0.0231070043090649036999908_RKC &
59 : , -0.0043742360975084077333218_RKC &
60 : , -0.0005765034226511854809364_RKC &
61 : , -0.0000109610223070923931242_RKC &
62 : , +0.0000251085470246442787982_RKC &
63 : , +0.0000105623360679477511955_RKC &
64 : , +0.0000027544123300306391503_RKC &
65 : , +0.0000004324844983283380689_RKC &
66 : , -0.0000000205303366552086916_RKC &
67 : , -0.0000000438915366654316784_RKC &
68 : , -0.0000000176840095080881795_RKC &
69 : , -0.0000000039912890280463420_RKC &
70 : , -0.0000000001869324124559212_RKC &
71 : , +0.0000000002729227396746077_RKC &
72 : , +0.0000000001328172131565497_RKC &
73 : , +0.0000000000318342484482286_RKC &
74 : , +0.0000000000016700607751926_RKC &
75 : , -0.0000000000020364649611537_RKC &
76 : , -0.0000000000009648468127965_RKC &
77 : , -0.0000000000002195672778128_RKC &
78 : , -0.0000000000000095689813014_RKC &
79 : , +0.0000000000000137032572230_RKC &
80 : , +0.0000000000000062538505417_RKC &
81 : , +0.0000000000000014584615266_RKC &
82 : , +0.0000000000000001078123993_RKC &
83 : , -0.0000000000000000709229988_RKC &
84 : , -0.0000000000000000391411775_RKC &
85 : , -0.0000000000000000111659209_RKC &
86 : , -0.0000000000000000015770366_RKC &
87 : , +0.0000000000000000002853149_RKC &
88 : , +0.0000000000000000002716662_RKC &
89 : , +0.0000000000000000000957770_RKC &
90 : , +0.0000000000000000000176835_RKC &
91 : , -0.0000000000000000000009828_RKC &
92 : , -0.0000000000000000000020464_RKC &
93 : , -0.0000000000000000000008020_RKC &
94 : , -0.0000000000000000000001650_RKC &
95 : , +0.9121588034175537733059200_RKC &
96 : , -0.0162662818676636958546661_RKC &
97 : , +0.0004335564729494453650589_RKC &
98 : , +0.0002144385700744592065205_RKC &
99 : , +0.0000026257510757648130176_RKC &
100 : , -0.0000030210910501037969912_RKC &
101 : , -0.0000000124060618367572157_RKC &
102 : , +0.0000000624066092999917380_RKC &
103 : , -0.0000000005401247900957858_RKC &
104 : , -0.0000000014232078975315910_RKC &
105 : , +0.0000000000343840281955305_RKC &
106 : , +0.0000000000335848703900138_RKC &
107 : , -0.0000000000014584288516512_RKC &
108 : , -0.0000000000008102174258833_RKC &
109 : , +0.0000000000000525324085874_RKC &
110 : , +0.0000000000000197115408612_RKC &
111 : , -0.0000000000000017494333828_RKC &
112 : , -0.0000000000000004800596619_RKC &
113 : , +0.0000000000000000557302987_RKC &
114 : , +0.0000000000000000116326054_RKC &
115 : , -0.0000000000000000017262489_RKC &
116 : , -0.0000000000000000002784973_RKC &
117 : , +0.0000000000000000000524481_RKC &
118 : , +0.0000000000000000000065270_RKC &
119 : , -0.0000000000000000000015707_RKC &
120 : , -0.0000000000000000000001475_RKC &
121 : , +0.0000000000000000000000450_RKC &
122 : , +0.9885750640661893136460358_RKC &
123 : , +0.0108577051845994776160281_RKC &
124 : , -0.0017511651027627952494825_RKC &
125 : , +0.0000211969932065633437984_RKC &
126 : , +0.0000156648714042435087911_RKC &
127 : , -0.0000005190416869103124261_RKC &
128 : , -0.0000000371357897426717780_RKC &
129 : , +0.0000000012174308662357429_RKC &
130 : , -0.0000000001768115526613442_RKC &
131 : , -0.0000000000119372182556161_RKC &
132 : , +0.0000000000003802505358299_RKC &
133 : , -0.0000000000000660188322362_RKC &
134 : , -0.0000000000000087917055170_RKC &
135 : , -0.0000000000000003506869329_RKC &
136 : , -0.0000000000000000697221497_RKC &
137 : , -0.0000000000000000109567941_RKC &
138 : , -0.0000000000000000011536390_RKC &
139 : , -0.0000000000000000001326235_RKC &
140 : , -0.0000000000000000000263938_RKC &
141 : , +0.0000000000000000000005341_RKC &
142 : , -0.0000000000000000000022610_RKC &
143 : , +0.0000000000000000000009552_RKC &
144 : , -0.0000000000000000000005250_RKC &
145 : , +0.0000000000000000000002487_RKC &
146 : , -0.0000000000000000000001134_RKC &
147 : , +0.0000000000000000000000420_RKC &
148 : , +0.9928853766189408231495800_RKC &
149 : , +0.1204675161431044864647846_RKC &
150 : , +0.0160781993420999447257039_RKC &
151 : , +0.0026867044371623158279591_RKC &
152 : , +0.0004996347302357262947170_RKC &
153 : , +0.0000988982185991204409911_RKC &
154 : , +0.0000203918127639944337340_RKC &
155 : , +0.0000043272716177354218758_RKC &
156 : , +0.0000009380814128593406758_RKC &
157 : , +0.0000002067347208683427411_RKC &
158 : , +0.0000000461596991054300078_RKC &
159 : , +0.0000000104166797027146217_RKC &
160 : , +0.0000000023715009995921222_RKC &
161 : , +0.0000000005439284068471390_RKC &
162 : , +0.0000000001255489864097987_RKC &
163 : , +0.0000000000291381803663201_RKC &
164 : , +0.0000000000067949421808797_RKC &
165 : , +0.0000000000015912343331469_RKC &
166 : , +0.0000000000003740250585245_RKC &
167 : , +0.0000000000000882087762421_RKC &
168 : , +0.0000000000000208650897725_RKC &
169 : , +0.0000000000000049488041039_RKC &
170 : , +0.0000000000000011766394740_RKC &
171 : , +0.0000000000000002803855725_RKC &
172 : , +0.0000000000000000669506638_RKC &
173 : , +0.0000000000000000160165495_RKC &
174 : , +0.0000000000000000038382583_RKC &
175 : , +0.0000000000000000009212851_RKC &
176 : , +0.0000000000000000002214615_RKC &
177 : , +0.0000000000000000000533091_RKC &
178 : , +0.0000000000000000000128488_RKC &
179 : , +0.0000000000000000000031006_RKC &
180 : , +0.0000000000000000000007491_RKC &
181 : , +0.0000000000000000000001812_RKC &
182 : , +0.0000000000000000000000439_RKC &
183 : , +0.0000000000000000000000106_RKC &
184 : , +0.0000000000000000000000026_RKC &
185 : , +0.0000000000000000000000006_RKC &
186 : , +0.0000000000000000000000002_RKC &
187 : ]
188 : real(RKC), parameter :: HALF_EPS = 0.5_RKC * epsilon(0._RKC)
189 : real(RKC) :: carg, cargTwice, w1, w2, w3
190 : ! COEF_LIM_FULL is an array containing MINXI, MAXXI, MINLAM, MAXLAM, MINDEL, MAXDEL, MINMU, MAXMU in locations -1..6 where,
191 : ! MAX... is the position in `coef` of the last coefficient of a Chebyshev polynomial expansion.
192 : ! MIN... is the position in `coef` of the first coefficient of a Chebyshev polynomial expansion.
193 : integer(IK) , parameter :: COEF_LIM_FULL(-1:6) = [MINXI, MAXXI, MINLAM, MAXLAM, MINDEL, MAXDEL, MINMU, MAXMU]
194 : ! Decide which approximation to use, and calculate the argument of the Chebyshev polynomial expansion.
195 : ! Compute the minimum index of a coefficient in the Chebyshev polynomial expansion to be used.
196 : ! That is, include only the coefficients whose precision is relevant to the requested output precision.
197 : ! \bug ifort 2021.8 : Cannot digest implied do loop below. Gfortran is fine.
198 : ! \remedy For now, the loop was unrolled.
199 : !integer(IK) , parameter :: COEF_LIM(-1:6) = [(COEF_LIM_FULL(i - 1), COEF_LIM_FULL(i) - count(abs(coef(COEF_LIM_FULL(i - 1) : COEF_LIM_FULL(i))) < HALF_EPS), i = 0, 6, 2)]
200 : integer(IK) , parameter :: COEF_LIM(-1:6) = [ COEF_LIM_FULL(-1), COEF_LIM_FULL(0) - count(abs(coef(COEF_LIM_FULL(-1) : COEF_LIM_FULL(0))) < HALF_EPS) &
201 : , COEF_LIM_FULL(+1), COEF_LIM_FULL(2) - count(abs(coef(COEF_LIM_FULL(+1) : COEF_LIM_FULL(2))) < HALF_EPS) &
202 : , COEF_LIM_FULL(+3), COEF_LIM_FULL(4) - count(abs(coef(COEF_LIM_FULL(+3) : COEF_LIM_FULL(4))) < HALF_EPS) &
203 : , COEF_LIM_FULL(+5), COEF_LIM_FULL(6) - count(abs(coef(COEF_LIM_FULL(+5) : COEF_LIM_FULL(6))) < HALF_EPS) &
204 : ]
205 : integer(IK) :: i, j, imin
206 24056 : if(-1._RKC < x .and. x < 1._RKC) then
207 24056 : if(1.e-7_RKC < abserr) then
208 12042 : temp = -log((1.0 - x)*(1.0 + x));
209 12042 : if (temp < 5._RKC) then
210 11992 : temp = temp - 2.5_RKC
211 : erfinv = +2.81022636e-08_RKC
212 11992 : erfinv = +3.43273939e-07_RKC + erfinv * temp
213 11992 : erfinv = -3.5233877e-06_RKC + erfinv * temp
214 11992 : erfinv = -4.39150654e-06_RKC + erfinv * temp
215 11992 : erfinv = +0.00021858087_RKC + erfinv * temp
216 11992 : erfinv = -0.00125372503_RKC + erfinv * temp
217 11992 : erfinv = -0.00417768164_RKC + erfinv * temp
218 11992 : erfinv = +0.246640727_RKC + erfinv * temp
219 11992 : erfinv = +1.50140941_RKC + erfinv * temp
220 : else
221 50 : temp = sqrt(temp) - 3._RKC
222 : erfinv = -0.000200214257_RKC
223 50 : erfinv = +0.000100950558_RKC + erfinv * temp
224 50 : erfinv = +0.00134934322_RKC + erfinv * temp
225 50 : erfinv = -0.00367342844_RKC + erfinv * temp
226 50 : erfinv = +0.00573950773_RKC + erfinv * temp
227 50 : erfinv = -0.0076224613_RKC + erfinv * temp
228 50 : erfinv = +0.00943887047_RKC + erfinv * temp
229 50 : erfinv = +1.00167406_RKC + erfinv * temp
230 50 : erfinv = +2.83297682_RKC + erfinv * temp
231 : end if
232 12042 : erfinv = erfinv * x
233 12014 : elseif (2.e-24_RKC < abserr) then
234 : ! This algorithm is based on the approximate formulae of Mathematics of Computation 22, (1968) PP144-158, which he claims to be accurate up to 23 decimal points.
235 : ! Amir Shahmoradi, Nov 10, 2009, 8:53 PM, Michigan
236 : ! Amir Shahmoradi, Wednesday 5:30 AM, XXX XX, 2012, Institute for Fusion Studies, The University of Texas Austin
237 : ! HAPPY BIRTHDAY AMIR! (Self-congrats syndrome caused by lack of sleep)
238 6001 : if (x == 0._RKC) then
239 1 : erfinv = 0._RKC
240 1 : return
241 : end if
242 6000 : absx = abs(x) ! The argument of the Chebyshev polynomial expansion.
243 6000 : if (absx <= 0.8_RKC) then
244 4800 : carg = 3.125_RKC * absx * absx - 1._RKC
245 : j = -1
246 : else
247 1200 : if (absx <= 0.9975_RKC) then
248 : j = 1
249 : else
250 : j = 3
251 : end if
252 1200 : absx = sqrt(-log((1._RKC - absx) * (1._RKC + absx)))
253 1200 : carg = chebyshScaler(j) * absx + chebyshScaler(j + 1)
254 : end if
255 : ! Evaluate the Chebyshev polynomial expansion.
256 6000 : cargTwice = carg + carg
257 : ! W1..W3 are the adjacent elements of the recurrence used to evaluate the Chebyshev polynomial expansion.
258 2000 : w1 = 0._RKC
259 2000 : w2 = 0._RKC
260 6000 : imin = COEF_LIM(j)
261 6000 : i = COEF_LIM(j + 1)
262 : do
263 71266 : w3 = w2
264 71266 : w2 = w1
265 158550 : w1 = (cargTwice * w2 - w3) + coef(i)
266 158550 : i = i - 1
267 158550 : if (imin < i) cycle
268 152550 : exit
269 : end do
270 6000 : erfinv = sign(absx * ((carg * w1 - w2) + coef(imin)), x)
271 : else
272 : ! Use Halley approximation.
273 : ! Verified accuracy up to at least 2 * 10**-26.
274 6013 : absx = 1._RKC - abs(x)
275 6013 : temp = sqrt(-2._RKC * log(0.5_RKC * absx))
276 6013 : erfinv = -0.70711_RKC * ((2.30753_RKC + temp * 0.27061_RKC) / (1._RKC + temp * (0.99229_RKC + temp * 0.04481_RKC)) - temp)
277 18039 : do i = 1, 2
278 12026 : temp = erfc(erfinv) - absx
279 18039 : erfinv = erfinv + temp / (TWO_OVER_SQRTPI * exp(-erfinv**2) - erfinv * temp)
280 : end do
281 6013 : if (x < 0._RKC) erfinv = -erfinv
282 : end if
283 : else
284 0 : CHECK_ASSERTION(__LINE__, -1._RKC < x .and. x < 1._RKC, SK_"@setErfInv(): The condition `-1. < x .and. x < 1.` must hold. x = "//getStr(x)) ! fpp
285 0 : erfinv = sign(huge(x), x)
286 : end if
287 : #else
288 : #error "Unrecognized interface."
289 : #endif
290 : #undef strecok_ENABLED
291 : #undef halley_ENABLED
|
