//
// ********************************************************************
// * License and Disclaimer                                           *
// *                                                                  *
// * The  Geant4 software  is  copyright of the Copyright Holders  of *
// * the Geant4 Collaboration.  It is provided  under  the terms  and *
// * conditions of the Geant4 Software License,  included in the file *
// * LICENSE and available at  http://cern.ch/geant4/license .  These *
// * include a list of copyright holders.                             *
// *                                                                  *
// * Neither the authors of this software system, nor their employing *
// * institutes,nor the agencies providing financial support for this *
// * work  make  any representation or  warranty, express or implied, *
// * regarding  this  software system or assume any liability for its *
// * use.  Please see the license in the file  LICENSE  and URL above *
// * for the full disclaimer and the limitation of liability.         *
// *                                                                  *
// * This  code  implementation is the result of  the  scientific and *
// * technical work of the GEANT4 collaboration.                      *
// * By using,  copying,  modifying or  distributing the software (or *
// * any work based  on the software)  you  agree  to acknowledge its *
// * use  in  resulting  scientific  publications,  and indicate your *
// * acceptance of all terms of the Geant4 Software license.          *
// ********************************************************************
//
// G4Integrator inline methods implementation
//
// Author: V.Grichine, 04.09.1999 - First implementation based on
//         G4SimpleIntegration class with H.P.Wellisch, G.Cosmo, and
//         E.TCherniaev advises
// --------------------------------------------------------------------

/////////////////////////////////////////////////////////////////////
//
// Sympson integration method
//
/////////////////////////////////////////////////////////////////////
//
// Integration of class member functions T::f by Simpson method.

template <class T, class F>
G4double G4Integrator<T, F>::Simpson(T& typeT, F f, G4double xInitial,
                                     G4double xFinal, G4int iterationNumber)
{
  G4int i;
  G4double step  = (xFinal - xInitial) / iterationNumber;
  G4double x     = xInitial;
  G4double xPlus = xInitial + 0.5 * step;
  G4double mean  = ((typeT.*f)(xInitial) + (typeT.*f)(xFinal)) * 0.5;
  G4double sum   = (typeT.*f)(xPlus);

  for(i = 1; i < iterationNumber; ++i)
  {
    x += step;
    xPlus += step;
    mean += (typeT.*f)(x);
    sum += (typeT.*f)(xPlus);
  }
  mean += 2.0 * sum;

  return mean * step / 3.0;
}

/////////////////////////////////////////////////////////////////////
//
// Integration of class member functions T::f by Simpson method.
// Convenient to use with 'this' pointer

template <class T, class F>
G4double G4Integrator<T, F>::Simpson(T* ptrT, F f, G4double xInitial,
                                     G4double xFinal, G4int iterationNumber)
{
  G4int i;
  G4double step  = (xFinal - xInitial) / iterationNumber;
  G4double x     = xInitial;
  G4double xPlus = xInitial + 0.5 * step;
  G4double mean  = ((ptrT->*f)(xInitial) + (ptrT->*f)(xFinal)) * 0.5;
  G4double sum   = (ptrT->*f)(xPlus);

  for(i = 1; i < iterationNumber; ++i)
  {
    x += step;
    xPlus += step;
    mean += (ptrT->*f)(x);
    sum += (ptrT->*f)(xPlus);
  }
  mean += 2.0 * sum;

  return mean * step / 3.0;
}

/////////////////////////////////////////////////////////////////////
//
// Integration of class member functions T::f by Simpson method.
// Convenient to use, when function f is defined in global scope, i.e. in main()
// program

template <class T, class F>
G4double G4Integrator<T, F>::Simpson(G4double (*f)(G4double), G4double xInitial,
                                     G4double xFinal, G4int iterationNumber)
{
  G4int i;
  G4double step  = (xFinal - xInitial) / iterationNumber;
  G4double x     = xInitial;
  G4double xPlus = xInitial + 0.5 * step;
  G4double mean  = ((*f)(xInitial) + (*f)(xFinal)) * 0.5;
  G4double sum   = (*f)(xPlus);

  for(i = 1; i < iterationNumber; ++i)
  {
    x += step;
    xPlus += step;
    mean += (*f)(x);
    sum += (*f)(xPlus);
  }
  mean += 2.0 * sum;

  return mean * step / 3.0;
}

//////////////////////////////////////////////////////////////////////////
//
// Adaptive Gauss method
//
//////////////////////////////////////////////////////////////////////////
//
//

template <class T, class F>
G4double G4Integrator<T, F>::Gauss(T& typeT, F f, G4double xInitial,
                                   G4double xFinal)
{
  static const G4double root = 1.0 / std::sqrt(3.0);

  G4double xMean = (xInitial + xFinal) / 2.0;
  G4double Step  = (xFinal - xInitial) / 2.0;
  G4double delta = Step * root;
  G4double sum   = ((typeT.*f)(xMean + delta) + (typeT.*f)(xMean - delta));

  return sum * Step;
}

//////////////////////////////////////////////////////////////////////
//
//

template <class T, class F>
G4double G4Integrator<T, F>::Gauss(T* ptrT, F f, G4double a, G4double b)
{
  return Gauss(*ptrT, f, a, b);
}

///////////////////////////////////////////////////////////////////////
//
//

template <class T, class F>
G4double G4Integrator<T, F>::Gauss(G4double (*f)(G4double), G4double xInitial,
                                   G4double xFinal)
{
  static const G4double root = 1.0 / std::sqrt(3.0);

  G4double xMean = (xInitial + xFinal) / 2.0;
  G4double Step  = (xFinal - xInitial) / 2.0;
  G4double delta = Step * root;
  G4double sum   = ((*f)(xMean + delta) + (*f)(xMean - delta));

  return sum * Step;
}

///////////////////////////////////////////////////////////////////////////
//
//

template <class T, class F>
void G4Integrator<T, F>::AdaptGauss(T& typeT, F f, G4double xInitial,
                                    G4double xFinal, G4double fTolerance,
                                    G4double& sum, G4int& depth)
{
  if(depth > 100)
  {
    G4cout << "G4Integrator<T,F>::AdaptGauss: WARNING !!!" << G4endl;
    G4cout << "Function varies too rapidly to get stated accuracy in 100 steps "
           << G4endl;

    return;
  }
  G4double xMean     = (xInitial + xFinal) / 2.0;
  G4double leftHalf  = Gauss(typeT, f, xInitial, xMean);
  G4double rightHalf = Gauss(typeT, f, xMean, xFinal);
  G4double full      = Gauss(typeT, f, xInitial, xFinal);
  if(std::fabs(leftHalf + rightHalf - full) < fTolerance)
  {
    sum += full;
  }
  else
  {
    ++depth;
    AdaptGauss(typeT, f, xInitial, xMean, fTolerance, sum, depth);
    AdaptGauss(typeT, f, xMean, xFinal, fTolerance, sum, depth);
  }
}

