// $Id: flatToGaussian.cc,v 1.4 2003/08/13 20:00:12 garren Exp $
// -*- C++ -*-
//
// -----------------------------------------------------------------------
//                             HEP Random
//                          --- flatToGaussian ---
//                      class implementation file
// -----------------------------------------------------------------------

// Contains the methods that depend on the 30K-footpring gaussTables.cdat.
//
// flatToGaussian (double x)
// inverseErf     (double x)
// erf		  (double x)

// =======================================================================
// M Fischler	  - Created 1/25/00.
//
// =======================================================================

#include "CLHEP/Random/Stat.h"
#include "CLHEP/Random/defs.h"
#include "CLHEP/Units/PhysicalConstants.h"
#include <iostream>
#include <cmath>

using namespace std;

namespace CLHEP {

double transformSmall (double r);


//
// Table of errInts, for use with transform(r) and quickTransform(r)
//

#ifdef Traces
#define Trace1
#define Trace2
#define Trace3
#endif

// Since all these are this is static to this compilation unit only, the 
// info is establised a priori and not at each invocation.

// The main data is of course the gaussTables table; the rest is all 
// bookkeeping to know what the tables mean.

#define Table0size   200
#define Table1size   250
#define Table2size   200
#define Table3size   250
#define Table4size  1000
#define TableSize   (Table0size+Table1size+Table2size+Table3size+Table4size)

static const int Tsizes[5] =   { Table0size,
				 Table1size,
				 Table2size,
				 Table3size,
				 Table4size };

#define Table0step  (2.0E-13)
#define Table1step  (4.0E-11)  
#define Table2step  (1.0E-8) 
#define Table3step  (2.0E-6) 
#define Table4step  (5.0E-4)

static const double Tsteps[5] = { Table0step,
				 Table1step,
				 Table2step,
				 Table3step,
				 Table4step };

#define Table0offset 0
#define Table1offset (2*(Table0size))
#define Table2offset (2*(Table0size + Table1size))
#define Table3offset (2*(Table0size + Table1size + Table2size))
#define Table4offset (2*(Table0size + Table1size + Table2size + Table3size))

static const int Toffsets[5] = { Table0offset,
				 Table1offset,
				 Table2offset,
				 Table3offset,
				 Table4offset };

