LibreOffice Module sc (master) 1
opinlinefun_statistical.cxx
Go to the documentation of this file.
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
12std::string MinDecl = "#define Min 2.22507e-308\n";
13std::string F_PIDecl="#define M_PI 3.1415926535897932384626433832795\n";
14std::string fBigInvDecl ="#define fBigInv 2.22045e-016\n";
15std::string fMachEpsDecl ="#define fMachEps 2.22045e-016\n";
16std::string fLogDblMaxDecl ="#define fLogDblMax log(1.79769e+308)\n";
17std::string fHalfMachEpsDecl ="#define fHalfMachEps 0.5*2.22045e-016\n";
19"#define fMaxGammaArgument 171.624376956302\n";
20std::string GetValueDecl=
21"double GetValue( double x,double fp,double fDF );\n";
22std::string GetValue=
23"double GetValue( double x,double fp,double fDF )\n"
24"{\n"
25" return fp - 2 * GetTDist(x, fDF);\n"
26"}\n";
28"double GetGammaSeries( double fA, double fX );\n";
29std::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";
48std::string GetGammaContFractionDecl = "double GetGammaContFraction( double "
49"fA, double fX );\n";
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";
98std::string GetLowRegIGammaDecl = "double GetLowRegIGamma( double "
99"fA, double fX );\n";
100std::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";
110std::string GetGammaDistDecl = "double GetGammaDist( double fX, double "
111"fAlpha, double fLambda );\n";
112std::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";
120std::string GetGammaDistPDFDecl = "double GetGammaDistPDF( double fX"
121", double fAlpha, double fLambda );\n";
122std::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";
174std::string GetBetaDistDecl =
175"double GetBetaDist(double fXin, double fAlpha, double fBeta);\n";
176std::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
220std::string GetFDistDecl =
221 "double GetFDist(double x, double fF1, double fF2);\n";
222std::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";
230std::string GetGammaInvValueDecl = "double"
231" GetGammaInvValue(double fAlpha,double fBeta,double fX1 );\n";
232std::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";
248std::string GetFInvValueDecl = "double GetFInvValue(double fF1,double fF2"
249",double fX1 );";
250std::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";
302 "double GetBinomDistPMF(double x, double n, double p);";
303std::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
331 "double lcl_GetBinomDistRange(double n, \n"
332"double xs, double xe, double fFactor, double p, double q);";
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
352std::string GetLogGammaDecl = "double GetLogGamma(double fZ);\n";
353std::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
365std::string GetChiDistDecl = "double GetChiDist(double fX, double fDF);\n";
366std::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
376"double GetChiSqDistCDF(double fX, double fDF);\n";
377std::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
387"double GetChiSqDistPDF(double fX, double fDF);\n";
388std::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*M_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
426"static double lcl_IterateInverseBetaInv(double fp, double fAlpha, \n"
427" double fBeta, double fAx, double fBx, bool *rConvError );\n";
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
523"static double lcl_IterateInverseChiInv"
524"(double fp, double fdf, double fAx, double fBx, bool *rConvError);\n";
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
620"static double lcl_IterateInverseChiSQInv( double fp, double fdf, \n"
621" double fAx, double fBx, bool *rConvError );\n";
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
717std::string gaussinvDecl = "double gaussinv(double x);\n";
718std::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
879"static double lcl_GetLogGammaHelper(double fZ);\n";
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";
888"static double lcl_GetGammaHelper(double fZ);\n";
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";
907"static double lcl_getLanczosSum(double fZ);\n";
908std::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
1031" double GetUpRegIGamma( double fA, double fX ) ;\n";
1032std::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
1044"static inline bool lcl_HasChangeOfSign( double u, double w );\n";
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
1051std::string GetTDistDecl=" double GetTDist(double T, double fDF);\n";
1052std::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
1058std::string GetBetaDecl=" double GetBeta(double fAlpha, double fBeta);\n";
1059std::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
1077std::string GetLogBetaDecl=
1078" double GetLogBeta(double fAlpha, double fBeta);\n";
1079std::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
1097"double GetBetaDistPDF(double fX, double fA, double fB);\n";
1098std::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
1163"double lcl_GetBetaHelperContFrac(double fX, double fA, double fB);\n";
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
1200"double lcl_IterateInverse("
1201"double fAx, double fBx, bool* rConvError,double fp,double fDF );\n";
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";
1294std::string phiDecl=
1295"double phi(double x);\n";
1296std::string phi =
1297"double phi(double x)\n"
1298"{\n"
1299" return 0.39894228040143268 * exp(-(x * x) / 2.0);\n"
1300"}\n";
1301std::string taylorDecl =
1302"double taylor(double* pPolynom, uint nMax, double x);\n";
1303std::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"}";
1313std::string gaussDecl = "double gauss(double x);\n";
1314std::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#endif
1365
1366/* vim:set shiftwidth=4 softtabstop=4 expandtab: */
std::string GetBetaDistDecl
std::string gaussDecl
std::string lcl_IterateInverseChiSQInvDecl
std::string GetLogGammaDecl
std::string GetValueDecl
std::string GetTDist
std::string taylorDecl
std::string GetFDist
std::string lcl_HasChangeOfSignDecl
std::string GetGammaDistPDF
std::string GetBetaDistPDF
std::string phi
std::string lcl_GetLogGammaHelper
std::string lcl_IterateInverseChiInv
std::string lcl_getLanczosSumDecl
std::string GetFInvValue
std::string GetUpRegIGammaDecl
std::string GetLogBeta
std::string lcl_getLanczosSum
std::string fHalfMachEpsDecl
std::string gauss
std::string GetValue
std::string GetBinomDistPMFDecl
std::string GetFInvValueDecl
std::string GetGammaSeriesDecl
std::string GetChiDistDecl
std::string GetGammaDist
std::string GetBetaDist
std::string lcl_IterateInverseBetaInvDecl
std::string GetBetaDistPDFDecl
std::string lcl_GetGammaHelper
std::string GetChiSqDistCDFDecl
std::string GetGammaContFraction
std::string gaussinvDecl
std::string fLogDblMaxDecl
std::string lcl_IterateInverse
std::string GetGammaInvValueDecl
std::string GetChiSqDistPDFDecl
std::string fMachEpsDecl
std::string lcl_IterateInverseChiSQInv
std::string GetLowRegIGammaDecl
std::string lcl_IterateInverseDecl
std::string lcl_GetGammaHelperDecl
std::string GetBetaDecl
std::string GetLowRegIGamma
std::string F_PIDecl
std::string GetUpRegIGamma
std::string GetChiSqDistPDF
std::string GetFDistDecl
std::string MinDecl
std::string GetGammaSeries
std::string lcl_GetBetaHelperContFracDecl
std::string GetBinomDistPMF
std::string lcl_GetBinomDistRangeDecl
std::string lcl_HasChangeOfSign
std::string GetChiDist
std::string phiDecl
std::string taylor
std::string lcl_GetBetaHelperContFrac
std::string GetGammaDistDecl
std::string lcl_GetLogGammaHelperDecl
std::string fBigInvDecl
std::string GetLogGamma
std::string GetChiSqDistCDF
std::string GetGammaDistPDFDecl
std::string GetLogBetaDecl
std::string fMaxGammaArgumentDecl
std::string GetBeta
std::string GetGammaInvValue
std::string lcl_GetBinomDistRange
std::string GetGammaContFractionDecl
std::string lcl_IterateInverseBetaInv
std::string lcl_IterateInverseChiInvDecl
std::string GetTDistDecl
std::string gaussinv