template <class T, class F>
void G4Integrator<T, F>::AdaptGauss(T* ptrT, F f, G4double xInitial,
                                    G4double xFinal, G4double fTolerance,
                                    G4double& sum, G4int& depth)
{
  AdaptGauss(*ptrT, f, xInitial, xFinal, fTolerance, sum, depth);
}

/////////////////////////////////////////////////////////////////////////
//
//
template <class T, class F>
void G4Integrator<T, F>::AdaptGauss(G4double (*f)(G4double), G4double xInitial,
                                    G4double xFinal, G4double fTolerance,
                                    G4double& sum, G4int& depth)
{
  if(depth > 100)
  {
    G4cout << "G4SimpleIntegration::AdaptGauss: WARNING !!!" << G4endl;
    G4cout << "Function varies too rapidly to get stated accuracy in 100 steps "
           << G4endl;

    return;
  }
  G4double xMean     = (xInitial + xFinal) / 2.0;
  G4double leftHalf  = Gauss(f, xInitial, xMean);
  G4double rightHalf = Gauss(f, xMean, xFinal);
  G4double full      = Gauss(f, xInitial, xFinal);
  if(std::fabs(leftHalf + rightHalf - full) < fTolerance)
  {
    sum += full;
  }
  else
  {
    ++depth;
    AdaptGauss(f, xInitial, xMean, fTolerance, sum, depth);
    AdaptGauss(f, xMean, xFinal, fTolerance, sum, depth);
  }
}

////////////////////////////////////////////////////////////////////////
//
// Adaptive Gauss integration with accuracy 'e'
// Convenient for using with class object typeT

template <class T, class F>
G4double G4Integrator<T, F>::AdaptiveGauss(T& typeT, F f, G4double xInitial,
                                           G4double xFinal, G4double e)
{
  G4int depth  = 0;
  G4double sum = 0.0;
  AdaptGauss(typeT, f, xInitial, xFinal, e, sum, depth);
  return sum;
}

////////////////////////////////////////////////////////////////////////
//
// Adaptive Gauss integration with accuracy 'e'
// Convenient for using with 'this' pointer

template <class T, class F>
G4double G4Integrator<T, F>::AdaptiveGauss(T* ptrT, F f, G4double xInitial,
                                           G4double xFinal, G4double e)
{
  return AdaptiveGauss(*ptrT, f, xInitial, xFinal, e);
}

////////////////////////////////////////////////////////////////////////
//
// Adaptive Gauss integration with accuracy 'e'
// Convenient for using with global scope function f

template <class T, class F>
G4double G4Integrator<T, F>::AdaptiveGauss(G4double (*f)(G4double),
                                           G4double xInitial, G4double xFinal,
                                           G4double e)
{
  G4int depth  = 0;
  G4double sum = 0.0;
  AdaptGauss(f, xInitial, xFinal, e, sum, depth);
  return sum;
}

////////////////////////////////////////////////////////////////////////////
// Gauss integration methods involving ortogonal polynomials
////////////////////////////////////////////////////////////////////////////
//
// Methods involving Legendre polynomials
//
/////////////////////////////////////////////////////////////////////////
//
// The value nLegendre set the accuracy required, i.e the number of points
// where the function pFunction will be evaluated during integration.
// The function creates the arrays for abscissas and weights that used
// in Gauss-Legendre quadrature method.
// The values a and b are the limits of integration of the function  f .
// nLegendre MUST BE EVEN !!!
// Returns the integral of the function f between a and b, by 2*fNumber point
// Gauss-Legendre integration: the function is evaluated exactly
// 2*fNumber times at interior points in the range of integration.
// Since the weights and abscissas are, in this case, symmetric around
// the midpoint of the range of integration, there are actually only
// fNumber distinct values of each.
// Convenient for using with some class object dataT

template <class T, class F>
G4double G4Integrator<T, F>::Legendre(T& typeT, F f, G4double a, G4double b,
                                      G4int nLegendre)
{
  G4double nwt, nwt1, temp1, temp2, temp3, temp;
  G4double xDiff, xMean, dx, integral;

  const G4double tolerance = 1.6e-10;
  G4int i, j, k = nLegendre;
  G4int fNumber = (nLegendre + 1) / 2;

  if(2 * fNumber != k)
  {
    G4Exception("G4Integrator<T,F>::Legendre(T&,F, ...)", "InvalidCall",
                FatalException, "Invalid (odd) nLegendre in constructor.");
  }

  G4double* fAbscissa = new G4double[fNumber];
  G4double* fWeight   = new G4double[fNumber];

  for(i = 1; i <= fNumber; ++i)  // Loop over the desired roots
  {
    nwt = std::cos(CLHEP::pi * (i - 0.25) /
                   (k + 0.5));  // Initial root approximation

    do  // loop of Newton's method
    {
      temp1 = 1.0;
      temp2 = 0.0;
      for(j = 1; j <= k; ++j)
      {
        temp3 = temp2;
        temp2 = temp1;
        temp1 = ((2.0 * j - 1.0) * nwt * temp2 - (j - 1.0) * temp3) / j;
      }
      temp = k * (nwt * temp1 - temp2) / (nwt * nwt - 1.0);
      nwt1 = nwt;
      nwt  = nwt1 - temp1 / temp;  // Newton's method
    } while(std::fabs(nwt - nwt1) > tolerance);

    fAbscissa[fNumber - i] = nwt;
    fWeight[fNumber - i]   = 2.0 / ((1.0 - nwt * nwt) * temp * temp);
  }

  //
  // Now we ready to get integral
  //

  xMean    = 0.5 * (a + b);
  xDiff    = 0.5 * (b - a);
  integral = 0.0;
  for(i = 0; i < fNumber; ++i)
  {
    dx = xDiff * fAbscissa[i];
    integral += fWeight[i] * ((typeT.*f)(xMean + dx) + (typeT.*f)(xMean - dx));
  }
  delete[] fAbscissa;
  delete[] fWeight;
  return integral *= xDiff;
}

///////////////////////////////////////////////////////////////////////
//
// Convenient for using with the pointer 'this'

template <class T, class F>
G4double G4Integrator<T, F>::Legendre(T* ptrT, F f, G4double a, G4double b,
                                      G4int nLegendre)
{
  return Legendre(*ptrT, f, a, b, nLegendre);
}

///////////////////////////////////////////////////////////////////////
//
// Convenient for using with global scope function f

