Utilities/Math/Math.cc

Go to the documentation of this file.

#include <sstream>
#include <vector>
#include <limits>

#include <utl/Math.h>
#include <utl/ErrorLogger.h>

using namespace std;


namespace utl {

  /*double
  LogGamma(const double x)
  {
    const double p[7] = {
        1.000000000190015,
       76.18009172947146,
      -86.50532032941677,
       24.01409824083091,
       -1.231739572450155,
        1.208650973866179e-3,
       -5.395239384953e-6
    };

    double series = p[0];
    for (int i = 1; i < 7; ++i)
      series += p[i]/(x+i);

    const double xx = x + 5.5;
    const double xxx = xx - log(xx) * (x + 0.5);
    const double sqrt2pi = 2.5066282746310005;

    return log(sqrt2pi/x * series) - xxx;
  }*/


  double
  IncompleteGammaPSeries(const double a, const double x)
  {
    if (x <= 0) {
      if (x < 0)
        ERROR("x less than 0 in function IncompleteGammaPSeries()");
      return 0;
    } else if (a <= 0) {
      ERROR("Invalid arguments in function IncompleteGammaP()");
      return 0;
    } else {
      const double eps = 1e-12;
      const double lnga = LogGamma(a);
      double ap = a;
      double del = 1/a;
      double sum = del;
      for (int n = 1; n <= 100; ++n) {
        ++ap;
        del *= x/ap;
        sum += del;
        if (fabs(del) < fabs(sum)*eps)
          return sum*exp(-x + a*log(x) - lnga);
      }
      ERROR("Too many iterations in function IncompleteGammaPSeries()");
      return 0;
    }
  }


  double
  IncompleteGammaPCF(const double a, const double x)
  {
    if (x < 0 || a <= 0) {
      ERROR("Invalid arguments in function IncompleteGammaP()");
      return 0;
    }

    // Continued fractions by modified Lentz’s method (section 5.2)
    const int maxIt = 100;
    const double eps = 3e-7;
    const double fpMin = 1e-30;

    double b = x + 1 - a;
    double c = 1 / fpMin;
    double d = 1 / b;
    double h = d;
    int i;
    for (i = 1; i <= maxIt; ++i) {
      double an = -i * (i-a);
      b += 2;
      d = an*d + b;
      if (fabs(d) < fpMin)
        d = fpMin;
      c = b + an/c;
      if (fabs(c) < fpMin)
        c = fpMin;
      d = 1 / d;
      double del = d*c;
      h *= del;
      if (fabs(del - 1) < eps)
        break;
    }
    if (i > maxIt)
      ERROR("Too many iterations in function IncompleteGammaPCF()");

    const double lnga = LogGamma(a);
    return h*exp(-x + a*log(x) - lnga);
  }


  double
  BetaCF(const double a, const double b, const double x)
  {
    const int maxIt = 100;
    const double eps = 3e-7;
    const double fpMin = 1e-30;

    const double qab = a + b;  //  These q's will be used in factors that occur
    const double qap = a + 1;  //  in the coefficients (6.4.6).
    const double qam = a - 1;
    double c = 1; // First step of Lentz's method.
    double d = 1 - qab*x/qap;
    if (fabs(d) < fpMin)
      d = fpMin;
    d = 1/d;
    double h = d;

    int m;
    for (m = 1; m <= maxIt; ++m) {
      const int m2 = 2*m;
      // One step (the even one) of the recurrence.
      const double aa1 = m * (b-m) * x / ((qam+m2) * (a+m2));
      d = 1 + aa1*d;
      if (fabs(d) < fpMin)
        d = fpMin;
      c = 1 + aa1/c;
      if (fabs(c) < fpMin)
        c = fpMin;
      d = 1/d;
      h *= d*c;
      // Next step of the recurrence (the odd one).
      const double aa2 = -(a+m) * (qab+m) * x / ((a+m2) * (qap+m2));
      d = 1 + aa2*d;
      if (fabs(d) < fpMin)
        d = fpMin;
      c = 1 + aa2/c;
      if (fabs(c) < fpMin)
        c = fpMin;
      d = 1/d;
      const double del = d*c;
      h *= del;
      if (fabs(del - 1) < eps)
        break; // Are we done?
    }
    if (m > maxIt)
      ERROR("Too many iterations in BetaCF()");

    return h;
  }


  double
  IncompleteBeta(const double a, const double b, const double x)
  {
    double bt;

    if (x < 0 || x > 1)
      ERROR("Bad x in routine incompleteBeta");

    if (!x || x == 1)
      bt = 0;
    else // Factors in front of the continued fraction.
      bt = exp(LogGamma(a + b) - LogGamma(a) - LogGamma(b) + a*log(x) + b*log(1-x));

    if (x < (a+1)/(a+b+2))
      // Use continued fraction directly.
      return bt * BetaCF(a, b, x) / a;
    else
      // Use continued fraction after making the symmetry transformation.
      return 1 - bt * BetaCF(b, a, 1-x) / b;
  }


  void
  PropagateAxisErrors(const Vector::Triple& u_v_w,
                      const Vector::Triple& sigma_u2_uv_v2,
                      double& thetaError, double& phiError, double& thetaPhiCorrelation)
  {
    double u;
    double v;
    double w;
    boost::tie(u, v, w) = u_v_w;
    double sigmaU2;
    double sigmaUV;
    double sigmaV2;
    boost::tie(sigmaU2, sigmaUV, sigmaV2) = sigma_u2_uv_v2;

    const double u2 = u*u;
    const double uv = u*v;
    const double v2 = v*v;
    const double duvsigmaUV = 2*uv*sigmaUV;
    const double u2pv2 = u2 + v2;

    thetaError = sqrt((u2*sigmaU2 + v2*sigmaV2 + duvsigmaUV) / u2pv2) / w;
    phiError = sqrt(v2*sigmaU2 + u2*sigmaV2 - duvsigmaUV) / u2pv2;
    thetaPhiCorrelation =
      (-u*v*sigmaU2 + (u2-v2)*sigmaUV + u*v*sigmaV2) / (w*pow(u2pv2, 1.5)*thetaError*phiError);
  }

}