sc/source/core/opencl/opinlinefun_statistical.cxx

   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 std::string MinDecl = "#define Min 2.22507e-308\n";
  13 std::string F_PIDecl="#define F_PI 3.1415926\n";
  14 std::string fBigInvDecl ="#define fBigInv  2.22045e-016\n";
  15 std::string fMachEpsDecl ="#define fMachEps  2.22045e-016\n";
  16 std::string fLogDblMaxDecl ="#define fLogDblMax  log(1.79769e+308)\n";
  17 std::string fHalfMachEpsDecl ="#define fHalfMachEps  0.5*2.22045e-016\n";
  18 std::string fMaxGammaArgumentDecl=
  19 "#define fMaxGammaArgument 171.624376956302\n";
  20 std::string GetValueDecl=
  21 "double  GetValue( double x,double fp,double fDF );\n";
  22 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 std::string GetGammaSeriesDecl=
  28 "double GetGammaSeries( double fA, double fX );\n";
  29 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 std::string GetGammaContFractionDecl =  "double GetGammaContFraction( double "
  49 "fA, double fX );\n";
  50 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 std::string GetLowRegIGammaDecl = "double GetLowRegIGamma( double "
  99 "fA, double fX );\n";
 100 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 std::string GetGammaDistDecl = "double GetGammaDist( double fX, double "
 111 "fAlpha, double fLambda );\n";
 112 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 std::string GetGammaDistPDFDecl =  "double GetGammaDistPDF( double fX"
 121 ", double fAlpha, double fLambda );\n";
 122 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 std::string GetBetaDistDecl =
 175 "double GetBetaDist(double fXin, double fAlpha, double fBeta);\n";
 176 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 std::string GetFDistDecl =
 221     "double GetFDist(double x, double fF1, double fF2);\n";
 222 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 std::string GetGammaInvValueDecl = "double"
 231 " GetGammaInvValue(double fAlpha,double fBeta,double fX1 );\n";
 232 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 std::string GetFInvValueDecl = "double GetFInvValue(double fF1,double fF2"
 249 ",double fX1 );";
 250 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 std::string GetBinomDistPMFDecl =
 302     "double GetBinomDistPMF(double x, double n, double p);";
 303 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 std::string lcl_GetBinomDistRangeDecl =
 331     "double lcl_GetBinomDistRange(double n, \n"
 332 "double xs, double xe, double fFactor, double p, double q);";
 333 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 std::string GetLogGammaDecl = "double GetLogGamma(double fZ);\n";
 353 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 std::string GetChiDistDecl = "double GetChiDist(double fX, double fDF);\n";
 366 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 std::string GetChiSqDistCDFDecl =
 376 "double GetChiSqDistCDF(double fX, double fDF);\n";
 377 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 std::string GetChiSqDistPDFDecl=
 387 "double GetChiSqDistPDF(double fX, double fDF);\n";
 388 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 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 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 std::string lcl_IterateInverseChiInvDecl =
 523 "static double lcl_IterateInverseChiInv"
 524 "(double fp, double fdf, double fAx, double fBx, bool *rConvError);\n";
 525 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 std::string lcl_IterateInverseChiSQInvDecl =
 620 "static double lcl_IterateInverseChiSQInv( double fp, double fdf, \n"
 621 "   double fAx, double fBx, bool *rConvError );\n";
 622 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 std::string gaussinvDecl = "double gaussinv(double x);\n";
 718 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 std::string lcl_GetLogGammaHelperDecl=
 879 "static double lcl_GetLogGammaHelper(double fZ);\n";
 880 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 std::string lcl_GetGammaHelperDecl=
 888 "static double lcl_GetGammaHelper(double fZ);\n";
 889 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 std::string lcl_getLanczosSumDecl=
 907 "static double lcl_getLanczosSum(double fZ);\n";
 908 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 std::string GetUpRegIGammaDecl=
1031 " double GetUpRegIGamma( double fA, double fX ) ;\n";
1032 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 std::string lcl_HasChangeOfSignDecl=
1044 "static inline bool lcl_HasChangeOfSign( double u, double w );\n";
1045 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 std::string GetTDistDecl=" double GetTDist(double T, double fDF);\n";
1052 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 std::string GetBetaDecl=" double GetBeta(double fAlpha, double fBeta);\n";
1059 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 std::string GetLogBetaDecl=
1078 " double GetLogBeta(double fAlpha, double fBeta);\n";
1079 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 std::string GetBetaDistPDFDecl=
1097 "double GetBetaDistPDF(double fX, double fA, double fB);\n";
1098 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 std::string lcl_GetBetaHelperContFracDecl=
1163 "double lcl_GetBetaHelperContFrac(double fX, double fA, double fB);\n";
1164 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 std::string lcl_IterateInverseDecl=
1200 "double lcl_IterateInverse("
1201 "double fAx, double fBx, bool* rConvError,double fp,double fDF );\n";
1202 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 std::string phiDecl=
1295 "double phi(double x);\n";
1296 std::string phi =
1297 "double phi(double x)\n"
1298 "{\n"
1299 "    return  0.39894228040143268 * exp(-(x * x) / 2.0);\n"
1300 "}\n";
1301 std::string taylorDecl =
1302 "double taylor(double* pPolynom, uint nMax, double x);\n";
1303 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 std::string gaussDecl = "double gauss(double x);\n";
1314 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 std::string lcl_Erfc0600Decl=
1365 "void lcl_Erfc0600( double x, double *fVal );\n";
1366 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 std::string lcl_Erfc2654Decl=
1429 "void lcl_Erfc2654( double x, double *fVal );\n";
1430 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 std::string lcl_Erf0065Decl=
1461 "void lcl_Erf0065( double x, double *fVal );\n";
1462 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 std::string rtl_math_erf_rdDecl=
1491 "double rtl_math_erf_rd( double x );\n";
1492 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 std::string rtl_math_erfc_rdDecl=
1513 "double rtl_math_erfc_rd( double x );\n";
1514 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: */