/* -*- 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 F_PI 3.1415926\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 );"; 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= 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/(fAlpha+fBeta);\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/fA;\n" " double fP = fA/(fA+fB);\n" " double fQ = fB/(fA+fB);\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 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/(fF2+fF1*x);\n" " double alpha = fF2/2.0;\n" " double beta = fF1/2.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/fBeta;\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/(fF2+fF1*fX1);\n" " double alpha = fF2/2.0;\n" " double beta = fF1/2.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/(fAlpha+fBeta);\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/fA;\n" " double fP = fA/(fA+fB);\n" " double fQ = fB/(fA+fB);\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)/(n + 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)/(i + 1)*q/p;\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)/(i + 1)*p/q;\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)/i*p/q;\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)/i*p/q;\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) / fZ);\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 /2.0, fX/2.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/2.0, fX/2.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 = 1/sqrt(fX*2*F_PI);\n" " fCount = 1.0;\n" " }\n" " while ( fCount < fDF)\n" " {\n" " fValue *= (fX / fCount);\n" " fCount += 2.0;\n" " }\n" " if (fX>=1425.0)\n" " fValue = exp(log(fValue)-fX/2);\n" " else\n" " fValue *= exp(-fX/2);\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 / (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 - 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 / (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 - 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" " /\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" " );\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" " /\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" " );\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" " /\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" " );\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 / 2 - 0.25);\n" " fGamma *= fHalfpower;\n" " fGamma /= exp(fZgHelp);\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" " for (nI = 11; nI >= 0; --nI)\n" " {\n" " fSumNum *= fZ;\n" " fSumNum += fNum[nI];\n" " fSumDenom *= fZ;\n" " fSumDenom += fDenom[nI];\n" " }\n" " }\n" " else\n" " {\n" " double fZInv = 1/fZ;\n" " fSumNum = fNum[0];\n" " fSumDenom = fDenom[0];\n" " for (nI = 1; nI <=12; ++nI)\n" " {\n" " fSumNum *= fZInv;\n" " fSumNum += fNum[nI];\n" " fSumDenom *= fZInv;\n" " fSumDenom += fDenom[nI];\n" " }\n" " }\n" " return fSumNum/fSumDenom;\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/(fDF+T*T), fDF/2.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" " if (fAlpha > fBeta)\n" " {\n" " fA = fAlpha; fB = fBeta;\n" " }\n" " else\n" " {\n" " fA = fBeta; fB = fAlpha;\n" " }\n" " if (fA+fB < fMaxGammaArgument)\n" " return tgamma(fA)/tgamma(fA + fB)*tgamma(fB);\n" " double fg = 6.024680040776729583740234375;\n" " double fgm = fg - 0.5;\n" " double fLanczos = lcl_getLanczosSum(fA);\n" " fLanczos /= lcl_getLanczosSum(fA + fB);\n" " fLanczos *= lcl_getLanczosSum(fB);\n" " double fABgm = fA + fB + fgm;\n" " fLanczos *= sqrt((fABgm/(fA + fgm))/(fB + fgm));\n" " double fTempA = fB/(fA + fgm);\n" " double fTempB = fA/(fB + fgm);\n" " double fResult = exp(-fA*log1p(fTempA) - fB*log1p(fTempB) - fgm);\n" " fResult *= fLanczos;\n" " return fResult;\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" " if (fAlpha > fBeta)\n" " {\n" " fA = fAlpha; fB = fBeta;\n" " }\n" " else\n" " {\n" " fA = fBeta; fB = fAlpha;\n" " }\n" " double fg = 6.024680040776729583740234375;\n" " double fgm = fg - 0.5;\n" " double fLanczos = lcl_getLanczosSum(fA);\n" " fLanczos /= lcl_getLanczosSum(fA + fB);\n" " fLanczos *= lcl_getLanczosSum(fB);\n" " double fLogLanczos = log(fLanczos);\n" " double fABgm = fA + fB + fgm;\n" " fLogLanczos += 0.5*(log(fABgm) - log(fA + fgm) - log(fB + fgm));\n" " double fTempA = fB/(fA + fgm);\n" " double fTempB = fA/(fB + fgm);\n" " double fResult = -fA * log1p(fTempA)\n" " -fB * log1p(fTempB)-fgm;\n" " fResult += fLogLanczos;\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) / GetBeta(fA,fB);\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)/(fA+1.0)*fX;\n" " if (b2 == 0.0)\n" " {\n" " a2 = 0.0;\n" " fnorm = 1.0;\n" " cf = 1.0;\n" " }\n" " else\n" " {\n" " a2 = 1.0;\n" " fnorm = 1.0/b2;\n" " cf = a2*fnorm;\n" " }\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/((apl2m-1.0)*apl2m);\n" " d2m1 = -(fA+rm)*(fA+fB+rm)*fX/(apl2m*(apl2m+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 = 1.0/b2;\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 0.0;//means wrong condition.\n" " double fAy = GetValue(fAx,fp,fDF);\n" " double fBy = GetValue(fBx,fp,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 = GetValue(fAx,fp,fDF);\n" " }\n" " else\n" " {\n" " fTemp = fBx;\n" " fBx += 2.0 * (fBx - fAx);\n" " fAx = fTemp;\n" " fAy = fBy;\n" " fBy = GetValue(fBx,fp,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) > 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 / (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" " \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" " 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"; #endif /* vim:set shiftwidth=4 softtabstop=4 expandtab: */