/* -*- 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/.
 */

#pragma once

const char MinDecl[] = "#define Min 2.22507e-308\n";
const char fBigInvDecl[] ="#define fBigInv  2.22045e-016\n";
const char fMachEpsDecl[] ="#define fMachEps  2.22045e-016\n";
const char fLogDblMaxDecl[] ="#define fLogDblMax  log(1.79769e+308)\n";
const char fHalfMachEpsDecl[] ="#define fHalfMachEps  0.5*2.22045e-016\n";
const char fMaxGammaArgumentDecl[] =
"#define fMaxGammaArgument 171.624376956302\n";
const char GetValueDecl[] =
"double  GetValue( double x,double fp,double fDF );\n";
const char GetValue[] =
"double  GetValue( double x,double fp,double fDF )\n"
"{\n"
"    return fp - 2 * GetTDist(x, fDF);\n"
"}\n";
const char fmin_countDecl[] = "double fmin_count(double a, double b, __private int *p);\n";
const char fmin_count[] =
"double fmin_count(double a, double b, __private int *p) {\n"
"    double result = fmin(a, b);\n"
"    bool t = isnan(result);\n"
"    (*p) += t?0:1;\n"
"    return result;\n"
"}\n";
const char fmax_countDecl[] =  "double fmax_count(double a, double b, __private int *p);\n";
const char fmax_count[] =
"double fmax_count(double a, double b, __private int *p) {\n"
"    double result = fmax(a, b);\n"
"    bool t = isnan(result);\n"
"    (*p) += t?0:1;\n"
"    return result;\n"
"}\n";
const char fsum_countDecl[] = "double fsum_count(double a, double b, __private int *p);\n";
const char fsum_count[] =
"double fsum_count(double a, double b, __private int *p) {\n"
"    bool t = isnan(a);\n"
"    (*p) += t?0:1;\n"
"    return t?b:a+b;\n"
"}\n";
const char GetGammaSeriesDecl[] =
"double GetGammaSeries( double fA, double fX );\n";
const char 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";
const char GetGammaContFractionDecl[] =  "double GetGammaContFraction( double "
"fA, double fX );\n";
const char 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";
const char GetLowRegIGammaDecl[] = "double GetLowRegIGamma( double "
"fA, double fX );\n";
const char 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";
const char GetGammaDistDecl[] = "double GetGammaDist( double fX, double "
"fAlpha, double fLambda );\n";
const char 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";
const char GetGammaDistPDFDecl[] =  "double GetGammaDistPDF( double fX"
", double fAlpha, double fLambda );\n";
const char 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";
const char GetBetaDistDecl[] =
"double GetBetaDist(double fXin, double fAlpha, double fBeta);\n";
const char 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)/fA;\n"
"    double fP = fA/(fA+fB);\n"
"    double fQ = fB/(fA+fB);\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";

const char GetFDistDecl[] =
    "double GetFDist(double x, double fF1, double fF2);\n";
const char 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";
const char GetGammaInvValueDecl[] = "double"
" GetGammaInvValue(double fAlpha,double fBeta,double fX1 );\n";
const char 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";
const char GetFInvValueDecl[] = "double GetFInvValue(double fF1,double fF2"
",double fX1 );";
const char 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";
const char GetBinomDistPMFDecl[] =
    "double GetBinomDistPMF(double x, double n, double p);";
const char 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)/(double)(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)/(double)(i + 1)*p/q;\n"
"       return fFactor;\n"
"   }\n"
"}\n";

const char lcl_GetBinomDistRangeDecl[] =
    "double lcl_GetBinomDistRange(double n, \n"
"double xs, double xe, double fFactor, double p, double q);";
const char 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)/(double)(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)/(double)(i)*p/q;\n"
"       fSum += fFactor;\n"
"   }\n"
"   return (fSum > 1.0) ? 1.0 : fSum;\n"
"}\n";

const char GetLogGammaDecl[] = "double GetLogGamma(double fZ);\n";
const char 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";

const char GetChiDistDecl[] = "double GetChiDist(double fX, double fDF);\n";
const char 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";

const char GetChiSqDistCDFDecl[] =
"double GetChiSqDistCDF(double fX, double fDF);\n";
const char 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";

