diff options
Diffstat (limited to 'sc/source/core/opencl/opinlinefun_statistical.cxx')
| -rw-r--r-- | sc/source/core/opencl/opinlinefun_statistical.cxx | 1366 |
1 files changed, 1366 insertions, 0 deletions
diff --git a/sc/source/core/opencl/opinlinefun_statistical.cxx b/sc/source/core/opencl/opinlinefun_statistical.cxx new file mode 100644 index 000000000..f981345df --- /dev/null +++ b/sc/source/core/opencl/opinlinefun_statistical.cxx @@ -0,0 +1,1366 @@ +/* -*- Mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */ +/* + * This file is part of the LibreOffice project. + * + * This Source Code Form is subject to the terms of the Mozilla Public + * License, v. 2.0. If a copy of the MPL was not distributed with this + * file, You can obtain one at http://mozilla.org/MPL/2.0/. + */ + +#ifndef SC_OPENCL_OPINLINFUN_statistical +#define SC_OPENCL_OPINLINFUN_statistical +std::string MinDecl = "#define Min 2.22507e-308\n"; +std::string F_PIDecl="#define M_PI 3.1415926535897932384626433832795\n"; +std::string fBigInvDecl ="#define fBigInv 2.22045e-016\n"; +std::string fMachEpsDecl ="#define fMachEps 2.22045e-016\n"; +std::string fLogDblMaxDecl ="#define fLogDblMax log(1.79769e+308)\n"; +std::string fHalfMachEpsDecl ="#define fHalfMachEps 0.5*2.22045e-016\n"; +std::string fMaxGammaArgumentDecl= +"#define fMaxGammaArgument 171.624376956302\n"; +std::string GetValueDecl= +"double GetValue( double x,double fp,double fDF );\n"; +std::string GetValue= +"double GetValue( double x,double fp,double fDF )\n" +"{\n" +" return fp - 2 * GetTDist(x, fDF);\n" +"}\n"; +std::string GetGammaSeriesDecl= +"double GetGammaSeries( double fA, double fX );\n"; +std::string GetGammaSeries = +"double GetGammaSeries( double fA, double fX )\n" +"{\n" +" double fDenomfactor = fA;\n" +" double fSummand = 1.0/fA;\n" +" double fSum = fSummand;\n" +" int nCount=1;\n" +" do\n" +" {\n" +" fDenomfactor = fDenomfactor + 1.0;\n" +" fSummand = fSummand * fX/fDenomfactor;\n" +" fSum = fSum + fSummand;\n" +" nCount = nCount+1;\n" +" } while ( fSummand/fSum > fHalfMachEps && nCount<=10000);\n" +" if (nCount>10000)\n" +" {\n" +" }\n" +" return fSum;\n" +"}\n"; +std::string GetGammaContFractionDecl = "double GetGammaContFraction( double " +"fA, double fX );\n"; +std::string GetGammaContFraction = +"double GetGammaContFraction( double fA, double fX )\n" +"{\n" +" double fBig = 1.0/fBigInv;\n" +" double fCount = 0.0;\n" +" double fNum = 0.0;\n" +" double fY = 1.0 - fA;\n" +" double fDenom = fX + 2.0-fA;\n" +" double fPk = 0.0;\n" +" double fPkm1 = fX + 1.0;\n" +" double fPkm2 = 1.0;\n" +" double fQk = 1.0;\n" +" double fQkm1 = fDenom * fX;\n" +" double fQkm2 = fX;\n" +" double fApprox = fPkm1/fQkm1;\n" +" bool bFinished = false;\n" +" double fR = 0.0;\n" +" do\n" +" {\n" +" fCount = fCount +1.0;\n" +" fY = fY+ 1.0;\n" +" fNum = fY * fCount;\n" +" fDenom = fDenom +2.0;\n" +" fPk = fPkm1 * fDenom - fPkm2 * fNum;\n" +" fQk = fQkm1 * fDenom - fQkm2 * fNum;\n" +" if (fQk != 0.0)\n" +" {\n" +" fR = fPk/fQk;\n" +" bFinished = (fabs( (fApprox - fR)/fR ) <= fHalfMachEps);\n" +" fApprox = fR;\n" +" }\n" +" fPkm2 = fPkm1;\n" +" fPkm1 = fPk;\n" +" fQkm2 = fQkm1;\n" +" fQkm1 = fQk;\n" +" if (fabs(fPk) > fBig)\n" +" {\n" +" fPkm2 = fPkm2 * fBigInv;\n" +" fPkm1 = fPkm1 * fBigInv;\n" +" fQkm2 = fQkm2 * fBigInv;\n" +" fQkm1 = fQkm1 * fBigInv;\n" +" }\n" +" } while (!bFinished && fCount<10000);\n" +" if (!bFinished)\n" +" {\n" +" }\n" +" return fApprox;\n" +"}\n"; +std::string GetLowRegIGammaDecl = "double GetLowRegIGamma( double " +"fA, double fX );\n"; +std::string GetLowRegIGamma = +"double GetLowRegIGamma( double fA, double fX )\n" +"{\n" +" double fLnFactor = fA * log(fX) - fX - lgamma(fA);\n" +" double fFactor = exp(fLnFactor);\n" +" if (fX>fA+1.0) \n" +" return 1.0 - fFactor * GetGammaContFraction(fA,fX);\n" +" else\n" +" return fFactor * GetGammaSeries(fA,fX);\n" +"}\n"; +std::string GetGammaDistDecl = "double GetGammaDist( double fX, double " +"fAlpha, double fLambda );\n"; +std::string GetGammaDist = +"double GetGammaDist( double fX, double fAlpha, double fLambda )\n" +"{\n" +" if (fX <= 0.0)\n" +" return 0.0;\n" +" else\n" +" return GetLowRegIGamma( fAlpha, fX / fLambda);\n" +"}\n"; +std::string GetGammaDistPDFDecl = "double GetGammaDistPDF( double fX" +", double fAlpha, double fLambda );\n"; +std::string GetGammaDistPDF = +"double GetGammaDistPDF( double fX, double fAlpha, double fLambda )\n" +"{\n" +" if (fX < 0.0)\n" +" return 0.0;\n" +" else if (fX == 0)\n" +" {\n" +" if (fAlpha < 1.0)\n" +" {\n" +" return HUGE_VAL;\n" +" }\n" +" else if (fAlpha == 1)\n" +" {\n" +" return (1.0 / fLambda);\n" +" }\n" +" else\n" +" {\n" +" return 0.0;\n" +" }\n" +" }\n" +" else\n" +" {\n" +" double fXr = fX / fLambda;\n" +" if (fXr > 1.0)\n" +" {\n" +" if (log(fXr) * (fAlpha-1.0) < fLogDblMax &&" +"fAlpha < fMaxGammaArgument)\n" +" {\n" +" return pow( fXr, fAlpha-1.0) * exp(-fXr) / " +"fLambda / tgamma(fAlpha);\n" +" }\n" +" else\n" +" {\n" +" return exp( (fAlpha-1.0) * log(fXr) - fXr - " +"log(fLambda) - lgamma(fAlpha));\n" +" }\n" +" }\n" +" else\n" +" {\n" +" if (fAlpha<fMaxGammaArgument)\n" +" {\n" +" return pow( fXr, fAlpha-1.0) * exp(-fXr) / " +"fLambda / tgamma(fAlpha);\n" +" }\n" +" else\n" +" {\n" +" return pow( fXr, fAlpha-1.0) * exp(-fXr) / " +"fLambda / exp( lgamma(fAlpha));\n" +" }\n" +" }\n" +" }\n" +"}\n"; +std::string GetBetaDistDecl = +"double GetBetaDist(double fXin, double fAlpha, double fBeta);\n"; +std::string GetBetaDist = +"double GetBetaDist(double fXin, double fAlpha, double fBeta)\n" +"{\n" +" if (fXin <= 0.0)\n" +" return 0.0;\n" +" if (fXin >= 1.0)\n" +" return 1.0;\n" +" if (fBeta == 1.0)\n" +" return pow(fXin, fAlpha);\n" +" if (fAlpha == 1.0)\n" +" return -expm1(fBeta * log1p(-fXin));\n" +" double fResult;\n" +" double fY = (0.5-fXin)+0.5;\n" +" double flnY = log1p(-fXin);\n" +" double fX = fXin;\n" +" double flnX = log(fXin);\n" +" double fA = fAlpha;\n" +" double fB = fBeta;\n" +" bool bReflect = fXin > fAlpha*pow((fAlpha+fBeta),-1.0);\n" +" if (bReflect)\n" +" {\n" +" fA = fBeta;\n" +" fB = fAlpha;\n" +" fX = fY;\n" +" fY = fXin;\n" +" flnX = flnY;\n" +" flnY = log(fXin);\n" +" }\n" +" fResult = lcl_GetBetaHelperContFrac(fX,fA,fB)*pow(fA,-1.0);\n" +" double fP = fA*pow((fA+fB),-1.0);\n" +" double fQ = fB*pow((fA+fB),-1.0);\n" +" if (fA > 1.0 && fB > 1.0 && fP < 0.97 && fQ < 0.97)\n" +" fResult *= GetBetaDistPDF(fX,fA,fB)*fX*fY;\n" +" else\n" +" fResult *= pow(exp(1.0),(fA*flnX + fB*flnY - GetLogBeta(fA,fB)));\n" +" if (bReflect)\n" +" fResult = 0.5 - fResult + 0.5;\n" +" if (fResult > 1.0)\n" +" fResult = 1.0;\n" +" if (fResult < 0.0)\n" +" fResult = 0.0;\n" +" return fResult;\n" +"}\n"; + +std::string GetFDistDecl = + "double GetFDist(double x, double fF1, double fF2);\n"; +std::string GetFDist = +"double GetFDist(double x, double fF1, double fF2)\n" +"{\n" +" double arg = fF2*pow((fF2+fF1*x),-1.0);\n" +" double alpha = fF2*pow(2.0,-1.0);\n" +" double beta = fF1*pow(2.0,-1.0);\n" +" return (GetBetaDist(arg, alpha, beta));\n" +"}\n"; +std::string GetGammaInvValueDecl = "double" +" GetGammaInvValue(double fAlpha,double fBeta,double fX1 );\n"; +std::string GetGammaInvValue = +"double GetGammaInvValue(double fAlpha,double fBeta,double fX1 )\n" +"{\n" +" if (fX1 <= 0.0)\n" +" return 0.0;\n" +" else\n" +" {\n" +" double fX=fX1*pow(fBeta,-1.0);\n" +" double fLnFactor = fAlpha * log(fX) - fX - lgamma(fAlpha);\n" +" double fFactor = exp(fLnFactor);\n" +" if (fX>fAlpha+1.0)\n" +" return 1.0 - fFactor * GetGammaContFraction(fAlpha,fX);\n" +" else\n" +" return fFactor * GetGammaSeries(fAlpha,fX);\n" +" }\n" +"}\n"; +std::string GetFInvValueDecl = "double GetFInvValue(double fF1,double fF2" +",double fX1 );"; +std::string GetFInvValue = +"double GetFInvValue(double fF1,double fF2,double fX1 )\n" +"{\n" +" double arg = fF2*pow((fF2+fF1*fX1),-1.0);\n" +" double alpha = fF2*pow(2.0,-1.0);\n" +" double beta = fF1*pow(2.0,-1.0);\n" +" double fXin,fAlpha,fBeta;\n" +" fXin=arg;fAlpha=alpha;fBeta=beta;\n" +" if (fXin <= 0.0)\n" +" return 0.0;\n" +" if (fXin >= 1.0)\n" +" return 1.0;\n" +" if (fBeta == 1.0)\n" +" return pow(fXin, fAlpha);\n" +" if (fAlpha == 1.0)\n" +" return -expm1(fBeta * log1p(-fXin));\n" +" double fResult;\n" +" double fY = (0.5-fXin)+0.5;\n" +" double flnY = log1p(-fXin);\n" +" double fX = fXin;\n" +" double flnX = log(fXin);\n" +" double fA = fAlpha;\n" +" double fB = fBeta;\n" +" bool bReflect = fXin > fAlpha*pow((fAlpha+fBeta),-1.0);\n" +" if (bReflect)\n" +" {\n" +" fA = fBeta;\n" +" fB = fAlpha;\n" +" fX = fY;\n" +" fY = fXin;\n" +" flnX = flnY;\n" +" flnY = log(fXin);\n" +" }\n" +" fResult = lcl_GetBetaHelperContFrac(fX,fA,fB);\n" +" fResult = fResult*pow(fA,-1.0);\n" +" double fP = fA*pow((fA+fB),-1.0);\n" +" double fQ = fB*pow((fA+fB),-1.0);\n" +" double fTemp;\n" +" if (fA > 1.0 && fB > 1.0 && fP < 0.97 && fQ < 0.97)\n" +" fTemp = GetBetaDistPDF(fX,fA,fB)*fX*fY;\n" +" else\n" +" fTemp = exp(fA*flnX + fB*flnY - GetLogBeta(fA,fB));\n" +" fResult *= fTemp;\n" +" if (bReflect)\n" +" fResult = 0.5 - fResult + 0.5;\n" +" if (fResult > 1.0)\n" +" fResult = 1.0;\n" +" if (fResult < 0.0)\n" +" fResult = 0.0;\n" +" return fResult;\n" +"}\n"; +std::string GetBinomDistPMFDecl = + "double GetBinomDistPMF(double x, double n, double p);"; +std::string GetBinomDistPMF = +"double GetBinomDistPMF(double x, double n, double p)\n" +"{\n" +" double q = (0.5 - p) + 0.5;\n" +" double fFactor = pow(q, n);\n" +" if (fFactor <= Min)\n" +" {\n" +" fFactor = pow(p, n);\n" +" if (fFactor <= Min)\n" +" return GetBetaDistPDF(p, x + 1.0, n - x + 1.0)*pow((n + 1.0),-1.0);\n" +" else\n" +" {\n" +" uint max = (uint)(n - x);\n" +" for (uint i = 0; i < max && fFactor > 0.0; ++i)\n" +" fFactor *= (n - i)*pow((i + 1),-1.0)*q*pow(p,-1.0);\n" +" return fFactor;\n" +" }\n" +" }\n" +" else\n" +" {\n" +" uint max = (uint)x;\n" +" for (uint i = 0; i < max && fFactor > 0.0; ++i)\n" +" fFactor *= (n - i)*pow((i + 1),-1.0)*p*pow(q,-1.0);\n" +" return fFactor;\n" +" }\n" +"}\n"; + +std::string lcl_GetBinomDistRangeDecl = + "double lcl_GetBinomDistRange(double n, \n" +"double xs, double xe, double fFactor, double p, double q);"; +std::string lcl_GetBinomDistRange= +"double lcl_GetBinomDistRange(double n, double xs, double xe,\n" +" double fFactor, double p, double q)\n" +"{\n" +" uint i;\n" +" double fSum;\n" +" uint nXs = (uint)xs;\n" +" for (i = 1; i <= nXs && fFactor > 0.0; ++i)\n" +" fFactor *= (n - i + 1)*pow(i,-1.0)*p*pow(q,-1.0);\n" +" fSum = fFactor;\n" +" uint nXe =(uint)xe;\n" +" for (i = nXs + 1; i <= nXe && fFactor > 0.0; ++i)\n" +" {\n" +" fFactor *= (n - i + 1)*pow(i,-1.0)*p*pow(q,-1.0);\n" +" fSum += fFactor;\n" +" }\n" +" return (fSum > 1.0) ? 1.0 : fSum;\n" +"}\n"; + +std::string GetLogGammaDecl = "double GetLogGamma(double fZ);\n"; +std::string GetLogGamma = +"double GetLogGamma(double fZ)\n" +"{\n" +" if (fZ >= fMaxGammaArgument)\n" +" return lcl_GetLogGammaHelper(fZ);\n" +" if (fZ >= 1.0)\n" +" return log(lcl_GetGammaHelper(fZ));\n" +" if (fZ >= 0.5)\n" +" return log( lcl_GetGammaHelper(fZ+1) *pow(fZ,-1.0));\n" +" return lcl_GetLogGammaHelper(fZ+2) - log(fZ+1) - log(fZ);\n" +"}\n"; + +std::string GetChiDistDecl = "double GetChiDist(double fX, double fDF);\n"; +std::string GetChiDist = +"double GetChiDist(double fX, double fDF)\n" +"{\n" +" if (fX <= 0.0)\n" +" return 1.0;\n" +" else\n" +" return GetUpRegIGamma( fDF*pow(2.0,-1.0), fX*pow(2.0,-1.0));\n" +"}\n"; + +std::string GetChiSqDistCDFDecl = +"double GetChiSqDistCDF(double fX, double fDF);\n"; +std::string GetChiSqDistCDF = +"double GetChiSqDistCDF(double fX, double fDF)\n" +"{\n" +" if (fX <= 0.0)\n" +" return 0.0;" +" else\n" +" return GetLowRegIGamma( fDF*pow(2.0,-1.0), fX*pow(2.0,-1.0));\n" +"}\n"; + +std::string GetChiSqDistPDFDecl= +"double GetChiSqDistPDF(double fX, double fDF);\n"; +std::string GetChiSqDistPDF = +"double GetChiSqDistPDF(double fX, double fDF)\n" +"{\n" +" double fValue;\n" +" if (fX <= 0.0)\n" +" return 0.0;\n" +" if (fDF*fX > 1391000.0)\n" +" {\n" +" fValue = exp((0.5*fDF - 1) * log(fX*0.5) - 0.5 * fX - log(2.0)" +" - lgamma(0.5*fDF));\n" +" }\n" +" else\n" +" {\n" +" double fCount;\n" +" if (fmod(fDF,2.0)<0.5)\n" +" {\n" +" fValue = 0.5;\n" +" fCount = 2.0;\n" +" }\n" +" else\n" +" {\n" +" fValue = pow(sqrt(fX*2*M_PI),-1.0);\n" +" fCount = 1.0;\n" +" }\n" +" while ( fCount < fDF)\n" +" {\n" +" fValue *= (fX *pow(fCount,-1.0));\n" +" fCount += 2.0;\n" +" }\n" +" if (fX>=1425.0)\n" +" fValue = exp(log(fValue)-fX*pow(2,-1.0));\n" +" else\n" +" fValue *= exp(-fX*pow(2,-1.0));\n" +" }\n" +" return fValue;\n" +"}\n"; + +std::string lcl_IterateInverseBetaInvDecl = +"static double lcl_IterateInverseBetaInv(double fp, double fAlpha, \n" +" double fBeta, double fAx, double fBx, bool *rConvError );\n"; +std::string lcl_IterateInverseBetaInv = +"static double lcl_IterateInverseBetaInv(double fp, double fAlpha, \n" +" double fBeta, double fAx, double fBx, bool *rConvError )\n" +"{\n" +" *rConvError = false;\n" +" double fYEps = 1.0E-307;\n" +" double fXEps = fMachEps;\n" +" if(!(fAx < fBx))\n" +" {\n" +" //print error\n" +" }\n" +" double fAy = fp - GetBetaDist(fAx, fAlpha, fBeta);\n" +" double fBy = fp - GetBetaDist(fBx, fAlpha, fBeta);\n" +" double fTemp;\n" +" unsigned short nCount;\n" +" for (nCount = 0; nCount < 1000 && !lcl_HasChangeOfSign(fAy,fBy);" +" nCount++)\n" +" {\n" +" if (fabs(fAy) <= fabs(fBy))\n" +" {\n" +" fTemp = fAx;\n" +" fAx += 2.0 * (fAx - fBx);\n" +" if (fAx < 0.0)\n" +" fAx = 0.0;\n" +" fBx = fTemp;\n" +" fBy = fAy;\n" +" fAy = fp - GetBetaDist(fAx, fAlpha, fBeta);\n" +" }\n" +" else\n" +" {\n" +" fTemp = fBx;\n" +" fBx += 2.0 * (fBx - fAx);\n" +" fAx = fTemp;\n" +" fAy = fBy;\n" +" fBy = fp - GetBetaDist(fBx, fAlpha, fBeta);\n" +" }\n" +" }\n" +" if (fAy == 0.0)\n" +" return fAx;\n" +" if (fBy == 0.0)\n" +" return fBx;\n" +" if (!lcl_HasChangeOfSign( fAy, fBy))\n" +" {\n" +" *rConvError = true;\n" +" return 0.0;\n" +" }\n" +" double fPx = fAx;\n" +" double fPy = fAy;\n" +" double fQx = fBx;\n" +" double fQy = fBy;\n" +" double fRx = fAx;\n" +" double fRy = fAy;\n" +" double fSx = 0.5 * (fAx + fBx);\n" +" bool bHasToInterpolate = true;\n" +" nCount = 0;\n" +" while ( nCount < 500 && fabs(fRy) > fYEps &&\n" +" (fBx-fAx) > fmax( fabs(fAx), fabs(fBx)) * fXEps )\n" +" {\n" +" if (bHasToInterpolate)\n" +" {\n" +" if (fPy!=fQy && fQy!=fRy && fRy!=fPy)\n" +" {\n" +" fSx = fPx*fRy*fQy*pow(fRy-fPy,-1.0)*pow(fQy-fPy,-1.0)\n" +" + fRx*fQy*fPy*pow(fQy-fRy,-1.0)*pow(fPy-fRy,-1.0)\n" +" + fQx*fPy*fRy*pow(fPy-fQy,-1.0)*pow(fRy-fQy,-1.0);\n" +" bHasToInterpolate = (fAx < fSx) && (fSx < fBx);\n" +" }\n" +" else\n" +" bHasToInterpolate = false;\n" +" }\n" +" if(!bHasToInterpolate)\n" +" {\n" +" fSx = 0.5 * (fAx + fBx);\n" +" fPx = fAx; fPy = fAy;\n" +" fQx = fBx; fQy = fBy;\n" +" bHasToInterpolate = true;\n" +" }\n" +" fPx = fQx; fQx = fRx; fRx = fSx;\n" +" fPy = fQy; fQy = fRy; fRy = fp - GetBetaDist(fSx, fAlpha, fBeta);\n" +" if (lcl_HasChangeOfSign( fAy, fRy))\n" +" {\n" +" fBx = fRx; fBy = fRy;\n" +" }\n" +" else\n" +" {\n" +" fAx = fRx; fAy = fRy;\n" +" }\n" +" bHasToInterpolate = bHasToInterpolate && (fabs(fRy) *" +" 2.0 <= fabs(fQy));\n" +" ++nCount;\n" +" }\n" +" return fRx;\n" +"}\n"; + +std::string lcl_IterateInverseChiInvDecl = +"static double lcl_IterateInverseChiInv" +"(double fp, double fdf, double fAx, double fBx, bool *rConvError);\n"; +std::string lcl_IterateInverseChiInv = +"static double lcl_IterateInverseChiInv" +"(double fp, double fdf, double fAx, double fBx, bool *rConvError)\n" +"{\n" +" *rConvError = false;\n" +" double fYEps = 1.0E-307;\n" +" double fXEps = fMachEps;\n" +" if(!(fAx < fBx))\n" +" {\n" +" //print error\n" +" }" +" double fAy = fp - GetChiDist(fAx, fdf);\n" +" double fBy = fp - GetChiDist(fBx, fdf);\n" +" double fTemp;\n" +" unsigned short nCount;\n" +" for (nCount = 0; nCount < 1000 && " +"!lcl_HasChangeOfSign(fAy,fBy); nCount++)\n" +" {\n" +" if (fabs(fAy) <= fabs(fBy))\n" +" {\n" +" fTemp = fAx;\n" +" fAx += 2.0 * (fAx - fBx);\n" +" if (fAx < 0.0)\n" +" fAx = 0.0;\n" +" fBx = fTemp;\n" +" fBy = fAy;\n" +" fAy = fp - GetChiDist(fAx, fdf);\n" +" }\n" +" else\n" +" {\n" +" fTemp = fBx;\n" +" fBx += 2.0 * (fBx - fAx);\n" +" fAx = fTemp;\n" +" fAy = fBy;\n" +" fBy = fp - GetChiDist(fBx, fdf);\n" +" }\n" +" }\n" +" if (fAy == 0.0)\n" +" return fAx;\n" +" if (fBy == 0.0)\n" +" return fBx;\n" +" if (!lcl_HasChangeOfSign( fAy, fBy))\n" +" {\n" +" *rConvError = true;\n" +" return 0.0;\n" +" }\n" +" double fPx = fAx;\n" +" double fPy = fAy;\n" +" double fQx = fBx;\n" +" double fQy = fBy;\n" +" double fRx = fAx;\n" +" double fRy = fAy;\n" +" double fSx = 0.5 * (fAx + fBx);\n" +" bool bHasToInterpolate = true;\n" +" nCount = 0;\n" +" while ( nCount < 500 && fabs(fRy) > fYEps &&\n" +" (fBx-fAx) > fmax( fabs(fAx), fabs(fBx)) * fXEps )\n" +" {\n" +" if (bHasToInterpolate)\n" +" {\n" +" if (fPy!=fQy && fQy!=fRy && fRy!=fPy)\n" +" {\n" +" fSx = fPx * fRy * fQy*pow(fRy-fPy,-1.0)*pow(fQy-fPy,-1.0)\n" +" + fRx * fQy * fPy*pow(fQy-fRy,-1.0)*pow(fPy-fRy,-1.0)\n" +" + fQx * fPy * fRy*pow(fPy-fQy,-1.0)*pow(fRy-fQy,-1.0);\n" +" bHasToInterpolate = (fAx < fSx) && (fSx < fBx);\n" +" }\n" +" else\n" +" bHasToInterpolate = false;\n" +" }\n" +" if(!bHasToInterpolate)\n" +" {\n" +" fSx = 0.5 * (fAx + fBx);\n" +" fPx = fAx; fPy = fAy;\n" +" fQx = fBx; fQy = fBy;\n" +" bHasToInterpolate = true;\n" +" }\n" +" fPx = fQx; fQx = fRx; fRx = fSx;\n" +" fPy = fQy; fQy = fRy; fRy = fp - GetChiDist(fSx, fdf);\n" +" if (lcl_HasChangeOfSign( fAy, fRy))\n" +" {\n" +" fBx = fRx; fBy = fRy;\n" +" }\n" +" else\n" +" {\n" +" fAx = fRx; fAy = fRy;\n" +" }\n" +" bHasToInterpolate = bHasToInterpolate && (fabs(fRy)" +" * 2.0 <= fabs(fQy));\n" +" ++nCount;\n" +" }\n" +" return fRx;\n" +"}\n"; + +std::string lcl_IterateInverseChiSQInvDecl = +"static double lcl_IterateInverseChiSQInv( double fp, double fdf, \n" +" double fAx, double fBx, bool *rConvError );\n"; +std::string lcl_IterateInverseChiSQInv = +"static double lcl_IterateInverseChiSQInv( double fp, double fdf, \n" +" double fAx, double fBx, bool *rConvError )\n" +"{\n" +" *rConvError = false;\n" +" double fYEps = 1.0E-307;\n" +" double fXEps = fMachEps;\n" + +" if(!(fAx < fBx))\n" +" {\n" +" //print error\n" +" }\n" +" double fAy = fp - GetChiSqDistCDF(fAx, fdf);\n" +" double fBy = fp - GetChiSqDistCDF(fBx, fdf);\n" +" double fTemp;\n" +" unsigned short nCount;\n" +" for (nCount = 0; nCount < 1000 && !lcl_HasChangeOfSign(fAy,fBy);" +" nCount++)\n" +" {\n" +" if (fabs(fAy) <= fabs(fBy))\n" +" {\n" +" fTemp = fAx;\n" +" fAx += 2.0 * (fAx - fBx);\n" +" if (fAx < 0.0)\n" +" fAx = 0.0;\n" +" fBx = fTemp;\n" +" fBy = fAy;\n" +" fAy = fp - GetChiSqDistCDF(fAx, fdf);\n" +" }\n" +" else\n" +" {\n" +" fTemp = fBx;\n" +" fBx += 2.0 * (fBx - fAx);\n" +" fAx = fTemp;\n" +" fAy = fBy;\n" +" fBy = fp - GetChiSqDistCDF(fBx, fdf);\n" +" }\n" +" }\n" +" if (fAy == 0.0)\n" +" return fAx;\n" +" if (fBy == 0.0)\n" +" return fBx;\n" +" if (!lcl_HasChangeOfSign( fAy, fBy))\n" +" {\n" +" *rConvError = true;\n" +" return 0.0;\n" +" }\n" +" double fPx = fAx;\n" +" double fPy = fAy;\n" +" double fQx = fBx;\n" +" double fQy = fBy;\n" +" double fRx = fAx;\n" +" double fRy = fAy;\n" +" double fSx = 0.5 * (fAx + fBx);\n" +" bool bHasToInterpolate = true;\n" +" nCount = 0;\n" +" while ( nCount < 500 && fabs(fRy) > fYEps &&\n" +" (fBx-fAx) > fmax( fabs(fAx), fabs(fBx)) * fXEps )\n" +" {\n" +" if (bHasToInterpolate)\n" +" {\n" +" if (fPy!=fQy && fQy!=fRy && fRy!=fPy)\n" +" {\n" +" fSx = fPx * fRy * fQy / (fRy-fPy) / (fQy-fPy)\n" +" + fRx * fQy * fPy / (fQy-fRy) / (fPy-fRy)\n" +" + fQx * fPy * fRy / (fPy-fQy) / (fRy-fQy);\n" +" bHasToInterpolate = (fAx < fSx) && (fSx < fBx);\n" +" }\n" +" else\n" +" bHasToInterpolate = false;\n" +" }\n" +" if(!bHasToInterpolate)\n" +" {\n" +" fSx = 0.5 * (fAx + fBx);\n" +" fPx = fAx; fPy = fAy;\n" +" fQx = fBx; fQy = fBy;\n" +" bHasToInterpolate = true;\n" +" }\n" +" fPx = fQx; fQx = fRx; fRx = fSx;\n" +" fPy = fQy; fQy = fRy; fRy = fp - GetChiSqDistCDF(fSx, fdf);\n" +" if (lcl_HasChangeOfSign( fAy, fRy))\n" +" {\n" +" fBx = fRx; fBy = fRy;\n" +" }\n" +" else\n" +" {\n" +" fAx = fRx; fAy = fRy;\n" +" }\n" +" bHasToInterpolate = bHasToInterpolate && (fabs(fRy) * 2.0" +" <= fabs(fQy));\n" +" ++nCount;\n" +" }\n" +" return fRx;\n" +"}\n"; + +std::string gaussinvDecl = "double gaussinv(double x);\n"; +std::string gaussinv = +"double gaussinv(double x)\n" +"{\n" +" double q,t,z;\n" +" q=x-0.5;\n" +" if(fabs(q)<=.425)\n" +" {\n" +" t=0.180625-q*q;\n" +" z=\n" +" q*\n" +" (\n" +" (\n" +" (\n" +" (\n" +" (\n" +" (\n" +" (\n" +" t*2509.0809287301226727+33430.575583588128105\n" +" )\n" +" *t+67265.770927008700853\n" +" )\n" +" *t+45921.953931549871457\n" +" )\n" +" *t+13731.693765509461125\n" +" )\n" +" *t+1971.5909503065514427\n" +" )\n" +" *t+133.14166789178437745\n" +" )\n" +" *t+3.387132872796366608\n" +" )\n" +" *pow\n" +" (\n" +" (\n" +" (\n" +" (\n" +" (\n" +" (\n" +" (\n" +" t*5226.495278852854561+28729.085735721942674\n" +" )\n" +" *t+39307.89580009271061\n" +" )\n" +" *t+21213.794301586595867\n" +" )\n" +" *t+5394.1960214247511077\n" +" )\n" +" *t+687.1870074920579083\n" +" )\n" +" *t+42.313330701600911252\n" +" )\n" +" *t+1.0\n" +" , -1.0);\n" +" }\n" +" else\n" +" {\n" +" if(q>0) t=1-x;\n" +" else t=x;\n" +" t=sqrt(-log(t));\n" +" if(t<=5.0)\n" +" {\n" +" t+=-1.6;\n" +" z=\n" +" (\n" +" (\n" +" (\n" +" (\n" +" (\n" +" (\n" +" (\n" +" t*7.7454501427834140764e-4+0.0227238449892691845833\n" +" )\n" +" *t+0.24178072517745061177\n" +" )\n" +" *t+1.27045825245236838258\n" +" )\n" +" *t+3.64784832476320460504\n" +" )\n" +" *t+5.7694972214606914055\n" +" )\n" +" *t+4.6303378461565452959\n" +" )\n" +" *t+1.42343711074968357734\n" +" )\n" +" *pow\n" +" (\n" +" (\n" +" (\n" +" (\n" +" (\n" +" (\n" +" (\n" +" t*1.05075007164441684324e-9+5.475938084995344946e-4\n" +" )\n" +" *t+0.0151986665636164571966\n" +" )\n" +" *t+0.14810397642748007459\n" +" )\n" +" *t+0.68976733498510000455\n" +" )\n" +" *t+1.6763848301838038494\n" +" )\n" +" *t+2.05319162663775882187\n" +" )\n" +" *t+1.0\n" +" , -1.0);\n" +" }\n" +" else\n" +" {\n" +" t+=-5.0;\n" +" z=\n" +" (\n" +" (\n" +" (\n" +" (\n" +" (\n" +" (\n" +" (\n" +" t*2.01033439929228813265e-7+2.71155556874348757815e-5\n" +" )\n" +" *t+0.0012426609473880784386\n" +" )\n" +" *t+0.026532189526576123093\n" +" )\n" +" *t+0.29656057182850489123\n" +" )\n" +" *t+1.7848265399172913358\n" +" )\n" +" *t+5.4637849111641143699\n" +" )\n" +" *t+6.6579046435011037772\n" +" )\n" +" *pow\n" +" (\n" +" (\n" +" (\n" +" (\n" +" (\n" +" (\n" +" (\n" +" t*2.04426310338993978564e-15+1.4215117583164458887e-7\n" +" )\n" +" *t+1.8463183175100546818e-5\n" +" )\n" +" *t+7.868691311456132591e-4\n" +" )\n" +" *t+0.0148753612908506148525\n" +" )\n" +" *t+0.13692988092273580531\n" +" )\n" +" *t+0.59983220655588793769\n" +" )\n" +" *t+1.0\n" +" , -1.0);\n" +" }\n" +" if(q<0.0) z=-z;\n" +" }\n" +" return z;\n" +"}\n"; + +std::string lcl_GetLogGammaHelperDecl= +"static double lcl_GetLogGammaHelper(double fZ);\n"; +std::string lcl_GetLogGammaHelper = +"static double lcl_GetLogGammaHelper(double fZ)\n" +"{\n" +" double fg = 6.024680040776729583740234375;\n" +" double fZgHelp = fZ + fg - 0.5;\n" +" return log( lcl_getLanczosSum(fZ)) + (fZ-0.5) * log(fZgHelp) - fZgHelp;\n" +"}\n"; +std::string lcl_GetGammaHelperDecl= +"static double lcl_GetGammaHelper(double fZ);\n"; +std::string lcl_GetGammaHelper = +"static double lcl_GetGammaHelper(double fZ)\n" +"{\n" +" double fGamma = lcl_getLanczosSum(fZ);\n" +" double fg = 6.024680040776729583740234375;\n" +" double fZgHelp = fZ + fg - 0.5;\n" +" double fHalfpower = pow( fZgHelp, fZ*pow(2,-1.0) - 0.25);\n" +" fGamma *= fHalfpower;\n" +" fGamma = fGamma*pow(exp(fZgHelp),-1.0);\n" +" fGamma *= fHalfpower;\n" +" fGamma = 120.4;\n" +" if (fZ <= 20.0 && fZ == (int)fZ)\n" +" {\n" +" fGamma = (int)(fGamma+0.5);\n" +" }\n" +" return fGamma;\n" +"}\n"; +std::string lcl_getLanczosSumDecl= +"static double lcl_getLanczosSum(double fZ);\n"; +std::string lcl_getLanczosSum = +"static double lcl_getLanczosSum(double fZ) \n" +"{ \n" +" double fNum[13] ={ \n" +" 23531376880.41075968857200767445163675473, \n" +" 42919803642.64909876895789904700198885093, \n" +" 35711959237.35566804944018545154716670596, \n" +" 17921034426.03720969991975575445893111267, \n" +" 6039542586.35202800506429164430729792107, \n" +" 1439720407.311721673663223072794912393972, \n" +" 248874557.8620541565114603864132294232163, \n" +" 31426415.58540019438061423162831820536287, \n" +" 2876370.628935372441225409051620849613599, \n" +" 186056.2653952234950402949897160456992822, \n" +" 8071.672002365816210638002902272250613822, \n" +" 210.8242777515793458725097339207133627117, \n" +" 2.506628274631000270164908177133837338626 \n" +" }; \n" +" double fDenom[13] = { \n" +" 0,\n" +" 39916800,\n" +" 120543840,\n" +" 150917976,\n" +" 105258076,\n" +" 45995730,\n" +" 13339535,\n" +" 2637558,\n" +" 357423,\n" +" 32670,\n" +" 1925,\n" +" 66,\n" +" 1\n" +" };\n" +" double fSumNum;\n" +" double fSumDenom;\n" +" int nI;\n" +" if (fZ<=1.0)\n" +" {\n" +" fSumNum = fNum[12];\n" +" fSumDenom = fDenom[12];\n" +" nI = 11;\n" +" fSumNum = fSumNum*fZ+fNum[nI];\n" +" fSumDenom = fSumDenom*fZ+fDenom[nI];\n" +" nI = 10;\n" +" fSumNum = fSumNum*fZ+fNum[nI];\n" +" fSumDenom = fSumDenom*fZ+fDenom[nI];\n" +" nI = 9;\n" +" fSumNum = fSumNum*fZ+fNum[nI];\n" +" fSumDenom = fSumDenom*fZ+fDenom[nI];\n" +" nI = 8;\n" +" fSumNum = fSumNum*fZ+fNum[nI];\n" +" fSumDenom = fSumDenom*fZ+fDenom[nI];\n" +" nI = 7;\n" +" fSumNum = fSumNum*fZ+fNum[nI];\n" +" fSumDenom = fSumDenom*fZ+fDenom[nI];\n" +" nI = 6;\n" +" fSumNum = fSumNum*fZ+fNum[nI];\n" +" fSumDenom = fSumDenom*fZ+fDenom[nI];\n" +" nI = 5;\n" +" fSumNum = fSumNum*fZ+fNum[nI];\n" +" fSumDenom = fSumDenom*fZ+fDenom[nI];\n" +" nI = 4;\n" +" fSumNum = fSumNum*fZ+fNum[nI];\n" +" fSumDenom = fSumDenom*fZ+fDenom[nI];\n" +" nI = 3;\n" +" fSumNum = fSumNum*fZ+fNum[nI];\n" +" fSumDenom = fSumDenom*fZ+fDenom[nI];\n" +" nI = 2;\n" +" fSumNum = fSumNum*fZ+fNum[nI];\n" +" fSumDenom = fSumDenom*fZ+fDenom[nI];\n" +" nI = 1;\n" +" fSumNum = fSumNum*fZ+fNum[nI];\n" +" fSumDenom = fSumDenom*fZ+fDenom[nI];\n" +" nI = 0;\n" +" fSumNum = fSumNum*fZ+fNum[nI];\n" +" fSumDenom = fSumDenom*fZ+fDenom[nI];\n" +" }\n" +" if (fZ>1.0)\n" +" {\n" +" double fZInv = pow(fZ,-1.0);\n" +" fSumNum = fNum[0];\n" +" fSumDenom = fDenom[0];\n" +" nI = 1;\n" +" fSumNum = fSumNum*fZInv+fNum[nI];\n" +" fSumDenom = fSumDenom*fZInv+fDenom[nI];\n" +" nI = 2;\n" +" fSumNum = fSumNum*fZInv+fNum[nI];\n" +" fSumDenom = fSumDenom*fZInv+fDenom[nI];\n" +" nI = 3;\n" +" fSumNum = fSumNum*fZInv+fNum[nI];\n" +" fSumDenom = fSumDenom*fZInv+fDenom[nI];\n" +" nI = 4;\n" +" fSumNum = fSumNum*fZInv+fNum[nI];\n" +" fSumDenom = fSumDenom*fZInv+fDenom[nI];\n" +" nI = 5;\n" +" fSumNum = fSumNum*fZInv+fNum[nI];\n" +" fSumDenom = fSumDenom*fZInv+fDenom[nI];\n" +" nI = 6;\n" +" fSumNum = fSumNum*fZInv+fNum[nI];\n" +" fSumDenom = fSumDenom*fZInv+fDenom[nI];\n" +" nI = 7;\n" +" fSumNum = fSumNum*fZInv+fNum[nI];\n" +" fSumDenom = fSumDenom*fZInv+fDenom[nI];\n" +" nI = 8;\n" +" fSumNum = fSumNum*fZInv+fNum[nI];\n" +" fSumDenom = fSumDenom*fZInv+fDenom[nI];\n" +" nI = 9;\n" +" fSumNum = fSumNum*fZInv+fNum[nI];\n" +" fSumDenom = fSumDenom*fZInv+fDenom[nI];\n" +" nI = 10;\n" +" fSumNum = fSumNum*fZInv+fNum[nI];\n" +" fSumDenom = fSumDenom*fZInv+fDenom[nI];\n" +" nI = 11;\n" +" fSumNum = fSumNum*fZInv+fNum[nI];\n" +" fSumDenom = fSumDenom*fZInv+fDenom[nI];\n" +" nI = 12;\n" +" fSumNum = fSumNum*fZInv+fNum[nI];\n" +" fSumDenom = fSumDenom*fZInv+fDenom[nI];\n" +" }\n" +" return fSumNum*pow(fSumDenom,-1.0);\n" +"}\n"; + +std::string GetUpRegIGammaDecl= +" double GetUpRegIGamma( double fA, double fX ) ;\n"; +std::string GetUpRegIGamma = +"double GetUpRegIGamma( double fA, double fX )\n" +"{\n" +" double fLnFactor= fA*log(fX)-fX-lgamma(fA);\n" +" double fFactor = exp(fLnFactor); \n" +" if (fX>fA+1.0) \n" +" return fFactor * GetGammaContFraction(fA,fX);\n" +" else \n" +" return 1.0 -fFactor * GetGammaSeries(fA,fX);\n" +"}\n"; + +std::string lcl_HasChangeOfSignDecl= +"static inline bool lcl_HasChangeOfSign( double u, double w );\n"; +std::string lcl_HasChangeOfSign = +"static inline bool lcl_HasChangeOfSign( double u, double w )\n" +"{\n" +" return (u < 0.0 && w > 0.0) || (u > 0.0 && w < 0.0);\n" +"}\n"; + +std::string GetTDistDecl=" double GetTDist(double T, double fDF);\n"; +std::string GetTDist = +"double GetTDist(double T, double fDF)\n" +"{\n" +" return 0.5 * GetBetaDist(fDF*pow(fDF+T*T,-1.0),fDF*pow(2.0,-1.0), 0.5);\n" +"}\n"; + +std::string GetBetaDecl=" double GetBeta(double fAlpha, double fBeta);\n"; +std::string GetBeta = +"double GetBeta(double fAlpha, double fBeta)\n" +"{\n" +" double fA;\n" +" double fB;\n" +" fAlpha>fBeta?(fA = fAlpha,fB = fBeta):(fA = fBeta,fB = fAlpha);\n" +" double fAB = fA + fB;\n" + +" if (fAB < fMaxGammaArgument)\n" +" return tgamma(fA)*pow(tgamma(fAB),-1.0)*tgamma(fB);\n" +" double fgm = 5.524680040776729583740234375;\n" +" double fLanczos = lcl_getLanczosSum(fA)*lcl_getLanczosSum(fB)\n" +" *pow(lcl_getLanczosSum(fAB),-1.0);\n" +" fLanczos *= sqrt(((fAB + fgm)*pow(fA + fgm,-1.0))*pow(fB + fgm,-1.0));\n" +" return fLanczos * pow(exp(1.0),(-fA*log1p(fB*pow(fA + fgm,-1.0)))" +" - fB*log1p(fA*pow(fB + fgm,-1.0)) - fgm);\n" +"}\n"; + +std::string GetLogBetaDecl= +" double GetLogBeta(double fAlpha, double fBeta);\n"; +std::string GetLogBeta = +"double GetLogBeta(double fAlpha, double fBeta)\n" +"{\n" +" double fA;\n" +" double fB;\n" +" fAlpha>fBeta?(fA = fAlpha,fB = fBeta):(fA = fBeta,fB = fAlpha);\n" +" double fgm = 5.524680040776729583740234375;\n" + +" double fLanczos = lcl_getLanczosSum(fA)*lcl_getLanczosSum(fB)*\n" +" pow(lcl_getLanczosSum(fA + fB),-1.0);\n" +" double fResult= -fA *log1p(fB*pow(fA + fgm,-1.0))" +"-fB *log1p(fA*pow(fB + fgm,-1.0))-fgm;\n" +" fResult += log(fLanczos)+0.5*(log(fA + fB + fgm) - log(fA + fgm)\n" +" - log(fB + fgm));\n" +" return fResult;\n" +"}\n"; + +std::string GetBetaDistPDFDecl= +"double GetBetaDistPDF(double fX, double fA, double fB);\n"; +std::string GetBetaDistPDF = +"double GetBetaDistPDF(double fX, double fA, double fB)\n" +"{\n" +" if (fA == 1.0) \n" +" {\n" +" if (fB == 1.0)\n" +" return 1.0;\n" +" if (fB == 2.0)\n" +" return -2.0*fX + 2.0;\n" +" if (fX == 1.0 && fB < 1.0)\n" +" {\n" +" return HUGE_VAL;\n" +" }\n" +" if (fX <= 0.01)\n" +" return fB + fB * expm1((fB-1.0) * log1p(-fX));\n" +" else \n" +" return fB * pow(0.5-fX+0.5,fB-1.0);\n" +" }\n" +" if (fB == 1.0) \n" +" {\n" +" if (fA == 2.0)\n" +" return fA * fX;\n" +" if (fX == 0.0 && fA < 1.0)\n" +" {\n" +" return HUGE_VAL;\n" +" }\n" +" return fA * pow(fX,fA-1);\n" +" }\n" +" if (fX <= 0.0)\n" +" {\n" +" if (fA < 1.0 && fX == 0.0)\n" +" {\n" +" return HUGE_VAL;\n" +" }\n" +" else\n" +" return 0.0;\n" +" }\n" +" if (fX >= 1.0)\n" +" {\n" +" if (fB < 1.0 && fX == 1.0)\n" +" {\n" +" return HUGE_VAL;\n" +" }\n" +" else \n" +" return 0.0;\n" +" }\n" +" double fLogDblMax = log( 1.79769e+308 );\n" +" double fLogDblMin = log( 2.22507e-308 );\n" +" double fLogY = (fX < 0.1) ? log1p(-fX) : log(0.5-fX+0.5);\n" +" double fLogX = log(fX);\n" +" double fAm1LogX = (fA-1.0) * fLogX;\n" +" double fBm1LogY = (fB-1.0) * fLogY;\n" +" double fLogBeta = GetLogBeta(fA,fB);\n" +" if ( fAm1LogX < fLogDblMax && fAm1LogX > fLogDblMin\n" +" && fBm1LogY < fLogDblMax && fBm1LogY > fLogDblMin\n" +" && fLogBeta < fLogDblMax && fLogBeta > fLogDblMin\n" +" && fAm1LogX + fBm1LogY < fLogDblMax && fAm1LogX + fBm1LogY > \n" +" fLogDblMin)\n" +" return pow(fX,fA-1.0)*pow(0.5-fX+0.5,fB-1.0)" +"*pow(GetBeta(fA,fB),-1.0);\n" +" else \n" +" return exp( fAm1LogX + fBm1LogY - fLogBeta);\n" +"}\n"; + +std::string lcl_GetBetaHelperContFracDecl= +"double lcl_GetBetaHelperContFrac(double fX, double fA, double fB);\n"; +std::string lcl_GetBetaHelperContFrac = +"double lcl_GetBetaHelperContFrac(double fX, double fA, double fB)\n" +"{ \n" + +" double a1, b1, a2, b2, fnorm, apl2m, d2m, d2m1, cfnew, cf;\n" +" a1 = 1.0; b1 = 1.0;\n" +" b2 = 1.0 - (fA+fB)*pow(fA+1.0,-1.0)*fX;\n" +" b2==0.0?(a2 = 0.0,fnorm = 1.0,cf = 1.0):\n" +" (a2 = 1.0,fnorm = pow(b2,-1.0),cf = a2*fnorm);\n" +" cfnew = 1.0;\n" +" double rm = 1.0;\n" +" double fMaxIter = 50000.0;\n" +" bool bfinished = false;\n" +" do\n" +" {\n" +" apl2m = fA + 2.0*rm;\n" +" d2m = (rm*(fB-rm))*fX*pow(apl2m*(apl2m-1.0),-1.0);\n" +" d2m1 = -((fA+rm)*(fA+rm+fB))*fX*pow(apl2m*(apl2m+1.0),-1.0);\n" +" a1 = (a2+d2m*a1)*fnorm;\n" +" b1 = (b2+d2m*b1)*fnorm;\n" +" a2 = a1 + d2m1*a2*fnorm;\n" +" b2 = b1 + d2m1*b2*fnorm;\n" +" if (b2 != 0.0) \n" +" {\n" +" fnorm = pow(b2,-1.0);\n" +" cfnew = a2*fnorm;\n" +" bfinished = (fabs(cf-cfnew) < fabs(cf)*fMachEps);\n" +" }\n" +" cf = cfnew;\n" +" rm += 1.0;\n" +" }\n" +" while (rm < fMaxIter && !bfinished);\n" +" return cf;\n" +"}\n"; + +std::string lcl_IterateInverseDecl= +"double lcl_IterateInverse(" +"double fAx, double fBx, bool* rConvError,double fp,double fDF );\n"; +std::string lcl_IterateInverse = +"double lcl_IterateInverse( " +"double fAx, double fBx, bool* rConvError,double fp,double fDF )\n" +"{\n" +" *rConvError = false;\n" +" double fYEps = 1.0E-307;\n" +" double fXEps =DBL_EPSILON;\n" +" if(fAx>fBx)\n" +" return DBL_MAX;\n" +" double fAy = GetValue(fAx,fp,fDF);\n" +" double fBy = GetValue(fBx,fp,fDF);\n" +" double fTemp;\n" +" unsigned short nCount;\n" +" double inter;\n" +" bool sign;\n" +" for (nCount =0;nCount<1000&&!lcl_HasChangeOfSign(fAy,fBy);nCount++)\n" +" {\n" +" inter = 2.0 * (fAx - fBx);\n" +" if (fabs(fAy) <= fabs(fBy)) \n" +" {\n" +" sign = true;\n" +" fTemp = fAx;\n" +" fAx += inter;\n" +" if (fAx < 0.0)\n" +" fAx = 0.0;\n" +" fBx = fTemp;\n" +" fBy = fAy;\n" +" fTemp = fAx;\n" +" }\n" +" else\n" +" {\n" +" sign = false;\n" +" fTemp = fBx;\n" +" fBx -= inter;\n" +" fAx = fTemp;\n" +" fAy = fBy;\n" +" fTemp = fBx;\n" +" }\n" +" fTemp = GetValue(fTemp,fp,fDF);\n" +" sign?(fAy = fTemp):(fBy = fTemp);\n" +" }\n" +" if (fAy == 0.0)\n" +" return fAx;\n" +" if (fBy == 0.0)\n" +" return fBx;\n" +" if (!lcl_HasChangeOfSign( fAy, fBy))\n" +" {\n" +" *rConvError = true;\n" +" return 0.0;\n" +" }\n" +" double fPx = fAx;\n" +" double fPy = fAy;\n" +" double fQx = fBx;\n" +" double fQy = fBy;\n" +" double fRx = fAx;\n" +" double fRy = fAy;\n" +" double fSx = 0.5 * (fAx + fBx); \n" +" bool bHasToInterpolate = true;\n" +" nCount = 0;\n" +" while ( nCount < 500 && fabs(fRy) > fYEps &&\n" +" (fBx-fAx) > max( fabs(fAx), fabs(fBx)) * fXEps )\n" +" {\n" +" if (bHasToInterpolate)\n" +" {\n" +" if (fPy!=fQy && fQy!=fRy && fRy!=fPy)\n" +" {\n" +" fSx = fPx * fRy * fQy * pow(fRy-fPy,-1.0)*pow(fQy-fPy,-1.0)\n" +" + fRx * fQy * fPy * pow(fQy-fRy,-1.0)*pow(fPy-fRy,-1.0)\n" +" + fQx * fPy * fRy * pow(fPy-fQy,-1.0)*pow(fRy-fQy,-1.0);\n" +" bHasToInterpolate = (fAx < fSx) && (fSx < fBx);\n" +" }\n" +" else\n" +" bHasToInterpolate = false;\n" +" }\n" +" if(!bHasToInterpolate)\n" +" {\n" +" fSx = 0.5 * (fAx + fBx);\n" +" \n" +" fPx = fAx; fPy = fAy;\n" +" fQx = fBx; fQy = fBy;\n" +" bHasToInterpolate = true;\n" +" }\n" +" fPx = fQx; fQx = fRx; fRx = fSx;\n" +" fPy = fQy; fQy = fRy; fRy = GetValue(fSx,fp,fDF);\n" +" lcl_HasChangeOfSign( fAy, fRy)?(fBx = fRx,fBy = fRy):\n" +" (fAx = fRx,fAy = fRy);\n" +" bHasToInterpolate =\n" +" bHasToInterpolate && (fabs(fRy) * 2.0 <= fabs(fQy));\n" +" ++nCount;\n" +" }\n" +" return fRx;\n" +"}\n"; +std::string phiDecl= +"double phi(double x);\n"; +std::string phi = +"double phi(double x)\n" +"{\n" +" return 0.39894228040143268 * exp(-(x * x) / 2.0);\n" +"}\n"; +std::string taylorDecl = +"double taylor(double* pPolynom, uint nMax, double x);\n"; +std::string taylor = +"double taylor(double* pPolynom, uint nMax, double x)\n" +"{\n" +" double nVal = pPolynom[nMax];\n" +" for (short i = nMax-1; i >= 0; i--)\n" +" {\n" +" nVal = pPolynom[i] + (nVal * x);\n" +" }\n" +" return nVal;\n" +"}"; +std::string gaussDecl = "double gauss(double x);\n"; +std::string gauss = +"double gauss(double x)\n" +"{\n" +" double xAbs = fabs(x);\n" +" uint xShort = (uint)(floor(xAbs));\n" +" double nVal = 0.0;\n" +" if (xShort == 0)\n" +" {\n" +" double t0[] =\n" +" { 0.39894228040143268, -0.06649038006690545, 0.00997355701003582,\n" +" -0.00118732821548045, 0.00011543468761616, -0.00000944465625950,\n" +" 0.00000066596935163, -0.00000004122667415, 0.00000000227352982,\n" +" 0.00000000011301172, 0.00000000000511243, -0.00000000000021218 };\n" +" nVal = taylor(t0, 11, (xAbs * xAbs)) * xAbs;\n" +" }\n" +" else if ((xShort >= 1) && (xShort <= 2))\n" +" {\n" +" double t2[] =\n" +" { 0.47724986805182079, 0.05399096651318805, -0.05399096651318805,\n" +" 0.02699548325659403, -0.00449924720943234, -0.00224962360471617,\n" +" 0.00134977416282970, -0.00011783742691370, -0.00011515930357476,\n" +" 0.00003704737285544, 0.00000282690796889, -0.00000354513195524,\n" +" 0.00000037669563126, 0.00000019202407921, -0.00000005226908590,\n" +" -0.00000000491799345, 0.00000000366377919, -0.00000000015981997,\n" +" -0.00000000017381238, 0.00000000002624031, 0.00000000000560919,\n" +" -0.00000000000172127, -0.00000000000008634, 0.00000000000007894 };\n" +" nVal = taylor(t2, 23, (xAbs - 2.0));\n" +" }\n" +" else if ((xShort >= 3) && (xShort <= 4))\n" +" {\n" +" double t4[] =\n" +" { 0.49996832875816688, 0.00013383022576489, -0.00026766045152977,\n" +" 0.00033457556441221, -0.00028996548915725, 0.00018178605666397,\n" +" -0.00008252863922168, 0.00002551802519049, -0.00000391665839292,\n" +" -0.00000074018205222, 0.00000064422023359, -0.00000017370155340,\n" +" 0.00000000909595465, 0.00000000944943118, -0.00000000329957075,\n" +" 0.00000000029492075, 0.00000000011874477, -0.00000000004420396,\n" +" 0.00000000000361422, 0.00000000000143638, -0.00000000000045848 };\n" +" nVal = taylor(t4, 20, (xAbs - 4.0));\n" +" }\n" +" else\n" +" {\n" +" double asympt[] = { -1.0, 1.0, -3.0, 15.0, -105.0 };\n" +" nVal = 0.5 + phi(xAbs) * taylor(asympt, 4, 1.0/(xAbs * xAbs))/xAbs;\n" +" }\n" +" if (x < 0.0)\n" +" return -nVal;\n" +" else\n" +" return nVal;\n" +"}\n"; +#endif + +/* vim:set shiftwidth=4 softtabstop=4 expandtab: */ |