template <class T, class F>
G4double G4Integrator<T, F>::Legendre(G4double (*f)(G4double), G4double a,
                                      G4double b, G4int nLegendre)
{
  G4double nwt, nwt1, temp1, temp2, temp3, temp;
  G4double xDiff, xMean, dx, integral;

  const G4double tolerance = 1.6e-10;
  G4int i, j, k = nLegendre;
  G4int fNumber = (nLegendre + 1) / 2;

  if(2 * fNumber != k)
  {
    G4Exception("G4Integrator<T,F>::Legendre(...)", "InvalidCall",
                FatalException, "Invalid (odd) nLegendre in constructor.");
  }

  G4double* fAbscissa = new G4double[fNumber];
  G4double* fWeight   = new G4double[fNumber];

  for(i = 1; i <= fNumber; i++)  // Loop over the desired roots
  {
    nwt = std::cos(CLHEP::pi * (i - 0.25) /
                   (k + 0.5));  // Initial root approximation

    do  // loop of Newton's method
    {
      temp1 = 1.0;
      temp2 = 0.0;
      for(j = 1; j <= k; ++j)
      {
        temp3 = temp2;
        temp2 = temp1;
        temp1 = ((2.0 * j - 1.0) * nwt * temp2 - (j - 1.0) * temp3) / j;
      }
      temp = k * (nwt * temp1 - temp2) / (nwt * nwt - 1.0);
      nwt1 = nwt;
      nwt  = nwt1 - temp1 / temp;  // Newton's method
    } while(std::fabs(nwt - nwt1) > tolerance);

    fAbscissa[fNumber - i] = nwt;
    fWeight[fNumber - i]   = 2.0 / ((1.0 - nwt * nwt) * temp * temp);
  }

  //
  // Now we ready to get integral
  //

  xMean    = 0.5 * (a + b);
  xDiff    = 0.5 * (b - a);
  integral = 0.0;
  for(i = 0; i < fNumber; ++i)
  {
    dx = xDiff * fAbscissa[i];
    integral += fWeight[i] * ((*f)(xMean + dx) + (*f)(xMean - dx));
  }
  delete[] fAbscissa;
  delete[] fWeight;

  return integral *= xDiff;
}

////////////////////////////////////////////////////////////////////////////
//
// Returns the integral of the function to be pointed by T::f between a and b,
// by ten point Gauss-Legendre integration: the function is evaluated exactly
// ten times at interior points in the range of integration. Since the weights
// and abscissas are, in this case, symmetric around the midpoint of the
// range of integration, there are actually only five distinct values of each
// Convenient for using with class object typeT

template <class T, class F>
G4double G4Integrator<T, F>::Legendre10(T& typeT, F f, G4double a, G4double b)
{
  G4int i;
  G4double xDiff, xMean, dx, integral;

  // From Abramowitz M., Stegan I.A. 1964 , Handbook of Math... , p. 916

  static const G4double abscissa[] = { 0.148874338981631, 0.433395394129247,
                                       0.679409568299024, 0.865063366688985,
                                       0.973906528517172 };

  static const G4double weight[] = { 0.295524224714753, 0.269266719309996,
                                     0.219086362515982, 0.149451349150581,
                                     0.066671344308688 };
  xMean                          = 0.5 * (a + b);
  xDiff                          = 0.5 * (b - a);
  integral                       = 0.0;
  for(i = 0; i < 5; ++i)
  {
    dx = xDiff * abscissa[i];
    integral += weight[i] * ((typeT.*f)(xMean + dx) + (typeT.*f)(xMean - dx));
  }
  return integral *= xDiff;
}

///////////////////////////////////////////////////////////////////////////
//
// Convenient for using with the pointer 'this'

template <class T, class F>
G4double G4Integrator<T, F>::Legendre10(T* ptrT, F f, G4double a, G4double b)
{
  return Legendre10(*ptrT, f, a, b);
}

//////////////////////////////////////////////////////////////////////////
//
// Convenient for using with global scope functions

template <class T, class F>
G4double G4Integrator<T, F>::Legendre10(G4double (*f)(G4double), G4double a,
                                        G4double b)
{
  G4int i;
  G4double xDiff, xMean, dx, integral;

  // From Abramowitz M., Stegan I.A. 1964 , Handbook of Math... , p. 916

  static const G4double abscissa[] = { 0.148874338981631, 0.433395394129247,
                                       0.679409568299024, 0.865063366688985,
                                       0.973906528517172 };

  static const G4double weight[] = { 0.295524224714753, 0.269266719309996,
                                     0.219086362515982, 0.149451349150581,
                                     0.066671344308688 };
  xMean                          = 0.5 * (a + b);
  xDiff                          = 0.5 * (b - a);
  integral                       = 0.0;
  for(i = 0; i < 5; ++i)
  {
    dx = xDiff * abscissa[i];
    integral += weight[i] * ((*f)(xMean + dx) + (*f)(xMean - dx));
  }
  return integral *= xDiff;
}

///////////////////////////////////////////////////////////////////////
//
// Returns the integral of the function to be pointed by T::f between a and b,
// by 96 point Gauss-Legendre integration: the function is evaluated exactly
// ten Times at interior points in the range of integration. Since the weights
// and abscissas are, in this case, symmetric around the midpoint of the
// range of integration, there are actually only five distinct values of each
// Convenient for using with some class object typeT

