namespace TMath

 The DiLogarithm function
 Code translated by R.Brun from CERNLIB DILOG function C332

 Computation of the error function erf(x).
 Erf(x) = (2/sqrt(pi)) Integral(exp(-t^2))dt between 0 and x

 Compute the complementary error function erfc(x).
 Erfc(x) = (2/sqrt(pi)) Integral(exp(-t^2))dt between x and infinity

 returns  the inverse error function
 x must be  <-1<x<1

 returns  the inverse of the complementary error function
 x must be  0<x<2
 implement using  the quantile of the normal distribution
 instead of ErfInverse for better numerical precision for large x

 Compute factorial(n).

 Computation of the normal frequency function freq(x).
 Freq(x) = (1/sqrt(2pi)) Integral(exp(-t^2/2))dt between -infinity and x.

 Translated from CERNLIB C300 by Rene Brun.

 Computation of gamma(z) for all z.

 C.Lanczos, SIAM Journal of Numerical Analysis B1 (1964), 86.

 Computation of the normalized lower incomplete gamma function P(a,x) as defined in the
 Handbook of Mathematical Functions by Abramowitz and Stegun, formula 6.5.1 on page 260 .
 Its normalization is such that TMath::Gamma(a,+infinity) = 1 .





--- Nve 14-nov-1998 UU-SAP Utrecht

 Calculate a Breit Wigner function with mean and gamma.

 Calculate a gaussian function with mean and sigma.
 If norm=kTRUE (default is kFALSE) the result is divided
 by sqrt(2*Pi)*sigma.

 The LANDAU function.
 mu is a location parameter and correspond approximatly to the most probable value
 and sigma is a scale parameter (not the sigma of the full distribution which is not defined)
 Note that for mu=0 and sigma=1 (default values) the exact location of the maximum of the distribution
 (most proble value) is at x = -0.22278
 This function has been adapted from the CERNLIB routine G110 denlan.
 If norm=kTRUE (default is kFALSE) the result is divided by sigma

 Computation of ln[gamma(z)] for all z.

 C.Lanczos, SIAM Journal of Numerical Analysis B1 (1964), 86.

 The accuracy of the result is better than 2e-10.

--- Nve 14-nov-1998 UU-SAP Utrecht

 Normalize a vector v in place.
 Returns the norm of the original vector.

 Normalize a vector v in place.
 Returns the norm of the original vector.
 This implementation (thanks Kevin Lynch <krlynch@bu.edu>) is protected
 against possible overflows.

 compute the Poisson distribution function for (x,par)
 The Poisson PDF is implemented by means of Euler's Gamma-function
 (for the factorial), so for any x integer argument it is correct.
 BUT for non-integer x values, it IS NOT equal to the Poisson distribution.
 see TMath::PoissonI to get a non-smooth function.
 Note that for large values of par, it is better to call
     TMath::Gaus(x,par,sqrt(par),kTRUE)

 compute the Poisson distribution function for (x,par)
 This is a non-smooth function.
 This function is equivalent to ROOT::Math::poisson_pdf

 Computation of the probability for a certain Chi-squared (chi2)
 and number of degrees of freedom (ndf).

 Calculations are based on the incomplete gamma function P(a,x),
 where a=ndf/2 and x=chi2/2.

 P(a,x) represents the probability that the observed Chi-squared
 for a correct model should be less than the value chi2.

 The returned probability corresponds to 1-P(a,x),
 which denotes the probability that an observed Chi-squared exceeds
 the value chi2 by chance, even for a correct model.

