LCOV - code coverage report
Current view: top level - sc/source/core/opencl - opinlinefun_statistical.cxx (source / functions) Hit Total Coverage
Test: commit 10e77ab3ff6f4314137acd6e2702a6e5c1ce1fae Lines: 85 85 100.0 %
Date: 2014-11-03 Functions: 0 0 -
Legend: Lines: hit not hit

          Line data    Source code
       1             : /* -*- Mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */
       2             : /*
       3             :  * This file is part of the LibreOffice project.
       4             :  *
       5             :  * This Source Code Form is subject to the terms of the Mozilla Public
       6             :  * License, v. 2.0. If a copy of the MPL was not distributed with this
       7             :  * file, You can obtain one at http://mozilla.org/MPL/2.0/.
       8             :  */
       9             : 
      10             : #ifndef SC_OPENCL_OPINLINFUN_statistical
      11             : #define SC_OPENCL_OPINLINFUN_statistical
      12          76 : std::string MinDecl = "#define Min 2.22507e-308\n";
      13          76 : std::string F_PIDecl="#define F_PI 3.1415926\n";
      14          76 : std::string fBigInvDecl ="#define fBigInv  2.22045e-016\n";
      15          76 : std::string fMachEpsDecl ="#define fMachEps  2.22045e-016\n";
      16          76 : std::string fLogDblMaxDecl ="#define fLogDblMax  log(1.79769e+308)\n";
      17          76 : std::string fHalfMachEpsDecl ="#define fHalfMachEps  0.5*2.22045e-016\n";
      18          76 : std::string fMaxGammaArgumentDecl=
      19             : "#define fMaxGammaArgument 171.624376956302\n";
      20          76 : std::string GetValueDecl=
      21             : "double  GetValue( double x,double fp,double fDF );\n";
      22          76 : std::string GetValue=
      23             : "double  GetValue( double x,double fp,double fDF )\n"
      24             : "{\n"
      25             : "    return fp - 2 * GetTDist(x, fDF);\n"
      26             : "}\n";
      27          76 : std::string GetGammaSeriesDecl=
      28             : "double GetGammaSeries( double fA, double fX );\n";
      29          76 : std::string GetGammaSeries =
      30             : "double GetGammaSeries( double fA, double fX )\n"
      31             : "{\n"
      32             : "    double fDenomfactor = fA;\n"
      33             : "     double fSummand = 1.0/fA;\n"
      34             : "    double fSum = fSummand;\n"
      35             : "    int nCount=1;\n"
      36             : "    do\n"
      37             : "    {\n"
      38             : "        fDenomfactor = fDenomfactor + 1.0;\n"
      39             : "        fSummand = fSummand * fX/fDenomfactor;\n"
      40             : "        fSum = fSum + fSummand;\n"
      41             : "        nCount = nCount+1;\n"
      42             : "    } while ( fSummand/fSum > fHalfMachEps && nCount<=10000);\n"
      43             : "    if (nCount>10000)\n"
      44             : "    {\n"
      45             : "    }\n"
      46             : "    return fSum;\n"
      47             : "}\n";
      48          76 : std::string GetGammaContFractionDecl =  "double GetGammaContFraction( double "
      49             : "fA, double fX );\n";
      50          76 : std::string GetGammaContFraction =
      51             : "double GetGammaContFraction( double fA, double fX )\n"
      52             : "{\n"
      53             : "    double fBig = 1.0/fBigInv;\n"
      54             : "    double fCount = 0.0;\n"
      55             : "    double fNum = 0.0;\n"
      56             : "    double fY = 1.0 - fA;\n"
      57             : "    double fDenom = fX + 2.0-fA;\n"
      58             : "    double fPk = 0.0;\n"
      59             : "    double fPkm1 = fX + 1.0;\n"
      60             : "    double fPkm2 = 1.0;\n"
      61             : "    double fQk = 1.0;\n"
      62             : "    double fQkm1 = fDenom * fX;\n"
      63             : "    double fQkm2 = fX;\n"
      64             : "    double fApprox = fPkm1/fQkm1;\n"
      65             : "    bool bFinished = false;\n"
      66             : "    double fR = 0.0;\n"
      67             : "    do\n"
      68             : "    {\n"
      69             : "        fCount = fCount +1.0;\n"
      70             : "        fY = fY+ 1.0;\n"
      71             : "        fNum = fY * fCount;\n"
      72             : "        fDenom = fDenom +2.0;\n"
      73             : "        fPk = fPkm1 * fDenom  -  fPkm2 * fNum;\n"
      74             : "        fQk = fQkm1 * fDenom  -  fQkm2 * fNum;\n"
      75             : "        if (fQk != 0.0)\n"
      76             : "        {\n"
      77             : "            fR = fPk/fQk;\n"
      78             : "            bFinished = (fabs( (fApprox - fR)/fR ) <= fHalfMachEps);\n"
      79             : "            fApprox = fR;\n"
      80             : "        }\n"
      81             : "        fPkm2 = fPkm1;\n"
      82             : "        fPkm1 = fPk;\n"
      83             : "        fQkm2 = fQkm1;\n"
      84             : "        fQkm1 = fQk;\n"
      85             : "        if (fabs(fPk) > fBig)\n"
      86             : "        {\n"
      87             : "            fPkm2 = fPkm2 * fBigInv;\n"
      88             : "            fPkm1 = fPkm1 * fBigInv;\n"
      89             : "            fQkm2 = fQkm2 * fBigInv;\n"
      90             : "            fQkm1 = fQkm1 * fBigInv;\n"
      91             : "        }\n"
      92             : "    } while (!bFinished && fCount<10000);\n"
      93             : "    if (!bFinished)\n"
      94             : "    {\n"
      95             : "    }\n"
      96             : "    return fApprox;\n"
      97             : "}\n";
      98          76 : std::string GetLowRegIGammaDecl = "double GetLowRegIGamma( double "
      99             : "fA, double fX );\n";
     100          76 : std::string GetLowRegIGamma =
     101             : "double GetLowRegIGamma( double fA, double fX )\n"
     102             : "{\n"
     103             : "    double fLnFactor = fA * log(fX) - fX - lgamma(fA);\n"
     104             : "    double fFactor = exp(fLnFactor);\n"
     105             : "    if (fX>fA+1.0) \n"
     106             : "        return 1.0 - fFactor * GetGammaContFraction(fA,fX);\n"
     107             : "    else\n"
     108             : "        return fFactor * GetGammaSeries(fA,fX);\n"
     109             : "}\n";
     110          76 : std::string GetGammaDistDecl = "double GetGammaDist( double fX, double "
     111             : "fAlpha, double fLambda );\n";
     112          76 : std::string GetGammaDist =
     113             : "double GetGammaDist( double fX, double fAlpha, double fLambda )\n"
     114             : "{\n"
     115             : "    if (fX <= 0.0)\n"
     116             : "        return 0.0;\n"
     117             : "    else\n"
     118             : "        return GetLowRegIGamma( fAlpha, fX / fLambda);\n"
     119             : "}\n";
     120          76 : std::string GetGammaDistPDFDecl =  "double GetGammaDistPDF( double fX"
     121             : ", double fAlpha, double fLambda );\n";
     122          76 : std::string GetGammaDistPDF =
     123             : "double GetGammaDistPDF( double fX, double fAlpha, double fLambda )\n"
     124             : "{\n"
     125             : "    if (fX < 0.0)\n"
     126             : "        return 0.0;\n"
     127             : "    else if (fX == 0)\n"
     128             : "    {\n"
     129             : "        if (fAlpha < 1.0)\n"
     130             : "        {\n"
     131             : "            return HUGE_VAL;\n"
     132             : "        }\n"
     133             : "        else if (fAlpha == 1)\n"
     134             : "        {\n"
     135             : "            return (1.0 / fLambda);\n"
     136             : "        }\n"
     137             : "        else\n"
     138             : "        {\n"
     139             : "            return 0.0;\n"
     140             : "        }\n"
     141             : "    }\n"
     142             : "    else\n"
     143             : "    {\n"
     144             : "        double fXr = fX / fLambda;\n"
     145             : "        if (fXr > 1.0)\n"
     146             : "        {\n"
     147             : "            if (log(fXr) * (fAlpha-1.0) < fLogDblMax &&"
     148             : "fAlpha < fMaxGammaArgument)\n"
     149             : "            {\n"
     150             : "                return pow( fXr, fAlpha-1.0) * exp(-fXr) / "
     151             : "fLambda / tgamma(fAlpha);\n"
     152             : "            }\n"
     153             : "            else\n"
     154             : "            {\n"
     155             : "                return exp( (fAlpha-1.0) * log(fXr) - fXr - "
     156             : "log(fLambda) - lgamma(fAlpha));\n"
     157             : "            }\n"
     158             : "        }\n"
     159             : "        else\n"
     160             : "        {\n"
     161             : "            if (fAlpha<fMaxGammaArgument)\n"
     162             : "            {\n"
     163             : "                return pow( fXr, fAlpha-1.0) * exp(-fXr) / "
     164             : "fLambda / tgamma(fAlpha);\n"
     165             : "            }\n"
     166             : "            else\n"
     167             : "            {\n"
     168             : "                return pow( fXr, fAlpha-1.0) * exp(-fXr) / "
     169             : "fLambda / exp( lgamma(fAlpha));\n"
     170             : "            }\n"
     171             : "        }\n"
     172             : "    }\n"
     173             : "}\n";
     174          76 : std::string GetBetaDistDecl =
     175             : "double GetBetaDist(double fXin, double fAlpha, double fBeta);\n";
     176          76 : std::string GetBetaDist =
     177             : "double GetBetaDist(double fXin, double fAlpha, double fBeta)\n"
     178             : "{\n"
     179             : "    if (fXin <= 0.0)\n"
     180             : "        return 0.0;\n"
     181             : "    if (fXin >= 1.0)\n"
     182             : "        return 1.0;\n"
     183             : "    if (fBeta == 1.0)\n"
     184             : "        return pow(fXin, fAlpha);\n"
     185             : "    if (fAlpha == 1.0)\n"
     186             : "        return -expm1(fBeta * log1p(-fXin));\n"
     187             : "    double fResult;\n"
     188             : "    double fY = (0.5-fXin)+0.5;\n"
     189             : "    double flnY = log1p(-fXin);\n"
     190             : "    double fX = fXin;\n"
     191             : "    double flnX = log(fXin);\n"
     192             : "    double fA = fAlpha;\n"
     193             : "    double fB = fBeta;\n"
     194             : "    bool bReflect = fXin > fAlpha*pow((fAlpha+fBeta),-1.0);\n"
     195             : "    if (bReflect)\n"
     196             : "    {\n"
     197             : "        fA = fBeta;\n"
     198             : "        fB = fAlpha;\n"
     199             : "        fX = fY;\n"
     200             : "        fY = fXin;\n"
     201             : "        flnX = flnY;\n"
     202             : "        flnY = log(fXin);\n"
     203             : "    }\n"
     204             : "    fResult = lcl_GetBetaHelperContFrac(fX,fA,fB)*pow(fA,-1.0);\n"
     205             : "    double fP = fA*pow((fA+fB),-1.0);\n"
     206             : "    double fQ = fB*pow((fA+fB),-1.0);\n"
     207             : "    if (fA > 1.0 && fB > 1.0 && fP < 0.97 && fQ < 0.97)\n"
     208             : "        fResult *= GetBetaDistPDF(fX,fA,fB)*fX*fY;\n"
     209             : "    else\n"
     210             : "        fResult *= pow(exp(1.0),(fA*flnX + fB*flnY - GetLogBeta(fA,fB)));\n"
     211             : "    if (bReflect)\n"
     212             : "        fResult = 0.5 - fResult + 0.5;\n"
     213             : "    if (fResult > 1.0)\n"
     214             : "        fResult = 1.0;\n"
     215             : "    if (fResult < 0.0)\n"
     216             : "        fResult = 0.0;\n"
     217             : "    return fResult;\n"
     218             : "}\n";
     219             : 
     220          76 : std::string GetFDistDecl =
     221             :     "double GetFDist(double x, double fF1, double fF2);\n";
     222          76 : std::string GetFDist =
     223             : "double GetFDist(double x, double fF1, double fF2)\n"
     224             : "{\n"
     225             : "    double arg = fF2*pow((fF2+fF1*x),-1.0);\n"
     226             : "    double alpha = fF2*pow(2.0,-1.0);\n"
     227             : "    double beta = fF1*pow(2.0,-1.0);\n"
     228             : "    return (GetBetaDist(arg, alpha, beta));\n"
     229             : "}\n";
     230          76 : std::string GetGammaInvValueDecl = "double"
     231             : " GetGammaInvValue(double fAlpha,double fBeta,double fX1 );\n";
     232          76 : std::string GetGammaInvValue =
     233             : "double GetGammaInvValue(double fAlpha,double fBeta,double fX1 )\n"
     234             : "{\n"
     235             : "    if (fX1 <= 0.0)\n"
     236             : "        return 0.0;\n"
     237             : "    else\n"
     238             : "    {\n"
     239             : "        double fX=fX1*pow(fBeta,-1.0);\n"
     240             : "        double fLnFactor = fAlpha * log(fX) - fX - lgamma(fAlpha);\n"
     241             : "        double fFactor = exp(fLnFactor);\n"
     242             : "        if (fX>fAlpha+1.0)\n"
     243             : "            return 1.0 - fFactor * GetGammaContFraction(fAlpha,fX);\n"
     244             : "        else\n"
     245             : "            return fFactor * GetGammaSeries(fAlpha,fX);\n"
     246             : "    }\n"
     247             : "}\n";
     248          76 : std::string GetFInvValueDecl = "double GetFInvValue(double fF1,double fF2"
     249             : ",double fX1 );";
     250          76 : std::string GetFInvValue =
     251             : "double GetFInvValue(double fF1,double fF2,double fX1 )\n"
     252             : "{\n"
     253             : "    double arg = fF2*pow((fF2+fF1*fX1),-1.0);\n"
     254             : "    double alpha = fF2*pow(2.0,-1.0);\n"
     255             : "    double beta = fF1*pow(2.0,-1.0);\n"
     256             : "    double fXin,fAlpha,fBeta;\n"
     257             : "        fXin=arg;fAlpha=alpha;fBeta=beta;\n"
     258             : "    if (fXin <= 0.0)\n"
     259             : "        return 0.0;\n"
     260             : "    if (fXin >= 1.0)\n"
     261             : "        return 1.0;\n"
     262             : "    if (fBeta == 1.0)\n"
     263             : "        return pow(fXin, fAlpha);\n"
     264             : "    if (fAlpha == 1.0)\n"
     265             : "        return -expm1(fBeta * log1p(-fXin));\n"
     266             : "    double fResult;\n"
     267             : "    double fY = (0.5-fXin)+0.5;\n"
     268             : "    double flnY = log1p(-fXin);\n"
     269             : "    double fX = fXin;\n"
     270             : "    double flnX = log(fXin);\n"
     271             : "    double fA = fAlpha;\n"
     272             : "    double fB = fBeta;\n"
     273             : "    bool bReflect = fXin > fAlpha*pow((fAlpha+fBeta),-1.0);\n"
     274             : "    if (bReflect)\n"
     275             : "    {\n"
     276             : "        fA = fBeta;\n"
     277             : "        fB = fAlpha;\n"
     278             : "        fX = fY;\n"
     279             : "        fY = fXin;\n"
     280             : "        flnX = flnY;\n"
     281             : "        flnY = log(fXin);\n"
     282             : "    }\n"
     283             : "    fResult = lcl_GetBetaHelperContFrac(fX,fA,fB);\n"
     284             : "    fResult = fResult*pow(fA,-1.0);\n"
     285             : "    double fP = fA*pow((fA+fB),-1.0);\n"
     286             : "    double fQ = fB*pow((fA+fB),-1.0);\n"
     287             : "    double fTemp;\n"
     288             : "    if (fA > 1.0 && fB > 1.0 && fP < 0.97 && fQ < 0.97)\n"
     289             : "        fTemp = GetBetaDistPDF(fX,fA,fB)*fX*fY;\n"
     290             : "    else\n"
     291             : "        fTemp = exp(fA*flnX + fB*flnY - GetLogBeta(fA,fB));\n"
     292             : "    fResult *= fTemp;\n"
     293             : "    if (bReflect)\n"
     294             : "        fResult = 0.5 - fResult + 0.5;\n"
     295             : "    if (fResult > 1.0)\n"
     296             : "        fResult = 1.0;\n"
     297             : "    if (fResult < 0.0)\n"
     298             : "        fResult = 0.0;\n"
     299             : "    return fResult;\n"
     300             : "}\n";
     301          76 : std::string GetBinomDistPMFDecl =
     302             :     "double GetBinomDistPMF(double x, double n, double p);";
     303          76 : std::string GetBinomDistPMF =
     304             : "double GetBinomDistPMF(double x, double n, double p)\n"
     305             : "{\n"
     306             : "   double q = (0.5 - p) + 0.5;\n"
     307             : "   double fFactor = pow(q, n);\n"
     308             : "   if (fFactor <= Min)\n"
     309             : "   {\n"
     310             : "       fFactor = pow(p, n);\n"
     311             : "       if (fFactor <= Min)\n"
     312             : "           return GetBetaDistPDF(p, x + 1.0, n - x + 1.0)*pow((n + 1.0),-1.0);\n"
     313             : "       else\n"
     314             : "       {\n"
     315             : "           uint max = (uint)(n - x);\n"
     316             : "           for (uint i = 0; i < max && fFactor > 0.0; ++i)\n"
     317             : "               fFactor *= (n - i)*pow((i + 1),-1.0)*q*pow(p,-1.0);\n"
     318             : "           return fFactor;\n"
     319             : "       }\n"
     320             : "   }\n"
     321             : "   else\n"
     322             : "   {\n"
     323             : "       uint max = (uint)x;\n"
     324             : "       for (uint i = 0; i < max && fFactor > 0.0; ++i)\n"
     325             : "           fFactor *= (n - i)*pow((i + 1),-1.0)*p*pow(q,-1.0);\n"
     326             : "       return fFactor;\n"
     327             : "   }\n"
     328             : "}\n";
     329             : 
     330          76 : std::string lcl_GetBinomDistRangeDecl =
     331             :     "double lcl_GetBinomDistRange(double n, \n"
     332             : "double xs, double xe, double fFactor, double p, double q);";
     333          76 : std::string lcl_GetBinomDistRange=
     334             : "double lcl_GetBinomDistRange(double n, double xs, double xe,\n"
     335             : "   double fFactor, double p, double q)\n"
     336             : "{\n"
     337             : "   uint i;\n"
     338             : "   double fSum;\n"
     339             : "   uint nXs = (uint)xs;\n"
     340             : "   for (i = 1; i <= nXs && fFactor > 0.0; ++i)\n"
     341             : "       fFactor *= (n - i + 1)*pow(i,-1.0)*p*pow(q,-1.0);\n"
     342             : "   fSum = fFactor;\n"
     343             : "   uint nXe =(uint)xe;\n"
     344             : "   for (i = nXs + 1; i <= nXe && fFactor > 0.0; ++i)\n"
     345             : "   {\n"
     346             : "       fFactor *= (n - i + 1)*pow(i,-1.0)*p*pow(q,-1.0);\n"
     347             : "       fSum += fFactor;\n"
     348             : "   }\n"
     349             : "   return (fSum > 1.0) ? 1.0 : fSum;\n"
     350             : "}\n";
     351             : 
     352          76 : std::string GetLogGammaDecl = "double GetLogGamma(double fZ);\n";
     353          76 : std::string GetLogGamma =
     354             : "double GetLogGamma(double fZ)\n"
     355             : "{\n"
     356             : "   if (fZ >= fMaxGammaArgument)\n"
     357             : "       return lcl_GetLogGammaHelper(fZ);\n"
     358             : "   if (fZ >= 1.0)\n"
     359             : "       return log(lcl_GetGammaHelper(fZ));\n"
     360             : "   if (fZ >= 0.5)\n"
     361             : "       return log( lcl_GetGammaHelper(fZ+1) *pow(fZ,-1.0));\n"
     362             : "   return lcl_GetLogGammaHelper(fZ+2) - log(fZ+1) - log(fZ);\n"
     363             : "}\n";
     364             : 
     365          76 : std::string GetChiDistDecl = "double GetChiDist(double fX, double fDF);\n";
     366          76 : std::string GetChiDist =
     367             : "double GetChiDist(double fX, double fDF)\n"
     368             : "{\n"
     369             : "   if (fX <= 0.0)\n"
     370             : "       return  1.0;\n"
     371             : "   else\n"
     372             : "       return GetUpRegIGamma( fDF*pow(2.0,-1.0), fX*pow(2.0,-1.0));\n"
     373             : "}\n";
     374             : 
     375          76 : std::string GetChiSqDistCDFDecl =
     376             : "double GetChiSqDistCDF(double fX, double fDF);\n";
     377          76 : std::string GetChiSqDistCDF =
     378             : "double GetChiSqDistCDF(double fX, double fDF)\n"
     379             : "{\n"
     380             : "   if (fX <= 0.0)\n"
     381             : "       return 0.0;"
     382             : "   else\n"
     383             : "       return GetLowRegIGamma( fDF*pow(2.0,-1.0), fX*pow(2.0,-1.0));\n"
     384             : "}\n";
     385             : 
     386          76 : std::string GetChiSqDistPDFDecl=
     387             : "double GetChiSqDistPDF(double fX, double fDF);\n";
     388          76 : std::string GetChiSqDistPDF =
     389             : "double GetChiSqDistPDF(double fX, double fDF)\n"
     390             : "{\n"
     391             : "   double fValue;\n"
     392             : "   if (fX <= 0.0)\n"
     393             : "       return 0.0;\n"
     394             : "   if (fDF*fX > 1391000.0)\n"
     395             : "   {\n"
     396             : "       fValue = exp((0.5*fDF - 1) * log(fX*0.5) - 0.5 * fX - log(2.0)"
     397             : " - lgamma(0.5*fDF));\n"
     398             : "   }\n"
     399             : "   else\n"
     400             : "   {\n"
     401             : "       double fCount;\n"
     402             : "       if (fmod(fDF,2.0)<0.5)\n"
     403             : "       {\n"
     404             : "           fValue = 0.5;\n"
     405             : "           fCount = 2.0;\n"
     406             : "       }\n"
     407             : "       else\n"
     408             : "       {\n"
     409             : "           fValue = pow(sqrt(fX*2*F_PI),-1.0);\n"
     410             : "           fCount = 1.0;\n"
     411             : "       }\n"
     412             : "       while ( fCount < fDF)\n"
     413             : "       {\n"
     414             : "           fValue *= (fX *pow(fCount,-1.0));\n"
     415             : "           fCount += 2.0;\n"
     416             : "       }\n"
     417             : "       if (fX>=1425.0)\n"
     418             : "           fValue = exp(log(fValue)-fX*pow(2,-1.0));\n"
     419             : "       else\n"
     420             : "           fValue *= exp(-fX*pow(2,-1.0));\n"
     421             : "   }\n"
     422             : "    return fValue;\n"
     423             : "}\n";
     424             : 
     425          76 : std::string lcl_IterateInverseBetaInvDecl =
     426             : "static double lcl_IterateInverseBetaInv(double fp, double fAlpha, \n"
     427             : "   double fBeta, double fAx, double fBx, bool *rConvError );\n";
     428          76 : std::string lcl_IterateInverseBetaInv =
     429             : "static double lcl_IterateInverseBetaInv(double fp, double fAlpha, \n"
     430             : "   double fBeta, double fAx, double fBx, bool *rConvError )\n"
     431             : "{\n"
     432             : "   *rConvError = false;\n"
     433             : "    double fYEps = 1.0E-307;\n"
     434             : "    double fXEps = fMachEps;\n"
     435             : "   if(!(fAx < fBx))\n"
     436             : "   {\n"
     437             : "       //print error\n"
     438             : "   }\n"
     439             : "   double fAy = fp - GetBetaDist(fAx, fAlpha, fBeta);\n"
     440             : "   double fBy = fp - GetBetaDist(fBx, fAlpha, fBeta);\n"
     441             : "   double fTemp;\n"
     442             : "   unsigned short nCount;\n"
     443             : "   for (nCount = 0; nCount < 1000 && !lcl_HasChangeOfSign(fAy,fBy);"
     444             : " nCount++)\n"
     445             : "   {\n"
     446             : "       if (fabs(fAy) <= fabs(fBy))\n"
     447             : "       {\n"
     448             : "           fTemp = fAx;\n"
     449             : "           fAx += 2.0 * (fAx - fBx);\n"
     450             : "           if (fAx < 0.0)\n"
     451             : "               fAx = 0.0;\n"
     452             : "           fBx = fTemp;\n"
     453             : "           fBy = fAy;\n"
     454             : "           fAy = fp - GetBetaDist(fAx, fAlpha, fBeta);\n"
     455             : "       }\n"
     456             : "       else\n"
     457             : "       {\n"
     458             : "           fTemp = fBx;\n"
     459             : "           fBx += 2.0 * (fBx - fAx);\n"
     460             : "           fAx = fTemp;\n"
     461             : "           fAy = fBy;\n"
     462             : "           fBy = fp - GetBetaDist(fBx, fAlpha, fBeta);\n"
     463             : "       }\n"
     464             : "   }\n"
     465             : "   if (fAy == 0.0)\n"
     466             : "       return fAx;\n"
     467             : "   if (fBy == 0.0)\n"
     468             : "       return fBx;\n"
     469             : "   if (!lcl_HasChangeOfSign( fAy, fBy))\n"
     470             : "   {\n"
     471             : "       *rConvError = true;\n"
     472             : "       return 0.0;\n"
     473             : "   }\n"
     474             : "   double fPx = fAx;\n"
     475             : "   double fPy = fAy;\n"
     476             : "   double fQx = fBx;\n"
     477             : "   double fQy = fBy;\n"
     478             : "   double fRx = fAx;\n"
     479             : "   double fRy = fAy;\n"
     480             : "   double fSx = 0.5 * (fAx + fBx);\n"
     481             : "   bool bHasToInterpolate = true;\n"
     482             : "   nCount = 0;\n"
     483             : "   while ( nCount < 500 && fabs(fRy) > fYEps &&\n"
     484             : "               (fBx-fAx) > fmax( fabs(fAx), fabs(fBx)) * fXEps )\n"
     485             : "   {\n"
     486             : "       if (bHasToInterpolate)\n"
     487             : "       {\n"
     488             : "           if (fPy!=fQy && fQy!=fRy && fRy!=fPy)\n"
     489             : "           {\n"
     490             : "               fSx = fPx*fRy*fQy*pow(fRy-fPy,-1.0)*pow(fQy-fPy,-1.0)\n"
     491             : "                   + fRx*fQy*fPy*pow(fQy-fRy,-1.0)*pow(fPy-fRy,-1.0)\n"
     492             : "                   + fQx*fPy*fRy*pow(fPy-fQy,-1.0)*pow(fRy-fQy,-1.0);\n"
     493             : "               bHasToInterpolate = (fAx < fSx) && (fSx < fBx);\n"
     494             : "           }\n"
     495             : "           else\n"
     496             : "               bHasToInterpolate = false;\n"
     497             : "       }\n"
     498             : "       if(!bHasToInterpolate)\n"
     499             : "       {\n"
     500             : "           fSx = 0.5 * (fAx + fBx);\n"
     501             : "           fPx = fAx; fPy = fAy;\n"
     502             : "           fQx = fBx; fQy = fBy;\n"
     503             : "           bHasToInterpolate = true;\n"
     504             : "       }\n"
     505             : "       fPx = fQx; fQx = fRx; fRx = fSx;\n"
     506             : "       fPy = fQy; fQy = fRy; fRy = fp - GetBetaDist(fSx, fAlpha, fBeta);\n"
     507             : "       if (lcl_HasChangeOfSign( fAy, fRy))\n"
     508             : "       {\n"
     509             : "           fBx = fRx; fBy = fRy;\n"
     510             : "       }\n"
     511             : "       else\n"
     512             : "       {\n"
     513             : "           fAx = fRx; fAy = fRy;\n"
     514             : "       }\n"
     515             : "       bHasToInterpolate = bHasToInterpolate && (fabs(fRy) *"
     516             : " 2.0 <= fabs(fQy));\n"
     517             : "       ++nCount;\n"
     518             : "   }\n"
     519             : "   return fRx;\n"
     520             : "}\n";
     521             : 
     522          76 : std::string lcl_IterateInverseChiInvDecl =
     523             : "static double lcl_IterateInverseChiInv"
     524             : "(double fp, double fdf, double fAx, double fBx, bool *rConvError);\n";
     525          76 : std::string lcl_IterateInverseChiInv =
     526             : "static double lcl_IterateInverseChiInv"
     527             : "(double fp, double fdf, double fAx, double fBx, bool *rConvError)\n"
     528             : "{\n"
     529             : "   *rConvError = false;\n"
     530             : "    double fYEps = 1.0E-307;\n"
     531             : "    double fXEps = fMachEps;\n"
     532             : "   if(!(fAx < fBx))\n"
     533             : "   {\n"
     534             : "       //print error\n"
     535             : "   }"
     536             : "   double fAy = fp - GetChiDist(fAx, fdf);\n"
     537             : "   double fBy = fp - GetChiDist(fBx, fdf);\n"
     538             : "   double fTemp;\n"
     539             : "   unsigned short nCount;\n"
     540             : "   for (nCount = 0; nCount < 1000 && "
     541             : "!lcl_HasChangeOfSign(fAy,fBy); nCount++)\n"
     542             : "   {\n"
     543             : "       if (fabs(fAy) <= fabs(fBy))\n"
     544             : "       {\n"
     545             : "           fTemp = fAx;\n"
     546             : "           fAx += 2.0 * (fAx - fBx);\n"
     547             : "           if (fAx < 0.0)\n"
     548             : "               fAx = 0.0;\n"
     549             : "           fBx = fTemp;\n"
     550             : "           fBy = fAy;\n"
     551             : "           fAy = fp - GetChiDist(fAx, fdf);\n"
     552             : "       }\n"
     553             : "       else\n"
     554             : "       {\n"
     555             : "           fTemp = fBx;\n"
     556             : "           fBx += 2.0 * (fBx - fAx);\n"
     557             : "           fAx = fTemp;\n"
     558             : "           fAy = fBy;\n"
     559             : "           fBy = fp - GetChiDist(fBx, fdf);\n"
     560             : "       }\n"
     561             : "   }\n"
     562             : "   if (fAy == 0.0)\n"
     563             : "       return fAx;\n"
     564             : "   if (fBy == 0.0)\n"
     565             : "       return fBx;\n"
     566             : "   if (!lcl_HasChangeOfSign( fAy, fBy))\n"
     567             : "   {\n"
     568             : "       *rConvError = true;\n"
     569             : "       return 0.0;\n"
     570             : "   }\n"
     571             : "   double fPx = fAx;\n"
     572             : "   double fPy = fAy;\n"
     573             : "   double fQx = fBx;\n"
     574             : "   double fQy = fBy;\n"
     575             : "   double fRx = fAx;\n"
     576             : "   double fRy = fAy;\n"
     577             : "   double fSx = 0.5 * (fAx + fBx);\n"
     578             : "   bool bHasToInterpolate = true;\n"
     579             : "   nCount = 0;\n"
     580             : "   while ( nCount < 500 && fabs(fRy) > fYEps &&\n"
     581             : "       (fBx-fAx) > fmax( fabs(fAx), fabs(fBx)) * fXEps )\n"
     582             : "   {\n"
     583             : "       if (bHasToInterpolate)\n"
     584             : "       {\n"
     585             : "           if (fPy!=fQy && fQy!=fRy && fRy!=fPy)\n"
     586             : "           {\n"
     587             : "               fSx = fPx * fRy * fQy*pow(fRy-fPy,-1.0)*pow(fQy-fPy,-1.0)\n"
     588             : "                   + fRx * fQy * fPy*pow(fQy-fRy,-1.0)*pow(fPy-fRy,-1.0)\n"
     589             : "                   + fQx * fPy * fRy*pow(fPy-fQy,-1.0)*pow(fRy-fQy,-1.0);\n"
     590             : "               bHasToInterpolate = (fAx < fSx) && (fSx < fBx);\n"
     591             : "           }\n"
     592             : "           else\n"
     593             : "               bHasToInterpolate = false;\n"
     594             : "       }\n"
     595             : "       if(!bHasToInterpolate)\n"
     596             : "       {\n"
     597             : "           fSx = 0.5 * (fAx + fBx);\n"
     598             : "           fPx = fAx; fPy = fAy;\n"
     599             : "           fQx = fBx; fQy = fBy;\n"
     600             : "           bHasToInterpolate = true;\n"
     601             : "       }\n"
     602             : "       fPx = fQx; fQx = fRx; fRx = fSx;\n"
     603             : "       fPy = fQy; fQy = fRy; fRy = fp - GetChiDist(fSx, fdf);\n"
     604             : "       if (lcl_HasChangeOfSign( fAy, fRy))\n"
     605             : "       {\n"
     606             : "           fBx = fRx; fBy = fRy;\n"
     607             : "       }\n"
     608             : "       else\n"
     609             : "       {\n"
     610             : "           fAx = fRx; fAy = fRy;\n"
     611             : "       }\n"
     612             : "       bHasToInterpolate = bHasToInterpolate && (fabs(fRy)"
     613             : " * 2.0 <= fabs(fQy));\n"
     614             : "       ++nCount;\n"
     615             : "   }\n"
     616             : "   return fRx;\n"
     617             : "}\n";
     618             : 
     619          76 : std::string lcl_IterateInverseChiSQInvDecl =
     620             : "static double lcl_IterateInverseChiSQInv( double fp, double fdf, \n"
     621             : "   double fAx, double fBx, bool *rConvError );\n";
     622          76 : std::string lcl_IterateInverseChiSQInv =
     623             : "static double lcl_IterateInverseChiSQInv( double fp, double fdf, \n"
     624             : "   double fAx, double fBx, bool *rConvError )\n"
     625             : "{\n"
     626             : "   *rConvError = false;\n"
     627             : "    double fYEps = 1.0E-307;\n"
     628             : "    double fXEps = fMachEps;\n"
     629             : 
     630             : "    if(!(fAx < fBx))\n"
     631             : "    {\n"
     632             : "        //print error\n"
     633             : "    }\n"
     634             : "    double fAy = fp - GetChiSqDistCDF(fAx, fdf);\n"
     635             : "    double fBy = fp - GetChiSqDistCDF(fBx, fdf);\n"
     636             : "    double fTemp;\n"
     637             : "    unsigned short nCount;\n"
     638             : "    for (nCount = 0; nCount < 1000 && !lcl_HasChangeOfSign(fAy,fBy);"
     639             : " nCount++)\n"
     640             : "    {\n"
     641             : "       if (fabs(fAy) <= fabs(fBy))\n"
     642             : "       {\n"
     643             : "           fTemp = fAx;\n"
     644             : "           fAx += 2.0 * (fAx - fBx);\n"
     645             : "           if (fAx < 0.0)\n"
     646             : "               fAx = 0.0;\n"
     647             : "           fBx = fTemp;\n"
     648             : "           fBy = fAy;\n"
     649             : "           fAy = fp - GetChiSqDistCDF(fAx, fdf);\n"
     650             : "       }\n"
     651             : "       else\n"
     652             : "       {\n"
     653             : "           fTemp = fBx;\n"
     654             : "           fBx += 2.0 * (fBx - fAx);\n"
     655             : "           fAx = fTemp;\n"
     656             : "           fAy = fBy;\n"
     657             : "           fBy = fp - GetChiSqDistCDF(fBx, fdf);\n"
     658             : "       }\n"
     659             : "   }\n"
     660             : "   if (fAy == 0.0)\n"
     661             : "       return fAx;\n"
     662             : "   if (fBy == 0.0)\n"
     663             : "       return fBx;\n"
     664             : "   if (!lcl_HasChangeOfSign( fAy, fBy))\n"
     665             : "   {\n"
     666             : "       *rConvError = true;\n"
     667             : "       return 0.0;\n"
     668             : "   }\n"
     669             : "   double fPx = fAx;\n"
     670             : "   double fPy = fAy;\n"
     671             : "   double fQx = fBx;\n"
     672             : "   double fQy = fBy;\n"
     673             : "   double fRx = fAx;\n"
     674             : "   double fRy = fAy;\n"
     675             : "   double fSx = 0.5 * (fAx + fBx);\n"
     676             : "   bool bHasToInterpolate = true;\n"
     677             : "   nCount = 0;\n"
     678             : "   while ( nCount < 500 && fabs(fRy) > fYEps &&\n"
     679             : "       (fBx-fAx) > fmax( fabs(fAx), fabs(fBx)) * fXEps )\n"
     680             : "   {\n"
     681             : "       if (bHasToInterpolate)\n"
     682             : "       {\n"
     683             : "           if (fPy!=fQy && fQy!=fRy && fRy!=fPy)\n"
     684             : "           {\n"
     685             : "               fSx = fPx * fRy * fQy / (fRy-fPy) / (fQy-fPy)\n"
     686             : "                   + fRx * fQy * fPy / (fQy-fRy) / (fPy-fRy)\n"
     687             : "                   + fQx * fPy * fRy / (fPy-fQy) / (fRy-fQy);\n"
     688             : "               bHasToInterpolate = (fAx < fSx) && (fSx < fBx);\n"
     689             : "           }\n"
     690             : "           else\n"
     691             : "               bHasToInterpolate = false;\n"
     692             : "       }\n"
     693             : "       if(!bHasToInterpolate)\n"
     694             : "       {\n"
     695             : "           fSx = 0.5 * (fAx + fBx);\n"
     696             : "           fPx = fAx; fPy = fAy;\n"
     697             : "           fQx = fBx; fQy = fBy;\n"
     698             : "           bHasToInterpolate = true;\n"
     699             : "       }\n"
     700             : "       fPx = fQx; fQx = fRx; fRx = fSx;\n"
     701             : "       fPy = fQy; fQy = fRy; fRy = fp - GetChiSqDistCDF(fSx, fdf);\n"
     702             : "       if (lcl_HasChangeOfSign( fAy, fRy))\n"
     703             : "       {\n"
     704             : "           fBx = fRx; fBy = fRy;\n"
     705             : "       }\n"
     706             : "       else\n"
     707             : "       {\n"
     708             : "           fAx = fRx; fAy = fRy;\n"
     709             : "       }\n"
     710             : "       bHasToInterpolate = bHasToInterpolate && (fabs(fRy) * 2.0"
     711             : " <= fabs(fQy));\n"
     712             : "       ++nCount;\n"
     713             : "   }\n"
     714             : "   return fRx;\n"
     715             : "}\n";
     716             : 
     717          76 : std::string gaussinvDecl = "double gaussinv(double x);\n";
     718          76 : std::string gaussinv =
     719             : "double gaussinv(double x)\n"
     720             : "{\n"
     721             : "    double q,t,z;\n"
     722             : "    q=x-0.5;\n"
     723             : "    if(fabs(q)<=.425)\n"
     724             : "    {\n"
     725             : "        t=0.180625-q*q;\n"
     726             : "        z=\n"
     727             : "            q*\n"
     728             : "            (\n"
     729             : "            (\n"
     730             : "            (\n"
     731             : "            (\n"
     732             : "            (\n"
     733             : "            (\n"
     734             : "            (\n"
     735             : "            t*2509.0809287301226727+33430.575583588128105\n"
     736             : "            )\n"
     737             : "            *t+67265.770927008700853\n"
     738             : "            )\n"
     739             : "            *t+45921.953931549871457\n"
     740             : "            )\n"
     741             : "            *t+13731.693765509461125\n"
     742             : "            )\n"
     743             : "            *t+1971.5909503065514427\n"
     744             : "            )\n"
     745             : "            *t+133.14166789178437745\n"
     746             : "            )\n"
     747             : "            *t+3.387132872796366608\n"
     748             : "            )\n"
     749             : "            *pow\n"
     750             : "            (\n"
     751             : "            (\n"
     752             : "            (\n"
     753             : "            (\n"
     754             : "            (\n"
     755             : "            (\n"
     756             : "            (\n"
     757             : "            t*5226.495278852854561+28729.085735721942674\n"
     758             : "            )\n"
     759             : "            *t+39307.89580009271061\n"
     760             : "            )\n"
     761             : "            *t+21213.794301586595867\n"
     762             : "            )\n"
     763             : "            *t+5394.1960214247511077\n"
     764             : "            )\n"
     765             : "            *t+687.1870074920579083\n"
     766             : "            )\n"
     767             : "            *t+42.313330701600911252\n"
     768             : "            )\n"
     769             : "            *t+1.0\n"
     770             : "            , -1.0);\n"
     771             : "    }\n"
     772             : "    else\n"
     773             : "    {\n"
     774             : "        if(q>0) t=1-x;\n"
     775             : "        else        t=x;\n"
     776             : "        t=sqrt(-log(t));\n"
     777             : "        if(t<=5.0)\n"
     778             : "        {\n"
     779             : "            t+=-1.6;\n"
     780             : "            z=\n"
     781             : "                (\n"
     782             : "                (\n"
     783             : "                (\n"
     784             : "                (\n"
     785             : "                (\n"
     786             : "                (\n"
     787             : "                (\n"
     788             : "                t*7.7454501427834140764e-4+0.0227238449892691845833\n"
     789             : "                )\n"
     790             : "                *t+0.24178072517745061177\n"
     791             : "                )\n"
     792             : "                *t+1.27045825245236838258\n"
     793             : "                )\n"
     794             : "                *t+3.64784832476320460504\n"
     795             : "                )\n"
     796             : "                *t+5.7694972214606914055\n"
     797             : "                )\n"
     798             : "                *t+4.6303378461565452959\n"
     799             : "                )\n"
     800             : "                *t+1.42343711074968357734\n"
     801             : "                )\n"
     802             : "                *pow\n"
     803             : "                (\n"
     804             : "                (\n"
     805             : "                (\n"
     806             : "                (\n"
     807             : "                (\n"
     808             : "                (\n"
     809             : "                (\n"
     810             : "                t*1.05075007164441684324e-9+5.475938084995344946e-4\n"
     811             : "                )\n"
     812             : "                *t+0.0151986665636164571966\n"
     813             : "                )\n"
     814             : "                *t+0.14810397642748007459\n"
     815             : "                )\n"
     816             : "                *t+0.68976733498510000455\n"
     817             : "                )\n"
     818             : "                *t+1.6763848301838038494\n"
     819             : "                )\n"
     820             : "                *t+2.05319162663775882187\n"
     821             : "                )\n"
     822             : "                *t+1.0\n"
     823             : "                , -1.0);\n"
     824             : "        }\n"
     825             : "        else\n"
     826             : "        {\n"
     827             : "            t+=-5.0;\n"
     828             : "            z=\n"
     829             : "                (\n"
     830             : "                (\n"
     831             : "                (\n"
     832             : "                (\n"
     833             : "                (\n"
     834             : "                (\n"
     835             : "                (\n"
     836             : "                t*2.01033439929228813265e-7+2.71155556874348757815e-5\n"
     837             : "                )\n"
     838             : "                *t+0.0012426609473880784386\n"
     839             : "                )\n"
     840             : "                *t+0.026532189526576123093\n"
     841             : "                )\n"
     842             : "                *t+0.29656057182850489123\n"
     843             : "                )\n"
     844             : "                *t+1.7848265399172913358\n"
     845             : "                )\n"
     846             : "                *t+5.4637849111641143699\n"
     847             : "                )\n"
     848             : "                *t+6.6579046435011037772\n"
     849             : "                )\n"
     850             : "                *pow\n"
     851             : "                (\n"
     852             : "                (\n"
     853             : "                (\n"
     854             : "                (\n"
     855             : "                (\n"
     856             : "                (\n"
     857             : "                (\n"
     858             : "                t*2.04426310338993978564e-15+1.4215117583164458887e-7\n"
     859             : "                )\n"
     860             : "                *t+1.8463183175100546818e-5\n"
     861             : "                )\n"
     862             : "                *t+7.868691311456132591e-4\n"
     863             : "                )\n"
     864             : "                *t+0.0148753612908506148525\n"
     865             : "                )\n"
     866             : "                *t+0.13692988092273580531\n"
     867             : "                )\n"
     868             : "                *t+0.59983220655588793769\n"
     869             : "                )\n"
     870             : "                *t+1.0\n"
     871             : "                , -1.0);\n"
     872             : "        }\n"
     873             : "        if(q<0.0) z=-z;\n"
     874             : "    }\n"
     875             : "    return z;\n"
     876             : "}\n";
     877             : 
     878          76 : std::string lcl_GetLogGammaHelperDecl=
     879             : "static double lcl_GetLogGammaHelper(double fZ);\n";
     880          76 : std::string lcl_GetLogGammaHelper =
     881             : "static double lcl_GetLogGammaHelper(double fZ)\n"
     882             : "{\n"
     883             : "    double fg = 6.024680040776729583740234375;\n"
     884             : "   double fZgHelp = fZ + fg - 0.5;\n"
     885             : "   return log( lcl_getLanczosSum(fZ)) + (fZ-0.5) * log(fZgHelp) - fZgHelp;\n"
     886             : "}\n";
     887          76 : std::string lcl_GetGammaHelperDecl=
     888             : "static double lcl_GetGammaHelper(double fZ);\n";
     889          76 : std::string lcl_GetGammaHelper =
     890             : "static double lcl_GetGammaHelper(double fZ)\n"
     891             : "{\n"
     892             : "   double fGamma = lcl_getLanczosSum(fZ);\n"
     893             : "   double fg = 6.024680040776729583740234375;\n"
     894             : "   double fZgHelp = fZ + fg - 0.5;\n"
     895             : "   double fHalfpower = pow( fZgHelp, fZ*pow(2,-1.0) - 0.25);\n"
     896             : "   fGamma *= fHalfpower;\n"
     897             : "   fGamma = fGamma*pow(exp(fZgHelp),-1.0);\n"
     898             : "   fGamma *= fHalfpower;\n"
     899             : "   fGamma = 120.4;\n"
     900             : "   if (fZ <= 20.0 && fZ == (int)fZ)\n"
     901             : "   {\n"
     902             : "     fGamma = (int)(fGamma+0.5);\n"
     903             : "   }\n"
     904             : "   return fGamma;\n"
     905             : "}\n";
     906          76 : std::string lcl_getLanczosSumDecl=
     907             : "static double lcl_getLanczosSum(double fZ);\n";
     908          76 : std::string lcl_getLanczosSum =
     909             : "static double lcl_getLanczosSum(double fZ)          \n"
     910             : "{                                                   \n"
     911             : "    double fNum[13] ={                        \n"
     912             : "        23531376880.41075968857200767445163675473,  \n"
     913             : "        42919803642.64909876895789904700198885093,  \n"
     914             : "        35711959237.35566804944018545154716670596,  \n"
     915             : "        17921034426.03720969991975575445893111267,  \n"
     916             : "        6039542586.35202800506429164430729792107,   \n"
     917             : "        1439720407.311721673663223072794912393972,  \n"
     918             : "        248874557.8620541565114603864132294232163,  \n"
     919             : "        31426415.58540019438061423162831820536287,  \n"
     920             : "        2876370.628935372441225409051620849613599,  \n"
     921             : "        186056.2653952234950402949897160456992822,  \n"
     922             : "        8071.672002365816210638002902272250613822,  \n"
     923             : "        210.8242777515793458725097339207133627117,  \n"
     924             : "        2.506628274631000270164908177133837338626   \n"
     925             : "        };                                          \n"
     926             : "    double fDenom[13] = {                     \n"
     927             : "        0,\n"
     928             : "        39916800,\n"
     929             : "        120543840,\n"
     930             : "        150917976,\n"
     931             : "        105258076,\n"
     932             : "        45995730,\n"
     933             : "        13339535,\n"
     934             : "        2637558,\n"
     935             : "        357423,\n"
     936             : "        32670,\n"
     937             : "        1925,\n"
     938             : "        66,\n"
     939             : "        1\n"
     940             : "        };\n"
     941             : "    double fSumNum;\n"
     942             : "    double fSumDenom;\n"
     943             : "    int nI;\n"
     944             : "    if (fZ<=1.0)\n"
     945             : "    {\n"
     946             : "        fSumNum = fNum[12];\n"
     947             : "        fSumDenom = fDenom[12];\n"
     948             : "        nI = 11;\n"
     949             : "        fSumNum = fSumNum*fZ+fNum[nI];\n"
     950             : "        fSumDenom = fSumDenom*fZ+fDenom[nI];\n"
     951             : "        nI = 10;\n"
     952             : "        fSumNum = fSumNum*fZ+fNum[nI];\n"
     953             : "        fSumDenom = fSumDenom*fZ+fDenom[nI];\n"
     954             : "        nI = 9;\n"
     955             : "        fSumNum = fSumNum*fZ+fNum[nI];\n"
     956             : "        fSumDenom = fSumDenom*fZ+fDenom[nI];\n"
     957             : "        nI = 8;\n"
     958             : "        fSumNum = fSumNum*fZ+fNum[nI];\n"
     959             : "        fSumDenom = fSumDenom*fZ+fDenom[nI];\n"
     960             : "        nI = 7;\n"
     961             : "        fSumNum = fSumNum*fZ+fNum[nI];\n"
     962             : "        fSumDenom = fSumDenom*fZ+fDenom[nI];\n"
     963             : "        nI = 6;\n"
     964             : "        fSumNum = fSumNum*fZ+fNum[nI];\n"
     965             : "        fSumDenom = fSumDenom*fZ+fDenom[nI];\n"
     966             : "        nI = 5;\n"
     967             : "        fSumNum = fSumNum*fZ+fNum[nI];\n"
     968             : "        fSumDenom = fSumDenom*fZ+fDenom[nI];\n"
     969             : "        nI = 4;\n"
     970             : "        fSumNum = fSumNum*fZ+fNum[nI];\n"
     971             : "        fSumDenom = fSumDenom*fZ+fDenom[nI];\n"
     972             : "        nI = 3;\n"
     973             : "        fSumNum = fSumNum*fZ+fNum[nI];\n"
     974             : "        fSumDenom = fSumDenom*fZ+fDenom[nI];\n"
     975             : "        nI = 2;\n"
     976             : "        fSumNum = fSumNum*fZ+fNum[nI];\n"
     977             : "        fSumDenom = fSumDenom*fZ+fDenom[nI];\n"
     978             : "        nI = 1;\n"
     979             : "        fSumNum = fSumNum*fZ+fNum[nI];\n"
     980             : "        fSumDenom = fSumDenom*fZ+fDenom[nI];\n"
     981             : "        nI = 0;\n"
     982             : "        fSumNum = fSumNum*fZ+fNum[nI];\n"
     983             : "        fSumDenom = fSumDenom*fZ+fDenom[nI];\n"
     984             : "    }\n"
     985             : "    if (fZ>1.0)\n"
     986             : "    {\n"
     987             : "        double fZInv = pow(fZ,-1.0);\n"
     988             : "        fSumNum = fNum[0];\n"
     989             : "        fSumDenom = fDenom[0];\n"
     990             : "        nI = 1;\n"
     991             : "        fSumNum = fSumNum*fZInv+fNum[nI];\n"
     992             : "        fSumDenom = fSumDenom*fZInv+fDenom[nI];\n"
     993             : "        nI = 2;\n"
     994             : "        fSumNum = fSumNum*fZInv+fNum[nI];\n"
     995             : "        fSumDenom = fSumDenom*fZInv+fDenom[nI];\n"
     996             : "        nI = 3;\n"
     997             : "        fSumNum = fSumNum*fZInv+fNum[nI];\n"
     998             : "        fSumDenom = fSumDenom*fZInv+fDenom[nI];\n"
     999             : "        nI = 4;\n"
    1000             : "        fSumNum = fSumNum*fZInv+fNum[nI];\n"
    1001             : "        fSumDenom = fSumDenom*fZInv+fDenom[nI];\n"
    1002             : "        nI = 5;\n"
    1003             : "        fSumNum = fSumNum*fZInv+fNum[nI];\n"
    1004             : "        fSumDenom = fSumDenom*fZInv+fDenom[nI];\n"
    1005             : "        nI = 6;\n"
    1006             : "        fSumNum = fSumNum*fZInv+fNum[nI];\n"
    1007             : "        fSumDenom = fSumDenom*fZInv+fDenom[nI];\n"
    1008             : "        nI = 7;\n"
    1009             : "        fSumNum = fSumNum*fZInv+fNum[nI];\n"
    1010             : "        fSumDenom = fSumDenom*fZInv+fDenom[nI];\n"
    1011             : "        nI = 8;\n"
    1012             : "        fSumNum = fSumNum*fZInv+fNum[nI];\n"
    1013             : "        fSumDenom = fSumDenom*fZInv+fDenom[nI];\n"
    1014             : "        nI = 9;\n"
    1015             : "        fSumNum = fSumNum*fZInv+fNum[nI];\n"
    1016             : "        fSumDenom = fSumDenom*fZInv+fDenom[nI];\n"
    1017             : "        nI = 10;\n"
    1018             : "        fSumNum = fSumNum*fZInv+fNum[nI];\n"
    1019             : "        fSumDenom = fSumDenom*fZInv+fDenom[nI];\n"
    1020             : "        nI = 11;\n"
    1021             : "        fSumNum = fSumNum*fZInv+fNum[nI];\n"
    1022             : "        fSumDenom = fSumDenom*fZInv+fDenom[nI];\n"
    1023             : "        nI = 12;\n"
    1024             : "        fSumNum = fSumNum*fZInv+fNum[nI];\n"
    1025             : "        fSumDenom = fSumDenom*fZInv+fDenom[nI];\n"
    1026             : "     }\n"
    1027             : "     return fSumNum*pow(fSumDenom,-1.0);\n"
    1028             : "}\n";
    1029             : 
    1030          76 : std::string GetUpRegIGammaDecl=
    1031             : " double GetUpRegIGamma( double fA, double fX ) ;\n";
    1032          76 : std::string GetUpRegIGamma =
    1033             : "double GetUpRegIGamma( double fA, double fX )\n"
    1034             : "{\n"
    1035             : "    double fLnFactor= fA*log(fX)-fX-lgamma(fA);\n"
    1036             : "    double fFactor = exp(fLnFactor); \n"
    1037             : "    if (fX>fA+1.0) \n"
    1038             : "            return fFactor * GetGammaContFraction(fA,fX);\n"
    1039             : "    else \n"
    1040             : "            return 1.0 -fFactor * GetGammaSeries(fA,fX);\n"
    1041             : "}\n";
    1042             : 
    1043          76 : std::string lcl_HasChangeOfSignDecl=
    1044             : "static inline bool lcl_HasChangeOfSign( double u, double w );\n";
    1045          76 : std::string lcl_HasChangeOfSign =
    1046             : "static inline bool lcl_HasChangeOfSign( double u, double w )\n"
    1047             : "{\n"
    1048             : "    return (u < 0.0 && w > 0.0) || (u > 0.0 && w < 0.0);\n"
    1049             : "}\n";
    1050             : 
    1051          76 : std::string GetTDistDecl=" double GetTDist(double T, double fDF);\n";
    1052          76 : std::string GetTDist =
    1053             : "double GetTDist(double T, double fDF)\n"
    1054             : "{\n"
    1055             : "    return 0.5 * GetBetaDist(fDF*pow(fDF+T*T,-1.0),fDF*pow(2.0,-1.0), 0.5);\n"
    1056             : "}\n";
    1057             : 
    1058          76 : std::string GetBetaDecl=" double GetBeta(double fAlpha, double fBeta);\n";
    1059          76 : std::string GetBeta =
    1060             : "double GetBeta(double fAlpha, double fBeta)\n"
    1061             : "{\n"
    1062             : "    double fA;\n"
    1063             : "    double fB;\n"
    1064             : "    fAlpha>fBeta?(fA = fAlpha,fB = fBeta):(fA = fBeta,fB = fAlpha);\n"
    1065             : "    double fAB = fA + fB;\n"
    1066             : 
    1067             : "    if (fAB < fMaxGammaArgument)\n"
    1068             : "        return tgamma(fA)*pow(tgamma(fAB),-1.0)*tgamma(fB);\n"
    1069             : "    double fgm = 5.524680040776729583740234375;\n"
    1070             : "    double fLanczos = lcl_getLanczosSum(fA)*lcl_getLanczosSum(fB)\n"
    1071             : "        *pow(lcl_getLanczosSum(fAB),-1.0);\n"
    1072             : "    fLanczos *= sqrt(((fAB + fgm)*pow(fA + fgm,-1.0))*pow(fB + fgm,-1.0));\n"
    1073             : "    return fLanczos * pow(exp(1.0),(-fA*log1p(fB*pow(fA + fgm,-1.0)))"
    1074             : "                    - fB*log1p(fA*pow(fB + fgm,-1.0)) - fgm);\n"
    1075             : "}\n";
    1076             : 
    1077          76 : std::string GetLogBetaDecl=
    1078             : " double GetLogBeta(double fAlpha, double fBeta);\n";
    1079          76 : std::string GetLogBeta =
    1080             : "double GetLogBeta(double fAlpha, double fBeta)\n"
    1081             : "{\n"
    1082             : "    double fA;\n"
    1083             : "    double fB;\n"
    1084             : "    fAlpha>fBeta?(fA = fAlpha,fB = fBeta):(fA = fBeta,fB = fAlpha);\n"
    1085             : "    double fgm = 5.524680040776729583740234375;\n"
    1086             : 
    1087             : "    double fLanczos = lcl_getLanczosSum(fA)*lcl_getLanczosSum(fB)*\n"
    1088             : "        pow(lcl_getLanczosSum(fA + fB),-1.0);\n"
    1089             : "    double fResult= -fA *log1p(fB*pow(fA + fgm,-1.0))"
    1090             : "-fB *log1p(fA*pow(fB + fgm,-1.0))-fgm;\n"
    1091             : "    fResult += log(fLanczos)+0.5*(log(fA + fB + fgm) - log(fA + fgm)\n"
    1092             : "        - log(fB + fgm));\n"
    1093             : "    return fResult;\n"
    1094             : "}\n";
    1095             : 
    1096          76 : std::string GetBetaDistPDFDecl=
    1097             : "double GetBetaDistPDF(double fX, double fA, double fB);\n";
    1098          76 : std::string GetBetaDistPDF =
    1099             : "double GetBetaDistPDF(double fX, double fA, double fB)\n"
    1100             : "{\n"
    1101             : "    if (fA == 1.0) \n"
    1102             : "    {\n"
    1103             : "        if (fB == 1.0)\n"
    1104             : "            return 1.0;\n"
    1105             : "        if (fB == 2.0)\n"
    1106             : "            return -2.0*fX + 2.0;\n"
    1107             : "        if (fX == 1.0 && fB < 1.0)\n"
    1108             : "        {\n"
    1109             : "            return HUGE_VAL;\n"
    1110             : "        }\n"
    1111             : "        if (fX <= 0.01)\n"
    1112             : "            return fB + fB * expm1((fB-1.0) * log1p(-fX));\n"
    1113             : "        else \n"
    1114             : "            return fB * pow(0.5-fX+0.5,fB-1.0);\n"
    1115             : "    }\n"
    1116             : "    if (fB == 1.0) \n"
    1117             : "    {\n"
    1118             : "    if (fA == 2.0)\n"
    1119             : "        return fA * fX;\n"
    1120             : "        if (fX == 0.0 && fA < 1.0)\n"
    1121             : "        {\n"
    1122             : "            return HUGE_VAL;\n"
    1123             : "        }\n"
    1124             : "        return fA * pow(fX,fA-1);\n"
    1125             : "    }\n"
    1126             : "    if (fX <= 0.0)\n"
    1127             : "    {\n"
    1128             : "        if (fA < 1.0 && fX == 0.0)\n"
    1129             : "        {\n"
    1130             : "            return HUGE_VAL;\n"
    1131             : "        }\n"
    1132             : "        else\n"
    1133             : "            return 0.0;\n"
    1134             : "    }\n"
    1135             : "    if (fX >= 1.0)\n"
    1136             : "    {\n"
    1137             : "        if (fB < 1.0 && fX == 1.0)\n"
    1138             : "        {\n"
    1139             : "            return HUGE_VAL;\n"
    1140             : "        }\n"
    1141             : "        else \n"
    1142             : "        return 0.0;\n"
    1143             : "    }\n"
    1144             : "    double fLogDblMax = log( 1.79769e+308 );\n"
    1145             : "    double fLogDblMin = log( 2.22507e-308 );\n"
    1146             : "    double fLogY = (fX < 0.1) ? log1p(-fX) : log(0.5-fX+0.5);\n"
    1147             : "    double fLogX = log(fX);\n"
    1148             : "    double fAm1LogX = (fA-1.0) * fLogX;\n"
    1149             : "    double fBm1LogY = (fB-1.0) * fLogY;\n"
    1150             : "    double fLogBeta = GetLogBeta(fA,fB);\n"
    1151             : "    if (   fAm1LogX < fLogDblMax  && fAm1LogX > fLogDblMin\n"
    1152             : "        && fBm1LogY < fLogDblMax  && fBm1LogY > fLogDblMin\n"
    1153             : "        && fLogBeta < fLogDblMax  && fLogBeta > fLogDblMin\n"
    1154             : "        && fAm1LogX + fBm1LogY < fLogDblMax && fAm1LogX + fBm1LogY > \n"
    1155             : "           fLogDblMin)\n"
    1156             : "        return pow(fX,fA-1.0)*pow(0.5-fX+0.5,fB-1.0)"
    1157             : "*pow(GetBeta(fA,fB),-1.0);\n"
    1158             : "    else \n"
    1159             : "         return exp( fAm1LogX + fBm1LogY - fLogBeta);\n"
    1160             : "}\n";
    1161             : 
    1162          76 : std::string lcl_GetBetaHelperContFracDecl=
    1163             : "double lcl_GetBetaHelperContFrac(double fX, double fA, double fB);\n";
    1164          76 : std::string lcl_GetBetaHelperContFrac =
    1165             : "double lcl_GetBetaHelperContFrac(double fX, double fA, double fB)\n"
    1166             : "{   \n"
    1167             : 
    1168             : "    double a1, b1, a2, b2, fnorm, apl2m, d2m, d2m1, cfnew, cf;\n"
    1169             : "    a1 = 1.0; b1 = 1.0;\n"
    1170             : "    b2 = 1.0 - (fA+fB)*pow(fA+1.0,-1.0)*fX;\n"
    1171             : "    b2==0.0?(a2 = 0.0,fnorm = 1.0,cf = 1.0):\n"
    1172             : "        (a2 = 1.0,fnorm = pow(b2,-1.0),cf = a2*fnorm);\n"
    1173             : "    cfnew = 1.0;\n"
    1174             : "    double rm = 1.0;\n"
    1175             : "    double fMaxIter = 50000.0;\n"
    1176             : "    bool bfinished = false;\n"
    1177             : "    do\n"
    1178             : "    {\n"
    1179             : "      apl2m = fA + 2.0*rm;\n"
    1180             : "      d2m = (rm*(fB-rm))*fX*pow(apl2m*(apl2m-1.0),-1.0);\n"
    1181             : "      d2m1 = -((fA+rm)*(fA+rm+fB))*fX*pow(apl2m*(apl2m+1.0),-1.0);\n"
    1182             : "      a1 = (a2+d2m*a1)*fnorm;\n"
    1183             : "      b1 = (b2+d2m*b1)*fnorm;\n"
    1184             : "      a2 = a1 + d2m1*a2*fnorm;\n"
    1185             : "      b2 = b1 + d2m1*b2*fnorm;\n"
    1186             : "      if (b2 != 0.0) \n"
    1187             : "      {\n"
    1188             : "        fnorm = pow(b2,-1.0);\n"
    1189             : "        cfnew = a2*fnorm;\n"
    1190             : "        bfinished = (fabs(cf-cfnew) < fabs(cf)*fMachEps);\n"
    1191             : "      }\n"
    1192             : "      cf = cfnew;\n"
    1193             : "      rm += 1.0;\n"
    1194             : "     }\n"
    1195             : "    while (rm < fMaxIter && !bfinished);\n"
    1196             : "    return cf;\n"
    1197             : "}\n";
    1198             : 
    1199          76 : std::string lcl_IterateInverseDecl=
    1200             : "double lcl_IterateInverse("
    1201             : "double fAx, double fBx, bool* rConvError,double fp,double fDF );\n";
    1202          76 : std::string lcl_IterateInverse =
    1203             : "double lcl_IterateInverse( "
    1204             : "double fAx, double fBx, bool* rConvError,double fp,double fDF )\n"
    1205             : "{\n"
    1206             : "    *rConvError = false;\n"
    1207             : "    double fYEps = 1.0E-307;\n"
    1208             : "    double fXEps =DBL_EPSILON;\n"
    1209             : "    if(fAx>fBx)\n"
    1210             : "      return DBL_MAX;\n"
    1211             : "    double fAy = GetValue(fAx,fp,fDF);\n"
    1212             : "    double fBy = GetValue(fBx,fp,fDF);\n"
    1213             : "    double fTemp;\n"
    1214             : "    unsigned short nCount;\n"
    1215             : "    double inter;\n"
    1216             : "    bool sign;\n"
    1217             : "    for (nCount =0;nCount<1000&&!lcl_HasChangeOfSign(fAy,fBy);nCount++)\n"
    1218             : "    {\n"
    1219             : "        inter = 2.0 * (fAx - fBx);\n"
    1220             : "        if (fabs(fAy) <= fabs(fBy)) \n"
    1221             : "        {\n"
    1222             : "            sign = true;\n"
    1223             : "            fTemp = fAx;\n"
    1224             : "            fAx += inter;\n"
    1225             : "            if (fAx < 0.0)\n"
    1226             : "                fAx = 0.0;\n"
    1227             : "            fBx = fTemp;\n"
    1228             : "            fBy = fAy;\n"
    1229             : "            fTemp = fAx;\n"
    1230             : "        }\n"
    1231             : "        else\n"
    1232             : "        {\n"
    1233             : "            sign = false;\n"
    1234             : "            fTemp = fBx;\n"
    1235             : "            fBx -= inter;\n"
    1236             : "            fAx = fTemp;\n"
    1237             : "            fAy = fBy;\n"
    1238             : "            fTemp = fBx;\n"
    1239             : "        }\n"
    1240             : "        fTemp = GetValue(fTemp,fp,fDF);\n"
    1241             : "        sign?(fAy = fTemp):(fBy = fTemp);\n"
    1242             : "    }\n"
    1243             : "    if (fAy == 0.0)\n"
    1244             : "        return fAx;\n"
    1245             : "    if (fBy == 0.0)\n"
    1246             : "        return fBx;\n"
    1247             : "    if (!lcl_HasChangeOfSign( fAy, fBy))\n"
    1248             : "    {\n"
    1249             : "        *rConvError = true;\n"
    1250             : "        return 0.0;\n"
    1251             : "    }\n"
    1252             : "    double fPx = fAx;\n"
    1253             : "    double fPy = fAy;\n"
    1254             : "    double fQx = fBx;\n"
    1255             : "    double fQy = fBy;\n"
    1256             : "    double fRx = fAx;\n"
    1257             : "    double fRy = fAy;\n"
    1258             : "    double fSx = 0.5 * (fAx + fBx); \n"
    1259             : "    bool bHasToInterpolate = true;\n"
    1260             : "    nCount = 0;\n"
    1261             : "    while ( nCount < 500 && fabs(fRy) > fYEps &&\n"
    1262             : "            (fBx-fAx) > max( fabs(fAx), fabs(fBx)) * fXEps )\n"
    1263             : "    {\n"
    1264             : "        if (bHasToInterpolate)\n"
    1265             : "        {\n"
    1266             : "           if (fPy!=fQy && fQy!=fRy && fRy!=fPy)\n"
    1267             : "           {\n"
    1268             : "               fSx = fPx * fRy * fQy * pow(fRy-fPy,-1.0)*pow(fQy-fPy,-1.0)\n"
    1269             : "                   + fRx * fQy * fPy * pow(fQy-fRy,-1.0)*pow(fPy-fRy,-1.0)\n"
    1270             : "                   + fQx * fPy * fRy * pow(fPy-fQy,-1.0)*pow(fRy-fQy,-1.0);\n"
    1271             : "               bHasToInterpolate = (fAx < fSx) && (fSx < fBx);\n"
    1272             : "           }\n"
    1273             : "            else\n"
    1274             : "                bHasToInterpolate = false;\n"
    1275             : "        }\n"
    1276             : "        if(!bHasToInterpolate)\n"
    1277             : "        {\n"
    1278             : "            fSx = 0.5 * (fAx + fBx);\n"
    1279             : "            \n"
    1280             : "            fPx = fAx; fPy = fAy;\n"
    1281             : "            fQx = fBx; fQy = fBy;\n"
    1282             : "            bHasToInterpolate = true;\n"
    1283             : "        }\n"
    1284             : "        fPx = fQx; fQx = fRx; fRx = fSx;\n"
    1285             : "        fPy = fQy; fQy = fRy; fRy = GetValue(fSx,fp,fDF);\n"
    1286             : "        lcl_HasChangeOfSign( fAy, fRy)?(fBx = fRx,fBy = fRy):\n"
    1287             : "            (fAx = fRx,fAy = fRy);\n"
    1288             : "        bHasToInterpolate =\n"
    1289             : "            bHasToInterpolate && (fabs(fRy) * 2.0 <= fabs(fQy));\n"
    1290             : "        ++nCount;\n"
    1291             : "    }\n"
    1292             : "    return fRx;\n"
    1293             : "}\n";
    1294          76 : std::string phiDecl=
    1295             : "double phi(double x);\n";
    1296          76 : std::string phi =
    1297             : "double phi(double x)\n"
    1298             : "{\n"
    1299             : "    return  0.39894228040143268 * exp(-(x * x) / 2.0);\n"
    1300             : "}\n";
    1301          76 : std::string taylorDecl =
    1302             : "double taylor(double* pPolynom, uint nMax, double x);\n";
    1303          76 : std::string taylor =
    1304             : "double taylor(double* pPolynom, uint nMax, double x)\n"
    1305             : "{\n"
    1306             : "    double nVal = pPolynom[nMax];\n"
    1307             : "    for (short i = nMax-1; i >= 0; i--)\n"
    1308             : "    {\n"
    1309             : "        nVal = pPolynom[i] + (nVal * x);\n"
    1310             : "    }\n"
    1311             : "    return nVal;\n"
    1312             : "}";
    1313          76 : std::string gaussDecl = "double gauss(double x);\n";
    1314          76 : std::string gauss =
    1315             : "double gauss(double x)\n"
    1316             : "{\n"
    1317             : "    double xAbs = fabs(x);\n"
    1318             : "    uint xShort = (uint)(floor(xAbs));\n"
    1319             : "    double nVal = 0.0;\n"
    1320             : "    if (xShort == 0)\n"
    1321             : "    {\n"
    1322             : "        double t0[] =\n"
    1323             : "        { 0.39894228040143268, -0.06649038006690545,  0.00997355701003582,\n"
    1324             : "         -0.00118732821548045,  0.00011543468761616, -0.00000944465625950,\n"
    1325             : "          0.00000066596935163, -0.00000004122667415,  0.00000000227352982,\n"
    1326             : "          0.00000000011301172,  0.00000000000511243, -0.00000000000021218 };\n"
    1327             : "        nVal = taylor(t0, 11, (xAbs * xAbs)) * xAbs;\n"
    1328             : "    }\n"
    1329             : "    else if ((xShort >= 1) && (xShort <= 2))\n"
    1330             : "    {\n"
    1331             : "        double t2[] =\n"
    1332             : "        { 0.47724986805182079,  0.05399096651318805, -0.05399096651318805,\n"
    1333             : "          0.02699548325659403, -0.00449924720943234, -0.00224962360471617,\n"
    1334             : "          0.00134977416282970, -0.00011783742691370, -0.00011515930357476,\n"
    1335             : "          0.00003704737285544,  0.00000282690796889, -0.00000354513195524,\n"
    1336             : "          0.00000037669563126,  0.00000019202407921, -0.00000005226908590,\n"
    1337             : "         -0.00000000491799345,  0.00000000366377919, -0.00000000015981997,\n"
    1338             : "         -0.00000000017381238,  0.00000000002624031,  0.00000000000560919,\n"
    1339             : "         -0.00000000000172127, -0.00000000000008634, 0.00000000000007894 };\n"
    1340             : "        nVal = taylor(t2, 23, (xAbs - 2.0));\n"
    1341             : "    }\n"
    1342             : "    else if ((xShort >= 3) && (xShort <= 4))\n"
    1343             : "    {\n"
    1344             : "       double t4[] =\n"
    1345             : "       { 0.49996832875816688,  0.00013383022576489, -0.00026766045152977,\n"
    1346             : "         0.00033457556441221, -0.00028996548915725,  0.00018178605666397,\n"
    1347             : "        -0.00008252863922168,  0.00002551802519049, -0.00000391665839292,\n"
    1348             : "        -0.00000074018205222,  0.00000064422023359, -0.00000017370155340,\n"
    1349             : "         0.00000000909595465,  0.00000000944943118, -0.00000000329957075,\n"
    1350             : "         0.00000000029492075,  0.00000000011874477, -0.00000000004420396,\n"
    1351             : "         0.00000000000361422,  0.00000000000143638, -0.00000000000045848 };\n"
    1352             : "        nVal = taylor(t4, 20, (xAbs - 4.0));\n"
    1353             : "    }\n"
    1354             : "    else\n"
    1355             : "    {\n"
    1356             : "        double asympt[] = { -1.0, 1.0, -3.0, 15.0, -105.0 };\n"
    1357             : "        nVal = 0.5 + phi(xAbs) * taylor(asympt, 4, 1.0/(xAbs * xAbs))/xAbs;\n"
    1358             : "    }\n"
    1359             : "    if (x < 0.0)\n"
    1360             : "        return -nVal;\n"
    1361             : "    else\n"
    1362             : "        return nVal;\n"
    1363             : "}\n";
    1364          76 : std::string lcl_Erfc0600Decl=
    1365             : "void lcl_Erfc0600( double x, double *fVal );\n";
    1366          76 : std::string lcl_Erfc0600 =
    1367             : "void lcl_Erfc0600( double  x, double *fVal )\n"
    1368             : "{\n"
    1369             : "    double fPSum = 0.0;\n"
    1370             : "    double fQSum = 0.0;\n"
    1371             : "    double fXPow = 1.0;\n"
    1372             : "    double *pn;\n"
    1373             : "    double *qn;\n"
    1374             : "    if ( x < 2.2 )\n"
    1375             : "    {\n"
    1376             : "        double pn22[] = {         \n"
    1377             : "        9.99999992049799098E-1,   \n"
    1378             : "        1.33154163936765307,      \n"
    1379             : "        8.78115804155881782E-1,   \n"
    1380             : "        3.31899559578213215E-1,   \n"
    1381             : "        7.14193832506776067E-2,   \n"
    1382             : "        7.06940843763253131E-3    \n"
    1383             : "      };                          \n"
    1384             : "        double qn22[] = {         \n"
    1385             : "        1.00000000000000000,      \n"
    1386             : "        2.45992070144245533,      \n"
    1387             : "        2.65383972869775752,      \n"
    1388             : "        1.61876655543871376,      \n"
    1389             : "        5.94651311286481502E-1,   \n"
    1390             : "        1.26579413030177940E-1,   \n"
    1391             : "        1.25304936549413393E-2    \n"
    1392             : "       };                         \n"
    1393             : "        pn = pn22;                \n"
    1394             : "        qn = qn22;                \n"
    1395             : "      }                           \n"
    1396             : "      else                        \n"
    1397             : "                                  \n"
    1398             : "       {                          \n"
    1399             : "            double pn60[] = {\n"
    1400             : "            9.99921140009714409E-1,\n"
    1401             : "            1.62356584489366647,\n"
    1402             : "            1.26739901455873222,\n"
    1403             : "            5.81528574177741135E-1,\n"
    1404             : "            1.57289620742838702E-1,\n"
    1405             : "            2.25716982919217555E-2\n"
    1406             : "            };\n"
    1407             : "            double qn60[] = {\n"
    1408             : "            1.00000000000000000,\n"
    1409             : "            2.75143870676376208,\n"
    1410             : "            3.37367334657284535,\n"
    1411             : "            2.38574194785344389,\n"
    1412             : "            1.05074004614827206,\n"
    1413             : "            2.78788439273628983E-1,\n"
    1414             : "            4.00072964526861362E-2\n"
    1415             : "            };\n"
    1416             : "            pn = pn60;\n"
    1417             : "            qn = qn60;\n"
    1418             : "       }\n"
    1419             : "    for ( unsigned int i = 0; i < 6; ++i )  \n"
    1420             : "    {\n"
    1421             : "        fPSum += pn[i]*fXPow;\n"
    1422             : "        fQSum += qn[i]*fXPow;\n"
    1423             : "        fXPow *= x;\n"
    1424             : "    }\n"
    1425             : "    fQSum += qn[6]*fXPow;\n"
    1426             : "    *fVal = exp((-1.0)*x*x)*fPSum*pow(fQSum, -1.0);\n"
    1427             : "   }\n";
    1428          76 : std::string lcl_Erfc2654Decl=
    1429             : "void lcl_Erfc2654( double x, double *fVal );\n";
    1430          76 : std::string lcl_Erfc2654 =
    1431             : "void lcl_Erfc2654( double x, double *fVal )\n"
    1432             : "{\n"
    1433             : "    double pn[] = {\n"
    1434             : "        5.64189583547756078E-1,\n"
    1435             : "        8.80253746105525775,\n"
    1436             : "        3.84683103716117320E1,\n"
    1437             : "        4.77209965874436377E1,\n"
    1438             : "        8.08040729052301677\n"
    1439             : "    };\n"
    1440             : "    double qn[] = {\n"
    1441             : "        1.00000000000000000,\n"
    1442             : "        1.61020914205869003E1,\n"
    1443             : "        7.54843505665954743E1,\n"
    1444             : "        1.12123870801026015E2,\n"
    1445             : "        3.73997570145040850E1\n"
    1446             : "    };\n"
    1447             : "\n"
    1448             : "    double fPSum = 0.0;\n"
    1449             : "    double fQSum = 0.0;\n"
    1450             : "    double fXPow = 1.0;\n"
    1451             : "\n"
    1452             : "    for ( unsigned int i = 0; i <= 4; ++i )\n"
    1453             : "    {\n"
    1454             : "        fPSum += pn[i]*fXPow;       \n"
    1455             : "        fQSum += qn[i]*fXPow;\n"
    1456             : "        fXPow *= pow(x*x, -1.0);\n"
    1457             : "    }\n"
    1458             : "    *fVal = exp((-1.0)*x*x)*fPSum*pow(x*fQSum, -1.0);\n"
    1459             : "}\n";
    1460          76 : std::string lcl_Erf0065Decl=
    1461             : "void lcl_Erf0065( double x, double *fVal );\n";
    1462          76 : std::string lcl_Erf0065 =
    1463             : "void lcl_Erf0065( double x, double *fVal )\n"
    1464             : "   {\n"
    1465             : "        double pn[] = {\n"
    1466             : "            1.12837916709551256,\n"
    1467             : "            1.35894887627277916E-1,\n"
    1468             : "            4.03259488531795274E-2,\n"
    1469             : "            1.20339380863079457E-3,\n"
    1470             : "            6.49254556481904354E-5\n"
    1471             : "            };\n"
    1472             : "        double qn[] = {\n"
    1473             : "            1.00000000000000000,\n"
    1474             : "            4.53767041780002545E-1,\n"
    1475             : "            8.69936222615385890E-2,\n"
    1476             : "            8.49717371168693357E-3,\n"
    1477             : "            3.64915280629351082E-4\n"
    1478             : "            };\n"
    1479             : "        double fPSum = 0.0;\n"
    1480             : "        double fQSum = 0.0;\n"
    1481             : "        double fXPow = 1.0;\n"
    1482             : "        for ( unsigned int i = 0; i <= 4; ++i )\n"
    1483             : "        {\n"
    1484             : "            fPSum += pn[i]*fXPow;\n"
    1485             : "            fQSum += qn[i]*fXPow;\n"
    1486             : "            fXPow *= x*x;\n"
    1487             : "        }\n"
    1488             : "        *fVal = x * fPSum * pow(fQSum, -1.0);\n"
    1489             : "   }\n";
    1490          76 : std::string rtl_math_erf_rdDecl=
    1491             : "double rtl_math_erf_rd( double x );\n";
    1492          76 : std::string rtl_math_erf_rd =
    1493             : " double rtl_math_erf_rd( double x )\n"
    1494             : " {\n"
    1495             : "     if( x == 0.0 )\n"
    1496             : "         return 0.0;\n"
    1497             : "     bool bNegative = false;\n"
    1498             : "     if ( x < 0.0 )\n"
    1499             : "     {\n"
    1500             : "         x = fabs( x );\n"
    1501             : "         bNegative = true;\n"
    1502             : "     }\n"
    1503             : "     double fErf = 1.0;\n"
    1504             : "     if ( x < 1.0e-10 )\n"
    1505             : "         fErf = (double) (x*1.1283791670955125738961589031215452);\n"
    1506             : "     else if ( x < 0.65 )\n"
    1507             : "         lcl_Erf0065( x, &fErf );\n"
    1508             : "     if ( bNegative )\n"
    1509             : "         fErf *= (-1.0);\n"
    1510             : "     return fErf;\n"
    1511             : " }\n";
    1512          76 : std::string rtl_math_erfc_rdDecl=
    1513             : "double rtl_math_erfc_rd( double x );\n";
    1514          76 : std::string rtl_math_erfc_rd =
    1515             : " double rtl_math_erfc_rd( double x )\n"
    1516             : " {\n"
    1517             : "     if ( x == 0.0 )\n"
    1518             : "        return 1.0;\n"
    1519             : "     bool bNegative = false;\n"
    1520             : "     if ( x < 0.0 )\n"
    1521             : "     {\n"
    1522             : "         x = fabs( x );\n"
    1523             : "         bNegative = true;\n"
    1524             : "     }\n"
    1525             : "     double fErfc = 0.0;\n"
    1526             : "     if ( x >= 0.65 )\n"
    1527             : "     {\n"
    1528             : "         if ( x < 6.0 )\n"
    1529             : "             lcl_Erfc0600( x, &fErfc );\n"
    1530             : "         else\n"
    1531             : "             lcl_Erfc2654( x, &fErfc );\n"
    1532             : "     }\n"
    1533             : "     else\n"
    1534             : "         fErfc = 1.0 - rtl_math_erf_rd( x );\n"
    1535             : "     if ( bNegative )\n"
    1536             : "         fErfc = 2.0 - fErfc;\n"
    1537             : "     return fErfc;\n"
    1538             : " }\n";
    1539             : #endif
    1540             : 
    1541             : /* vim:set shiftwidth=4 softtabstop=4 expandtab: */

Generated by: LCOV version 1.10