template <class T, class F>
G4double G4Integrator<T, F>::Legendre96(T& typeT, F f, G4double a, G4double b)
{
  G4int i;
  G4double xDiff, xMean, dx, integral;

  // From Abramowitz M., Stegan I.A. 1964 , Handbook of Math... , p. 919

  static const G4double abscissa[] = {
    0.016276744849602969579, 0.048812985136049731112,
    0.081297495464425558994, 0.113695850110665920911,
    0.145973714654896941989, 0.178096882367618602759,  // 6

    0.210031310460567203603, 0.241743156163840012328,
    0.273198812591049141487, 0.304364944354496353024,
    0.335208522892625422616, 0.365696861472313635031,  // 12

    0.395797649828908603285, 0.425478988407300545365,
    0.454709422167743008636, 0.483457973920596359768,
    0.511694177154667673586, 0.539388108324357436227,  // 18

    0.566510418561397168404, 0.593032364777572080684,
    0.618925840125468570386, 0.644163403784967106798,
    0.668718310043916153953, 0.692564536642171561344,  // 24

    0.715676812348967626225, 0.738030643744400132851,
    0.759602341176647498703, 0.780369043867433217604,
    0.800308744139140817229, 0.819400310737931675539,  // 30

    0.837623511228187121494, 0.854959033434601455463,
    0.871388505909296502874, 0.886894517402420416057,
    0.901460635315852341319, 0.915071423120898074206,  // 36

    0.927712456722308690965, 0.939370339752755216932,
    0.950032717784437635756, 0.959688291448742539300,
    0.968326828463264212174, 0.975939174585136466453,  // 42

    0.982517263563014677447, 0.988054126329623799481,
    0.992543900323762624572, 0.995981842987209290650,
    0.998364375863181677724, 0.999689503883230766828  // 48
  };

  static const G4double weight[] = {
    0.032550614492363166242, 0.032516118713868835987,
    0.032447163714064269364, 0.032343822568575928429,
    0.032206204794030250669, 0.032034456231992663218,  // 6

    0.031828758894411006535, 0.031589330770727168558,
    0.031316425596862355813, 0.031010332586313837423,
    0.030671376123669149014, 0.030299915420827593794,  // 12

    0.029896344136328385984, 0.029461089958167905970,
    0.028994614150555236543, 0.028497411065085385646,
    0.027970007616848334440, 0.027412962726029242823,  // 18

    0.026826866725591762198, 0.026212340735672413913,
    0.025570036005349361499, 0.024900633222483610288,
    0.024204841792364691282, 0.023483399085926219842,  // 24

    0.022737069658329374001, 0.021966644438744349195,
    0.021172939892191298988, 0.020356797154333324595,
    0.019519081140145022410, 0.018660679627411467385,  // 30

    0.017782502316045260838, 0.016885479864245172450,
    0.015970562902562291381, 0.015038721026994938006,
    0.014090941772314860916, 0.013128229566961572637,  // 36

    0.012151604671088319635, 0.011162102099838498591,
    0.010160770535008415758, 0.009148671230783386633,
    0.008126876925698759217, 0.007096470791153865269,  // 42

    0.006058545504235961683, 0.005014202742927517693,
    0.003964554338444686674, 0.002910731817934946408,
    0.001853960788946921732, 0.000796792065552012429  // 48
  };
  xMean    = 0.5 * (a + b);
  xDiff    = 0.5 * (b - a);
  integral = 0.0;
  for(i = 0; i < 48; ++i)
  {
    dx = xDiff * abscissa[i];
    integral += weight[i] * ((typeT.*f)(xMean + dx) + (typeT.*f)(xMean - dx));
  }
  return integral *= xDiff;
}

///////////////////////////////////////////////////////////////////////
//
// Convenient for using with the pointer 'this'

template <class T, class F>
G4double G4Integrator<T, F>::Legendre96(T* ptrT, F f, G4double a, G4double b)
{
  return Legendre96(*ptrT, f, a, b);
}

///////////////////////////////////////////////////////////////////////
//
// Convenient for using with global scope function f

template <class T, class F>
G4double G4Integrator<T, F>::Legendre96(G4double (*f)(G4double), G4double a,
                                        G4double b)
{
  G4int i;
  G4double xDiff, xMean, dx, integral;

  // From Abramowitz M., Stegan I.A. 1964 , Handbook of Math... , p. 919

  static const G4double abscissa[] = {
    0.016276744849602969579, 0.048812985136049731112,
    0.081297495464425558994, 0.113695850110665920911,
    0.145973714654896941989, 0.178096882367618602759,  // 6

    0.210031310460567203603, 0.241743156163840012328,
    0.273198812591049141487, 0.304364944354496353024,
    0.335208522892625422616, 0.365696861472313635031,  // 12

    0.395797649828908603285, 0.425478988407300545365,
    0.454709422167743008636, 0.483457973920596359768,
    0.511694177154667673586, 0.539388108324357436227,  // 18

    0.566510418561397168404, 0.593032364777572080684,
    0.618925840125468570386, 0.644163403784967106798,
    0.668718310043916153953, 0.692564536642171561344,  // 24

    0.715676812348967626225, 0.738030643744400132851,
    0.759602341176647498703, 0.780369043867433217604,
    0.800308744139140817229, 0.819400310737931675539,  // 30

    0.837623511228187121494, 0.854959033434601455463,
    0.871388505909296502874, 0.886894517402420416057,
    0.901460635315852341319, 0.915071423120898074206,  // 36

    0.927712456722308690965, 0.939370339752755216932,
    0.950032717784437635756, 0.959688291448742539300,
    0.968326828463264212174, 0.975939174585136466453,  // 42

    0.982517263563014677447, 0.988054126329623799481,
    0.992543900323762624572, 0.995981842987209290650,
    0.998364375863181677724, 0.999689503883230766828  // 48
  };

  static const G4double weight[] = {
    0.032550614492363166242, 0.032516118713868835987,
    0.032447163714064269364, 0.032343822568575928429,
    0.032206204794030250669, 0.032034456231992663218,  // 6

    0.031828758894411006535, 0.031589330770727168558,
    0.031316425596862355813, 0.031010332586313837423,
    0.030671376123669149014, 0.030299915420827593794,  // 12

    0.029896344136328385984, 0.029461089958167905970,
    0.028994614150555236543, 0.028497411065085385646,
    0.027970007616848334440, 0.027412962726029242823,  // 18

    0.026826866725591762198, 0.026212340735672413913,
    0.025570036005349361499, 0.024900633222483610288,
    0.024204841792364691282, 0.023483399085926219842,  // 24

    0.022737069658329374001, 0.021966644438744349195,
    0.021172939892191298988, 0.020356797154333324595,
    0.019519081140145022410, 0.018660679627411467385,  // 30

    0.017782502316045260838, 0.016885479864245172450,
    0.015970562902562291381, 0.015038721026994938006,
    0.014090941772314860916, 0.013128229566961572637,  // 36

    0.012151604671088319635, 0.011162102099838498591,
    0.010160770535008415758, 0.009148671230783386633,
    0.008126876925698759217, 0.007096470791153865269,  // 42

    0.006058545504235961683, 0.005014202742927517693,
    0.003964554338444686674, 0.002910731817934946408,
    0.001853960788946921732, 0.000796792065552012429  // 48
  };
  xMean    = 0.5 * (a + b);
  xDiff    = 0.5 * (b - a);
  integral = 0.0;
  for(i = 0; i < 48; ++i)
  {
    dx = xDiff * abscissa[i];
    integral += weight[i] * ((*f)(xMean + dx) + (*f)(xMean - dx));
  }
  return integral *= xDiff;
}

//////////////////////////////////////////////////////////////////////////////
//
// Methods involving Chebyshev polynomials
//
///////////////////////////////////////////////////////////////////////////
//
// Integrates function pointed by T::f from a to b by Gauss-Chebyshev
// quadrature method.
// Convenient for using with class object typeT