--- NvE 14-nov-1998 UU-SAP Utrecht

 Calculates the Kolmogorov distribution function,

 which gives the probability that Kolmogorov's test statistic will exceed
 the value z assuming the null hypothesis. This gives a very powerful
 test for comparing two one-dimensional distributions.
 see, for example, Eadie et al, "statistocal Methods in Experimental
 Physics', pp 269-270).

 This function returns the confidence level for the null hypothesis, where:
   z = dn*sqrt(n), and
   dn  is the maximum deviation between a hypothetical distribution
       function and an experimental distribution with
   n    events

 NOTE: To compare two experimental distributions with m and n events,
       use z = sqrt(m*n/(m+n))*dn

 Accuracy: The function is far too accurate for any imaginable application.
           Probabilities less than 10^-15 are returned as zero.
           However, remember that the formula is only valid for "large" n.
 Theta function inversion formula is used for z <= 1

 This function was translated by Rene Brun from PROBKL in CERNLIB.

  Statistical test whether two one-dimensional sets of points are compatible
  with coming from the same parent distribution, using the Kolmogorov test.
  That is, it is used to compare two experimental distributions of unbinned data.

  Input:
  a,b: One-dimensional arrays of length na, nb, respectively.
       The elements of a and b must be given in ascending order.
  option is a character string to specify options
         "D" Put out a line of "Debug" printout
         "M" Return the Maximum Kolmogorov distance instead of prob

  Output:
 The returned value prob is a calculated confidence level which gives a
 statistical test for compatibility of a and b.
 Values of prob close to zero are taken as indicating a small probability
 of compatibility. For two point sets drawn randomly from the same parent
 distribution, the value of prob should be uniformly distributed between
 zero and one.
   in case of error the function return -1
   If the 2 sets have a different number of points, the minimum of
   the two sets is used.

 Method:
 The Kolmogorov test is used. The test statistic is the maximum deviation
 between the two integrated distribution functions, multiplied by the
 normalizing factor (rdmax*sqrt(na*nb/(na+nb)).

  Code adapted by Rene Brun from CERNLIB routine TKOLMO (Fred James)
   (W.T. Eadie, D. Drijard, F.E. James, M. Roos and B. Sadoulet,
      Statistical Methods in Experimental Physics, (North-Holland,
      Amsterdam 1971) 269-271)

  Method Improvement by Jason A Detwiler (JADetwiler@lbl.gov)

   The nuts-and-bolts of the TMath::KolmogorovTest() algorithm is a for-loop
   over the two sorted arrays a and b representing empirical distribution
   functions. The for-loop handles 3 cases: when the next points to be
   evaluated satisfy a>b, a<b, or a=b:

      for (Int_t i=0;i<na+nb;i++) {
         if (a[ia-1] < b[ib-1]) {
            rdiff -= sa;
            ia++;
            if (ia > na) {ok = kTRUE; break;}
         } else if (a[ia-1] > b[ib-1]) {
            rdiff += sb;
            ib++;
            if (ib > nb) {ok = kTRUE; break;}
         } else {
            rdiff += sb - sa;
            ia++;
            ib++;
            if (ia > na) {ok = kTRUE; break;}
            if (ib > nb) {ok = kTRUE; break;}
        }
         rdmax = TMath::Max(rdmax,TMath::Abs(rdiff));
      }

   For the last case, a=b, the algorithm advances each array by one index in an
   attempt to move through the equality. However, this is incorrect when one or
   the other of a or b (or both) have a repeated value, call it x. For the KS
   statistic to be computed properly, rdiff needs to be calculated after all of
   the a and b at x have been tallied (this is due to the definition of the
   empirical distribution function; another way to convince yourself that the
   old CERNLIB method is wrong is that it implies that the function defined as the
   difference between a and b is multi-valued at x -- besides being ugly, this
   would invalidate Kolmogorov's theorem).

   The solution is to just add while-loops into the equality-case handling to
   perform the tally:

         } else {
            double x = a[ia-1];
            while(a[ia-1] == x && ia <= na) {
              rdiff -= sa;
              ia++;
            }
            while(b[ib-1] == x && ib <= nb) {
              rdiff += sb;
              ib++;
            }
            if (ia > na) {ok = kTRUE; break;}
            if (ib > nb) {ok = kTRUE; break;}
         }


  NOTE1
  A good description of the Kolmogorov test can be seen at:
    http://www.itl.nist.gov/div898/handbook/eda/section3/eda35g.htm

 Computation of Voigt function (normalised).
 Voigt is a convolution of
 gauss(xx) = 1/(sqrt(2*pi)*sigma) * exp(xx*xx/(2*sigma*sigma)
 and
 lorentz(xx) = (1/pi) * (lg/2) / (xx*xx + lg*lg/4)
 functions.

 The Voigt function is known to be the real part of Faddeeva function also
 called complex error function [2].

 The algoritm was developed by J. Humlicek [1].
 This code is based on fortran code presented by R. J. Wells [2].
 Translated and adapted by Miha D. Puc

 To calculate the Faddeeva function with relative error less than 10^(-r).
 r can be set by the the user subject to the constraints 2 <= r <= 5.

 [1] J. Humlicek, JQSRT, 21, 437 (1982).
 [2] R.J. Wells "Rapid Approximation to the Voigt/Faddeeva Function and its
 Derivatives" JQSRT 62 (1999), pp 29-48.
 http://www-atm.physics.ox.ac.uk/user/wells/voigt.html

 Calculates roots of polynomial of 3rd order a*x^3 + b*x^2 + c*x + d, where
 a == coef[3], b == coef[2], c == coef[1], d == coef[0]
coef[3] must be different from 0
 If the boolean returned by the method is false:
    ==> there are 3 real roots a,b,c
 If the boolean returned by the method is true:
    ==> there is one real root a and 2 complex conjugates roots (b+i*c,b-i*c)
 Author: Francois-Xavier Gentit

Computes sample quantiles, corresponding to the given probabilities
Parameters:
  x -the data sample
  n - its size
  quantiles - computed quantiles are returned in there
  prob - probabilities where to compute quantiles
  nprob - size of prob array
  isSorted - is the input array x sorted?
  NOTE, that when the input is not sorted, an array of integers of size n needs
        to be allocated. It can be passed by the user in parameter index,
        or, if not passed, it will be allocated inside the function

  type - method to compute (from 1 to 9). Following types are provided:
  Discontinuous:
    type=1 - inverse of the empirical distribution function
    type=2 - like type 1, but with averaging at discontinuities
    type=3 - SAS definition: nearest even order statistic
  Piecwise linear continuous:
    In this case, sample quantiles can be obtained by linear interpolation
    between the k-th order statistic and p(k).
    type=4 - linear interpolation of empirical cdf, p(k)=k/n;
    type=5 - a very popular definition, p(k) = (k-0.5)/n;
    type=6 - used by Minitab and SPSS, p(k) = k/(n+1);
    type=7 - used by S-Plus and R, p(k) = (k-1)/(n-1);
    type=8 - resulting sample quantiles are approximately median unbiased
             regardless of the distribution of x. p(k) = (k-1/3)/(n+1/3);
    type=9 - resulting sample quantiles are approximately unbiased, when
             the sample comes from Normal distribution. p(k)=(k-3/8)/(n+1/4);

    default type = 7

 References:
 1) Hyndman, R.J and Fan, Y, (1996) "Sample quantiles in statistical packages"
                                     American Statistician, 50, 361-365
 2) R Project documentation for the function quantile of package {stats}

 Bubble sort variant to obtain the order of an array's elements into
 an index in order to do more useful things than the standard built
 in functions.
 *arr1 is unchanged;
 *arr2 is the array of indicies corresponding to the decending value
 of arr1 with arr2[0] corresponding to the largest arr1 value and
 arr2[Narr] the smallest.

  Author: Adrian Bevan (bevan@slac.stanford.edu)

 Opposite ordering of the array arr2[] to that of BubbleHigh.

  Author: Adrian Bevan (bevan@slac.stanford.edu)

 Calculates hash index from any char string.
 Based on precalculated table of 256 specially selected numbers.
 These numbers are selected in such a way, that for string
 length == 4 (integer number) the hash is unambigous, i.e.
 from hash value we can recalculate input (no degeneration).

 The quality of hash method is good enough, that
 "random" numbers made as R = Hash(1), Hash(2), ...Hash(N)
 tested by <R>, <R*R>, <Ri*Ri+1> gives the same result
 as for libc rand().

 For string:  i = TMath::Hash(string,nstring);
 For int:     i = TMath::Hash(&intword,sizeof(int));
 For pointer: i = TMath::Hash(&pointer,sizeof(void*));

              V.Perev
 This function is kept for back compatibility. The code previously in this function
 has been moved to the static function TString::Hash

 Compute the modified Bessel function I_0(x) for any real x.

