Line data Source code
1 : /* -*- Mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */
2 : /*
3 : * This file is part of the LibreOffice project.
4 : *
5 : * This Source Code Form is subject to the terms of the Mozilla Public
6 : * License, v. 2.0. If a copy of the MPL was not distributed with this
7 : * file, You can obtain one at http://mozilla.org/MPL/2.0/.
8 : */
9 :
10 : #ifndef SC_OPENCL_OPINLINFUN_statistical
11 : #define SC_OPENCL_OPINLINFUN_statistical
12 52 : std::string MinDecl = "#define Min 2.22507e-308\n";
13 52 : std::string F_PIDecl="#define F_PI 3.1415926\n";
14 52 : std::string fBigInvDecl ="#define fBigInv 2.22045e-016\n";
15 52 : std::string fMachEpsDecl ="#define fMachEps 2.22045e-016\n";
16 52 : std::string fLogDblMaxDecl ="#define fLogDblMax log(1.79769e+308)\n";
17 52 : std::string fHalfMachEpsDecl ="#define fHalfMachEps 0.5*2.22045e-016\n";
18 52 : std::string fMaxGammaArgumentDecl=
19 : "#define fMaxGammaArgument 171.624376956302\n";
20 52 : std::string GetValueDecl=
21 : "double GetValue( double x,double fp,double fDF );\n";
22 52 : std::string GetValue=
23 : "double GetValue( double x,double fp,double fDF )\n"
24 : "{\n"
25 : " return fp - 2 * GetTDist(x, fDF);\n"
26 : "}\n";
27 52 : std::string GetGammaSeriesDecl=
28 : "double GetGammaSeries( double fA, double fX );\n";
29 52 : std::string GetGammaSeries =
30 : "double GetGammaSeries( double fA, double fX )\n"
31 : "{\n"
32 : " double fDenomfactor = fA;\n"
33 : " double fSummand = 1.0/fA;\n"
34 : " double fSum = fSummand;\n"
35 : " int nCount=1;\n"
36 : " do\n"
37 : " {\n"
38 : " fDenomfactor = fDenomfactor + 1.0;\n"
39 : " fSummand = fSummand * fX/fDenomfactor;\n"
40 : " fSum = fSum + fSummand;\n"
41 : " nCount = nCount+1;\n"
42 : " } while ( fSummand/fSum > fHalfMachEps && nCount<=10000);\n"
43 : " if (nCount>10000)\n"
44 : " {\n"
45 : " }\n"
46 : " return fSum;\n"
47 : "}\n";
48 52 : std::string GetGammaContFractionDecl = "double GetGammaContFraction( double "
49 : "fA, double fX );\n";
50 52 : std::string GetGammaContFraction =
51 : "double GetGammaContFraction( double fA, double fX )\n"
52 : "{\n"
53 : " double fBig = 1.0/fBigInv;\n"
54 : " double fCount = 0.0;\n"
55 : " double fNum = 0.0;\n"
56 : " double fY = 1.0 - fA;\n"
57 : " double fDenom = fX + 2.0-fA;\n"
58 : " double fPk = 0.0;\n"
59 : " double fPkm1 = fX + 1.0;\n"
60 : " double fPkm2 = 1.0;\n"
61 : " double fQk = 1.0;\n"
62 : " double fQkm1 = fDenom * fX;\n"
63 : " double fQkm2 = fX;\n"
64 : " double fApprox = fPkm1/fQkm1;\n"
65 : " bool bFinished = false;\n"
66 : " double fR = 0.0;\n"
67 : " do\n"
68 : " {\n"
69 : " fCount = fCount +1.0;\n"
70 : " fY = fY+ 1.0;\n"
71 : " fNum = fY * fCount;\n"
72 : " fDenom = fDenom +2.0;\n"
73 : " fPk = fPkm1 * fDenom - fPkm2 * fNum;\n"
74 : " fQk = fQkm1 * fDenom - fQkm2 * fNum;\n"
75 : " if (fQk != 0.0)\n"
76 : " {\n"
77 : " fR = fPk/fQk;\n"
78 : " bFinished = (fabs( (fApprox - fR)/fR ) <= fHalfMachEps);\n"
79 : " fApprox = fR;\n"
80 : " }\n"
81 : " fPkm2 = fPkm1;\n"
82 : " fPkm1 = fPk;\n"
83 : " fQkm2 = fQkm1;\n"
84 : " fQkm1 = fQk;\n"
85 : " if (fabs(fPk) > fBig)\n"
86 : " {\n"
87 : " fPkm2 = fPkm2 * fBigInv;\n"
88 : " fPkm1 = fPkm1 * fBigInv;\n"
89 : " fQkm2 = fQkm2 * fBigInv;\n"
90 : " fQkm1 = fQkm1 * fBigInv;\n"
91 : " }\n"
92 : " } while (!bFinished && fCount<10000);\n"
93 : " if (!bFinished)\n"
94 : " {\n"
95 : " }\n"
96 : " return fApprox;\n"
97 : "}\n";
98 52 : std::string GetLowRegIGammaDecl = "double GetLowRegIGamma( double "
99 : "fA, double fX );\n";
100 52 : std::string GetLowRegIGamma =
101 : "double GetLowRegIGamma( double fA, double fX )\n"
102 : "{\n"
103 : " double fLnFactor = fA * log(fX) - fX - lgamma(fA);\n"
104 : " double fFactor = exp(fLnFactor);\n"
105 : " if (fX>fA+1.0) \n"
106 : " return 1.0 - fFactor * GetGammaContFraction(fA,fX);\n"
107 : " else\n"
108 : " return fFactor * GetGammaSeries(fA,fX);\n"
109 : "}\n";
110 52 : std::string GetGammaDistDecl = "double GetGammaDist( double fX, double "
111 : "fAlpha, double fLambda );\n";
112 52 : std::string GetGammaDist =
113 : "double GetGammaDist( double fX, double fAlpha, double fLambda )\n"
114 : "{\n"
115 : " if (fX <= 0.0)\n"
116 : " return 0.0;\n"
117 : " else\n"
118 : " return GetLowRegIGamma( fAlpha, fX / fLambda);\n"
119 : "}\n";
120 52 : std::string GetGammaDistPDFDecl = "double GetGammaDistPDF( double fX"
121 : ", double fAlpha, double fLambda );\n";
122 52 : std::string GetGammaDistPDF =
123 : "double GetGammaDistPDF( double fX, double fAlpha, double fLambda )\n"
124 : "{\n"
125 : " if (fX < 0.0)\n"
126 : " return 0.0;\n"
127 : " else if (fX == 0)\n"
128 : " {\n"
129 : " if (fAlpha < 1.0)\n"
130 : " {\n"
131 : " return HUGE_VAL;\n"
132 : " }\n"
133 : " else if (fAlpha == 1)\n"
134 : " {\n"
135 : " return (1.0 / fLambda);\n"
136 : " }\n"
137 : " else\n"
138 : " {\n"
139 : " return 0.0;\n"
140 : " }\n"
141 : " }\n"
142 : " else\n"
143 : " {\n"
144 : " double fXr = fX / fLambda;\n"
145 : " if (fXr > 1.0)\n"
146 : " {\n"
147 : " if (log(fXr) * (fAlpha-1.0) < fLogDblMax &&"
148 : "fAlpha < fMaxGammaArgument)\n"
149 : " {\n"
150 : " return pow( fXr, fAlpha-1.0) * exp(-fXr) / "
151 : "fLambda / tgamma(fAlpha);\n"
152 : " }\n"
153 : " else\n"
154 : " {\n"
155 : " return exp( (fAlpha-1.0) * log(fXr) - fXr - "
156 : "log(fLambda) - lgamma(fAlpha));\n"
157 : " }\n"
158 : " }\n"
159 : " else\n"
160 : " {\n"
161 : " if (fAlpha<fMaxGammaArgument)\n"
162 : " {\n"
163 : " return pow( fXr, fAlpha-1.0) * exp(-fXr) / "
164 : "fLambda / tgamma(fAlpha);\n"
165 : " }\n"
166 : " else\n"
167 : " {\n"
168 : " return pow( fXr, fAlpha-1.0) * exp(-fXr) / "
169 : "fLambda / exp( lgamma(fAlpha));\n"
170 : " }\n"
171 : " }\n"
172 : " }\n"
173 : "}\n";
174 52 : std::string GetBetaDistDecl =
175 : "double GetBetaDist(double fXin, double fAlpha, double fBeta);\n";
176 52 : std::string GetBetaDist =
177 : "double GetBetaDist(double fXin, double fAlpha, double fBeta)\n"
178 : "{\n"
179 : " if (fXin <= 0.0)\n"
180 : " return 0.0;\n"
181 : " if (fXin >= 1.0)\n"
182 : " return 1.0;\n"
183 : " if (fBeta == 1.0)\n"
184 : " return pow(fXin, fAlpha);\n"
185 : " if (fAlpha == 1.0)\n"
186 : " return -expm1(fBeta * log1p(-fXin));\n"
187 : " double fResult;\n"
188 : " double fY = (0.5-fXin)+0.5;\n"
189 : " double flnY = log1p(-fXin);\n"
190 : " double fX = fXin;\n"
191 : " double flnX = log(fXin);\n"
192 : " double fA = fAlpha;\n"
193 : " double fB = fBeta;\n"
194 : " bool bReflect = fXin > fAlpha*pow((fAlpha+fBeta),-1.0);\n"
195 : " if (bReflect)\n"
196 : " {\n"
197 : " fA = fBeta;\n"
198 : " fB = fAlpha;\n"
199 : " fX = fY;\n"
200 : " fY = fXin;\n"
201 : " flnX = flnY;\n"
202 : " flnY = log(fXin);\n"
203 : " }\n"
204 : " fResult = lcl_GetBetaHelperContFrac(fX,fA,fB)*pow(fA,-1.0);\n"
205 : " double fP = fA*pow((fA+fB),-1.0);\n"
206 : " double fQ = fB*pow((fA+fB),-1.0);\n"
207 : " if (fA > 1.0 && fB > 1.0 && fP < 0.97 && fQ < 0.97)\n"
208 : " fResult *= GetBetaDistPDF(fX,fA,fB)*fX*fY;\n"
209 : " else\n"
210 : " fResult *= pow(exp(1.0),(fA*flnX + fB*flnY - GetLogBeta(fA,fB)));\n"
211 : " if (bReflect)\n"
212 : " fResult = 0.5 - fResult + 0.5;\n"
213 : " if (fResult > 1.0)\n"
214 : " fResult = 1.0;\n"
215 : " if (fResult < 0.0)\n"
216 : " fResult = 0.0;\n"
217 : " return fResult;\n"
218 : "}\n";
219 :
220 52 : std::string GetFDistDecl =
221 : "double GetFDist(double x, double fF1, double fF2);\n";
222 52 : std::string GetFDist =
223 : "double GetFDist(double x, double fF1, double fF2)\n"
224 : "{\n"
225 : " double arg = fF2*pow((fF2+fF1*x),-1.0);\n"
226 : " double alpha = fF2*pow(2.0,-1.0);\n"
227 : " double beta = fF1*pow(2.0,-1.0);\n"
228 : " return (GetBetaDist(arg, alpha, beta));\n"
229 : "}\n";
230 52 : std::string GetGammaInvValueDecl = "double"
231 : " GetGammaInvValue(double fAlpha,double fBeta,double fX1 );\n";
232 52 : std::string GetGammaInvValue =
233 : "double GetGammaInvValue(double fAlpha,double fBeta,double fX1 )\n"
234 : "{\n"
235 : " if (fX1 <= 0.0)\n"
236 : " return 0.0;\n"
237 : " else\n"
238 : " {\n"
239 : " double fX=fX1*pow(fBeta,-1.0);\n"
240 : " double fLnFactor = fAlpha * log(fX) - fX - lgamma(fAlpha);\n"
241 : " double fFactor = exp(fLnFactor);\n"
242 : " if (fX>fAlpha+1.0)\n"
243 : " return 1.0 - fFactor * GetGammaContFraction(fAlpha,fX);\n"
244 : " else\n"
245 : " return fFactor * GetGammaSeries(fAlpha,fX);\n"
246 : " }\n"
247 : "}\n";
248 52 : std::string GetFInvValueDecl = "double GetFInvValue(double fF1,double fF2"
249 : ",double fX1 );";
250 52 : std::string GetFInvValue =
251 : "double GetFInvValue(double fF1,double fF2,double fX1 )\n"
252 : "{\n"
253 : " double arg = fF2*pow((fF2+fF1*fX1),-1.0);\n"
254 : " double alpha = fF2*pow(2.0,-1.0);\n"
255 : " double beta = fF1*pow(2.0,-1.0);\n"
256 : " double fXin,fAlpha,fBeta;\n"
257 : " fXin=arg;fAlpha=alpha;fBeta=beta;\n"
258 : " if (fXin <= 0.0)\n"
259 : " return 0.0;\n"
260 : " if (fXin >= 1.0)\n"
261 : " return 1.0;\n"
262 : " if (fBeta == 1.0)\n"
263 : " return pow(fXin, fAlpha);\n"
264 : " if (fAlpha == 1.0)\n"
265 : " return -expm1(fBeta * log1p(-fXin));\n"
266 : " double fResult;\n"
267 : " double fY = (0.5-fXin)+0.5;\n"
268 : " double flnY = log1p(-fXin);\n"
269 : " double fX = fXin;\n"
270 : " double flnX = log(fXin);\n"
271 : " double fA = fAlpha;\n"
272 : " double fB = fBeta;\n"
273 : " bool bReflect = fXin > fAlpha*pow((fAlpha+fBeta),-1.0);\n"
274 : " if (bReflect)\n"
275 : " {\n"
276 : " fA = fBeta;\n"
277 : " fB = fAlpha;\n"
278 : " fX = fY;\n"
279 : " fY = fXin;\n"
280 : " flnX = flnY;\n"
281 : " flnY = log(fXin);\n"
282 : " }\n"
283 : " fResult = lcl_GetBetaHelperContFrac(fX,fA,fB);\n"
284 : " fResult = fResult*pow(fA,-1.0);\n"
285 : " double fP = fA*pow((fA+fB),-1.0);\n"
286 : " double fQ = fB*pow((fA+fB),-1.0);\n"
287 : " double fTemp;\n"
288 : " if (fA > 1.0 && fB > 1.0 && fP < 0.97 && fQ < 0.97)\n"
289 : " fTemp = GetBetaDistPDF(fX,fA,fB)*fX*fY;\n"
290 : " else\n"
291 : " fTemp = exp(fA*flnX + fB*flnY - GetLogBeta(fA,fB));\n"
292 : " fResult *= fTemp;\n"
293 : " if (bReflect)\n"
294 : " fResult = 0.5 - fResult + 0.5;\n"
295 : " if (fResult > 1.0)\n"
296 : " fResult = 1.0;\n"
297 : " if (fResult < 0.0)\n"
298 : " fResult = 0.0;\n"
299 : " return fResult;\n"
300 : "}\n";
301 52 : std::string GetBinomDistPMFDecl =
302 : "double GetBinomDistPMF(double x, double n, double p);";
303 52 : std::string GetBinomDistPMF =
304 : "double GetBinomDistPMF(double x, double n, double p)\n"
305 : "{\n"
306 : " double q = (0.5 - p) + 0.5;\n"
307 : " double fFactor = pow(q, n);\n"
308 : " if (fFactor <= Min)\n"
309 : " {\n"
310 : " fFactor = pow(p, n);\n"
311 : " if (fFactor <= Min)\n"
312 : " return GetBetaDistPDF(p, x + 1.0, n - x + 1.0)*pow((n + 1.0),-1.0);\n"
313 : " else\n"
314 : " {\n"
315 : " uint max = (uint)(n - x);\n"
316 : " for (uint i = 0; i < max && fFactor > 0.0; ++i)\n"
317 : " fFactor *= (n - i)*pow((i + 1),-1.0)*q*pow(p,-1.0);\n"
318 : " return fFactor;\n"
319 : " }\n"
320 : " }\n"
321 : " else\n"
322 : " {\n"
323 : " uint max = (uint)x;\n"
324 : " for (uint i = 0; i < max && fFactor > 0.0; ++i)\n"
325 : " fFactor *= (n - i)*pow((i + 1),-1.0)*p*pow(q,-1.0);\n"
326 : " return fFactor;\n"
327 : " }\n"
328 : "}\n";
329 :
330 52 : std::string lcl_GetBinomDistRangeDecl =
331 : "double lcl_GetBinomDistRange(double n, \n"
332 : "double xs, double xe, double fFactor, double p, double q);";
333 52 : std::string lcl_GetBinomDistRange=
334 : "double lcl_GetBinomDistRange(double n, double xs, double xe,\n"
335 : " double fFactor, double p, double q)\n"
336 : "{\n"
337 : " uint i;\n"
338 : " double fSum;\n"
339 : " uint nXs = (uint)xs;\n"
340 : " for (i = 1; i <= nXs && fFactor > 0.0; ++i)\n"
341 : " fFactor *= (n - i + 1)*pow(i,-1.0)*p*pow(q,-1.0);\n"
342 : " fSum = fFactor;\n"
343 : " uint nXe =(uint)xe;\n"
344 : " for (i = nXs + 1; i <= nXe && fFactor > 0.0; ++i)\n"
345 : " {\n"
346 : " fFactor *= (n - i + 1)*pow(i,-1.0)*p*pow(q,-1.0);\n"
347 : " fSum += fFactor;\n"
348 : " }\n"
349 : " return (fSum > 1.0) ? 1.0 : fSum;\n"
350 : "}\n";
351 :
352 52 : std::string GetLogGammaDecl = "double GetLogGamma(double fZ);\n";
353 52 : std::string GetLogGamma =
354 : "double GetLogGamma(double fZ)\n"
355 : "{\n"
356 : " if (fZ >= fMaxGammaArgument)\n"
357 : " return lcl_GetLogGammaHelper(fZ);\n"
358 : " if (fZ >= 1.0)\n"
359 : " return log(lcl_GetGammaHelper(fZ));\n"
360 : " if (fZ >= 0.5)\n"
361 : " return log( lcl_GetGammaHelper(fZ+1) *pow(fZ,-1.0));\n"
362 : " return lcl_GetLogGammaHelper(fZ+2) - log(fZ+1) - log(fZ);\n"
363 : "}\n";
364 :
365 52 : std::string GetChiDistDecl = "double GetChiDist(double fX, double fDF);\n";
366 52 : std::string GetChiDist =
367 : "double GetChiDist(double fX, double fDF)\n"
368 : "{\n"
369 : " if (fX <= 0.0)\n"
370 : " return 1.0;\n"
371 : " else\n"
372 : " return GetUpRegIGamma( fDF*pow(2.0,-1.0), fX*pow(2.0,-1.0));\n"
373 : "}\n";
374 :
375 52 : std::string GetChiSqDistCDFDecl =
376 : "double GetChiSqDistCDF(double fX, double fDF);\n";
377 52 : std::string GetChiSqDistCDF =
378 : "double GetChiSqDistCDF(double fX, double fDF)\n"
379 : "{\n"
380 : " if (fX <= 0.0)\n"
381 : " return 0.0;"
382 : " else\n"
383 : " return GetLowRegIGamma( fDF*pow(2.0,-1.0), fX*pow(2.0,-1.0));\n"
384 : "}\n";
385 :
386 52 : std::string GetChiSqDistPDFDecl=
387 : "double GetChiSqDistPDF(double fX, double fDF);\n";
388 52 : std::string GetChiSqDistPDF =
389 : "double GetChiSqDistPDF(double fX, double fDF)\n"
390 : "{\n"
391 : " double fValue;\n"
392 : " if (fX <= 0.0)\n"
393 : " return 0.0;\n"
394 : " if (fDF*fX > 1391000.0)\n"
395 : " {\n"
396 : " fValue = exp((0.5*fDF - 1) * log(fX*0.5) - 0.5 * fX - log(2.0)"
397 : " - lgamma(0.5*fDF));\n"
398 : " }\n"
399 : " else\n"
400 : " {\n"
401 : " double fCount;\n"
402 : " if (fmod(fDF,2.0)<0.5)\n"
403 : " {\n"
404 : " fValue = 0.5;\n"
405 : " fCount = 2.0;\n"
406 : " }\n"
407 : " else\n"
408 : " {\n"
409 : " fValue = pow(sqrt(fX*2*F_PI),-1.0);\n"
410 : " fCount = 1.0;\n"
411 : " }\n"
412 : " while ( fCount < fDF)\n"
413 : " {\n"
414 : " fValue *= (fX *pow(fCount,-1.0));\n"
415 : " fCount += 2.0;\n"
416 : " }\n"
417 : " if (fX>=1425.0)\n"
418 : " fValue = exp(log(fValue)-fX*pow(2,-1.0));\n"
419 : " else\n"
420 : " fValue *= exp(-fX*pow(2,-1.0));\n"
421 : " }\n"
422 : " return fValue;\n"
423 : "}\n";
424 :
425 52 : std::string lcl_IterateInverseBetaInvDecl =
426 : "static double lcl_IterateInverseBetaInv(double fp, double fAlpha, \n"
427 : " double fBeta, double fAx, double fBx, bool *rConvError );\n";
428 52 : std::string lcl_IterateInverseBetaInv =
429 : "static double lcl_IterateInverseBetaInv(double fp, double fAlpha, \n"
430 : " double fBeta, double fAx, double fBx, bool *rConvError )\n"
431 : "{\n"
432 : " *rConvError = false;\n"
433 : " double fYEps = 1.0E-307;\n"
434 : " double fXEps = fMachEps;\n"
435 : " if(!(fAx < fBx))\n"
436 : " {\n"
437 : " //print error\n"
438 : " }\n"
439 : " double fAy = fp - GetBetaDist(fAx, fAlpha, fBeta);\n"
440 : " double fBy = fp - GetBetaDist(fBx, fAlpha, fBeta);\n"
441 : " double fTemp;\n"
442 : " unsigned short nCount;\n"
443 : " for (nCount = 0; nCount < 1000 && !lcl_HasChangeOfSign(fAy,fBy);"
444 : " nCount++)\n"
445 : " {\n"
446 : " if (fabs(fAy) <= fabs(fBy))\n"
447 : " {\n"
448 : " fTemp = fAx;\n"
449 : " fAx += 2.0 * (fAx - fBx);\n"
450 : " if (fAx < 0.0)\n"
451 : " fAx = 0.0;\n"
452 : " fBx = fTemp;\n"
453 : " fBy = fAy;\n"
454 : " fAy = fp - GetBetaDist(fAx, fAlpha, fBeta);\n"
455 : " }\n"
456 : " else\n"
457 : " {\n"
458 : " fTemp = fBx;\n"
459 : " fBx += 2.0 * (fBx - fAx);\n"
460 : " fAx = fTemp;\n"
461 : " fAy = fBy;\n"
462 : " fBy = fp - GetBetaDist(fBx, fAlpha, fBeta);\n"
463 : " }\n"
464 : " }\n"
465 : " if (fAy == 0.0)\n"
466 : " return fAx;\n"
467 : " if (fBy == 0.0)\n"
468 : " return fBx;\n"
469 : " if (!lcl_HasChangeOfSign( fAy, fBy))\n"
470 : " {\n"
471 : " *rConvError = true;\n"
472 : " return 0.0;\n"
473 : " }\n"
474 : " double fPx = fAx;\n"
475 : " double fPy = fAy;\n"
476 : " double fQx = fBx;\n"
477 : " double fQy = fBy;\n"
478 : " double fRx = fAx;\n"
479 : " double fRy = fAy;\n"
480 : " double fSx = 0.5 * (fAx + fBx);\n"
481 : " bool bHasToInterpolate = true;\n"
482 : " nCount = 0;\n"
483 : " while ( nCount < 500 && fabs(fRy) > fYEps &&\n"
484 : " (fBx-fAx) > fmax( fabs(fAx), fabs(fBx)) * fXEps )\n"
485 : " {\n"
486 : " if (bHasToInterpolate)\n"
487 : " {\n"
488 : " if (fPy!=fQy && fQy!=fRy && fRy!=fPy)\n"
489 : " {\n"
490 : " fSx = fPx*fRy*fQy*pow(fRy-fPy,-1.0)*pow(fQy-fPy,-1.0)\n"
491 : " + fRx*fQy*fPy*pow(fQy-fRy,-1.0)*pow(fPy-fRy,-1.0)\n"
492 : " + fQx*fPy*fRy*pow(fPy-fQy,-1.0)*pow(fRy-fQy,-1.0);\n"
493 : " bHasToInterpolate = (fAx < fSx) && (fSx < fBx);\n"
494 : " }\n"
495 : " else\n"
496 : " bHasToInterpolate = false;\n"
497 : " }\n"
498 : " if(!bHasToInterpolate)\n"
499 : " {\n"
500 : " fSx = 0.5 * (fAx + fBx);\n"
501 : " fPx = fAx; fPy = fAy;\n"
502 : " fQx = fBx; fQy = fBy;\n"
503 : " bHasToInterpolate = true;\n"
504 : " }\n"
505 : " fPx = fQx; fQx = fRx; fRx = fSx;\n"
506 : " fPy = fQy; fQy = fRy; fRy = fp - GetBetaDist(fSx, fAlpha, fBeta);\n"
507 : " if (lcl_HasChangeOfSign( fAy, fRy))\n"
508 : " {\n"
509 : " fBx = fRx; fBy = fRy;\n"
510 : " }\n"
511 : " else\n"
512 : " {\n"
513 : " fAx = fRx; fAy = fRy;\n"
514 : " }\n"
515 : " bHasToInterpolate = bHasToInterpolate && (fabs(fRy) *"
516 : " 2.0 <= fabs(fQy));\n"
517 : " ++nCount;\n"
518 : " }\n"
519 : " return fRx;\n"
520 : "}\n";
521 :
522 52 : std::string lcl_IterateInverseChiInvDecl =
523 : "static double lcl_IterateInverseChiInv"
524 : "(double fp, double fdf, double fAx, double fBx, bool *rConvError);\n";
525 52 : std::string lcl_IterateInverseChiInv =
526 : "static double lcl_IterateInverseChiInv"
527 : "(double fp, double fdf, double fAx, double fBx, bool *rConvError)\n"
528 : "{\n"
529 : " *rConvError = false;\n"
530 : " double fYEps = 1.0E-307;\n"
531 : " double fXEps = fMachEps;\n"
532 : " if(!(fAx < fBx))\n"
533 : " {\n"
534 : " //print error\n"
535 : " }"
536 : " double fAy = fp - GetChiDist(fAx, fdf);\n"
537 : " double fBy = fp - GetChiDist(fBx, fdf);\n"
538 : " double fTemp;\n"
539 : " unsigned short nCount;\n"
540 : " for (nCount = 0; nCount < 1000 && "
541 : "!lcl_HasChangeOfSign(fAy,fBy); nCount++)\n"
542 : " {\n"
543 : " if (fabs(fAy) <= fabs(fBy))\n"
544 : " {\n"
545 : " fTemp = fAx;\n"
546 : " fAx += 2.0 * (fAx - fBx);\n"
547 : " if (fAx < 0.0)\n"
548 : " fAx = 0.0;\n"
549 : " fBx = fTemp;\n"
550 : " fBy = fAy;\n"
551 : " fAy = fp - GetChiDist(fAx, fdf);\n"
552 : " }\n"
553 : " else\n"
554 : " {\n"
555 : " fTemp = fBx;\n"
556 : " fBx += 2.0 * (fBx - fAx);\n"
557 : " fAx = fTemp;\n"
558 : " fAy = fBy;\n"
559 : " fBy = fp - GetChiDist(fBx, fdf);\n"
560 : " }\n"
561 : " }\n"
562 : " if (fAy == 0.0)\n"
563 : " return fAx;\n"
564 : " if (fBy == 0.0)\n"
565 : " return fBx;\n"
566 : " if (!lcl_HasChangeOfSign( fAy, fBy))\n"
567 : " {\n"
568 : " *rConvError = true;\n"
569 : " return 0.0;\n"
570 : " }\n"
571 : " double fPx = fAx;\n"
572 : " double fPy = fAy;\n"
573 : " double fQx = fBx;\n"
574 : " double fQy = fBy;\n"
575 : " double fRx = fAx;\n"
576 : " double fRy = fAy;\n"
577 : " double fSx = 0.5 * (fAx + fBx);\n"
578 : " bool bHasToInterpolate = true;\n"
579 : " nCount = 0;\n"
580 : " while ( nCount < 500 && fabs(fRy) > fYEps &&\n"
581 : " (fBx-fAx) > fmax( fabs(fAx), fabs(fBx)) * fXEps )\n"
582 : " {\n"
583 : " if (bHasToInterpolate)\n"
584 : " {\n"
585 : " if (fPy!=fQy && fQy!=fRy && fRy!=fPy)\n"
586 : " {\n"
587 : " fSx = fPx * fRy * fQy*pow(fRy-fPy,-1.0)*pow(fQy-fPy,-1.0)\n"
588 : " + fRx * fQy * fPy*pow(fQy-fRy,-1.0)*pow(fPy-fRy,-1.0)\n"
589 : " + fQx * fPy * fRy*pow(fPy-fQy,-1.0)*pow(fRy-fQy,-1.0);\n"
590 : " bHasToInterpolate = (fAx < fSx) && (fSx < fBx);\n"
591 : " }\n"
592 : " else\n"
593 : " bHasToInterpolate = false;\n"
594 : " }\n"
595 : " if(!bHasToInterpolate)\n"
596 : " {\n"
597 : " fSx = 0.5 * (fAx + fBx);\n"
598 : " fPx = fAx; fPy = fAy;\n"
599 : " fQx = fBx; fQy = fBy;\n"
600 : " bHasToInterpolate = true;\n"
601 : " }\n"
602 : " fPx = fQx; fQx = fRx; fRx = fSx;\n"
603 : " fPy = fQy; fQy = fRy; fRy = fp - GetChiDist(fSx, fdf);\n"
604 : " if (lcl_HasChangeOfSign( fAy, fRy))\n"
605 : " {\n"
606 : " fBx = fRx; fBy = fRy;\n"
607 : " }\n"
608 : " else\n"
609 : " {\n"
610 : " fAx = fRx; fAy = fRy;\n"
611 : " }\n"
612 : " bHasToInterpolate = bHasToInterpolate && (fabs(fRy)"
613 : " * 2.0 <= fabs(fQy));\n"
614 : " ++nCount;\n"
615 : " }\n"
616 : " return fRx;\n"
617 : "}\n";
618 :
619 52 : std::string lcl_IterateInverseChiSQInvDecl =
620 : "static double lcl_IterateInverseChiSQInv( double fp, double fdf, \n"
621 : " double fAx, double fBx, bool *rConvError );\n";
622 52 : std::string lcl_IterateInverseChiSQInv =
623 : "static double lcl_IterateInverseChiSQInv( double fp, double fdf, \n"
624 : " double fAx, double fBx, bool *rConvError )\n"
625 : "{\n"
626 : " *rConvError = false;\n"
627 : " double fYEps = 1.0E-307;\n"
628 : " double fXEps = fMachEps;\n"
629 :
630 : " if(!(fAx < fBx))\n"
631 : " {\n"
632 : " //print error\n"
633 : " }\n"
634 : " double fAy = fp - GetChiSqDistCDF(fAx, fdf);\n"
635 : " double fBy = fp - GetChiSqDistCDF(fBx, fdf);\n"
636 : " double fTemp;\n"
637 : " unsigned short nCount;\n"
638 : " for (nCount = 0; nCount < 1000 && !lcl_HasChangeOfSign(fAy,fBy);"
639 : " nCount++)\n"
640 : " {\n"
641 : " if (fabs(fAy) <= fabs(fBy))\n"
642 : " {\n"
643 : " fTemp = fAx;\n"
644 : " fAx += 2.0 * (fAx - fBx);\n"
645 : " if (fAx < 0.0)\n"
646 : " fAx = 0.0;\n"
647 : " fBx = fTemp;\n"
648 : " fBy = fAy;\n"
649 : " fAy = fp - GetChiSqDistCDF(fAx, fdf);\n"
650 : " }\n"
651 : " else\n"
652 : " {\n"
653 : " fTemp = fBx;\n"
654 : " fBx += 2.0 * (fBx - fAx);\n"
655 : " fAx = fTemp;\n"
656 : " fAy = fBy;\n"
657 : " fBy = fp - GetChiSqDistCDF(fBx, fdf);\n"
658 : " }\n"
659 : " }\n"
660 : " if (fAy == 0.0)\n"
661 : " return fAx;\n"
662 : " if (fBy == 0.0)\n"
663 : " return fBx;\n"
664 : " if (!lcl_HasChangeOfSign( fAy, fBy))\n"
665 : " {\n"
666 : " *rConvError = true;\n"
667 : " return 0.0;\n"
668 : " }\n"
669 : " double fPx = fAx;\n"
670 : " double fPy = fAy;\n"
671 : " double fQx = fBx;\n"
672 : " double fQy = fBy;\n"
673 : " double fRx = fAx;\n"
674 : " double fRy = fAy;\n"
675 : " double fSx = 0.5 * (fAx + fBx);\n"
676 : " bool bHasToInterpolate = true;\n"
677 : " nCount = 0;\n"
678 : " while ( nCount < 500 && fabs(fRy) > fYEps &&\n"
679 : " (fBx-fAx) > fmax( fabs(fAx), fabs(fBx)) * fXEps )\n"
680 : " {\n"
681 : " if (bHasToInterpolate)\n"
682 : " {\n"
683 : " if (fPy!=fQy && fQy!=fRy && fRy!=fPy)\n"
684 : " {\n"
685 : " fSx = fPx * fRy * fQy / (fRy-fPy) / (fQy-fPy)\n"
686 : " + fRx * fQy * fPy / (fQy-fRy) / (fPy-fRy)\n"
687 : " + fQx * fPy * fRy / (fPy-fQy) / (fRy-fQy);\n"
688 : " bHasToInterpolate = (fAx < fSx) && (fSx < fBx);\n"
689 : " }\n"
690 : " else\n"
691 : " bHasToInterpolate = false;\n"
692 : " }\n"
693 : " if(!bHasToInterpolate)\n"
694 : " {\n"
695 : " fSx = 0.5 * (fAx + fBx);\n"
696 : " fPx = fAx; fPy = fAy;\n"
697 : " fQx = fBx; fQy = fBy;\n"
698 : " bHasToInterpolate = true;\n"
699 : " }\n"
700 : " fPx = fQx; fQx = fRx; fRx = fSx;\n"
701 : " fPy = fQy; fQy = fRy; fRy = fp - GetChiSqDistCDF(fSx, fdf);\n"
702 : " if (lcl_HasChangeOfSign( fAy, fRy))\n"
703 : " {\n"
704 : " fBx = fRx; fBy = fRy;\n"
705 : " }\n"
706 : " else\n"
707 : " {\n"
708 : " fAx = fRx; fAy = fRy;\n"
709 : " }\n"
710 : " bHasToInterpolate = bHasToInterpolate && (fabs(fRy) * 2.0"
711 : " <= fabs(fQy));\n"
712 : " ++nCount;\n"
713 : " }\n"
714 : " return fRx;\n"
715 : "}\n";
716 :
717 52 : std::string gaussinvDecl = "double gaussinv(double x);\n";
718 52 : std::string gaussinv =
719 : "double gaussinv(double x)\n"
720 : "{\n"
721 : " double q,t,z;\n"
722 : " q=x-0.5;\n"
723 : " if(fabs(q)<=.425)\n"
724 : " {\n"
725 : " t=0.180625-q*q;\n"
726 : " z=\n"
727 : " q*\n"
728 : " (\n"
729 : " (\n"
730 : " (\n"
731 : " (\n"
732 : " (\n"
733 : " (\n"
734 : " (\n"
735 : " t*2509.0809287301226727+33430.575583588128105\n"
736 : " )\n"
737 : " *t+67265.770927008700853\n"
738 : " )\n"
739 : " *t+45921.953931549871457\n"
740 : " )\n"
741 : " *t+13731.693765509461125\n"
742 : " )\n"
743 : " *t+1971.5909503065514427\n"
744 : " )\n"
745 : " *t+133.14166789178437745\n"
746 : " )\n"
747 : " *t+3.387132872796366608\n"
748 : " )\n"
749 : " *pow\n"
750 : " (\n"
751 : " (\n"
752 : " (\n"
753 : " (\n"
754 : " (\n"
755 : " (\n"
756 : " (\n"
757 : " t*5226.495278852854561+28729.085735721942674\n"
758 : " )\n"
759 : " *t+39307.89580009271061\n"
760 : " )\n"
761 : " *t+21213.794301586595867\n"
762 : " )\n"
763 : " *t+5394.1960214247511077\n"
764 : " )\n"
765 : " *t+687.1870074920579083\n"
766 : " )\n"
767 : " *t+42.313330701600911252\n"
768 : " )\n"
769 : " *t+1.0\n"
770 : " , -1.0);\n"
771 : " }\n"
772 : " else\n"
773 : " {\n"
774 : " if(q>0) t=1-x;\n"
775 : " else t=x;\n"
776 : " t=sqrt(-log(t));\n"
777 : " if(t<=5.0)\n"
778 : " {\n"
779 : " t+=-1.6;\n"
780 : " z=\n"
781 : " (\n"
782 : " (\n"
783 : " (\n"
784 : " (\n"
785 : " (\n"
786 : " (\n"
787 : " (\n"
788 : " t*7.7454501427834140764e-4+0.0227238449892691845833\n"
789 : " )\n"
790 : " *t+0.24178072517745061177\n"
791 : " )\n"
792 : " *t+1.27045825245236838258\n"
793 : " )\n"
794 : " *t+3.64784832476320460504\n"
795 : " )\n"
796 : " *t+5.7694972214606914055\n"
797 : " )\n"
798 : " *t+4.6303378461565452959\n"
799 : " )\n"
800 : " *t+1.42343711074968357734\n"
801 : " )\n"
802 : " *pow\n"
803 : " (\n"
804 : " (\n"
805 : " (\n"
806 : " (\n"
807 : " (\n"
808 : " (\n"
809 : " (\n"
810 : " t*1.05075007164441684324e-9+5.475938084995344946e-4\n"
811 : " )\n"
812 : " *t+0.0151986665636164571966\n"
813 : " )\n"
814 : " *t+0.14810397642748007459\n"
815 : " )\n"
816 : " *t+0.68976733498510000455\n"
817 : " )\n"
818 : " *t+1.6763848301838038494\n"
819 : " )\n"
820 : " *t+2.05319162663775882187\n"
821 : " )\n"
822 : " *t+1.0\n"
823 : " , -1.0);\n"
824 : " }\n"
825 : " else\n"
826 : " {\n"
827 : " t+=-5.0;\n"
828 : " z=\n"
829 : " (\n"
830 : " (\n"
831 : " (\n"
832 : " (\n"
833 : " (\n"
834 : " (\n"
835 : " (\n"
836 : " t*2.01033439929228813265e-7+2.71155556874348757815e-5\n"
837 : " )\n"
838 : " *t+0.0012426609473880784386\n"
839 : " )\n"
840 : " *t+0.026532189526576123093\n"
841 : " )\n"
842 : " *t+0.29656057182850489123\n"
843 : " )\n"
844 : " *t+1.7848265399172913358\n"
845 : " )\n"
846 : " *t+5.4637849111641143699\n"
847 : " )\n"
848 : " *t+6.6579046435011037772\n"
849 : " )\n"
850 : " *pow\n"
851 : " (\n"
852 : " (\n"
853 : " (\n"
854 : " (\n"
855 : " (\n"
856 : " (\n"
857 : " (\n"
858 : " t*2.04426310338993978564e-15+1.4215117583164458887e-7\n"
859 : " )\n"
860 : " *t+1.8463183175100546818e-5\n"
861 : " )\n"
862 : " *t+7.868691311456132591e-4\n"
863 : " )\n"
864 : " *t+0.0148753612908506148525\n"
865 : " )\n"
866 : " *t+0.13692988092273580531\n"
867 : " )\n"
868 : " *t+0.59983220655588793769\n"
869 : " )\n"
870 : " *t+1.0\n"
871 : " , -1.0);\n"
872 : " }\n"
873 : " if(q<0.0) z=-z;\n"
874 : " }\n"
875 : " return z;\n"
876 : "}\n";
877 :
878 52 : std::string lcl_GetLogGammaHelperDecl=
879 : "static double lcl_GetLogGammaHelper(double fZ);\n";
880 52 : std::string lcl_GetLogGammaHelper =
881 : "static double lcl_GetLogGammaHelper(double fZ)\n"
882 : "{\n"
883 : " double fg = 6.024680040776729583740234375;\n"
884 : " double fZgHelp = fZ + fg - 0.5;\n"
885 : " return log( lcl_getLanczosSum(fZ)) + (fZ-0.5) * log(fZgHelp) - fZgHelp;\n"
886 : "}\n";
887 52 : std::string lcl_GetGammaHelperDecl=
888 : "static double lcl_GetGammaHelper(double fZ);\n";
889 52 : std::string lcl_GetGammaHelper =
890 : "static double lcl_GetGammaHelper(double fZ)\n"
891 : "{\n"
892 : " double fGamma = lcl_getLanczosSum(fZ);\n"
893 : " double fg = 6.024680040776729583740234375;\n"
894 : " double fZgHelp = fZ + fg - 0.5;\n"
895 : " double fHalfpower = pow( fZgHelp, fZ*pow(2,-1.0) - 0.25);\n"
896 : " fGamma *= fHalfpower;\n"
897 : " fGamma = fGamma*pow(exp(fZgHelp),-1.0);\n"
898 : " fGamma *= fHalfpower;\n"
899 : " fGamma = 120.4;\n"
900 : " if (fZ <= 20.0 && fZ == (int)fZ)\n"
901 : " {\n"
902 : " fGamma = (int)(fGamma+0.5);\n"
903 : " }\n"
904 : " return fGamma;\n"
905 : "}\n";
906 52 : std::string lcl_getLanczosSumDecl=
907 : "static double lcl_getLanczosSum(double fZ);\n";
908 52 : std::string lcl_getLanczosSum =
909 : "static double lcl_getLanczosSum(double fZ) \n"
910 : "{ \n"
911 : " double fNum[13] ={ \n"
912 : " 23531376880.41075968857200767445163675473, \n"
913 : " 42919803642.64909876895789904700198885093, \n"
914 : " 35711959237.35566804944018545154716670596, \n"
915 : " 17921034426.03720969991975575445893111267, \n"
916 : " 6039542586.35202800506429164430729792107, \n"
917 : " 1439720407.311721673663223072794912393972, \n"
918 : " 248874557.8620541565114603864132294232163, \n"
919 : " 31426415.58540019438061423162831820536287, \n"
920 : " 2876370.628935372441225409051620849613599, \n"
921 : " 186056.2653952234950402949897160456992822, \n"
922 : " 8071.672002365816210638002902272250613822, \n"
923 : " 210.8242777515793458725097339207133627117, \n"
924 : " 2.506628274631000270164908177133837338626 \n"
925 : " }; \n"
926 : " double fDenom[13] = { \n"
927 : " 0,\n"
928 : " 39916800,\n"
929 : " 120543840,\n"
930 : " 150917976,\n"
931 : " 105258076,\n"
932 : " 45995730,\n"
933 : " 13339535,\n"
934 : " 2637558,\n"
935 : " 357423,\n"
936 : " 32670,\n"
937 : " 1925,\n"
938 : " 66,\n"
939 : " 1\n"
940 : " };\n"
941 : " double fSumNum;\n"
942 : " double fSumDenom;\n"
943 : " int nI;\n"
944 : " if (fZ<=1.0)\n"
945 : " {\n"
946 : " fSumNum = fNum[12];\n"
947 : " fSumDenom = fDenom[12];\n"
948 : " nI = 11;\n"
949 : " fSumNum = fSumNum*fZ+fNum[nI];\n"
950 : " fSumDenom = fSumDenom*fZ+fDenom[nI];\n"
951 : " nI = 10;\n"
952 : " fSumNum = fSumNum*fZ+fNum[nI];\n"
953 : " fSumDenom = fSumDenom*fZ+fDenom[nI];\n"
954 : " nI = 9;\n"
955 : " fSumNum = fSumNum*fZ+fNum[nI];\n"
956 : " fSumDenom = fSumDenom*fZ+fDenom[nI];\n"
957 : " nI = 8;\n"
958 : " fSumNum = fSumNum*fZ+fNum[nI];\n"
959 : " fSumDenom = fSumDenom*fZ+fDenom[nI];\n"
960 : " nI = 7;\n"
961 : " fSumNum = fSumNum*fZ+fNum[nI];\n"
962 : " fSumDenom = fSumDenom*fZ+fDenom[nI];\n"
963 : " nI = 6;\n"
964 : " fSumNum = fSumNum*fZ+fNum[nI];\n"
965 : " fSumDenom = fSumDenom*fZ+fDenom[nI];\n"
966 : " nI = 5;\n"
967 : " fSumNum = fSumNum*fZ+fNum[nI];\n"
968 : " fSumDenom = fSumDenom*fZ+fDenom[nI];\n"
969 : " nI = 4;\n"
970 : " fSumNum = fSumNum*fZ+fNum[nI];\n"
971 : " fSumDenom = fSumDenom*fZ+fDenom[nI];\n"
972 : " nI = 3;\n"
973 : " fSumNum = fSumNum*fZ+fNum[nI];\n"
974 : " fSumDenom = fSumDenom*fZ+fDenom[nI];\n"
975 : " nI = 2;\n"
976 : " fSumNum = fSumNum*fZ+fNum[nI];\n"
977 : " fSumDenom = fSumDenom*fZ+fDenom[nI];\n"
978 : " nI = 1;\n"
979 : " fSumNum = fSumNum*fZ+fNum[nI];\n"
980 : " fSumDenom = fSumDenom*fZ+fDenom[nI];\n"
981 : " nI = 0;\n"
982 : " fSumNum = fSumNum*fZ+fNum[nI];\n"
983 : " fSumDenom = fSumDenom*fZ+fDenom[nI];\n"
984 : " }\n"
985 : " if (fZ>1.0)\n"
986 : " {\n"
987 : " double fZInv = pow(fZ,-1.0);\n"
988 : " fSumNum = fNum[0];\n"
989 : " fSumDenom = fDenom[0];\n"
990 : " nI = 1;\n"
991 : " fSumNum = fSumNum*fZInv+fNum[nI];\n"
992 : " fSumDenom = fSumDenom*fZInv+fDenom[nI];\n"
993 : " nI = 2;\n"
994 : " fSumNum = fSumNum*fZInv+fNum[nI];\n"
995 : " fSumDenom = fSumDenom*fZInv+fDenom[nI];\n"
996 : " nI = 3;\n"
997 : " fSumNum = fSumNum*fZInv+fNum[nI];\n"
998 : " fSumDenom = fSumDenom*fZInv+fDenom[nI];\n"
999 : " nI = 4;\n"
1000 : " fSumNum = fSumNum*fZInv+fNum[nI];\n"
1001 : " fSumDenom = fSumDenom*fZInv+fDenom[nI];\n"
1002 : " nI = 5;\n"
1003 : " fSumNum = fSumNum*fZInv+fNum[nI];\n"
1004 : " fSumDenom = fSumDenom*fZInv+fDenom[nI];\n"
1005 : " nI = 6;\n"
1006 : " fSumNum = fSumNum*fZInv+fNum[nI];\n"
1007 : " fSumDenom = fSumDenom*fZInv+fDenom[nI];\n"
1008 : " nI = 7;\n"
1009 : " fSumNum = fSumNum*fZInv+fNum[nI];\n"
1010 : " fSumDenom = fSumDenom*fZInv+fDenom[nI];\n"
1011 : " nI = 8;\n"
1012 : " fSumNum = fSumNum*fZInv+fNum[nI];\n"
1013 : " fSumDenom = fSumDenom*fZInv+fDenom[nI];\n"
1014 : " nI = 9;\n"
1015 : " fSumNum = fSumNum*fZInv+fNum[nI];\n"
1016 : " fSumDenom = fSumDenom*fZInv+fDenom[nI];\n"
1017 : " nI = 10;\n"
1018 : " fSumNum = fSumNum*fZInv+fNum[nI];\n"
1019 : " fSumDenom = fSumDenom*fZInv+fDenom[nI];\n"
1020 : " nI = 11;\n"
1021 : " fSumNum = fSumNum*fZInv+fNum[nI];\n"
1022 : " fSumDenom = fSumDenom*fZInv+fDenom[nI];\n"
1023 : " nI = 12;\n"
1024 : " fSumNum = fSumNum*fZInv+fNum[nI];\n"
1025 : " fSumDenom = fSumDenom*fZInv+fDenom[nI];\n"
1026 : " }\n"
1027 : " return fSumNum*pow(fSumDenom,-1.0);\n"
1028 : "}\n";
1029 :
1030 52 : std::string GetUpRegIGammaDecl=
1031 : " double GetUpRegIGamma( double fA, double fX ) ;\n";
1032 52 : std::string GetUpRegIGamma =
1033 : "double GetUpRegIGamma( double fA, double fX )\n"
1034 : "{\n"
1035 : " double fLnFactor= fA*log(fX)-fX-lgamma(fA);\n"
1036 : " double fFactor = exp(fLnFactor); \n"
1037 : " if (fX>fA+1.0) \n"
1038 : " return fFactor * GetGammaContFraction(fA,fX);\n"
1039 : " else \n"
1040 : " return 1.0 -fFactor * GetGammaSeries(fA,fX);\n"
1041 : "}\n";
1042 :
1043 52 : std::string lcl_HasChangeOfSignDecl=
1044 : "static inline bool lcl_HasChangeOfSign( double u, double w );\n";
1045 52 : std::string lcl_HasChangeOfSign =
1046 : "static inline bool lcl_HasChangeOfSign( double u, double w )\n"
1047 : "{\n"
1048 : " return (u < 0.0 && w > 0.0) || (u > 0.0 && w < 0.0);\n"
1049 : "}\n";
1050 :
1051 52 : std::string GetTDistDecl=" double GetTDist(double T, double fDF);\n";
1052 52 : std::string GetTDist =
1053 : "double GetTDist(double T, double fDF)\n"
1054 : "{\n"
1055 : " return 0.5 * GetBetaDist(fDF*pow(fDF+T*T,-1.0),fDF*pow(2.0,-1.0), 0.5);\n"
1056 : "}\n";
1057 :
1058 52 : std::string GetBetaDecl=" double GetBeta(double fAlpha, double fBeta);\n";
1059 52 : std::string GetBeta =
1060 : "double GetBeta(double fAlpha, double fBeta)\n"
1061 : "{\n"
1062 : " double fA;\n"
1063 : " double fB;\n"
1064 : " fAlpha>fBeta?(fA = fAlpha,fB = fBeta):(fA = fBeta,fB = fAlpha);\n"
1065 : " double fAB = fA + fB;\n"
1066 :
1067 : " if (fAB < fMaxGammaArgument)\n"
1068 : " return tgamma(fA)*pow(tgamma(fAB),-1.0)*tgamma(fB);\n"
1069 : " double fgm = 5.524680040776729583740234375;\n"
1070 : " double fLanczos = lcl_getLanczosSum(fA)*lcl_getLanczosSum(fB)\n"
1071 : " *pow(lcl_getLanczosSum(fAB),-1.0);\n"
1072 : " fLanczos *= sqrt(((fAB + fgm)*pow(fA + fgm,-1.0))*pow(fB + fgm,-1.0));\n"
1073 : " return fLanczos * pow(exp(1.0),(-fA*log1p(fB*pow(fA + fgm,-1.0)))"
1074 : " - fB*log1p(fA*pow(fB + fgm,-1.0)) - fgm);\n"
1075 : "}\n";
1076 :
1077 52 : std::string GetLogBetaDecl=
1078 : " double GetLogBeta(double fAlpha, double fBeta);\n";
1079 52 : std::string GetLogBeta =
1080 : "double GetLogBeta(double fAlpha, double fBeta)\n"
1081 : "{\n"
1082 : " double fA;\n"
1083 : " double fB;\n"
1084 : " fAlpha>fBeta?(fA = fAlpha,fB = fBeta):(fA = fBeta,fB = fAlpha);\n"
1085 : " double fgm = 5.524680040776729583740234375;\n"
1086 :
1087 : " double fLanczos = lcl_getLanczosSum(fA)*lcl_getLanczosSum(fB)*\n"
1088 : " pow(lcl_getLanczosSum(fA + fB),-1.0);\n"
1089 : " double fResult= -fA *log1p(fB*pow(fA + fgm,-1.0))"
1090 : "-fB *log1p(fA*pow(fB + fgm,-1.0))-fgm;\n"
1091 : " fResult += log(fLanczos)+0.5*(log(fA + fB + fgm) - log(fA + fgm)\n"
1092 : " - log(fB + fgm));\n"
1093 : " return fResult;\n"
1094 : "}\n";
1095 :
1096 52 : std::string GetBetaDistPDFDecl=
1097 : "double GetBetaDistPDF(double fX, double fA, double fB);\n";
1098 52 : std::string GetBetaDistPDF =
1099 : "double GetBetaDistPDF(double fX, double fA, double fB)\n"
1100 : "{\n"
1101 : " if (fA == 1.0) \n"
1102 : " {\n"
1103 : " if (fB == 1.0)\n"
1104 : " return 1.0;\n"
1105 : " if (fB == 2.0)\n"
1106 : " return -2.0*fX + 2.0;\n"
1107 : " if (fX == 1.0 && fB < 1.0)\n"
1108 : " {\n"
1109 : " return HUGE_VAL;\n"
1110 : " }\n"
1111 : " if (fX <= 0.01)\n"
1112 : " return fB + fB * expm1((fB-1.0) * log1p(-fX));\n"
1113 : " else \n"
1114 : " return fB * pow(0.5-fX+0.5,fB-1.0);\n"
1115 : " }\n"
1116 : " if (fB == 1.0) \n"
1117 : " {\n"
1118 : " if (fA == 2.0)\n"
1119 : " return fA * fX;\n"
1120 : " if (fX == 0.0 && fA < 1.0)\n"
1121 : " {\n"
1122 : " return HUGE_VAL;\n"
1123 : " }\n"
1124 : " return fA * pow(fX,fA-1);\n"
1125 : " }\n"
1126 : " if (fX <= 0.0)\n"
1127 : " {\n"
1128 : " if (fA < 1.0 && fX == 0.0)\n"
1129 : " {\n"
1130 : " return HUGE_VAL;\n"
1131 : " }\n"
1132 : " else\n"
1133 : " return 0.0;\n"
1134 : " }\n"
1135 : " if (fX >= 1.0)\n"
1136 : " {\n"
1137 : " if (fB < 1.0 && fX == 1.0)\n"
1138 : " {\n"
1139 : " return HUGE_VAL;\n"
1140 : " }\n"
1141 : " else \n"
1142 : " return 0.0;\n"
1143 : " }\n"
1144 : " double fLogDblMax = log( 1.79769e+308 );\n"
1145 : " double fLogDblMin = log( 2.22507e-308 );\n"
1146 : " double fLogY = (fX < 0.1) ? log1p(-fX) : log(0.5-fX+0.5);\n"
1147 : " double fLogX = log(fX);\n"
1148 : " double fAm1LogX = (fA-1.0) * fLogX;\n"
1149 : " double fBm1LogY = (fB-1.0) * fLogY;\n"
1150 : " double fLogBeta = GetLogBeta(fA,fB);\n"
1151 : " if ( fAm1LogX < fLogDblMax && fAm1LogX > fLogDblMin\n"
1152 : " && fBm1LogY < fLogDblMax && fBm1LogY > fLogDblMin\n"
1153 : " && fLogBeta < fLogDblMax && fLogBeta > fLogDblMin\n"
1154 : " && fAm1LogX + fBm1LogY < fLogDblMax && fAm1LogX + fBm1LogY > \n"
1155 : " fLogDblMin)\n"
1156 : " return pow(fX,fA-1.0)*pow(0.5-fX+0.5,fB-1.0)"
1157 : "*pow(GetBeta(fA,fB),-1.0);\n"
1158 : " else \n"
1159 : " return exp( fAm1LogX + fBm1LogY - fLogBeta);\n"
1160 : "}\n";
1161 :
1162 52 : std::string lcl_GetBetaHelperContFracDecl=
1163 : "double lcl_GetBetaHelperContFrac(double fX, double fA, double fB);\n";
1164 52 : std::string lcl_GetBetaHelperContFrac =
1165 : "double lcl_GetBetaHelperContFrac(double fX, double fA, double fB)\n"
1166 : "{ \n"
1167 :
1168 : " double a1, b1, a2, b2, fnorm, apl2m, d2m, d2m1, cfnew, cf;\n"
1169 : " a1 = 1.0; b1 = 1.0;\n"
1170 : " b2 = 1.0 - (fA+fB)*pow(fA+1.0,-1.0)*fX;\n"
1171 : " b2==0.0?(a2 = 0.0,fnorm = 1.0,cf = 1.0):\n"
1172 : " (a2 = 1.0,fnorm = pow(b2,-1.0),cf = a2*fnorm);\n"
1173 : " cfnew = 1.0;\n"
1174 : " double rm = 1.0;\n"
1175 : " double fMaxIter = 50000.0;\n"
1176 : " bool bfinished = false;\n"
1177 : " do\n"
1178 : " {\n"
1179 : " apl2m = fA + 2.0*rm;\n"
1180 : " d2m = (rm*(fB-rm))*fX*pow(apl2m*(apl2m-1.0),-1.0);\n"
1181 : " d2m1 = -((fA+rm)*(fA+rm+fB))*fX*pow(apl2m*(apl2m+1.0),-1.0);\n"
1182 : " a1 = (a2+d2m*a1)*fnorm;\n"
1183 : " b1 = (b2+d2m*b1)*fnorm;\n"
1184 : " a2 = a1 + d2m1*a2*fnorm;\n"
1185 : " b2 = b1 + d2m1*b2*fnorm;\n"
1186 : " if (b2 != 0.0) \n"
1187 : " {\n"
1188 : " fnorm = pow(b2,-1.0);\n"
1189 : " cfnew = a2*fnorm;\n"
1190 : " bfinished = (fabs(cf-cfnew) < fabs(cf)*fMachEps);\n"
1191 : " }\n"
1192 : " cf = cfnew;\n"
1193 : " rm += 1.0;\n"
1194 : " }\n"
1195 : " while (rm < fMaxIter && !bfinished);\n"
1196 : " return cf;\n"
1197 : "}\n";
1198 :
1199 52 : std::string lcl_IterateInverseDecl=
1200 : "double lcl_IterateInverse("
1201 : "double fAx, double fBx, bool* rConvError,double fp,double fDF );\n";
1202 52 : std::string lcl_IterateInverse =
1203 : "double lcl_IterateInverse( "
1204 : "double fAx, double fBx, bool* rConvError,double fp,double fDF )\n"
1205 : "{\n"
1206 : " *rConvError = false;\n"
1207 : " double fYEps = 1.0E-307;\n"
1208 : " double fXEps =DBL_EPSILON;\n"
1209 : " if(fAx>fBx)\n"
1210 : " return DBL_MAX;\n"
1211 : " double fAy = GetValue(fAx,fp,fDF);\n"
1212 : " double fBy = GetValue(fBx,fp,fDF);\n"
1213 : " double fTemp;\n"
1214 : " unsigned short nCount;\n"
1215 : " double inter;\n"
1216 : " bool sign;\n"
1217 : " for (nCount =0;nCount<1000&&!lcl_HasChangeOfSign(fAy,fBy);nCount++)\n"
1218 : " {\n"
1219 : " inter = 2.0 * (fAx - fBx);\n"
1220 : " if (fabs(fAy) <= fabs(fBy)) \n"
1221 : " {\n"
1222 : " sign = true;\n"
1223 : " fTemp = fAx;\n"
1224 : " fAx += inter;\n"
1225 : " if (fAx < 0.0)\n"
1226 : " fAx = 0.0;\n"
1227 : " fBx = fTemp;\n"
1228 : " fBy = fAy;\n"
1229 : " fTemp = fAx;\n"
1230 : " }\n"
1231 : " else\n"
1232 : " {\n"
1233 : " sign = false;\n"
1234 : " fTemp = fBx;\n"
1235 : " fBx -= inter;\n"
1236 : " fAx = fTemp;\n"
1237 : " fAy = fBy;\n"
1238 : " fTemp = fBx;\n"
1239 : " }\n"
1240 : " fTemp = GetValue(fTemp,fp,fDF);\n"
1241 : " sign?(fAy = fTemp):(fBy = fTemp);\n"
1242 : " }\n"
1243 : " if (fAy == 0.0)\n"
1244 : " return fAx;\n"
1245 : " if (fBy == 0.0)\n"
1246 : " return fBx;\n"
1247 : " if (!lcl_HasChangeOfSign( fAy, fBy))\n"
1248 : " {\n"
1249 : " *rConvError = true;\n"
1250 : " return 0.0;\n"
1251 : " }\n"
1252 : " double fPx = fAx;\n"
1253 : " double fPy = fAy;\n"
1254 : " double fQx = fBx;\n"
1255 : " double fQy = fBy;\n"
1256 : " double fRx = fAx;\n"
1257 : " double fRy = fAy;\n"
1258 : " double fSx = 0.5 * (fAx + fBx); \n"
1259 : " bool bHasToInterpolate = true;\n"
1260 : " nCount = 0;\n"
1261 : " while ( nCount < 500 && fabs(fRy) > fYEps &&\n"
1262 : " (fBx-fAx) > max( fabs(fAx), fabs(fBx)) * fXEps )\n"
1263 : " {\n"
1264 : " if (bHasToInterpolate)\n"
1265 : " {\n"
1266 : " if (fPy!=fQy && fQy!=fRy && fRy!=fPy)\n"
1267 : " {\n"
1268 : " fSx = fPx * fRy * fQy * pow(fRy-fPy,-1.0)*pow(fQy-fPy,-1.0)\n"
1269 : " + fRx * fQy * fPy * pow(fQy-fRy,-1.0)*pow(fPy-fRy,-1.0)\n"
1270 : " + fQx * fPy * fRy * pow(fPy-fQy,-1.0)*pow(fRy-fQy,-1.0);\n"
1271 : " bHasToInterpolate = (fAx < fSx) && (fSx < fBx);\n"
1272 : " }\n"
1273 : " else\n"
1274 : " bHasToInterpolate = false;\n"
1275 : " }\n"
1276 : " if(!bHasToInterpolate)\n"
1277 : " {\n"
1278 : " fSx = 0.5 * (fAx + fBx);\n"
1279 : " \n"
1280 : " fPx = fAx; fPy = fAy;\n"
1281 : " fQx = fBx; fQy = fBy;\n"
1282 : " bHasToInterpolate = true;\n"
1283 : " }\n"
1284 : " fPx = fQx; fQx = fRx; fRx = fSx;\n"
1285 : " fPy = fQy; fQy = fRy; fRy = GetValue(fSx,fp,fDF);\n"
1286 : " lcl_HasChangeOfSign( fAy, fRy)?(fBx = fRx,fBy = fRy):\n"
1287 : " (fAx = fRx,fAy = fRy);\n"
1288 : " bHasToInterpolate =\n"
1289 : " bHasToInterpolate && (fabs(fRy) * 2.0 <= fabs(fQy));\n"
1290 : " ++nCount;\n"
1291 : " }\n"
1292 : " return fRx;\n"
1293 : "}\n";
1294 52 : std::string phiDecl=
1295 : "double phi(double x);\n";
1296 52 : std::string phi =
1297 : "double phi(double x)\n"
1298 : "{\n"
1299 : " return 0.39894228040143268 * exp(-(x * x) / 2.0);\n"
1300 : "}\n";
1301 52 : std::string taylorDecl =
1302 : "double taylor(double* pPolynom, uint nMax, double x);\n";
1303 52 : std::string taylor =
1304 : "double taylor(double* pPolynom, uint nMax, double x)\n"
1305 : "{\n"
1306 : " double nVal = pPolynom[nMax];\n"
1307 : " for (short i = nMax-1; i >= 0; i--)\n"
1308 : " {\n"
1309 : " nVal = pPolynom[i] + (nVal * x);\n"
1310 : " }\n"
1311 : " return nVal;\n"
1312 : "}";
1313 52 : std::string gaussDecl = "double gauss(double x);\n";
1314 52 : std::string gauss =
1315 : "double gauss(double x)\n"
1316 : "{\n"
1317 : " double xAbs = fabs(x);\n"
1318 : " uint xShort = (uint)(floor(xAbs));\n"
1319 : " double nVal = 0.0;\n"
1320 : " if (xShort == 0)\n"
1321 : " {\n"
1322 : " double t0[] =\n"
1323 : " { 0.39894228040143268, -0.06649038006690545, 0.00997355701003582,\n"
1324 : " -0.00118732821548045, 0.00011543468761616, -0.00000944465625950,\n"
1325 : " 0.00000066596935163, -0.00000004122667415, 0.00000000227352982,\n"
1326 : " 0.00000000011301172, 0.00000000000511243, -0.00000000000021218 };\n"
1327 : " nVal = taylor(t0, 11, (xAbs * xAbs)) * xAbs;\n"
1328 : " }\n"
1329 : " else if ((xShort >= 1) && (xShort <= 2))\n"
1330 : " {\n"
1331 : " double t2[] =\n"
1332 : " { 0.47724986805182079, 0.05399096651318805, -0.05399096651318805,\n"
1333 : " 0.02699548325659403, -0.00449924720943234, -0.00224962360471617,\n"
1334 : " 0.00134977416282970, -0.00011783742691370, -0.00011515930357476,\n"
1335 : " 0.00003704737285544, 0.00000282690796889, -0.00000354513195524,\n"
1336 : " 0.00000037669563126, 0.00000019202407921, -0.00000005226908590,\n"
1337 : " -0.00000000491799345, 0.00000000366377919, -0.00000000015981997,\n"
1338 : " -0.00000000017381238, 0.00000000002624031, 0.00000000000560919,\n"
1339 : " -0.00000000000172127, -0.00000000000008634, 0.00000000000007894 };\n"
1340 : " nVal = taylor(t2, 23, (xAbs - 2.0));\n"
1341 : " }\n"
1342 : " else if ((xShort >= 3) && (xShort <= 4))\n"
1343 : " {\n"
1344 : " double t4[] =\n"
1345 : " { 0.49996832875816688, 0.00013383022576489, -0.00026766045152977,\n"
1346 : " 0.00033457556441221, -0.00028996548915725, 0.00018178605666397,\n"
1347 : " -0.00008252863922168, 0.00002551802519049, -0.00000391665839292,\n"
1348 : " -0.00000074018205222, 0.00000064422023359, -0.00000017370155340,\n"
1349 : " 0.00000000909595465, 0.00000000944943118, -0.00000000329957075,\n"
1350 : " 0.00000000029492075, 0.00000000011874477, -0.00000000004420396,\n"
1351 : " 0.00000000000361422, 0.00000000000143638, -0.00000000000045848 };\n"
1352 : " nVal = taylor(t4, 20, (xAbs - 4.0));\n"
1353 : " }\n"
1354 : " else\n"
1355 : " {\n"
1356 : " double asympt[] = { -1.0, 1.0, -3.0, 15.0, -105.0 };\n"
1357 : " nVal = 0.5 + phi(xAbs) * taylor(asympt, 4, 1.0/(xAbs * xAbs))/xAbs;\n"
1358 : " }\n"
1359 : " if (x < 0.0)\n"
1360 : " return -nVal;\n"
1361 : " else\n"
1362 : " return nVal;\n"
1363 : "}\n";
1364 52 : std::string lcl_Erfc0600Decl=
1365 : "void lcl_Erfc0600( double x, double *fVal );\n";
1366 52 : std::string lcl_Erfc0600 =
1367 : "void lcl_Erfc0600( double x, double *fVal )\n"
1368 : "{\n"
1369 : " double fPSum = 0.0;\n"
1370 : " double fQSum = 0.0;\n"
1371 : " double fXPow = 1.0;\n"
1372 : " double *pn;\n"
1373 : " double *qn;\n"
1374 : " if ( x < 2.2 )\n"
1375 : " {\n"
1376 : " double pn22[] = { \n"
1377 : " 9.99999992049799098E-1, \n"
1378 : " 1.33154163936765307, \n"
1379 : " 8.78115804155881782E-1, \n"
1380 : " 3.31899559578213215E-1, \n"
1381 : " 7.14193832506776067E-2, \n"
1382 : " 7.06940843763253131E-3 \n"
1383 : " }; \n"
1384 : " double qn22[] = { \n"
1385 : " 1.00000000000000000, \n"
1386 : " 2.45992070144245533, \n"
1387 : " 2.65383972869775752, \n"
1388 : " 1.61876655543871376, \n"
1389 : " 5.94651311286481502E-1, \n"
1390 : " 1.26579413030177940E-1, \n"
1391 : " 1.25304936549413393E-2 \n"
1392 : " }; \n"
1393 : " pn = pn22; \n"
1394 : " qn = qn22; \n"
1395 : " } \n"
1396 : " else \n"
1397 : " \n"
1398 : " { \n"
1399 : " double pn60[] = {\n"
1400 : " 9.99921140009714409E-1,\n"
1401 : " 1.62356584489366647,\n"
1402 : " 1.26739901455873222,\n"
1403 : " 5.81528574177741135E-1,\n"
1404 : " 1.57289620742838702E-1,\n"
1405 : " 2.25716982919217555E-2\n"
1406 : " };\n"
1407 : " double qn60[] = {\n"
1408 : " 1.00000000000000000,\n"
1409 : " 2.75143870676376208,\n"
1410 : " 3.37367334657284535,\n"
1411 : " 2.38574194785344389,\n"
1412 : " 1.05074004614827206,\n"
1413 : " 2.78788439273628983E-1,\n"
1414 : " 4.00072964526861362E-2\n"
1415 : " };\n"
1416 : " pn = pn60;\n"
1417 : " qn = qn60;\n"
1418 : " }\n"
1419 : " for ( unsigned int i = 0; i < 6; ++i ) \n"
1420 : " {\n"
1421 : " fPSum += pn[i]*fXPow;\n"
1422 : " fQSum += qn[i]*fXPow;\n"
1423 : " fXPow *= x;\n"
1424 : " }\n"
1425 : " fQSum += qn[6]*fXPow;\n"
1426 : " *fVal = exp((-1.0)*x*x)*fPSum*pow(fQSum, -1.0);\n"
1427 : " }\n";
1428 52 : std::string lcl_Erfc2654Decl=
1429 : "void lcl_Erfc2654( double x, double *fVal );\n";
1430 52 : std::string lcl_Erfc2654 =
1431 : "void lcl_Erfc2654( double x, double *fVal )\n"
1432 : "{\n"
1433 : " double pn[] = {\n"
1434 : " 5.64189583547756078E-1,\n"
1435 : " 8.80253746105525775,\n"
1436 : " 3.84683103716117320E1,\n"
1437 : " 4.77209965874436377E1,\n"
1438 : " 8.08040729052301677\n"
1439 : " };\n"
1440 : " double qn[] = {\n"
1441 : " 1.00000000000000000,\n"
1442 : " 1.61020914205869003E1,\n"
1443 : " 7.54843505665954743E1,\n"
1444 : " 1.12123870801026015E2,\n"
1445 : " 3.73997570145040850E1\n"
1446 : " };\n"
1447 : "\n"
1448 : " double fPSum = 0.0;\n"
1449 : " double fQSum = 0.0;\n"
1450 : " double fXPow = 1.0;\n"
1451 : "\n"
1452 : " for ( unsigned int i = 0; i <= 4; ++i )\n"
1453 : " {\n"
1454 : " fPSum += pn[i]*fXPow; \n"
1455 : " fQSum += qn[i]*fXPow;\n"
1456 : " fXPow *= pow(x*x, -1.0);\n"
1457 : " }\n"
1458 : " *fVal = exp((-1.0)*x*x)*fPSum*pow(x*fQSum, -1.0);\n"
1459 : "}\n";
1460 52 : std::string lcl_Erf0065Decl=
1461 : "void lcl_Erf0065( double x, double *fVal );\n";
1462 52 : std::string lcl_Erf0065 =
1463 : "void lcl_Erf0065( double x, double *fVal )\n"
1464 : " {\n"
1465 : " double pn[] = {\n"
1466 : " 1.12837916709551256,\n"
1467 : " 1.35894887627277916E-1,\n"
1468 : " 4.03259488531795274E-2,\n"
1469 : " 1.20339380863079457E-3,\n"
1470 : " 6.49254556481904354E-5\n"
1471 : " };\n"
1472 : " double qn[] = {\n"
1473 : " 1.00000000000000000,\n"
1474 : " 4.53767041780002545E-1,\n"
1475 : " 8.69936222615385890E-2,\n"
1476 : " 8.49717371168693357E-3,\n"
1477 : " 3.64915280629351082E-4\n"
1478 : " };\n"
1479 : " double fPSum = 0.0;\n"
1480 : " double fQSum = 0.0;\n"
1481 : " double fXPow = 1.0;\n"
1482 : " for ( unsigned int i = 0; i <= 4; ++i )\n"
1483 : " {\n"
1484 : " fPSum += pn[i]*fXPow;\n"
1485 : " fQSum += qn[i]*fXPow;\n"
1486 : " fXPow *= x*x;\n"
1487 : " }\n"
1488 : " *fVal = x * fPSum * pow(fQSum, -1.0);\n"
1489 : " }\n";
1490 52 : std::string rtl_math_erf_rdDecl=
1491 : "double rtl_math_erf_rd( double x );\n";
1492 52 : std::string rtl_math_erf_rd =
1493 : " double rtl_math_erf_rd( double x )\n"
1494 : " {\n"
1495 : " if( x == 0.0 )\n"
1496 : " return 0.0;\n"
1497 : " bool bNegative = false;\n"
1498 : " if ( x < 0.0 )\n"
1499 : " {\n"
1500 : " x = fabs( x );\n"
1501 : " bNegative = true;\n"
1502 : " }\n"
1503 : " double fErf = 1.0;\n"
1504 : " if ( x < 1.0e-10 )\n"
1505 : " fErf = (double) (x*1.1283791670955125738961589031215452);\n"
1506 : " else if ( x < 0.65 )\n"
1507 : " lcl_Erf0065( x, &fErf );\n"
1508 : " if ( bNegative )\n"
1509 : " fErf *= (-1.0);\n"
1510 : " return fErf;\n"
1511 : " }\n";
1512 52 : std::string rtl_math_erfc_rdDecl=
1513 : "double rtl_math_erfc_rd( double x );\n";
1514 52 : std::string rtl_math_erfc_rd =
1515 : " double rtl_math_erfc_rd( double x )\n"
1516 : " {\n"
1517 : " if ( x == 0.0 )\n"
1518 : " return 1.0;\n"
1519 : " bool bNegative = false;\n"
1520 : " if ( x < 0.0 )\n"
1521 : " {\n"
1522 : " x = fabs( x );\n"
1523 : " bNegative = true;\n"
1524 : " }\n"
1525 : " double fErfc = 0.0;\n"
1526 : " if ( x >= 0.65 )\n"
1527 : " {\n"
1528 : " if ( x < 6.0 )\n"
1529 : " lcl_Erfc0600( x, &fErfc );\n"
1530 : " else\n"
1531 : " lcl_Erfc2654( x, &fErfc );\n"
1532 : " }\n"
1533 : " else\n"
1534 : " fErfc = 1.0 - rtl_math_erf_rd( x );\n"
1535 : " if ( bNegative )\n"
1536 : " fErfc = 2.0 - fErfc;\n"
1537 : " return fErfc;\n"
1538 : " }\n";
1539 : #endif
1540 :
1541 : /* vim:set shiftwidth=4 softtabstop=4 expandtab: */
|