template <class T, class F>
G4double G4Integrator<T, F>::Chebyshev(T& typeT, F f, G4double a, G4double b,
                                       G4int nChebyshev)
{
  G4int i;
  G4double xDiff, xMean, dx, integral = 0.0;

  G4int fNumber       = nChebyshev;  // Try to reduce fNumber twice ??
  G4double cof        = CLHEP::pi / fNumber;
  G4double* fAbscissa = new G4double[fNumber];
  G4double* fWeight   = new G4double[fNumber];
  for(i = 0; i < fNumber; ++i)
  {
    fAbscissa[i] = std::cos(cof * (i + 0.5));
    fWeight[i]   = cof * std::sqrt(1 - fAbscissa[i] * fAbscissa[i]);
  }

  //
  // Now we ready to estimate the integral
  //

  xMean = 0.5 * (a + b);
  xDiff = 0.5 * (b - a);
  for(i = 0; i < fNumber; ++i)
  {
    dx = xDiff * fAbscissa[i];
    integral += fWeight[i] * (typeT.*f)(xMean + dx);
  }
  delete[] fAbscissa;
  delete[] fWeight;
  return integral *= xDiff;
}

///////////////////////////////////////////////////////////////////////
//
// Convenient for using with 'this' pointer

template <class T, class F>
G4double G4Integrator<T, F>::Chebyshev(T* ptrT, F f, G4double a, G4double b,
                                       G4int n)
{
  return Chebyshev(*ptrT, f, a, b, n);
}

////////////////////////////////////////////////////////////////////////
//
// For use with global scope functions f

template <class T, class F>
G4double G4Integrator<T, F>::Chebyshev(G4double (*f)(G4double), G4double a,
                                       G4double b, G4int nChebyshev)
{
  G4int i;
  G4double xDiff, xMean, dx, integral = 0.0;

  G4int fNumber       = nChebyshev;  // Try to reduce fNumber twice ??
  G4double cof        = CLHEP::pi / fNumber;
  G4double* fAbscissa = new G4double[fNumber];
  G4double* fWeight   = new G4double[fNumber];
  for(i = 0; i < fNumber; ++i)
  {
    fAbscissa[i] = std::cos(cof * (i + 0.5));
    fWeight[i]   = cof * std::sqrt(1 - fAbscissa[i] * fAbscissa[i]);
  }

  //
  // Now we ready to estimate the integral
  //

  xMean = 0.5 * (a + b);
  xDiff = 0.5 * (b - a);
  for(i = 0; i < fNumber; ++i)
  {
    dx = xDiff * fAbscissa[i];
    integral += fWeight[i] * (*f)(xMean + dx);
  }
  delete[] fAbscissa;
  delete[] fWeight;
  return integral *= xDiff;
}

//////////////////////////////////////////////////////////////////////
//
// Method involving Laguerre polynomials
//
//////////////////////////////////////////////////////////////////////
//
// Integral from zero to infinity of std::pow(x,alpha)*std::exp(-x)*f(x).
// The value of nLaguerre sets the accuracy.
// The function creates arrays fAbscissa[0,..,nLaguerre-1] and
// fWeight[0,..,nLaguerre-1] .
// Convenient for using with class object 'typeT' and (typeT.*f) function
// (T::f)

template <class T, class F>
G4double G4Integrator<T, F>::Laguerre(T& typeT, F f, G4double alpha,
                                      G4int nLaguerre)
{
  const G4double tolerance = 1.0e-10;
  const G4int maxNumber    = 12;
  G4int i, j, k;
  G4double nwt      = 0., nwt1, temp1, temp2, temp3, temp, cofi;
  G4double integral = 0.0;

  G4int fNumber       = nLaguerre;
  G4double* fAbscissa = new G4double[fNumber];
  G4double* fWeight   = new G4double[fNumber];

  for(i = 1; i <= fNumber; ++i)  // Loop over the desired roots
  {
    if(i == 1)
    {
      nwt = (1.0 + alpha) * (3.0 + 0.92 * alpha) /
            (1.0 + 2.4 * fNumber + 1.8 * alpha);
    }
    else if(i == 2)
    {
      nwt += (15.0 + 6.25 * alpha) / (1.0 + 0.9 * alpha + 2.5 * fNumber);
    }
    else
    {
      cofi = i - 2;
      nwt += ((1.0 + 2.55 * cofi) / (1.9 * cofi) +
              1.26 * cofi * alpha / (1.0 + 3.5 * cofi)) *
             (nwt - fAbscissa[i - 3]) / (1.0 + 0.3 * alpha);
    }
    for(k = 1; k <= maxNumber; ++k)
    {
      temp1 = 1.0;
      temp2 = 0.0;

      for(j = 1; j <= fNumber; ++j)
      {
        temp3 = temp2;
        temp2 = temp1;
        temp1 =
          ((2 * j - 1 + alpha - nwt) * temp2 - (j - 1 + alpha) * temp3) / j;
      }
      temp = (fNumber * temp1 - (fNumber + alpha) * temp2) / nwt;
      nwt1 = nwt;
      nwt  = nwt1 - temp1 / temp;

      if(std::fabs(nwt - nwt1) <= tolerance)
      {
        break;
      }
    }
    if(k > maxNumber)
    {
      G4Exception("G4Integrator<T,F>::Laguerre(T,F, ...)", "Error",
                  FatalException, "Too many (>12) iterations.");
    }

    fAbscissa[i - 1] = nwt;
    fWeight[i - 1]   = -std::exp(GammaLogarithm(alpha + fNumber) -
                               GammaLogarithm((G4double) fNumber)) /
                     (temp * fNumber * temp2);
  }

  //
  // Integral evaluation
  //

  for(i = 0; i < fNumber; ++i)
  {
    integral += fWeight[i] * (typeT.*f)(fAbscissa[i]);
  }
  delete[] fAbscissa;
  delete[] fWeight;
  return integral;
}

//////////////////////////////////////////////////////////////////////
//
//

template <class T, class F>
G4double G4Integrator<T, F>::Laguerre(T* ptrT, F f, G4double alpha,
                                      G4int nLaguerre)
{
  return Laguerre(*ptrT, f, alpha, nLaguerre);
}

////////////////////////////////////////////////////////////////////////
//
// For use with global scope functions f