  // Here comes the big (30K bytes) table, kept in a file ---

static const double gaussTables [2*TableSize] = {
#include "gaussTables.cdat"
};


double HepStat::flatToGaussian (double r) {

  double sign = +1.0;	// We always compute a negative number of 
				// sigmas.  For r > 0 we will multiply by
				// sign = -1 to return a positive number.
#ifdef Trace1
std::cout << "r = " << r << "\n";
#endif

  if ( r > .5 ) {
    r = 1-r;
    sign = -1.0;
#ifdef Trace1
std::cout << "r = " << r << "sign negative \n";
#endif
  } else if ( r == .5 ) {
    return 0.0;
  }  

  // Find a pointer to the proper table entries, along with the fraction 
  // of the way in the relevant bin dx and the bin size h.
  
  // Optimize for the case of table 4 by testing for that first.  
  // By removing that case from the for loop below, we save not only
  // several table lookups, but also several index calculations that
  // now become (compile-time) constants.
  //
  // Past the case of table 4, we need not be as concerned about speed since
  // this will happen only .1% of the time.

  const double* tptr = 0;
  double  dx = 0;
  double  h = 0;

  // The following big if block will locate tptr, h, and dx from whichever
  // table is applicable:

  register int index;

  if ( r >= Table4step ) {

    index = int((Table4size<<1) * r);	// 1 to Table4size-1 
    if (index <= 0) index = 1; 			// in case of rounding problem
    if (index >= Table4size) index = Table4size-1;
    dx = (Table4size<<1) * r - index; 		// fraction of way to next bin
    h = Table4step;
#ifdef Trace2 
std::cout << "index = " << index << " dx = " << dx << " h = " << h << "\n";
#endif
    index = (index<<1) + (Table4offset-2);	
	// at r = table4step+eps, index refers to the start of table 4 
	// and at r = .5 - eps, index refers to the next-to-last entry:
	// remember, the table has two numbers per actual entry.
#ifdef Trace2 
std::cout << "offset index = " << index << "\n";
#endif

    tptr = &gaussTables [index];
    
  } else if (r < Tsteps[0])  {

    // If r is so small none of the tables apply, use the asymptotic formula.
    return (sign * transformSmall (r));

  } else {
    
    for ( int tableN = 3; tableN >= 0; tableN-- ) {
      if ( r < Tsteps[tableN] ) {continue;} 	// can't happen when tableN==0
#ifdef Trace2 
std::cout << "Using table " << tableN << "\n";
#endif
      double step = Tsteps[tableN];
      index = int(r/step);			// 1 to TableNsize-1 
        // The following two tests should probably never be true, but
        // roundoff might cause index to be outside its proper range.
        // In such a case, the interpolation still makes sense, but we
        // need to take care that tptr points to valid data in the right table.
      if (index == 0) index = 1; 			
      if (index >= Tsizes[tableN]) index = Tsizes[tableN] - 1;
      dx =  r/step - index; 			// fraction of way to next bin
      h  =  step;
#ifdef Trace2 
std::cout << "index = " << index << " dx = " << dx << " h = " << h << "\n";
#endif
      index = (index<<1) + Toffsets[tableN] - 2;
      tptr = &gaussTables [index];
      break;
    } // end of the for loop which finds tptr, dx and h in one of the tables

    // The code can only get to here by a break statement, having set dx etc.
    // It not get to here without going through one of the breaks,
    // because before the for loop we test for the case of r < Tsteps[0].

  } // End of the big if block.

  // At the end of this if block, we have tptr, dx and h -- and if r is less 
  // than the smallest step, we have already returned the proper answer.  

  double  y0 = *tptr++;
  double  d0 = *tptr++;
  double  y1 = *tptr++;
  double  d1 = *tptr;

#ifdef Trace3
std::cout << "y0: " << y0 << " y1: " << y1 << " d0: " << d0 << " d1: " << d1 << "\n";
#endif

  double  x2 = dx * dx;
  double  oneMinusX = 1 - dx;
  double  oneMinusX2 = oneMinusX * oneMinusX;

  double  f0 = (2. * dx + 1.) * oneMinusX2;
  double  f1 = (3. - 2. * dx) * x2;
  double  g0 =  h * dx * oneMinusX2;
  double  g1 =  - h * oneMinusX * x2;

#ifdef Trace3
std::cout << "f0: " << f0 << " f1: " << f1 << " g0: " << g0 << " g1: " << g1 << "\n";
#endif

  double answer = f0 * y0 + f1 * y1 + g0 * d0 + g1 * d1;

#ifdef Trace1
std::cout << "variate is: " << sign*answer << "\n";
#endif

  return sign * answer;

} // flatToGaussian


double transformSmall (double r) {

  // Solve for -v in the asymtotic formula 
  //
  // errInt (-v) =  exp(-v*v/2)         1     1*3    1*3*5
  //		   ------------ * (1 - ---- + ---- - ----- + ... )
  //		   v*sqrt(2*pi)        v**2   v**4   v**6

  // The value of r (=errInt(-v)) supplied is going to less than 2.0E-13,
  // which is such that v < -7.25.  Since the value of r is meaningful only
  // to an absolute error of 1E-16 (double precision accuracy for a number 
  // which on the high side could be of the form 1-epsilon), computing
  // v to more than 3-4 digits of accuracy is suspect; however, to ensure 
  // smoothness with the table generator (which uses quite a few terms) we
  // also use terms up to 1*3*5* ... *13/v**14, and insist on accuracy of
  // solution at the level of 1.0e-7.

  // This routine is called less than one time in a trillion firings, so
  // speed is of no concern.  As a matter of technique, we terminate the
  // iterations in case they would be infinite, but this should not happen.

  double eps = 1.0e-7;
  double guess = 7.5;
  double v;
  
  for ( int i = 1; i < 50; i++ ) {
    double vn2 = 1.0/(guess*guess);
    double s = -13*11*9*7*5*3 * vn2*vn2*vn2*vn2*vn2*vn2*vn2;
    	      s +=    11*9*7*5*3 * vn2*vn2*vn2*vn2*vn2*vn2;
    	      s +=      -9*7*5*3 * vn2*vn2*vn2*vn2*vn2;
	      s +=         7*5*3 * vn2*vn2*vn2*vn2;
    	      s +=          -5*3 * vn2*vn2*vn2;
	      s += 	       3 * vn2*vn2    - vn2  +    1.0;
    v = sqrt ( 2.0 * log ( s / (r*guess*sqrt(CLHEP::twopi)) ) );
    if ( abs(v-guess) < eps ) break;
    guess = v;
  }
 
  return -v;

} // transformSmall()


double HepStat::inverseErf (double t) {

  // This uses erf(x) = 2*gaussCDF(sqrt(2)*x) - 1

  return sqrt(0.5) * flatToGaussian(.5*(t+1));

}


double HepStat::erf (double x) {

// Math:
//
// For any given x we can "quickly" find t0 = erfQ (x) = erf(x) + epsilon.
//
// Then we can find x1 = inverseErf (t0) = inverseErf (erf(x)+epsilon)
//
// Expanding in the small epsion, 
// 
//  x1 = x + epsilon * [deriv(inverseErf(u),u) at u = t0] + O(epsilon**2)
//
// so epsilon is (x1-x) / [deriv(inverseErf(u),u) at u = t0] + O(epsilon**2)
//
// But the derivative of an inverse function is one over the derivative of the
// function, so 
// epsilon  = (x1-x) * [deriv(erf(v),v) at v = x] + O(epsilon**2)
//
// And the definition of the erf integral makes that derivative trivial.
// Ultimately,
//
// erf(x) = erfQ(x) - (inverseErf(erfQ(x))-x) * exp(-x**2) * 2/sqrt(CLHEP::pi)
//

  double t0 = erfQ(x);
  double deriv = exp(-x*x) * (2.0 / sqrt(CLHEP::pi));

  return t0 - (inverseErf (t0) - x) * deriv;

}


}  // namespace CLHEP