--- NvE 12-mar-2000 UU-SAP Utrecht

 Compute the modified Bessel function K_0(x) for positive real x.

  M.Abramowitz and I.A.Stegun, Handbook of Mathematical Functions,
     Applied Mathematics Series vol. 55 (1964), Washington.

--- NvE 12-mar-2000 UU-SAP Utrecht

 Compute the modified Bessel function I_1(x) for any real x.

  M.Abramowitz and I.A.Stegun, Handbook of Mathematical Functions,
     Applied Mathematics Series vol. 55 (1964), Washington.

--- NvE 12-mar-2000 UU-SAP Utrecht

 Compute the modified Bessel function K_1(x) for positive real x.

  M.Abramowitz and I.A.Stegun, Handbook of Mathematical Functions,
     Applied Mathematics Series vol. 55 (1964), Washington.

--- NvE 12-mar-2000 UU-SAP Utrecht

 Compute the Integer Order Modified Bessel function K_n(x)
 for n=0,1,2,... and positive real x.

--- NvE 12-mar-2000 UU-SAP Utrecht

 Compute the Integer Order Modified Bessel function I_n(x)
 for n=0,1,2,... and any real x.

--- NvE 12-mar-2000 UU-SAP Utrecht

 Returns the Bessel function J0(x) for any real x.

 Returns the Bessel function J1(x) for any real x.

 Returns the Bessel function Y0(x) for positive x.

 Returns the Bessel function Y1(x) for positive x.

 Struve Functions of Order 0

 Converted from CERNLIB M342 by Rene Brun.

 Struve Functions of Order 1

 Converted from CERNLIB M342 by Rene Brun.

 Modified Struve Function of Order 0.
 By Kirill Filimonov.

 Modified Struve Function of Order 1.
 By Kirill Filimonov.

 Calculates Beta-function Gamma(p)*Gamma(q)/Gamma(p+q).

 Continued fraction evaluation by modified Lentz's method
 used in calculation of incomplete Beta function.

 Computes the probability density function of the Beta distribution
 (the distribution function is computed in BetaDistI).
 The first argument is the point, where the function will be
 computed, second and third are the function parameters.
 Since the Beta distribution is bounded on both sides, it's often
 used to represent processes with natural lower and upper limits.

 Computes the distribution function of the Beta distribution.
 The first argument is the point, where the function will be
 computed, second and third are the function parameters.
 Since the Beta distribution is bounded on both sides, it's often
 used to represent processes with natural lower and upper limits.

 Calculates the incomplete Beta-function.

 Calculate the binomial coefficient n over k.

 Suppose an event occurs with probability _p_ per trial
 Then the probability P of its occuring _k_ or more times
 in _n_ trials is termed a cumulative binomial probability
 the formula is P = sum_from_j=k_to_n(TMath::Binomial(n, j)*
 *TMath::Power(p, j)*TMath::Power(1-p, n-j)
 For _n_ larger than 12 BetaIncomplete is a much better way
 to evaluate the sum than would be the straightforward sum calculation
 for _n_ smaller than 12 either method is acceptable
 ("Numerical Recipes")
     --implementation by Anna Kreshuk

 Computes the density of Cauchy distribution at point x
 by default, standard Cauchy distribution is used (t=0, s=1)
    t is the location parameter
    s is the scale parameter
 The Cauchy distribution, also called Lorentzian distribution,
 is a continuous distribution describing resonance behavior
 The mean and standard deviation of the Cauchy distribution are undefined.
 The practical meaning of this is that collecting 1,000 data points gives
 no more accurate an estimate of the mean and standard deviation than
 does a single point.
 The formula was taken from "Engineering Statistics Handbook" on site
 http://www.itl.nist.gov/div898/handbook/eda/section3/eda3663.htm
 Implementation by Anna Kreshuk.
 Example:
    TF1* fc = new TF1("fc", "TMath::CauchyDist(x, [0], [1])", -5, 5);
    fc->SetParameters(0, 1);
    fc->Draw();

 Evaluate the quantiles of the chi-squared probability distribution function.
 Algorithm AS 91   Appl. Statist. (1975) Vol.24, P.35
 implemented by Anna Kreshuk.
 Incorporates the suggested changes in AS R85 (vol.40(1), pp.233-5, 1991)
 Parameters:
   p   - the probability value, at which the quantile is computed
   ndf - number of degrees of freedom

 Computes the density function of F-distribution
 (probability function, integral of density, is computed in FDistI).

 Parameters N and M stand for degrees of freedom of chi-squares
 mentioned above parameter F is the actual variable x of the
 density function p(x) and the point at which the density function
 is calculated.

 About F distribution:
 F-distribution arises in testing whether two random samples
 have the same variance. It is the ratio of two chi-square
 distributions, with N and M degrees of freedom respectively,
 where each chi-square is first divided by it's number of degrees
 of freedom.
 Implementation by Anna Kreshuk.

 Calculates the cumulative distribution function of F-distribution,
 this function occurs in the statistical test of whether two observed
 samples have the same variance. For this test a certain statistic F,
 the ratio of observed dispersion of the first sample to that of the
 second sample, is calculated. N and M stand for numbers of degrees
 of freedom in the samples 1-FDistI() is the significance level at
 which the hypothesis "1 has smaller variance than 2" can be rejected.
 A small numerical value of 1 - FDistI() implies a very significant
 rejection, in turn implying high confidence in the hypothesis
 "1 has variance greater than 2".
 Implementation by Anna Kreshuk.

 Computes the density function of Gamma distribution at point x.
   gamma - shape parameter
   mu    - location parameter
   beta  - scale parameter

 The definition can be found in "Engineering Statistics Handbook" on site
 http://www.itl.nist.gov/div898/handbook/eda/section3/eda366b.htm
 use now implementation in ROOT::Math::gamma_pdf

 Computes the probability density function of Laplace distribution
 at point x, with location parameter alpha and shape parameter beta.
 By default, alpha=0, beta=1
 This distribution is known under different names, most common is
 double exponential distribution, but it also appears as
 the two-tailed exponential or the bilateral exponential distribution

 Computes the distribution function of Laplace distribution
 at point x, with location parameter alpha and shape parameter beta.
 By default, alpha=0, beta=1
 This distribution is known under different names, most common is
 double exponential distribution, but it also appears as
 the two-tailed exponential or the bilateral exponential distribution

 Computes the density of LogNormal distribution at point x.
 Variable X has lognormal distribution if Y=Ln(X) has normal distribution
 sigma is the shape parameter
 theta is the location parameter
 m is the scale parameter
 The formula was taken from "Engineering Statistics Handbook" on site
 http://www.itl.nist.gov/div898/handbook/eda/section3/eda3669.htm
 Implementation using ROOT::Math::lognormal_pdf

 Computes quantiles for standard normal distribution N(0, 1)
 at probability p
 ALGORITHM AS241  APPL. STATIST. (1988) VOL. 37, NO. 3, 477-484.

 Simple recursive algorithm to find the permutations of
 n natural numbers, not necessarily all distinct
 adapted from CERNLIB routine PERMU.
 The input array has to be initialised with a non descending
 sequence. The method returns kFALSE when
 all combinations are exhausted.

 Computes density function for Student's t- distribution
 (the probability function (integral of density) is computed in StudentI).

 First parameter stands for x - the actual variable of the
 density function p(x) and the point at which the density is calculated.
 Second parameter stands for number of degrees of freedom.

 About Student distribution:
 Student's t-distribution is used for many significance tests, for example,
 for the Student's t-tests for the statistical significance of difference
 between two sample means and for confidence intervals for the difference
 between two population means.

 Example: suppose we have a random sample of size n drawn from normal
 distribution with mean Mu and st.deviation Sigma. Then the variable

   t = (sample_mean - Mu)/(sample_deviation / sqrt(n))

 has Student's t-distribution with n-1 degrees of freedom.

 NOTE that this function's second argument is number of degrees of freedom,
 not the sample size.

 As the number of degrees of freedom grows, t-distribution approaches
 Normal(0,1) distribution.
 Implementation by Anna Kreshuk.

 Calculates the cumulative distribution function of Student's
 t-distribution second parameter stands for number of degrees of freedom,
 not for the number of samples
 if x has Student's t-distribution, the function returns the probability of
 x being less than T.
 Implementation by Anna Kreshuk.

 Computes quantiles of the Student's t-distribution
 1st argument is the probability, at which the quantile is computed
 2nd argument - the number of degrees of freedom of the
 Student distribution
 When the 3rd argument lower_tail is kTRUE (default)-
 the algorithm returns such x0, that
   P(x < x0)=p
 upper tail (lower_tail is kFALSE)- the algorithm returns such x0, that
   P(x > x0)=p
 the algorithm was taken from
   G.W.Hill, "Algorithm 396, Student's t-quantiles"
             "Communications of the ACM", 13(10), October 1970

Returns the value of the Vavilov density function
Parameters: 1st - the point were the density function is evaluated
            2nd - value of kappa (distribution parameter)
            3rd - value of beta2 (distribution parameter)
The algorithm was taken from the CernLib function vavden(G115)
Reference: A.Rotondi and P.Montagna, Fast Calculation of Vavilov distribution
Nucl.Instr. and Meth. B47(1990), 215-224
Accuracy: quote from the reference above:
"The resuls of our code have been compared with the values of the Vavilov
density function computed numerically in an accurate way: our approximation
shows a difference of less than 3% around the peak of the density function, slowly
increasing going towards the extreme tails to the right and to the left"

Returns the value of the Vavilov distribution function
Parameters: 1st - the point were the density function is evaluated
            2nd - value of kappa (distribution parameter)
            3rd - value of beta2 (distribution parameter)
The algorithm was taken from the CernLib function vavden(G115)
Reference: A.Rotondi and P.Montagna, Fast Calculation of Vavilov distribution
Nucl.Instr. and Meth. B47(1990), 215-224
Accuracy: quote from the reference above:
"The resuls of our code have been compared with the values of the Vavilov
density function computed numerically in an accurate way: our approximation
shows a difference of less than 3% around the peak of the density function, slowly
increasing going towards the extreme tails to the right and to the left"

Returns the value of the Landau distribution function at point x.
The algorithm was taken from the Cernlib function dislan(G110)
Reference: K.S.Kolbig and B.Schorr, "A program package for the Landau
distribution", Computer Phys.Comm., 31(1984), 97-111

NextPrime is used by the Core classes.