template <class T, class F>
G4double G4Integrator<T, F>::Laguerre(G4double (*f)(G4double), G4double alpha,
                                      G4int nLaguerre)
{
  const G4double tolerance = 1.0e-10;
  const G4int maxNumber    = 12;
  G4int i, j, k;
  G4double nwt      = 0., nwt1, temp1, temp2, temp3, temp, cofi;
  G4double integral = 0.0;

  G4int fNumber       = nLaguerre;
  G4double* fAbscissa = new G4double[fNumber];
  G4double* fWeight   = new G4double[fNumber];

  for(i = 1; i <= fNumber; ++i)  // Loop over the desired roots
  {
    if(i == 1)
    {
      nwt = (1.0 + alpha) * (3.0 + 0.92 * alpha) /
            (1.0 + 2.4 * fNumber + 1.8 * alpha);
    }
    else if(i == 2)
    {
      nwt += (15.0 + 6.25 * alpha) / (1.0 + 0.9 * alpha + 2.5 * fNumber);
    }
    else
    {
      cofi = i - 2;
      nwt += ((1.0 + 2.55 * cofi) / (1.9 * cofi) +
              1.26 * cofi * alpha / (1.0 + 3.5 * cofi)) *
             (nwt - fAbscissa[i - 3]) / (1.0 + 0.3 * alpha);
    }
    for(k = 1; k <= maxNumber; ++k)
    {
      temp1 = 1.0;
      temp2 = 0.0;

      for(j = 1; j <= fNumber; ++j)
      {
        temp3 = temp2;
        temp2 = temp1;
        temp1 =
          ((2 * j - 1 + alpha - nwt) * temp2 - (j - 1 + alpha) * temp3) / j;
      }
      temp = (fNumber * temp1 - (fNumber + alpha) * temp2) / nwt;
      nwt1 = nwt;
      nwt  = nwt1 - temp1 / temp;

      if(std::fabs(nwt - nwt1) <= tolerance)
      {
        break;
      }
    }
    if(k > maxNumber)
    {
      G4Exception("G4Integrator<T,F>::Laguerre( ...)", "Error", FatalException,
                  "Too many (>12) iterations.");
    }

    fAbscissa[i - 1] = nwt;
    fWeight[i - 1]   = -std::exp(GammaLogarithm(alpha + fNumber) -
                               GammaLogarithm((G4double) fNumber)) /
                     (temp * fNumber * temp2);
  }

  //
  // Integral evaluation
  //

  for(i = 0; i < fNumber; i++)
  {
    integral += fWeight[i] * (*f)(fAbscissa[i]);
  }
  delete[] fAbscissa;
  delete[] fWeight;
  return integral;
}

///////////////////////////////////////////////////////////////////////
//
// Auxiliary function which returns the value of std::log(gamma-function(x))
// Returns the value ln(Gamma(xx) for xx > 0.  Full accuracy is obtained for
// xx > 1. For 0 < xx < 1. the reflection formula (6.1.4) can be used first.
// (Adapted from Numerical Recipes in C)
//

template <class T, class F>
G4double G4Integrator<T, F>::GammaLogarithm(G4double xx)
{
  static const G4double cof[6] = { 76.18009172947146,     -86.50532032941677,
                                   24.01409824083091,     -1.231739572450155,
                                   0.1208650973866179e-2, -0.5395239384953e-5 };
  G4int j;
  G4double x   = xx - 1.0;
  G4double tmp = x + 5.5;
  tmp -= (x + 0.5) * std::log(tmp);
  G4double ser = 1.000000000190015;

  for(j = 0; j <= 5; ++j)
  {
    x += 1.0;
    ser += cof[j] / x;
  }
  return -tmp + std::log(2.5066282746310005 * ser);
}

///////////////////////////////////////////////////////////////////////
//
// Method involving Hermite polynomials
//
///////////////////////////////////////////////////////////////////////
//
//
// Gauss-Hermite method for integration of std::exp(-x*x)*f(x)
// from minus infinity to plus infinity .
//

template <class T, class F>
G4double G4Integrator<T, F>::Hermite(T& typeT, F f, G4int nHermite)
{
  const G4double tolerance = 1.0e-12;
  const G4int maxNumber    = 12;

  G4int i, j, k;
  G4double integral = 0.0;
  G4double nwt      = 0., nwt1, temp1, temp2, temp3, temp;

  G4double piInMinusQ =
    std::pow(CLHEP::pi, -0.25);  // 1.0/std::sqrt(std::sqrt(pi)) ??

  G4int fNumber       = (nHermite + 1) / 2;
  G4double* fAbscissa = new G4double[fNumber];
  G4double* fWeight   = new G4double[fNumber];

  for(i = 1; i <= fNumber; ++i)
  {
    if(i == 1)
    {
      nwt = std::sqrt((G4double)(2 * nHermite + 1)) -
            1.85575001 * std::pow((G4double)(2 * nHermite + 1), -0.16666999);
    }
    else if(i == 2)
    {
      nwt -= 1.14001 * std::pow((G4double) nHermite, 0.425999) / nwt;
    }
    else if(i == 3)
    {
      nwt = 1.86002 * nwt - 0.86002 * fAbscissa[0];
    }
    else if(i == 4)
    {
      nwt = 1.91001 * nwt - 0.91001 * fAbscissa[1];
    }
    else
    {
      nwt = 2.0 * nwt - fAbscissa[i - 3];
    }
    for(k = 1; k <= maxNumber; ++k)
    {
      temp1 = piInMinusQ;
      temp2 = 0.0;

      for(j = 1; j <= nHermite; ++j)
      {
        temp3 = temp2;
        temp2 = temp1;
        temp1 = nwt * std::sqrt(2.0 / j) * temp2 -
                std::sqrt(((G4double)(j - 1)) / j) * temp3;
      }
      temp = std::sqrt((G4double) 2 * nHermite) * temp2;
      nwt1 = nwt;
      nwt  = nwt1 - temp1 / temp;

      if(std::fabs(nwt - nwt1) <= tolerance)
      {
        break;
      }
    }
    if(k > maxNumber)
    {
      G4Exception("G4Integrator<T,F>::Hermite(T,F, ...)", "Error",
                  FatalException, "Too many (>12) iterations.");
    }
    fAbscissa[i - 1] = nwt;
    fWeight[i - 1]   = 2.0 / (temp * temp);
  }

  //
  // Integral calculation
  //

  for(i = 0; i < fNumber; ++i)
  {
    integral +=
      fWeight[i] * ((typeT.*f)(fAbscissa[i]) + (typeT.*f)(-fAbscissa[i]));
  }
  delete[] fAbscissa;
  delete[] fWeight;
  return integral;
}

////////////////////////////////////////////////////////////////////////
//
// For use with 'this' pointer

template <class T, class F>
G4double G4Integrator<T, F>::Hermite(T* ptrT, F f, G4int n)
{
  return Hermite(*ptrT, f, n);
}

////////////////////////////////////////////////////////////////////////
//
// For use with global scope f