const char GetChiSqDistPDFDecl[] =
"double GetChiSqDistPDF(double fX, double fDF);\n";
const char 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.0/sqrt(fX*2*M_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";

const char lcl_IterateInverseBetaInvDecl[] =
"static double lcl_IterateInverseBetaInv(double fp, double fAlpha, \n"
"   double fBeta, double fAx, double fBx, bool *rConvError );\n";
const char 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";

const char lcl_IterateInverseChiInvDecl[] =
"static double lcl_IterateInverseChiInv"
"(double fp, double fdf, double fAx, double fBx, bool *rConvError);\n";
const char 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";

const char lcl_IterateInverseChiSQInvDecl[] =
"static double lcl_IterateInverseChiSQInv( double fp, double fdf, \n"
"   double fAx, double fBx, bool *rConvError );\n";
const char 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";

const char gaussinvDecl[] = "double gaussinv(double x);\n";
const char 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"
"    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"
"        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"
"        if(q<0.0) z=-z;\n"
"    }\n"
"    return z;\n"
"}\n";

const char lcl_GetLogGammaHelperDecl[] =
"static double lcl_GetLogGammaHelper(double fZ);\n";
const char 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";
const char lcl_GetGammaHelperDecl[] =
"static double lcl_GetGammaHelper(double fZ);\n";
const char 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 = 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";
const char lcl_getLanczosSumDecl[] =
"static double lcl_getLanczosSum(double fZ);\n";
const char 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 = 1.0/fZ;\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/fSumDenom;\n"
"}\n";

const char GetUpRegIGammaDecl[] =
" double GetUpRegIGamma( double fA, double fX ) ;\n";
const char 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";

const char lcl_HasChangeOfSignDecl[] =
"static inline bool lcl_HasChangeOfSign( double u, double w );\n";
const char 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";

const char GetTDistDecl[] =" double GetTDist(double T, double fDF);\n";
const char GetTDist[] =
"double GetTDist(double T, double fDF)\n"
"{\n"
"    return 0.5 * GetBetaDist(fDF/(fDF+T*T),fDF/2.0, 0.5);\n"
"}\n";

const char GetBetaDecl[] =" double GetBeta(double fAlpha, double fBeta);\n";
const char 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)/tgamma(fAB)*tgamma(fB);\n"
"    double fgm = 5.524680040776729583740234375;\n"
"    double fLanczos = lcl_getLanczosSum(fA)*lcl_getLanczosSum(fB)\n"
"        /lcl_getLanczosSum(fAB);\n"
"    fLanczos *= sqrt(((fAB + fgm)/(fA + fgm))/(fB + fgm));\n"
"    return fLanczos * pow(exp(1.0),(-fA*log1p(fB/(fA + fgm)))"
"                    - fB*log1p(fA/(fB + fgm)) - fgm);\n"
"}\n";

const char GetLogBetaDecl[] =
" double GetLogBeta(double fAlpha, double fBeta);\n";
const char 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"
"        /lcl_getLanczosSum(fA + fB);\n"
"    double fResult= -fA *log1p(fB/(fA + fgm))"
"-fB *log1p(fA/(fB + fgm))-fgm;\n"
"    fResult += log(fLanczos)+0.5*(log(fA + fB + fgm) - log(fA + fgm)\n"
"        - log(fB + fgm));\n"
"    return fResult;\n"
"}\n";

const char GetBetaDistPDFDecl[] =
"double GetBetaDistPDF(double fX, double fA, double fB);\n";
const char 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";

const char lcl_GetBetaHelperContFracDecl[] =
"double lcl_GetBetaHelperContFrac(double fX, double fA, double fB);\n";
const char 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"
"    b2==0.0?(a2 = 0.0,fnorm = 1.0,cf = 1.0):\n"
"        (a2 = 1.0,fnorm = 1.0/b2,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/(apl2m*(apl2m-1.0));\n"
"      d2m1 = -((fA+rm)*(fA+rm+fB))*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";

const char lcl_IterateInverseDecl[] =
"double lcl_IterateInverse("
"double fAx, double fBx, bool* rConvError,double fp,double fDF );\n";
const char 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 / (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"
"        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";
const char phiDecl[] =
"double phi(double x);\n";
const char phi[] =
"double phi(double x)\n"
"{\n"
"    return  0.39894228040143268 * exp(-(x * x) / 2.0);\n"
"}\n";
const char taylorDecl[] =
"double taylor(double* pPolynom, uint nMax, double x);\n";
const char 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"
"}";
const char gaussDecl[] = "double gauss(double x);\n";
const char 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";

/* vim:set shiftwidth=4 softtabstop=4 expandtab: */
