GPU Calc: implemented GAMMADIST

AMLOEXT-157 FIX

Change-Id: I806b6522e7481af87c03fbd1b1d500eb5cd1b9bf
Signed-off-by: haochen <haochen@multicorewareinc.com>
Signed-off-by: I-Jui (Ray) Sung <ray@multicorewareinc.com>

4 files changed, 1181 insertions, 1 deletions

diff --git a/sc/source/core/opencl/formulagroupcl.cxx b/sc/source/core/opencl/formulagroupcl.cxx
index 68f90eeff458..f191af07b305 100644
--- a/sc/source/core/opencl/formulagroupcl.cxx
+++ b/sc/source/core/opencl/formulagroupcl.cxx
@@ -1188,6 +1188,10 @@ DynamicKernelSoPArguments::DynamicKernelSoPArguments(
                 mvSubArguments.push_back(SoPHelper(ts,
                          ft->Children[i], new OpLogNormDist));
                 break;
+            case ocGammaDist:
+                mvSubArguments.push_back(SoPHelper(ts,
+                         ft->Children[i], new OpGammaDist));
+                break;
             case ocExternal:
                 if ( !(pChild->GetExternal().compareTo(OUString(
                     "com.sun.star.sheet.addin.Analysis.getEffect"))))
diff --git a/sc/source/core/opencl/op_statistical.cxx b/sc/source/core/opencl/op_statistical.cxx
index a58eee8a91f3..8e39c363b5f7 100644
--- a/sc/source/core/opencl/op_statistical.cxx
+++ b/sc/source/core/opencl/op_statistical.cxx
@@ -2471,7 +2471,88 @@ void OpLogNormDist::GenSlidingWindowFunction(std::stringstream &ss,
     ss << "    return tmp;\n";
     ss << "}\n";
 }
+void OpGammaDist::BinInlineFun(std::set<std::string>& decls,
+    std::set<std::string>& funs)
+{
+    decls.insert(fBigInvDecl);decls.insert(fLogDblMaxDecl);
+    decls.insert(fHalfMachEpsDecl);decls.insert(fMaxGammaArgumentDecl);
+    decls.insert(GetGammaSeriesDecl);decls.insert(GetGammaContFractionDecl);
+    decls.insert(GetLowRegIGammaDecl);decls.insert(GetGammaDistDecl);
+    decls.insert(GetGammaDistPDFDecl);
+    funs.insert(GetGammaSeries);funs.insert(GetGammaContFraction);
+    funs.insert(GetLowRegIGamma);funs.insert(GetGammaDist);
+    funs.insert(GetGammaDistPDF);
+}
 
+void OpGammaDist::GenSlidingWindowFunction(std::stringstream &ss,
+            const std::string sSymName, SubArguments &vSubArguments)
+{
+        FormulaToken *tmpCur0 = vSubArguments[0]->GetFormulaToken();
+        const formula::SingleVectorRefToken*tmpCurDVR0= dynamic_cast<const
+formula::SingleVectorRefToken *>(tmpCur0);
+        FormulaToken *tmpCur1 = vSubArguments[1]->GetFormulaToken();
+        const formula::SingleVectorRefToken*tmpCurDVR1= dynamic_cast<const
+formula::SingleVectorRefToken *>(tmpCur1);
+        FormulaToken *tmpCur2 = vSubArguments[2]->GetFormulaToken();
+        const formula::SingleVectorRefToken*tmpCurDVR2= dynamic_cast<const
+formula::SingleVectorRefToken *>(tmpCur2);
+        FormulaToken *tmpCur3 = vSubArguments[3]->GetFormulaToken();
+        const formula::SingleVectorRefToken*tmpCurDVR3= dynamic_cast<const
+formula::SingleVectorRefToken *>(tmpCur3);
+        ss << "\ndouble " << sSymName;
+        ss << "_"<< BinFuncName() <<"(";
+        for (unsigned i = 0; i < vSubArguments.size(); i++)
+        {
+            if (i)
+                ss << ",";
+            vSubArguments[i]->GenSlidingWindowDecl(ss);
+        }
+        ss << ") {\n";
+        ss << "    int gid0=get_global_id(0);\n";
+        ss << "    double arg0 = ";
+        ss << vSubArguments[0]->GenSlidingWindowDeclRef();
+        ss << ";\n";
+        ss << "    double arg1 = ";
+        ss << vSubArguments[1]->GenSlidingWindowDeclRef();
+        ss << ";\n";
+        ss << "    double arg2 = ";
+        ss << vSubArguments[2]->GenSlidingWindowDeclRef();
+        ss << ";\n";
+        ss << "    double arg3 = ";
+        ss << vSubArguments[3]->GenSlidingWindowDeclRef();
+        ss << ";\n";
+        ss << "    double tmp;\n";
+#ifdef ISNAN
+        ss << "    if(isNan(arg0)||(gid0>=";
+        ss << tmpCurDVR0->GetArrayLength();
+        ss << "))\n";
+        ss << "        arg0 = 0;\n";
+#endif
+#ifdef ISNAN
+        ss << "    if(isNan(arg1)||(gid0>=";
+        ss << tmpCurDVR1->GetArrayLength();
+        ss << "))\n";
+        ss << "        arg1 = 0;\n";
+#endif
+#ifdef ISNAN
+        ss << "    if(isNan(arg2)||(gid0>=";
+        ss << tmpCurDVR2->GetArrayLength();
+        ss << "))\n";
+        ss << "        arg2 = 0;\n";
+#endif
+#ifdef ISNAN
+        ss << "    if(isNan(arg3)||(gid0>=";
+        ss << tmpCurDVR3->GetArrayLength();
+        ss << "))\n";
+        ss << "        arg3 = 0;\n";
+#endif
+        ss << "    if (arg3)\n";
+        ss << "        tmp=GetGammaDist( arg0, arg1, arg2);\n";
+        ss << "    else\n";
+        ss << "        tmp=GetGammaDistPDF( arg0, arg1, arg2);\n";
+        ss << "    return tmp;\n";
+        ss << "}\n";
+}
 }}
 
 /* vim:set shiftwidth=4 softtabstop=4 expandtab: */