template <class T, class F>
G4double G4Integrator<T, F>::Hermite(G4double (*f)(G4double), G4int nHermite)
{
  const G4double tolerance = 1.0e-12;
  const G4int maxNumber    = 12;

  G4int i, j, k;
  G4double integral = 0.0;
  G4double nwt      = 0., nwt1, temp1, temp2, temp3, temp;

  G4double piInMinusQ =
    std::pow(CLHEP::pi, -0.25);  // 1.0/std::sqrt(std::sqrt(pi)) ??

  G4int fNumber       = (nHermite + 1) / 2;
  G4double* fAbscissa = new G4double[fNumber];
  G4double* fWeight   = new G4double[fNumber];

  for(i = 1; i <= fNumber; ++i)
  {
    if(i == 1)
    {
      nwt = std::sqrt((G4double)(2 * nHermite + 1)) -
            1.85575001 * std::pow((G4double)(2 * nHermite + 1), -0.16666999);
    }
    else if(i == 2)
    {
      nwt -= 1.14001 * std::pow((G4double) nHermite, 0.425999) / nwt;
    }
    else if(i == 3)
    {
      nwt = 1.86002 * nwt - 0.86002 * fAbscissa[0];
    }
    else if(i == 4)
    {
      nwt = 1.91001 * nwt - 0.91001 * fAbscissa[1];
    }
    else
    {
      nwt = 2.0 * nwt - fAbscissa[i - 3];
    }
    for(k = 1; k <= maxNumber; ++k)
    {
      temp1 = piInMinusQ;
      temp2 = 0.0;

      for(j = 1; j <= nHermite; ++j)
      {
        temp3 = temp2;
        temp2 = temp1;
        temp1 = nwt * std::sqrt(2.0 / j) * temp2 -
                std::sqrt(((G4double)(j - 1)) / j) * temp3;
      }
      temp = std::sqrt((G4double) 2 * nHermite) * temp2;
      nwt1 = nwt;
      nwt  = nwt1 - temp1 / temp;

      if(std::fabs(nwt - nwt1) <= tolerance)
      {
        break;
      }
    }
    if(k > maxNumber)
    {
      G4Exception("G4Integrator<T,F>::Hermite(...)", "Error", FatalException,
                  "Too many (>12) iterations.");
    }
    fAbscissa[i - 1] = nwt;
    fWeight[i - 1]   = 2.0 / (temp * temp);
  }

  //
  // Integral calculation
  //

  for(i = 0; i < fNumber; ++i)
  {
    integral += fWeight[i] * ((*f)(fAbscissa[i]) + (*f)(-fAbscissa[i]));
  }
  delete[] fAbscissa;
  delete[] fWeight;
  return integral;
}

////////////////////////////////////////////////////////////////////////////
//
// Method involving Jacobi polynomials
//
////////////////////////////////////////////////////////////////////////////
//
// Gauss-Jacobi method for integration of ((1-x)^alpha)*((1+x)^beta)*f(x)
// from minus unit to plus unit .
//

template <class T, class F>
G4double G4Integrator<T, F>::Jacobi(T& typeT, F f, G4double alpha,
                                    G4double beta, G4int nJacobi)
{
  const G4double tolerance = 1.0e-12;
  const G4double maxNumber = 12;
  G4int i, k, j;
  G4double alphaBeta, alphaReduced, betaReduced, root1 = 0., root2 = 0.,
                                                 root3 = 0.;
  G4double a, b, c, nwt1, nwt2, nwt3, nwt, temp, root = 0., rootTemp;

  G4int fNumber       = nJacobi;
  G4double* fAbscissa = new G4double[fNumber];
  G4double* fWeight   = new G4double[fNumber];

  for(i = 1; i <= nJacobi; ++i)
  {
    if(i == 1)
    {
      alphaReduced = alpha / nJacobi;
      betaReduced  = beta / nJacobi;
      root1        = (1.0 + alpha) * (2.78002 / (4.0 + nJacobi * nJacobi) +
                               0.767999 * alphaReduced / nJacobi);
      root2        = 1.0 + 1.48 * alphaReduced + 0.96002 * betaReduced +
              0.451998 * alphaReduced * alphaReduced +
              0.83001 * alphaReduced * betaReduced;
      root = 1.0 - root1 / root2;
    }
    else if(i == 2)
    {
      root1 = (4.1002 + alpha) / ((1.0 + alpha) * (1.0 + 0.155998 * alpha));
      root2 = 1.0 + 0.06 * (nJacobi - 8.0) * (1.0 + 0.12 * alpha) / nJacobi;
      root3 =
        1.0 + 0.012002 * beta * (1.0 + 0.24997 * std::fabs(alpha)) / nJacobi;
      root -= (1.0 - root) * root1 * root2 * root3;
    }
    else if(i == 3)
    {
      root1 = (1.67001 + 0.27998 * alpha) / (1.0 + 0.37002 * alpha);
      root2 = 1.0 + 0.22 * (nJacobi - 8.0) / nJacobi;
      root3 = 1.0 + 8.0 * beta / ((6.28001 + beta) * nJacobi * nJacobi);
      root -= (fAbscissa[0] - root) * root1 * root2 * root3;
    }
    else if(i == nJacobi - 1)
    {
      root1 = (1.0 + 0.235002 * beta) / (0.766001 + 0.118998 * beta);
      root2 = 1.0 / (1.0 + 0.639002 * (nJacobi - 4.0) /
                             (1.0 + 0.71001 * (nJacobi - 4.0)));
      root3 = 1.0 / (1.0 + 20.0 * alpha / ((7.5 + alpha) * nJacobi * nJacobi));
      root += (root - fAbscissa[nJacobi - 4]) * root1 * root2 * root3;
    }
    else if(i == nJacobi)
    {
      root1 = (1.0 + 0.37002 * beta) / (1.67001 + 0.27998 * beta);
      root2 = 1.0 / (1.0 + 0.22 * (nJacobi - 8.0) / nJacobi);
      root3 =
        1.0 / (1.0 + 8.0 * alpha / ((6.28002 + alpha) * nJacobi * nJacobi));
      root += (root - fAbscissa[nJacobi - 3]) * root1 * root2 * root3;
    }
    else
    {
      root = 3.0 * fAbscissa[i - 2] - 3.0 * fAbscissa[i - 3] + fAbscissa[i - 4];
    }
    alphaBeta = alpha + beta;
    for(k = 1; k <= maxNumber; ++k)
    {
      temp = 2.0 + alphaBeta;
      nwt1 = (alpha - beta + temp * root) / 2.0;
      nwt2 = 1.0;
      for(j = 2; j <= nJacobi; ++j)
      {
        nwt3 = nwt2;
        nwt2 = nwt1;
        temp = 2 * j + alphaBeta;
        a    = 2 * j * (j + alphaBeta) * (temp - 2.0);
        b    = (temp - 1.0) *
            (alpha * alpha - beta * beta + temp * (temp - 2.0) * root);
        c    = 2.0 * (j - 1 + alpha) * (j - 1 + beta) * temp;
        nwt1 = (b * nwt2 - c * nwt3) / a;
      }
      nwt = (nJacobi * (alpha - beta - temp * root) * nwt1 +
             2.0 * (nJacobi + alpha) * (nJacobi + beta) * nwt2) /
            (temp * (1.0 - root * root));
      rootTemp = root;
      root     = rootTemp - nwt1 / nwt;
      if(std::fabs(root - rootTemp) <= tolerance)
      {
        break;
      }
    }
    if(k > maxNumber)
    {
      G4Exception("G4Integrator<T,F>::Jacobi(T,F, ...)", "Error",
                  FatalException, "Too many (>12) iterations.");
    }
    fAbscissa[i - 1] = root;
    fWeight[i - 1] =
      std::exp(GammaLogarithm((G4double)(alpha + nJacobi)) +
               GammaLogarithm((G4double)(beta + nJacobi)) -
               GammaLogarithm((G4double)(nJacobi + 1.0)) -
               GammaLogarithm((G4double)(nJacobi + alphaBeta + 1.0))) *
      temp * std::pow(2.0, alphaBeta) / (nwt * nwt2);
  }

  //
  // Calculation of the integral
  //

  G4double integral = 0.0;
  for(i = 0; i < fNumber; ++i)
  {
    integral += fWeight[i] * (typeT.*f)(fAbscissa[i]);
  }
  delete[] fAbscissa;
  delete[] fWeight;
  return integral;
}

