Moyal.cc

Go to the documentation of this file.

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

using namespace std;


namespace utl {

  template<int branch>
  inline
  double
  InverseMoyalBranchPointExpansion(const double y)
  {
    const double p = branch * sqrt(2*(y-1));
    const double p2 = Sqr(p);
    return p + (1./6)*p2 + (1./36)*p*p2 + (1./270)*p2*p2;
  }


  template<int branch>
  inline
  double
  InverseMoyalApproximation(const double y)
  {
    if (branch == 1)
      return (y > 4.2) ? y : InverseMoyalBranchPointExpansion<branch>(y);
    else
      return (y > 32.) ? -log(y) : InverseMoyalBranchPointExpansion<branch>(y);
  }


  template<int branch>
  double
  InverseMoyal(const double y, const double epsilon)
  {
    if (y <= 0)
      return numeric_limits<double>::quiet_NaN();
    const double yy = -2*log(y);
    if (yy < 1)
      return numeric_limits<double>::quiet_NaN();
    double x = InverseMoyalApproximation<branch>(yy);
    const double eps = branch*epsilon;
    for (int i = 0; i < 30; ++i) {
      const double expmx = exp(-x);
      const double delta = (x + expmx - yy)/(1 - expmx);
      x -= delta;
      if (fabs(delta) < eps*x)
        return x;
    }
    WARNING("convergence not reached.");
    return x;
  }


  // instantiate branch +1 and -1

  template double InverseMoyal<+1>(const double y, const double eps);
  template double InverseMoyal<-1>(const double y, const double eps);

}