diff --git a/sc/source/core/opencl/op_statistical.hxx b/sc/source/core/opencl/op_statistical.hxx
index 92b67dc4b5d9..7d8fa5ebb931 100644
--- a/sc/source/core/opencl/op_statistical.hxx
+++ b/sc/source/core/opencl/op_statistical.hxx
@@ -225,6 +225,14 @@ public:
             const std::string sSymName, SubArguments &vSubArguments);
     virtual std::string BinFuncName(void) const { return "LogNormdist"; }
 };
+class OpGammaDist: public Normal
+{
+public:
+    virtual void GenSlidingWindowFunction(std::stringstream &ss,
+            const std::string sSymName, SubArguments &vSubArguments);
+    void BinInlineFun(std::set<std::string>& decls,std::set<std::string>& funs);
+    virtual std::string BinFuncName(void) const { return "GammaDist"; }
+};
 }}
 
 #endif
diff --git a/sc/source/core/opencl/opinlinefun_statistical.cxx b/sc/source/core/opencl/opinlinefun_statistical.cxx
index 31e9dc251206..f9d357995f62 100644
--- a/sc/source/core/opencl/opinlinefun_statistical.cxx
+++ b/sc/source/core/opencl/opinlinefun_statistical.cxx
@@ -9,13 +9,711 @@
 
 #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 MinDecl = "#define Min 2.22507e-308\n";
 
 std::string fMaxGammaArgumentDecl =
 "#define fMaxGammaArgument  171.624376956302\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/(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"
+"   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 =
@@ -178,6 +876,395 @@ std::string gaussinv =
 "    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