/////////////////////////////////////////////////////////////////////////
//
// For use with 'this' pointer

template <class T, class F>
G4double G4Integrator<T, F>::Jacobi(T* ptrT, F f, G4double alpha, G4double beta,
                                    G4int n)
{
  return Jacobi(*ptrT, f, alpha, beta, n);
}

/////////////////////////////////////////////////////////////////////////
//
// For use with global scope f

template <class T, class F>
G4double G4Integrator<T, F>::Jacobi(G4double (*f)(G4double), G4double alpha,
                                    G4double beta, G4int nJacobi)
{
  const G4double tolerance = 1.0e-12;
  const G4double maxNumber = 12;
  G4int i, k, j;
  G4double alphaBeta, alphaReduced, betaReduced, root1 = 0., root2 = 0.,
                                                 root3 = 0.;
  G4double a, b, c, nwt1, nwt2, nwt3, nwt, temp, root = 0., rootTemp;

  G4int fNumber       = nJacobi;
  G4double* fAbscissa = new G4double[fNumber];
  G4double* fWeight   = new G4double[fNumber];

  for(i = 1; i <= nJacobi; ++i)
  {
    if(i == 1)
    {
      alphaReduced = alpha / nJacobi;
      betaReduced  = beta / nJacobi;
      root1        = (1.0 + alpha) * (2.78002 / (4.0 + nJacobi * nJacobi) +
                               0.767999 * alphaReduced / nJacobi);
      root2        = 1.0 + 1.48 * alphaReduced + 0.96002 * betaReduced +
              0.451998 * alphaReduced * alphaReduced +
              0.83001 * alphaReduced * betaReduced;
      root = 1.0 - root1 / root2;
    }
    else if(i == 2)
    {
      root1 = (4.1002 + alpha) / ((1.0 + alpha) * (1.0 + 0.155998 * alpha));
      root2 = 1.0 + 0.06 * (nJacobi - 8.0) * (1.0 + 0.12 * alpha) / nJacobi;
      root3 =
        1.0 + 0.012002 * beta * (1.0 + 0.24997 * std::fabs(alpha)) / nJacobi;
      root -= (1.0 - root) * root1 * root2 * root3;
    }
    else if(i == 3)
    {
      root1 = (1.67001 + 0.27998 * alpha) / (1.0 + 0.37002 * alpha);
      root2 = 1.0 + 0.22 * (nJacobi - 8.0) / nJacobi;
      root3 = 1.0 + 8.0 * beta / ((6.28001 + beta) * nJacobi * nJacobi);
      root -= (fAbscissa[0] - root) * root1 * root2 * root3;
    }
    else if(i == nJacobi - 1)
    {
      root1 = (1.0 + 0.235002 * beta) / (0.766001 + 0.118998 * beta);
      root2 = 1.0 / (1.0 + 0.639002 * (nJacobi - 4.0) /
                             (1.0 + 0.71001 * (nJacobi - 4.0)));
      root3 = 1.0 / (1.0 + 20.0 * alpha / ((7.5 + alpha) * nJacobi * nJacobi));
      root += (root - fAbscissa[nJacobi - 4]) * root1 * root2 * root3;
    }
    else if(i == nJacobi)
    {
      root1 = (1.0 + 0.37002 * beta) / (1.67001 + 0.27998 * beta);
      root2 = 1.0 / (1.0 + 0.22 * (nJacobi - 8.0) / nJacobi);
      root3 =
        1.0 / (1.0 + 8.0 * alpha / ((6.28002 + alpha) * nJacobi * nJacobi));
      root += (root - fAbscissa[nJacobi - 3]) * root1 * root2 * root3;
    }
    else
    {
      root = 3.0 * fAbscissa[i - 2] - 3.0 * fAbscissa[i - 3] + fAbscissa[i - 4];
    }
    alphaBeta = alpha + beta;
    for(k = 1; k <= maxNumber; ++k)
    {
      temp = 2.0 + alphaBeta;
      nwt1 = (alpha - beta + temp * root) / 2.0;
      nwt2 = 1.0;
      for(j = 2; j <= nJacobi; ++j)
      {
        nwt3 = nwt2;
        nwt2 = nwt1;
        temp = 2 * j + alphaBeta;
        a    = 2 * j * (j + alphaBeta) * (temp - 2.0);
        b    = (temp - 1.0) *
            (alpha * alpha - beta * beta + temp * (temp - 2.0) * root);
        c    = 2.0 * (j - 1 + alpha) * (j - 1 + beta) * temp;
        nwt1 = (b * nwt2 - c * nwt3) / a;
      }
      nwt = (nJacobi * (alpha - beta - temp * root) * nwt1 +
             2.0 * (nJacobi + alpha) * (nJacobi + beta) * nwt2) /
            (temp * (1.0 - root * root));
      rootTemp = root;
      root     = rootTemp - nwt1 / nwt;
      if(std::fabs(root - rootTemp) <= tolerance)
      {
        break;
      }
    }
    if(k > maxNumber)
    {
      G4Exception("G4Integrator<T,F>::Jacobi(...)", "Error", FatalException,
                  "Too many (>12) iterations.");
    }
    fAbscissa[i - 1] = root;
    fWeight[i - 1] =
      std::exp(GammaLogarithm((G4double)(alpha + nJacobi)) +
               GammaLogarithm((G4double)(beta + nJacobi)) -
               GammaLogarithm((G4double)(nJacobi + 1.0)) -
               GammaLogarithm((G4double)(nJacobi + alphaBeta + 1.0))) *
      temp * std::pow(2.0, alphaBeta) / (nwt * nwt2);
  }

  //
  // Calculation of the integral
  //

  G4double integral = 0.0;
  for(i = 0; i < fNumber; ++i)
  {
    integral += fWeight[i] * (*f)(fAbscissa[i]);
  }
  delete[] fAbscissa;
  delete[] fWeight;
  return integral;
}

//
//
///////////////////////////////////////////////////////////////////