context: /*************************************************************************
  ALGLIB 3.11.0 (source code generated 2017-05-11)
  Copyright (c) Sergey Bochkanov (ALGLIB project).
  
  >>> SOURCE LICENSE >>>
  This program is free software; you can redistribute it and/or modify
  it under the terms of the GNU General Public License as published by
  the Free Software Foundation (www.fsf.org); either version 2 of the 
  License, or (at your option) any later version.
  
  This program is distributed in the hope that it will be useful,
  but WITHOUT ANY WARRANTY; without even the implied warranty of
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  GNU General Public License for more details.
  
  A copy of the GNU General Public License is available at
  http://www.fsf.org/licensing/licenses
  >>> END OF LICENSE >>>
  *************************************************************************/
  #pragma warning disable 162
  #pragma warning disable 164
  #pragma warning disable 219
  using System;
  
  public partial class alglib
  {
  
      
      /*************************************************************************
      Gamma function
  
      Input parameters:
          X   -   argument
  
      Domain:
          0 < X < 171.6
          -170 < X < 0, X is not an integer.
  
      Relative error:
       arithmetic   domain     # trials      peak         rms
          IEEE    -170,-33      20000       2.3e-15     3.3e-16
          IEEE     -33,  33     20000       9.4e-16     2.2e-16
          IEEE      33, 171.6   20000       2.3e-15     3.2e-16
  
      Cephes Math Library Release 2.8:  June, 2000
      Original copyright 1984, 1987, 1989, 1992, 2000 by Stephen L. Moshier
      Translated to AlgoPascal by Bochkanov Sergey (2005, 2006, 2007).
      *************************************************************************/
      public static double gammafunction(double x)
      {
      
          double result = gammafunc.gammafunction(x);
          return result;
      }
      
      /*************************************************************************
      Natural logarithm of gamma function
  
      Input parameters:
          X       -   argument
  
      Result:
          logarithm of the absolute value of the Gamma(X).
  
      Output parameters:
          SgnGam  -   sign(Gamma(X))
  
      Domain:
          0 < X < 2.55e305
          -2.55e305 < X < 0, X is not an integer.
  
      ACCURACY:
      arithmetic      domain        # trials     peak         rms
         IEEE    0, 3                 28000     5.4e-16     1.1e-16
         IEEE    2.718, 2.556e305     40000     3.5e-16     8.3e-17
      The error criterion was relative when the function magnitude
      was greater than one but absolute when it was less than one.
  
      The following test used the relative error criterion, though
      at certain points the relative error could be much higher than
      indicated.
         IEEE    -200, -4             10000     4.8e-16     1.3e-16
  
      Cephes Math Library Release 2.8:  June, 2000
      Copyright 1984, 1987, 1989, 1992, 2000 by Stephen L. Moshier
      Translated to AlgoPascal by Bochkanov Sergey (2005, 2006, 2007).
      *************************************************************************/
      public static double lngamma(double x, out double sgngam)
      {
          sgngam = 0;
          double result = gammafunc.lngamma(x, ref sgngam);
          return result;
      }
  
  }
  public partial class alglib
  {
  
      
      /*************************************************************************
      Error function
  
      The integral is
  
                                x
                                 -
                      2         | |          2
        erf(x)  =  --------     |    exp( - t  ) dt.
                   sqrt(pi)   | |
                               -
                                0
  
      For 0 <= |x| < 1, erf(x) = x * P4(x**2)/Q5(x**2); otherwise
      erf(x) = 1 - erfc(x).
  
  
      ACCURACY:
  
                           Relative error:
      arithmetic   domain     # trials      peak         rms
         IEEE      0,1         30000       3.7e-16     1.0e-16
  
      Cephes Math Library Release 2.8:  June, 2000
      Copyright 1984, 1987, 1988, 1992, 2000 by Stephen L. Moshier
      *************************************************************************/
      public static double errorfunction(double x)
      {
      
          double result = normaldistr.errorfunction(x);
          return result;
      }
      
      /*************************************************************************
      Complementary error function
  
       1 - erf(x) =
  
                                inf.
                                  -
                       2         | |          2
        erfc(x)  =  --------     |    exp( - t  ) dt
                    sqrt(pi)   | |
                                -
                                 x
  
  
      For small x, erfc(x) = 1 - erf(x); otherwise rational
      approximations are computed.
  
  
      ACCURACY:
  
                           Relative error:
      arithmetic   domain     # trials      peak         rms
         IEEE      0,26.6417   30000       5.7e-14     1.5e-14
  
      Cephes Math Library Release 2.8:  June, 2000
      Copyright 1984, 1987, 1988, 1992, 2000 by Stephen L. Moshier
      *************************************************************************/
      public static double errorfunctionc(double x)
      {
      
          double result = normaldistr.errorfunctionc(x);
          return result;
      }
      
      /*************************************************************************
      Normal distribution function
  
      Returns the area under the Gaussian probability density
      function, integrated from minus infinity to x:
  
                                 x
                                  -
                        1        | |          2
         ndtr(x)  = ---------    |    exp( - t /2 ) dt
                    sqrt(2pi)  | |
                                -
                               -inf.
  
                  =  ( 1 + erf(z) ) / 2
                  =  erfc(z) / 2
  
      where z = x/sqrt(2). Computation is via the functions
      erf and erfc.
  
  
      ACCURACY:
  
                           Relative error:
      arithmetic   domain     # trials      peak         rms
         IEEE     -13,0        30000       3.4e-14     6.7e-15
  
      Cephes Math Library Release 2.8:  June, 2000
      Copyright 1984, 1987, 1988, 1992, 2000 by Stephen L. Moshier
      *************************************************************************/
      public static double normaldistribution(double x)
      {
      
          double result = normaldistr.normaldistribution(x);
          return result;
      }
      
      /*************************************************************************
      Inverse of the error function
  
      Cephes Math Library Release 2.8:  June, 2000
      Copyright 1984, 1987, 1988, 1992, 2000 by Stephen L. Moshier
      *************************************************************************/
      public static double inverf(double e)
      {
      
          double result = normaldistr.inverf(e);
          return result;
      }
      
      /*************************************************************************
      Inverse of Normal distribution function
  
      Returns the argument, x, for which the area under the
      Gaussian probability density function (integrated from
      minus infinity to x) is equal to y.
  
  
      For small arguments 0 < y < exp(-2), the program computes
      z = sqrt( -2.0 * log(y) );  then the approximation is
      x = z - log(z)/z  - (1/z) P(1/z) / Q(1/z).
      There are two rational functions P/Q, one for 0 < y < exp(-32)
      and the other for y up to exp(-2).  For larger arguments,
      w = y - 0.5, and  x/sqrt(2pi) = w + w**3 R(w**2)/S(w**2)).
  
      ACCURACY:
  
                           Relative error:
      arithmetic   domain        # trials      peak         rms
         IEEE     0.125, 1        20000       7.2e-16     1.3e-16
         IEEE     3e-308, 0.135   50000       4.6e-16     9.8e-17
  
      Cephes Math Library Release 2.8:  June, 2000
      Copyright 1984, 1987, 1988, 1992, 2000 by Stephen L. Moshier
      *************************************************************************/
      public static double invnormaldistribution(double y0)
      {
      
          double result = normaldistr.invnormaldistribution(y0);
          return result;
      }
  
  }
  public partial class alglib
  {
  
      
      /*************************************************************************
      Incomplete gamma integral
  
      The function is defined by
  
                                x
                                 -
                        1       | |  -t  a-1
       igam(a,x)  =   -----     |   e   t   dt.
                       -      | |
                      | (a)    -
                                0
  
  
      In this implementation both arguments must be positive.
      The integral is evaluated by either a power series or
      continued fraction expansion, depending on the relative
      values of a and x.
  
      ACCURACY:
  
                           Relative error:
      arithmetic   domain     # trials      peak         rms
         IEEE      0,30       200000       3.6e-14     2.9e-15
         IEEE      0,100      300000       9.9e-14     1.5e-14
  
      Cephes Math Library Release 2.8:  June, 2000
      Copyright 1985, 1987, 2000 by Stephen L. Moshier
      *************************************************************************/
      public static double incompletegamma(double a, double x)
      {
      
          double result = igammaf.incompletegamma(a, x);
          return result;
      }
      
      /*************************************************************************
      Complemented incomplete gamma integral
  
      The function is defined by
  
  
       igamc(a,x)   =   1 - igam(a,x)
  
                                 inf.
                                   -
                          1       | |  -t  a-1
                    =   -----     |   e   t   dt.
                         -      | |
                        | (a)    -
                                  x
  
  
      In this implementation both arguments must be positive.
      The integral is evaluated by either a power series or
      continued fraction expansion, depending on the relative
      values of a and x.
  
      ACCURACY:
  
      Tested at random a, x.
                     a         x                      Relative error:
      arithmetic   domain   domain     # trials      peak         rms
         IEEE     0.5,100   0,100      200000       1.9e-14     1.7e-15
         IEEE     0.01,0.5  0,100      200000       1.4e-13     1.6e-15
  
      Cephes Math Library Release 2.8:  June, 2000
      Copyright 1985, 1987, 2000 by Stephen L. Moshier
      *************************************************************************/
      public static double incompletegammac(double a, double x)
      {
      
          double result = igammaf.incompletegammac(a, x);
          return result;
      }
      
      /*************************************************************************
      Inverse of complemented imcomplete gamma integral
  
      Given p, the function finds x such that
  
       igamc( a, x ) = p.
  
      Starting with the approximate value
  
              3
       x = a t
  
       where
  
       t = 1 - d - ndtri(p) sqrt(d)
  
      and
  
       d = 1/9a,
  
      the routine performs up to 10 Newton iterations to find the
      root of igamc(a,x) - p = 0.
  
      ACCURACY:
  
      Tested at random a, p in the intervals indicated.
  
                     a        p                      Relative error:
      arithmetic   domain   domain     # trials      peak         rms
         IEEE     0.5,100   0,0.5       100000       1.0e-14     1.7e-15
         IEEE     0.01,0.5  0,0.5       100000       9.0e-14     3.4e-15
         IEEE    0.5,10000  0,0.5        20000       2.3e-13     3.8e-14
  
      Cephes Math Library Release 2.8:  June, 2000
      Copyright 1984, 1987, 1995, 2000 by Stephen L. Moshier
      *************************************************************************/
      public static double invincompletegammac(double a, double y0)
      {
      
          double result = igammaf.invincompletegammac(a, y0);
          return result;
      }
  
  }
  public partial class alglib
  {
  
      
      /*************************************************************************
      Complete elliptic integral of the first kind
  
      Approximates the integral
  
  
  
                 pi/2
                  -
                 | |
                 |           dt
      K(m)  =    |    ------------------
                 |                   2
               | |    sqrt( 1 - m sin t )
                -
                 0
  
      using the approximation
  
          P(x)  -  log x Q(x).
  
      ACCURACY:
  
                           Relative error:
      arithmetic   domain     # trials      peak         rms
         IEEE       0,1        30000       2.5e-16     6.8e-17
  
      Cephes Math Library, Release 2.8:  June, 2000
      Copyright 1984, 1987, 2000 by Stephen L. Moshier
      *************************************************************************/
      public static double ellipticintegralk(double m)
      {
      
          double result = elliptic.ellipticintegralk(m);
          return result;
      }
      
      /*************************************************************************
      Complete elliptic integral of the first kind
  
      Approximates the integral
  
  
  
                 pi/2
                  -
                 | |
                 |           dt
      K(m)  =    |    ------------------
                 |                   2
               | |    sqrt( 1 - m sin t )
                -
                 0
  
      where m = 1 - m1, using the approximation
  
          P(x)  -  log x Q(x).
  
      The argument m1 is used rather than m so that the logarithmic
      singularity at m = 1 will be shifted to the origin; this
      preserves maximum accuracy.
  
      K(0) = pi/2.
  
      ACCURACY:
  
                           Relative error:
      arithmetic   domain     # trials      peak         rms
         IEEE       0,1        30000       2.5e-16     6.8e-17
  
      Cephes Math Library, Release 2.8:  June, 2000
      Copyright 1984, 1987, 2000 by Stephen L. Moshier
      *************************************************************************/
      public static double ellipticintegralkhighprecision(double m1)
      {
      
          double result = elliptic.ellipticintegralkhighprecision(m1);
          return result;
      }
      
      /*************************************************************************
      Incomplete elliptic integral of the first kind F(phi|m)
  
      Approximates the integral
  
  
  
                     phi
                      -
                     | |
                     |           dt
      F(phi_\m)  =    |    ------------------
                     |                   2
                   | |    sqrt( 1 - m sin t )
                    -
                     0
  
      of amplitude phi and modulus m, using the arithmetic -
      geometric mean algorithm.
  
  
  
  
      ACCURACY:
  
      Tested at random points with m in [0, 1] and phi as indicated.
  
                           Relative error:
      arithmetic   domain     # trials      peak         rms
         IEEE     -10,10       200000      7.4e-16     1.0e-16
  
      Cephes Math Library Release 2.8:  June, 2000
      Copyright 1984, 1987, 2000 by Stephen L. Moshier
      *************************************************************************/
      public static double incompleteellipticintegralk(double phi, double m)
      {
      
          double result = elliptic.incompleteellipticintegralk(phi, m);
          return result;
      }
      
      /*************************************************************************
      Complete elliptic integral of the second kind
  
      Approximates the integral
  
  
                 pi/2
                  -
                 | |                 2
      E(m)  =    |    sqrt( 1 - m sin t ) dt
               | |
                -
                 0
  
      using the approximation
  
           P(x)  -  x log x Q(x).
  
      ACCURACY:
  
                           Relative error:
      arithmetic   domain     # trials      peak         rms
         IEEE       0, 1       10000       2.1e-16     7.3e-17
  
      Cephes Math Library, Release 2.8: June, 2000
      Copyright 1984, 1987, 1989, 2000 by Stephen L. Moshier
      *************************************************************************/
      public static double ellipticintegrale(double m)
      {
      
          double result = elliptic.ellipticintegrale(m);
          return result;
      }
      
      /*************************************************************************
      Incomplete elliptic integral of the second kind
  
      Approximates the integral
  
  
                     phi
                      -
                     | |
                     |                   2
      E(phi_\m)  =    |    sqrt( 1 - m sin t ) dt
                     |
                   | |
                    -
                     0
  
      of amplitude phi and modulus m, using the arithmetic -
      geometric mean algorithm.
  
      ACCURACY:
  
      Tested at random arguments with phi in [-10, 10] and m in
      [0, 1].
                           Relative error:
      arithmetic   domain     # trials      peak         rms
         IEEE     -10,10      150000       3.3e-15     1.4e-16
  
      Cephes Math Library Release 2.8:  June, 2000
      Copyright 1984, 1987, 1993, 2000 by Stephen L. Moshier
      *************************************************************************/
      public static double incompleteellipticintegrale(double phi, double m)
      {
      
          double result = elliptic.incompleteellipticintegrale(phi, m);
          return result;
      }
  
  }
  public partial class alglib
  {
  
      
      /*************************************************************************
      Calculation of the value of the Hermite polynomial.
  
      Parameters:
          n   -   degree, n>=0
          x   -   argument
  
      Result:
          the value of the Hermite polynomial Hn at x
      *************************************************************************/
      public static double hermitecalculate(int n, double x)
      {
      
          double result = hermite.hermitecalculate(n, x);
          return result;
      }
      
      /*************************************************************************
      Summation of Hermite polynomials using Clenshaw's recurrence formula.
  
      This routine calculates
          c[0]*H0(x) + c[1]*H1(x) + ... + c[N]*HN(x)
  
      Parameters:
          n   -   degree, n>=0
          x   -   argument
  
      Result:
          the value of the Hermite polynomial at x
      *************************************************************************/
      public static double hermitesum(double[] c, int n, double x)
      {
      
          double result = hermite.hermitesum(c, n, x);
          return result;
      }
      
      /*************************************************************************
      Representation of Hn as C[0] + C[1]*X + ... + C[N]*X^N
  
      Input parameters:
          N   -   polynomial degree, n>=0
  
      Output parameters:
          C   -   coefficients
      *************************************************************************/
      public static void hermitecoefficients(int n, out double[] c)
      {
          c = new double[0];
          hermite.hermitecoefficients(n, ref c);
          return;
      }
  
  }
  public partial class alglib
  {
  
      
      /*************************************************************************
      Dawson's Integral
  
      Approximates the integral
  
                                  x
                                  -
                           2     | |        2
       dawsn(x)  =  exp( -x  )   |    exp( t  ) dt
                               | |
                                -
                                0
  
      Three different rational approximations are employed, for
      the intervals 0 to 3.25; 3.25 to 6.25; and 6.25 up.
  
      ACCURACY:
  
                           Relative error:
      arithmetic   domain     # trials      peak         rms
         IEEE      0,10        10000       6.9e-16     1.0e-16
  
      Cephes Math Library Release 2.8:  June, 2000
      Copyright 1984, 1987, 1989, 2000 by Stephen L. Moshier
      *************************************************************************/
      public static double dawsonintegral(double x)
      {
      
          double result = dawson.dawsonintegral(x);
          return result;
      }
  
  }
  public partial class alglib
  {
  
      
      /*************************************************************************
      Sine and cosine integrals
  
      Evaluates the integrals
  
                               x
                               -
                              |  cos t - 1
        Ci(x) = eul + ln x +  |  --------- dt,
                              |      t
                             -
                              0
                  x
                  -
                 |  sin t
        Si(x) =  |  ----- dt
                 |    t
                -
                 0
  
      where eul = 0.57721566490153286061 is Euler's constant.
      The integrals are approximated by rational functions.
      For x > 8 auxiliary functions f(x) and g(x) are employed
      such that
  
      Ci(x) = f(x) sin(x) - g(x) cos(x)
      Si(x) = pi/2 - f(x) cos(x) - g(x) sin(x)
  
  
      ACCURACY:
         Test interval = [0,50].
      Absolute error, except relative when > 1:
      arithmetic   function   # trials      peak         rms
         IEEE        Si        30000       4.4e-16     7.3e-17
         IEEE        Ci        30000       6.9e-16     5.1e-17
  
      Cephes Math Library Release 2.1:  January, 1989
      Copyright 1984, 1987, 1989 by Stephen L. Moshier
      *************************************************************************/
      public static void sinecosineintegrals(double x, out double si, out double ci)
      {
          si = 0;
          ci = 0;
          trigintegrals.sinecosineintegrals(x, ref si, ref ci);
          return;
      }
      
      /*************************************************************************
      Hyperbolic sine and cosine integrals
  
      Approximates the integrals
  
                                 x
                                 -
                                | |   cosh t - 1
        Chi(x) = eul + ln x +   |    -----------  dt,
                              | |          t
                               -
                               0
  
                    x
                    -
                   | |  sinh t
        Shi(x) =   |    ------  dt
                 | |       t
                  -
                  0
  
      where eul = 0.57721566490153286061 is Euler's constant.
      The integrals are evaluated by power series for x < 8
      and by Chebyshev expansions for x between 8 and 88.
      For large x, both functions approach exp(x)/2x.
      Arguments greater than 88 in magnitude return MAXNUM.
  
  
      ACCURACY:
  
      Test interval 0 to 88.
                           Relative error:
      arithmetic   function  # trials      peak         rms
         IEEE         Shi      30000       6.9e-16     1.6e-16
             Absolute error, except relative when |Chi| > 1:
         IEEE         Chi      30000       8.4e-16     1.4e-16
  
      Cephes Math Library Release 2.8:  June, 2000
      Copyright 1984, 1987, 2000 by Stephen L. Moshier
      *************************************************************************/
      public static void hyperbolicsinecosineintegrals(double x, out double shi, out double chi)
      {
          shi = 0;
          chi = 0;
          trigintegrals.hyperbolicsinecosineintegrals(x, ref shi, ref chi);
          return;
      }
  
  }
  public partial class alglib
  {
  
      
      /*************************************************************************
      Poisson distribution
  
      Returns the sum of the first k+1 terms of the Poisson
      distribution:
  
        k         j
        --   -m  m
        >   e    --
        --       j!
       j=0
  
      The terms are not summed directly; instead the incomplete
      gamma integral is employed, according to the relation
  
      y = pdtr( k, m ) = igamc( k+1, m ).
  
      The arguments must both be positive.
      ACCURACY:
  
      See incomplete gamma function
  
      Cephes Math Library Release 2.8:  June, 2000
      Copyright 1984, 1987, 1995, 2000 by Stephen L. Moshier
      *************************************************************************/
      public static double poissondistribution(int k, double m)
      {
      
          double result = poissondistr.poissondistribution(k, m);
          return result;
      }
      
      /*************************************************************************
      Complemented Poisson distribution
  
      Returns the sum of the terms k+1 to infinity of the Poisson
      distribution:
  
       inf.       j
        --   -m  m
        >   e    --
        --       j!
       j=k+1
  
      The terms are not summed directly; instead the incomplete
      gamma integral is employed, according to the formula
  
      y = pdtrc( k, m ) = igam( k+1, m ).
  
      The arguments must both be positive.
  
      ACCURACY:
  
      See incomplete gamma function
  
      Cephes Math Library Release 2.8:  June, 2000
      Copyright 1984, 1987, 1995, 2000 by Stephen L. Moshier
      *************************************************************************/
      public static double poissoncdistribution(int k, double m)
      {
      
          double result = poissondistr.poissoncdistribution(k, m);
          return result;
      }
      
      /*************************************************************************
      Inverse Poisson distribution
  
      Finds the Poisson variable x such that the integral
      from 0 to x of the Poisson density is equal to the
      given probability y.
  
      This is accomplished using the inverse gamma integral
      function and the relation
  
         m = igami( k+1, y ).
  
      ACCURACY:
  
      See inverse incomplete gamma function
  
      Cephes Math Library Release 2.8:  June, 2000
      Copyright 1984, 1987, 1995, 2000 by Stephen L. Moshier
      *************************************************************************/
      public static double invpoissondistribution(int k, double y)
      {
      
          double result = poissondistr.invpoissondistribution(k, y);
          return result;
      }
  
  }
  public partial class alglib
  {
  
      
      /*************************************************************************
      Bessel function of order zero
  
      Returns Bessel function of order zero of the argument.
  
      The domain is divided into the intervals [0, 5] and
      (5, infinity). In the first interval the following rational
      approximation is used:
  
  
             2         2
      (w - r  ) (w - r  ) P (w) / Q (w)
            1         2    3       8
  
                 2
      where w = x  and the two r's are zeros of the function.
  
      In the second interval, the Hankel asymptotic expansion
      is employed with two rational functions of degree 6/6
      and 7/7.
  
      ACCURACY:
  
                           Absolute error:
      arithmetic   domain     # trials      peak         rms
         IEEE      0, 30       60000       4.2e-16     1.1e-16
  
      Cephes Math Library Release 2.8:  June, 2000
      Copyright 1984, 1987, 1989, 2000 by Stephen L. Moshier
      *************************************************************************/
      public static double besselj0(double x)
      {
      
          double result = bessel.besselj0(x);
          return result;
      }
      
      /*************************************************************************
      Bessel function of order one
  
      Returns Bessel function of order one of the argument.
  
      The domain is divided into the intervals [0, 8] and
      (8, infinity). In the first interval a 24 term Chebyshev
      expansion is used. In the second, the asymptotic
      trigonometric representation is employed using two
      rational functions of degree 5/5.
  
      ACCURACY:
  
                           Absolute error:
      arithmetic   domain      # trials      peak         rms
         IEEE      0, 30       30000       2.6e-16     1.1e-16
  
      Cephes Math Library Release 2.8:  June, 2000
      Copyright 1984, 1987, 1989, 2000 by Stephen L. Moshier
      *************************************************************************/
      public static double besselj1(double x)
      {
      
          double result = bessel.besselj1(x);
          return result;
      }
      
      /*************************************************************************
      Bessel function of integer order
  
      Returns Bessel function of order n, where n is a
      (possibly negative) integer.
  
      The ratio of jn(x) to j0(x) is computed by backward
      recurrence.  First the ratio jn/jn-1 is found by a
      continued fraction expansion.  Then the recurrence
      relating successive orders is applied until j0 or j1 is
      reached.
  
      If n = 0 or 1 the routine for j0 or j1 is called
      directly.
  
      ACCURACY:
  
                           Absolute error:
      arithmetic   range      # trials      peak         rms
         IEEE      0, 30        5000       4.4e-16     7.9e-17
  
  
      Not suitable for large n or x. Use jv() (fractional order) instead.
  
      Cephes Math Library Release 2.8:  June, 2000
      Copyright 1984, 1987, 2000 by Stephen L. Moshier
      *************************************************************************/
      public static double besseljn(int n, double x)
      {
      
          double result = bessel.besseljn(n, x);
          return result;
      }
      
      /*************************************************************************
      Bessel function of the second kind, order zero
  
      Returns Bessel function of the second kind, of order
      zero, of the argument.
  
      The domain is divided into the intervals [0, 5] and
      (5, infinity). In the first interval a rational approximation
      R(x) is employed to compute
        y0(x)  = R(x)  +   2 * log(x) * j0(x) / PI.
      Thus a call to j0() is required.
  
      In the second interval, the Hankel asymptotic expansion
      is employed with two rational functions of degree 6/6
      and 7/7.
  
  
  
      ACCURACY:
  
       Absolute error, when y0(x) < 1; else relative error:
  
      arithmetic   domain     # trials      peak         rms
         IEEE      0, 30       30000       1.3e-15     1.6e-16
  
      Cephes Math Library Release 2.8:  June, 2000
      Copyright 1984, 1987, 1989, 2000 by Stephen L. Moshier
      *************************************************************************/
      public static double bessely0(double x)
      {
      
          double result = bessel.bessely0(x);
          return result;
      }
      
      /*************************************************************************
      Bessel function of second kind of order one
  
      Returns Bessel function of the second kind of order one
      of the argument.
  
      The domain is divided into the intervals [0, 8] and
      (8, infinity). In the first interval a 25 term Chebyshev
      expansion is used, and a call to j1() is required.
      In the second, the asymptotic trigonometric representation
      is employed using two rational functions of degree 5/5.
  
      ACCURACY:
  
                           Absolute error:
      arithmetic   domain      # trials      peak         rms
         IEEE      0, 30       30000       1.0e-15     1.3e-16
  
      Cephes Math Library Release 2.8:  June, 2000
      Copyright 1984, 1987, 1989, 2000 by Stephen L. Moshier
      *************************************************************************/
      public static double bessely1(double x)
      {
      
          double result = bessel.bessely1(x);
          return result;
      }
      
      /*************************************************************************
      Bessel function of second kind of integer order
  
      Returns Bessel function of order n, where n is a
      (possibly negative) integer.
  
      The function is evaluated by forward recurrence on
      n, starting with values computed by the routines
      y0() and y1().
  
      If n = 0 or 1 the routine for y0 or y1 is called
      directly.
  
      ACCURACY:
                           Absolute error, except relative
                           when y > 1:
      arithmetic   domain     # trials      peak         rms
         IEEE      0, 30       30000       3.4e-15     4.3e-16
  
      Cephes Math Library Release 2.8:  June, 2000
      Copyright 1984, 1987, 2000 by Stephen L. Moshier
      *************************************************************************/
      public static double besselyn(int n, double x)
      {
      
          double result = bessel.besselyn(n, x);
          return result;
      }
      
      /*************************************************************************
      Modified Bessel function of order zero
  
      Returns modified Bessel function of order zero of the
      argument.
  
      The function is defined as i0(x) = j0( ix ).
  
      The range is partitioned into the two intervals [0,8] and
      (8, infinity).  Chebyshev polynomial expansions are employed
      in each interval.
  
      ACCURACY:
  
                           Relative error:
      arithmetic   domain     # trials      peak         rms
         IEEE      0,30        30000       5.8e-16     1.4e-16
  
      Cephes Math Library Release 2.8:  June, 2000
      Copyright 1984, 1987, 2000 by Stephen L. Moshier
      *************************************************************************/
      public static double besseli0(double x)
      {
      
          double result = bessel.besseli0(x);
          return result;
      }
      
      /*************************************************************************
      Modified Bessel function of order one
  
      Returns modified Bessel function of order one of the
      argument.
  
      The function is defined as i1(x) = -i j1( ix ).
  
      The range is partitioned into the two intervals [0,8] and
      (8, infinity).  Chebyshev polynomial expansions are employed
      in each interval.
  
      ACCURACY:
  
                           Relative error:
      arithmetic   domain     # trials      peak         rms
         IEEE      0, 30       30000       1.9e-15     2.1e-16
  
      Cephes Math Library Release 2.8:  June, 2000
      Copyright 1985, 1987, 2000 by Stephen L. Moshier
      *************************************************************************/
      public static double besseli1(double x)
      {
      
          double result = bessel.besseli1(x);
          return result;
      }
      
      /*************************************************************************
      Modified Bessel function, second kind, order zero
  
      Returns modified Bessel function of the second kind
      of order zero of the argument.
  
      The range is partitioned into the two intervals [0,8] and
      (8, infinity).  Chebyshev polynomial expansions are employed
      in each interval.
  
      ACCURACY:
  
      Tested at 2000 random points between 0 and 8.  Peak absolute
      error (relative when K0 > 1) was 1.46e-14; rms, 4.26e-15.
                           Relative error:
      arithmetic   domain     # trials      peak         rms
         IEEE      0, 30       30000       1.2e-15     1.6e-16
  
      Cephes Math Library Release 2.8:  June, 2000
      Copyright 1984, 1987, 2000 by Stephen L. Moshier
      *************************************************************************/
      public static double besselk0(double x)
      {
      
          double result = bessel.besselk0(x);
          return result;
      }
      
      /*************************************************************************
      Modified Bessel function, second kind, order one
  
      Computes the modified Bessel function of the second kind
      of order one of the argument.
  
      The range is partitioned into the two intervals [0,2] and
      (2, infinity).  Chebyshev polynomial expansions are employed
      in each interval.
  
      ACCURACY:
  
                           Relative error:
      arithmetic   domain     # trials      peak         rms
         IEEE      0, 30       30000       1.2e-15     1.6e-16
  
      Cephes Math Library Release 2.8:  June, 2000
      Copyright 1984, 1987, 2000 by Stephen L. Moshier
      *************************************************************************/
      public static double besselk1(double x)
      {
      
          double result = bessel.besselk1(x);
          return result;
      }
      
      /*************************************************************************
      Modified Bessel function, second kind, integer order
  
      Returns modified Bessel function of the second kind
      of order n of the argument.
  
      The range is partitioned into the two intervals [0,9.55] and
      (9.55, infinity).  An ascending power series is used in the
      low range, and an asymptotic expansion in the high range.
  
      ACCURACY:
  
                           Relative error:
      arithmetic   domain     # trials      peak         rms
         IEEE      0,30        90000       1.8e-8      3.0e-10
  
      Error is high only near the crossover point x = 9.55
      between the two expansions used.
  
      Cephes Math Library Release 2.8:  June, 2000
      Copyright 1984, 1987, 1988, 2000 by Stephen L. Moshier
      *************************************************************************/
      public static double besselkn(int nn, double x)
      {
      
          double result = bessel.besselkn(nn, x);
          return result;
      }
  
  }
  public partial class alglib
  {
  
      
      /*************************************************************************
      Incomplete beta integral
  
      Returns incomplete beta integral of the arguments, evaluated
      from zero to x.  The function is defined as
  
                       x
          -            -
         | (a+b)      | |  a-1     b-1
       -----------    |   t   (1-t)   dt.
        -     -     | |
       | (a) | (b)   -
                      0
  
      The domain of definition is 0 <= x <= 1.  In this
      implementation a and b are restricted to positive values.
      The integral from x to 1 may be obtained by the symmetry
      relation
  
         1 - incbet( a, b, x )  =  incbet( b, a, 1-x ).
  
      The integral is evaluated by a continued fraction expansion
      or, when b*x is small, by a power series.
  
      ACCURACY:
  
      Tested at uniformly distributed random points (a,b,x) with a and b
      in "domain" and x between 0 and 1.
                                             Relative error
      arithmetic   domain     # trials      peak         rms
         IEEE      0,5         10000       6.9e-15     4.5e-16
         IEEE      0,85       250000       2.2e-13     1.7e-14
         IEEE      0,1000      30000       5.3e-12     6.3e-13
         IEEE      0,10000    250000       9.3e-11     7.1e-12
         IEEE      0,100000    10000       8.7e-10     4.8e-11
      Outputs smaller than the IEEE gradual underflow threshold
      were excluded from these statistics.
  
      Cephes Math Library, Release 2.8:  June, 2000
      Copyright 1984, 1995, 2000 by Stephen L. Moshier
      *************************************************************************/
      public static double incompletebeta(double a, double b, double x)
      {
      
          double result = ibetaf.incompletebeta(a, b, x);
          return result;
      }
      
      /*************************************************************************
      Inverse of imcomplete beta integral
  
      Given y, the function finds x such that
  
       incbet( a, b, x ) = y .
  
      The routine performs interval halving or Newton iterations to find the
      root of incbet(a,b,x) - y = 0.
  
  
      ACCURACY:
  
                           Relative error:
                     x     a,b
      arithmetic   domain  domain  # trials    peak       rms
         IEEE      0,1    .5,10000   50000    5.8e-12   1.3e-13
         IEEE      0,1   .25,100    100000    1.8e-13   3.9e-15
         IEEE      0,1     0,5       50000    1.1e-12   5.5e-15
      With a and b constrained to half-integer or integer values:
         IEEE      0,1    .5,10000   50000    5.8e-12   1.1e-13
         IEEE      0,1    .5,100    100000    1.7e-14   7.9e-16
      With a = .5, b constrained to half-integer or integer values:
         IEEE      0,1    .5,10000   10000    8.3e-11   1.0e-11
  
      Cephes Math Library Release 2.8:  June, 2000
      Copyright 1984, 1996, 2000 by Stephen L. Moshier
      *************************************************************************/
      public static double invincompletebeta(double a, double b, double y)
      {
      
          double result = ibetaf.invincompletebeta(a, b, y);
          return result;
      }
  
  }
  public partial class alglib
  {
  
      
      /*************************************************************************
      F distribution
  
      Returns the area from zero to x under the F density
      function (also known as Snedcor's density or the
      variance ratio density).  This is the density
      of x = (u1/df1)/(u2/df2), where u1 and u2 are random
      variables having Chi square distributions with df1
      and df2 degrees of freedom, respectively.
      The incomplete beta integral is used, according to the
      formula
  
      P(x) = incbet( df1/2, df2/2, (df1*x/(df2 + df1*x) ).
  
  
      The arguments a and b are greater than zero, and x is
      nonnegative.
  
      ACCURACY:
  
      Tested at random points (a,b,x).
  
                     x     a,b                     Relative error:
      arithmetic  domain  domain     # trials      peak         rms
         IEEE      0,1    0,100       100000      9.8e-15     1.7e-15
         IEEE      1,5    0,100       100000      6.5e-15     3.5e-16
         IEEE      0,1    1,10000     100000      2.2e-11     3.3e-12
         IEEE      1,5    1,10000     100000      1.1e-11     1.7e-13
  
      Cephes Math Library Release 2.8:  June, 2000
      Copyright 1984, 1987, 1995, 2000 by Stephen L. Moshier
      *************************************************************************/
      public static double fdistribution(int a, int b, double x)
      {
      
          double result = fdistr.fdistribution(a, b, x);
          return result;
      }
      
      /*************************************************************************
      Complemented F distribution
  
      Returns the area from x to infinity under the F density
      function (also known as Snedcor's density or the
      variance ratio density).
  
  
                           inf.
                            -
                   1       | |  a-1      b-1
      1-P(x)  =  ------    |   t    (1-t)    dt
                 B(a,b)  | |
                          -
                           x
  
  
      The incomplete beta integral is used, according to the
      formula
  
      P(x) = incbet( df2/2, df1/2, (df2/(df2 + df1*x) ).
  
  
      ACCURACY:
  
      Tested at random points (a,b,x) in the indicated intervals.
                     x     a,b                     Relative error:
      arithmetic  domain  domain     # trials      peak         rms
         IEEE      0,1    1,100       100000      3.7e-14     5.9e-16
         IEEE      1,5    1,100       100000      8.0e-15     1.6e-15
         IEEE      0,1    1,10000     100000      1.8e-11     3.5e-13
         IEEE      1,5    1,10000     100000      2.0e-11     3.0e-12
  
      Cephes Math Library Release 2.8:  June, 2000
      Copyright 1984, 1987, 1995, 2000 by Stephen L. Moshier
      *************************************************************************/
      public static double fcdistribution(int a, int b, double x)
      {
      
          double result = fdistr.fcdistribution(a, b, x);
          return result;
      }
      
      /*************************************************************************
      Inverse of complemented F distribution
  
      Finds the F density argument x such that the integral
      from x to infinity of the F density is equal to the
      given probability p.
  
      This is accomplished using the inverse beta integral
      function and the relations
  
           z = incbi( df2/2, df1/2, p )
           x = df2 (1-z) / (df1 z).
  
      Note: the following relations hold for the inverse of
      the uncomplemented F distribution:
  
           z = incbi( df1/2, df2/2, p )
           x = df2 z / (df1 (1-z)).
  
      ACCURACY:
  
      Tested at random points (a,b,p).
  
                   a,b                     Relative error:
      arithmetic  domain     # trials      peak         rms
       For p between .001 and 1:
         IEEE     1,100       100000      8.3e-15     4.7e-16
         IEEE     1,10000     100000      2.1e-11     1.4e-13
       For p between 10^-6 and 10^-3:
         IEEE     1,100        50000      1.3e-12     8.4e-15
         IEEE     1,10000      50000      3.0e-12     4.8e-14
  
      Cephes Math Library Release 2.8:  June, 2000
      Copyright 1984, 1987, 1995, 2000 by Stephen L. Moshier
      *************************************************************************/
      public static double invfdistribution(int a, int b, double y)
      {
      
          double result = fdistr.invfdistribution(a, b, y);
          return result;
      }
  
  }
  public partial class alglib
  {
  
      
      /*************************************************************************
      Fresnel integral
  
      Evaluates the Fresnel integrals
  
                x
                -
               | |
      C(x) =   |   cos(pi/2 t**2) dt,
             | |
              -
               0
  
                x
                -
               | |
      S(x) =   |   sin(pi/2 t**2) dt.
             | |
              -
               0
  
  
      The integrals are evaluated by a power series for x < 1.
      For x >= 1 auxiliary functions f(x) and g(x) are employed
      such that
  
      C(x) = 0.5 + f(x) sin( pi/2 x**2 ) - g(x) cos( pi/2 x**2 )
      S(x) = 0.5 - f(x) cos( pi/2 x**2 ) - g(x) sin( pi/2 x**2 )
  
  
  
      ACCURACY:
  
       Relative error.
  
      Arithmetic  function   domain     # trials      peak         rms
        IEEE       S(x)      0, 10       10000       2.0e-15     3.2e-16
        IEEE       C(x)      0, 10       10000       1.8e-15     3.3e-16
  
      Cephes Math Library Release 2.8:  June, 2000
      Copyright 1984, 1987, 1989, 2000 by Stephen L. Moshier
      *************************************************************************/
      public static void fresnelintegral(double x, ref double c, ref double s)
      {
      
          fresnel.fresnelintegral(x, ref c, ref s);
          return;
      }
  
  }
  public partial class alglib
  {
  
      
      /*************************************************************************
      Jacobian Elliptic Functions
  
      Evaluates the Jacobian elliptic functions sn(u|m), cn(u|m),
      and dn(u|m) of parameter m between 0 and 1, and real
      argument u.
  
      These functions are periodic, with quarter-period on the
      real axis equal to the complete elliptic integral
      ellpk(1.0-m).
  
      Relation to incomplete elliptic integral:
      If u = ellik(phi,m), then sn(u|m) = sin(phi),
      and cn(u|m) = cos(phi).  Phi is called the amplitude of u.
  
      Computation is by means of the arithmetic-geometric mean
      algorithm, except when m is within 1e-9 of 0 or 1.  In the
      latter case with m close to 1, the approximation applies
      only for phi < pi/2.
  
      ACCURACY:
  
      Tested at random points with u between 0 and 10, m between
      0 and 1.
  
                 Absolute error (* = relative error):
      arithmetic   function   # trials      peak         rms
         IEEE      phi         10000       9.2e-16*    1.4e-16*
         IEEE      sn          50000       4.1e-15     4.6e-16
         IEEE      cn          40000       3.6e-15     4.4e-16
         IEEE      dn          10000       1.3e-12     1.8e-14
  
       Peak error observed in consistency check using addition
      theorem for sn(u+v) was 4e-16 (absolute).  Also tested by
      the above relation to the incomplete elliptic integral.
      Accuracy deteriorates when u is large.
  
      Cephes Math Library Release 2.8:  June, 2000
      Copyright 1984, 1987, 2000 by Stephen L. Moshier
      *************************************************************************/
      public static void jacobianellipticfunctions(double u, double m, out double sn, out double cn, out double dn, out double ph)
      {
          sn = 0;
          cn = 0;
          dn = 0;
          ph = 0;
          jacobianelliptic.jacobianellipticfunctions(u, m, ref sn, ref cn, ref dn, ref ph);
          return;
      }
  
  }
  public partial class alglib
  {
  
      
      /*************************************************************************
      Psi (digamma) function
  
                   d      -
        psi(x)  =  -- ln | (x)
                   dx
  
      is the logarithmic derivative of the gamma function.
      For integer x,
                        n-1
                         -
      psi(n) = -EUL  +   >  1/k.
                         -
                        k=1
  
      This formula is used for 0 < n <= 10.  If x is negative, it
      is transformed to a positive argument by the reflection
      formula  psi(1-x) = psi(x) + pi cot(pi x).
      For general positive x, the argument is made greater than 10
      using the recurrence  psi(x+1) = psi(x) + 1/x.
      Then the following asymptotic expansion is applied:
  
                                inf.   B
                                 -      2k
      psi(x) = log(x) - 1/2x -   >   -------
                                 -        2k
                                k=1   2k x
  
      where the B2k are Bernoulli numbers.
  
      ACCURACY:
         Relative error (except absolute when |psi| < 1):
      arithmetic   domain     # trials      peak         rms
         IEEE      0,30        30000       1.3e-15     1.4e-16
         IEEE      -30,0       40000       1.5e-15     2.2e-16
  
      Cephes Math Library Release 2.8:  June, 2000
      Copyright 1984, 1987, 1992, 2000 by Stephen L. Moshier
      *************************************************************************/
      public static double psi(double x)
      {
      
          double result = psif.psi(x);
          return result;
      }
  
  }
  public partial class alglib
  {
  
      
      /*************************************************************************
      Exponential integral Ei(x)
  
                    x
                     -     t
                    | |   e
         Ei(x) =   -|-   ---  dt .
                  | |     t
                   -
                  -inf
  
      Not defined for x <= 0.
      See also expn.c.
  
  
  
      ACCURACY:
  
                           Relative error:
      arithmetic   domain     # trials      peak         rms
         IEEE       0,100       50000      8.6e-16     1.3e-16
  
      Cephes Math Library Release 2.8:  May, 1999
      Copyright 1999 by Stephen L. Moshier
      *************************************************************************/
      public static double exponentialintegralei(double x)
      {
      
          double result = expintegrals.exponentialintegralei(x);
          return result;
      }
      
      /*************************************************************************
      Exponential integral En(x)
  
      Evaluates the exponential integral
  
                      inf.
                        -
                       | |   -xt
                       |    e
           E (x)  =    |    ----  dt.
            n          |      n
                     | |     t
                      -
                       1
  
  
      Both n and x must be nonnegative.
  
      The routine employs either a power series, a continued
      fraction, or an asymptotic formula depending on the
      relative values of n and x.
  
      ACCURACY:
  
                           Relative error:
      arithmetic   domain     # trials      peak         rms
         IEEE      0, 30       10000       1.7e-15     3.6e-16
  
      Cephes Math Library Release 2.8:  June, 2000
      Copyright 1985, 2000 by Stephen L. Moshier
      *************************************************************************/
      public static double exponentialintegralen(double x, int n)
      {
      
          double result = expintegrals.exponentialintegralen(x, n);
          return result;
      }
  
  }
  public partial class alglib
  {
  
      
      /*************************************************************************
      Calculation of the value of the Laguerre polynomial.
  
      Parameters:
          n   -   degree, n>=0
          x   -   argument
  
      Result:
          the value of the Laguerre polynomial Ln at x
      *************************************************************************/
      public static double laguerrecalculate(int n, double x)
      {
      
          double result = laguerre.laguerrecalculate(n, x);
          return result;
      }
      
      /*************************************************************************
      Summation of Laguerre polynomials using Clenshaw's recurrence formula.
  
      This routine calculates c[0]*L0(x) + c[1]*L1(x) + ... + c[N]*LN(x)
  
      Parameters:
          n   -   degree, n>=0
          x   -   argument
  
      Result:
          the value of the Laguerre polynomial at x
      *************************************************************************/
      public static double laguerresum(double[] c, int n, double x)
      {
      
          double result = laguerre.laguerresum(c, n, x);
          return result;
      }
      
      /*************************************************************************
      Representation of Ln as C[0] + C[1]*X + ... + C[N]*X^N
  
      Input parameters:
          N   -   polynomial degree, n>=0
  
      Output parameters:
          C   -   coefficients
      *************************************************************************/
      public static void laguerrecoefficients(int n, out double[] c)
      {
          c = new double[0];
          laguerre.laguerrecoefficients(n, ref c);
          return;
      }
  
  }
  public partial class alglib
  {
  
      
      /*************************************************************************
      Chi-square distribution
  
      Returns the area under the left hand tail (from 0 to x)
      of the Chi square probability density function with
      v degrees of freedom.
  
  
                                        x
                                         -
                             1          | |  v/2-1  -t/2
       P( x | v )   =   -----------     |   t      e     dt
                         v/2  -       | |
                        2    | (v/2)   -
                                        0
  
      where x is the Chi-square variable.
  
      The incomplete gamma integral is used, according to the
      formula
  
      y = chdtr( v, x ) = igam( v/2.0, x/2.0 ).
  
      The arguments must both be positive.
  
      ACCURACY:
  
      See incomplete gamma function
  
  
      Cephes Math Library Release 2.8:  June, 2000
      Copyright 1984, 1987, 2000 by Stephen L. Moshier
      *************************************************************************/
      public static double chisquaredistribution(double v, double x)
      {
      
          double result = chisquaredistr.chisquaredistribution(v, x);
          return result;
      }
      
      /*************************************************************************
      Complemented Chi-square distribution
  
      Returns the area under the right hand tail (from x to
      infinity) of the Chi square probability density function
      with v degrees of freedom:
  
                                       inf.
                                         -
                             1          | |  v/2-1  -t/2
       P( x | v )   =   -----------     |   t      e     dt
                         v/2  -       | |
                        2    | (v/2)   -
                                        x
  
      where x is the Chi-square variable.
  
      The incomplete gamma integral is used, according to the
      formula
  
      y = chdtr( v, x ) = igamc( v/2.0, x/2.0 ).
  
      The arguments must both be positive.
  
      ACCURACY:
  
      See incomplete gamma function
  
      Cephes Math Library Release 2.8:  June, 2000
      Copyright 1984, 1987, 2000 by Stephen L. Moshier
      *************************************************************************/
      public static double chisquarecdistribution(double v, double x)
      {
      
          double result = chisquaredistr.chisquarecdistribution(v, x);
          return result;
      }
      
      /*************************************************************************
      Inverse of complemented Chi-square distribution
  
      Finds the Chi-square argument x such that the integral
      from x to infinity of the Chi-square density is equal
      to the given cumulative probability y.
  
      This is accomplished using the inverse gamma integral
      function and the relation
  
         x/2 = igami( df/2, y );
  
      ACCURACY:
  
      See inverse incomplete gamma function
  
  
      Cephes Math Library Release 2.8:  June, 2000
      Copyright 1984, 1987, 2000 by Stephen L. Moshier
      *************************************************************************/
      public static double invchisquaredistribution(double v, double y)
      {
      
          double result = chisquaredistr.invchisquaredistribution(v, y);
          return result;
      }
  
  }
  public partial class alglib
  {
  
      
      /*************************************************************************
      Calculation of the value of the Legendre polynomial Pn.
  
      Parameters:
          n   -   degree, n>=0
          x   -   argument
  
      Result:
          the value of the Legendre polynomial Pn at x
      *************************************************************************/
      public static double legendrecalculate(int n, double x)
      {
      
          double result = legendre.legendrecalculate(n, x);
          return result;
      }
      
      /*************************************************************************
      Summation of Legendre polynomials using Clenshaw's recurrence formula.
  
      This routine calculates
          c[0]*P0(x) + c[1]*P1(x) + ... + c[N]*PN(x)
  
      Parameters:
          n   -   degree, n>=0
          x   -   argument
  
      Result:
          the value of the Legendre polynomial at x
      *************************************************************************/
      public static double legendresum(double[] c, int n, double x)
      {
      
          double result = legendre.legendresum(c, n, x);
          return result;
      }
      
      /*************************************************************************
      Representation of Pn as C[0] + C[1]*X + ... + C[N]*X^N
  
      Input parameters:
          N   -   polynomial degree, n>=0
  
      Output parameters:
          C   -   coefficients
      *************************************************************************/
      public static void legendrecoefficients(int n, out double[] c)
      {
          c = new double[0];
          legendre.legendrecoefficients(n, ref c);
          return;
      }
  
  }
  public partial class alglib
  {
  
      
      /*************************************************************************
      Beta function
  
  
                        -     -
                       | (a) | (b)
      beta( a, b )  =  -----------.
                          -
                         | (a+b)
  
      For large arguments the logarithm of the function is
      evaluated using lgam(), then exponentiated.
  
      ACCURACY:
  
                           Relative error:
      arithmetic   domain     # trials      peak         rms
         IEEE       0,30       30000       8.1e-14     1.1e-14
  
      Cephes Math Library Release 2.0:  April, 1987
      Copyright 1984, 1987 by Stephen L. Moshier
      *************************************************************************/
      public static double beta(double a, double b)
      {
      
          double result = betaf.beta(a, b);
          return result;
      }
  
  }
  public partial class alglib
  {
  
      
      /*************************************************************************
      Calculation of the value of the Chebyshev polynomials of the
      first and second kinds.
  
      Parameters:
          r   -   polynomial kind, either 1 or 2.
          n   -   degree, n>=0
          x   -   argument, -1 <= x <= 1
  
      Result:
          the value of the Chebyshev polynomial at x
      *************************************************************************/
      public static double chebyshevcalculate(int r, int n, double x)
      {
      
          double result = chebyshev.chebyshevcalculate(r, n, x);
          return result;
      }
      
      /*************************************************************************
      Summation of Chebyshev polynomials using Clenshaw's recurrence formula.
  
      This routine calculates
          c[0]*T0(x) + c[1]*T1(x) + ... + c[N]*TN(x)
      or
          c[0]*U0(x) + c[1]*U1(x) + ... + c[N]*UN(x)
      depending on the R.
  
      Parameters:
          r   -   polynomial kind, either 1 or 2.
          n   -   degree, n>=0
          x   -   argument
  
      Result:
          the value of the Chebyshev polynomial at x
      *************************************************************************/
      public static double chebyshevsum(double[] c, int r, int n, double x)
      {
      
          double result = chebyshev.chebyshevsum(c, r, n, x);
          return result;
      }
      
      /*************************************************************************
      Representation of Tn as C[0] + C[1]*X + ... + C[N]*X^N
  
      Input parameters:
          N   -   polynomial degree, n>=0
  
      Output parameters:
          C   -   coefficients
      *************************************************************************/
      public static void chebyshevcoefficients(int n, out double[] c)
      {
          c = new double[0];
          chebyshev.chebyshevcoefficients(n, ref c);
          return;
      }
      
      /*************************************************************************
      Conversion of a series of Chebyshev polynomials to a power series.
  
      Represents A[0]*T0(x) + A[1]*T1(x) + ... + A[N]*Tn(x) as
      B[0] + B[1]*X + ... + B[N]*X^N.
  
      Input parameters:
          A   -   Chebyshev series coefficients
          N   -   degree, N>=0
  
      Output parameters
          B   -   power series coefficients
      *************************************************************************/
      public static void fromchebyshev(double[] a, int n, out double[] b)
      {
          b = new double[0];
          chebyshev.fromchebyshev(a, n, ref b);
          return;
      }
  
  }
  public partial class alglib
  {
  
      
      /*************************************************************************
      Student's t distribution
  
      Computes the integral from minus infinity to t of the Student
      t distribution with integer k > 0 degrees of freedom:
  
                                           t
                                           -
                                          | |
                   -                      |         2   -(k+1)/2
                  | ( (k+1)/2 )           |  (     x   )
            ----------------------        |  ( 1 + --- )        dx
                          -               |  (      k  )
            sqrt( k pi ) | ( k/2 )        |
                                        | |
                                         -
                                        -inf.
  
      Relation to incomplete beta integral:
  
             1 - stdtr(k,t) = 0.5 * incbet( k/2, 1/2, z )
      where
             z = k/(k + t**2).
  
      For t < -2, this is the method of computation.  For higher t,
      a direct method is derived from integration by parts.
      Since the function is symmetric about t=0, the area under the
      right tail of the density is found by calling the function
      with -t instead of t.
  
      ACCURACY:
  
      Tested at random 1 <= k <= 25.  The "domain" refers to t.
                           Relative error:
      arithmetic   domain     # trials      peak         rms
         IEEE     -100,-2      50000       5.9e-15     1.4e-15
         IEEE     -2,100      500000       2.7e-15     4.9e-17
  
      Cephes Math Library Release 2.8:  June, 2000
      Copyright 1984, 1987, 1995, 2000 by Stephen L. Moshier
      *************************************************************************/
      public static double studenttdistribution(int k, double t)
      {
      
          double result = studenttdistr.studenttdistribution(k, t);
          return result;
      }
      
      /*************************************************************************
      Functional inverse of Student's t distribution
  
      Given probability p, finds the argument t such that stdtr(k,t)
      is equal to p.
  
      ACCURACY:
  
      Tested at random 1 <= k <= 100.  The "domain" refers to p:
                           Relative error:
      arithmetic   domain     # trials      peak         rms
         IEEE    .001,.999     25000       5.7e-15     8.0e-16
         IEEE    10^-6,.001    25000       2.0e-12     2.9e-14
  
      Cephes Math Library Release 2.8:  June, 2000
      Copyright 1984, 1987, 1995, 2000 by Stephen L. Moshier
      *************************************************************************/
      public static double invstudenttdistribution(int k, double p)
      {
      
          double result = studenttdistr.invstudenttdistribution(k, p);
          return result;
      }
  
  }
  public partial class alglib
  {
  
      
      /*************************************************************************
      Binomial distribution
  
      Returns the sum of the terms 0 through k of the Binomial
      probability density:
  
        k
        --  ( n )   j      n-j
        >   (   )  p  (1-p)
        --  ( j )
       j=0
  
      The terms are not summed directly; instead the incomplete
      beta integral is employed, according to the formula
  
      y = bdtr( k, n, p ) = incbet( n-k, k+1, 1-p ).
  
      The arguments must be positive, with p ranging from 0 to 1.
  
      ACCURACY:
  
      Tested at random points (a,b,p), with p between 0 and 1.
  
                    a,b                     Relative error:
      arithmetic  domain     # trials      peak         rms
       For p between 0.001 and 1:
         IEEE     0,100       100000      4.3e-15     2.6e-16
  
      Cephes Math Library Release 2.8:  June, 2000
      Copyright 1984, 1987, 1995, 2000 by Stephen L. Moshier
      *************************************************************************/
      public static double binomialdistribution(int k, int n, double p)
      {
      
          double result = binomialdistr.binomialdistribution(k, n, p);
          return result;
      }
      
      /*************************************************************************
      Complemented binomial distribution
  
      Returns the sum of the terms k+1 through n of the Binomial
      probability density:
  
        n
        --  ( n )   j      n-j
        >   (   )  p  (1-p)
        --  ( j )
       j=k+1
  
      The terms are not summed directly; instead the incomplete
      beta integral is employed, according to the formula
  
      y = bdtrc( k, n, p ) = incbet( k+1, n-k, p ).
  
      The arguments must be positive, with p ranging from 0 to 1.
  
      ACCURACY:
  
      Tested at random points (a,b,p).
  
                    a,b                     Relative error:
      arithmetic  domain     # trials      peak         rms
       For p between 0.001 and 1:
         IEEE     0,100       100000      6.7e-15     8.2e-16
       For p between 0 and .001:
         IEEE     0,100       100000      1.5e-13     2.7e-15
  
      Cephes Math Library Release 2.8:  June, 2000
      Copyright 1984, 1987, 1995, 2000 by Stephen L. Moshier
      *************************************************************************/
      public static double binomialcdistribution(int k, int n, double p)
      {
      
          double result = binomialdistr.binomialcdistribution(k, n, p);
          return result;
      }
      
      /*************************************************************************
      Inverse binomial distribution
  
      Finds the event probability p such that the sum of the
      terms 0 through k of the Binomial probability density
      is equal to the given cumulative probability y.
  
      This is accomplished using the inverse beta integral
      function and the relation
  
      1 - p = incbi( n-k, k+1, y ).
  
      ACCURACY:
  
      Tested at random points (a,b,p).
  
                    a,b                     Relative error:
      arithmetic  domain     # trials      peak         rms
       For p between 0.001 and 1:
         IEEE     0,100       100000      2.3e-14     6.4e-16
         IEEE     0,10000     100000      6.6e-12     1.2e-13
       For p between 10^-6 and 0.001:
         IEEE     0,100       100000      2.0e-12     1.3e-14
         IEEE     0,10000     100000      1.5e-12     3.2e-14
  
      Cephes Math Library Release 2.8:  June, 2000
      Copyright 1984, 1987, 1995, 2000 by Stephen L. Moshier
      *************************************************************************/
      public static double invbinomialdistribution(int k, int n, double y)
      {
      
          double result = binomialdistr.invbinomialdistribution(k, n, y);
          return result;
      }
  
  }
  public partial class alglib
  {
  
      
      /*************************************************************************
      Airy function
  
      Solution of the differential equation
  
      y"(x) = xy.
  
      The function returns the two independent solutions Ai, Bi
      and their first derivatives Ai'(x), Bi'(x).
  
      Evaluation is by power series summation for small x,
      by rational minimax approximations for large x.
  
  
  
      ACCURACY:
      Error criterion is absolute when function <= 1, relative
      when function > 1, except * denotes relative error criterion.
      For large negative x, the absolute error increases as x^1.5.
      For large positive x, the relative error increases as x^1.5.
  
      Arithmetic  domain   function  # trials      peak         rms
      IEEE        -10, 0     Ai        10000       1.6e-15     2.7e-16
      IEEE          0, 10    Ai        10000       2.3e-14*    1.8e-15*
      IEEE        -10, 0     Ai'       10000       4.6e-15     7.6e-16
      IEEE          0, 10    Ai'       10000       1.8e-14*    1.5e-15*
      IEEE        -10, 10    Bi        30000       4.2e-15     5.3e-16
      IEEE        -10, 10    Bi'       30000       4.9e-15     7.3e-16
  
      Cephes Math Library Release 2.8:  June, 2000
      Copyright 1984, 1987, 1989, 2000 by Stephen L. Moshier
      *************************************************************************/
      public static void airy(double x, out double ai, out double aip, out double bi, out double bip)
      {
          ai = 0;
          aip = 0;
          bi = 0;
          bip = 0;
          airyf.airy(x, ref ai, ref aip, ref bi, ref bip);
          return;
      }
  
  }
  public partial class alglib
  {
      public class gammafunc
      {
          /*************************************************************************
          Gamma function
  
          Input parameters:
              X   -   argument
  
          Domain:
              0 < X < 171.6
              -170 < X < 0, X is not an integer.
  
          Relative error:
           arithmetic   domain     # trials      peak         rms
              IEEE    -170,-33      20000       2.3e-15     3.3e-16
              IEEE     -33,  33     20000       9.4e-16     2.2e-16
              IEEE      33, 171.6   20000       2.3e-15     3.2e-16
  
          Cephes Math Library Release 2.8:  June, 2000
          Original copyright 1984, 1987, 1989, 1992, 2000 by Stephen L. Moshier
          Translated to AlgoPascal by Bochkanov Sergey (2005, 2006, 2007).
          *************************************************************************/
          public static double gammafunction(double x)
          {
              double result = 0;
              double p = 0;
              double pp = 0;
              double q = 0;
              double qq = 0;
              double z = 0;
              int i = 0;
              double sgngam = 0;
  
              sgngam = 1;
              q = Math.Abs(x);
              if( (double)(q)>(double)(33.0) )
              {
                  if( (double)(x)<(double)(0.0) )
                  {
                      p = (int)Math.Floor(q);
                      i = (int)Math.Round(p);
                      if( i%2==0 )
                      {
                          sgngam = -1;
                      }
                      z = q-p;
                      if( (double)(z)>(double)(0.5) )
                      {
                          p = p+1;
                          z = q-p;
                      }
                      z = q*Math.Sin(Math.PI*z);
                      z = Math.Abs(z);
                      z = Math.PI/(z*gammastirf(q));
                  }
                  else
                  {
                      z = gammastirf(x);
                  }
                  result = sgngam*z;
                  return result;
              }
              z = 1;
              while( (double)(x)>=(double)(3) )
              {
                  x = x-1;
                  z = z*x;
              }
              while( (double)(x)<(double)(0) )
              {
                  if( (double)(x)>(double)(-0.000000001) )
                  {
                      result = z/((1+0.5772156649015329*x)*x);
                      return result;
                  }
                  z = z/x;
                  x = x+1;
              }
              while( (double)(x)<(double)(2) )
              {
                  if( (double)(x)<(double)(0.000000001) )
                  {
                      result = z/((1+0.5772156649015329*x)*x);
                      return result;
                  }
                  z = z/x;
                  x = x+1.0;
              }
              if( (double)(x)==(double)(2) )
              {
                  result = z;
                  return result;
              }
              x = x-2.0;
              pp = 1.60119522476751861407E-4;
              pp = 1.19135147006586384913E-3+x*pp;
              pp = 1.04213797561761569935E-2+x*pp;
              pp = 4.76367800457137231464E-2+x*pp;
              pp = 2.07448227648435975150E-1+x*pp;
              pp = 4.94214826801497100753E-1+x*pp;
              pp = 9.99999999999999996796E-1+x*pp;
              qq = -2.31581873324120129819E-5;
              qq = 5.39605580493303397842E-4+x*qq;
              qq = -4.45641913851797240494E-3+x*qq;
              qq = 1.18139785222060435552E-2+x*qq;
              qq = 3.58236398605498653373E-2+x*qq;
              qq = -2.34591795718243348568E-1+x*qq;
              qq = 7.14304917030273074085E-2+x*qq;
              qq = 1.00000000000000000320+x*qq;
              result = z*pp/qq;
              return result;
          }
  
  
          /*************************************************************************
          Natural logarithm of gamma function
  
          Input parameters:
              X       -   argument
  
          Result:
              logarithm of the absolute value of the Gamma(X).
  
          Output parameters:
              SgnGam  -   sign(Gamma(X))
  
          Domain:
              0 < X < 2.55e305
              -2.55e305 < X < 0, X is not an integer.
  
          ACCURACY:
          arithmetic      domain        # trials     peak         rms
             IEEE    0, 3                 28000     5.4e-16     1.1e-16
             IEEE    2.718, 2.556e305     40000     3.5e-16     8.3e-17
          The error criterion was relative when the function magnitude
          was greater than one but absolute when it was less than one.
  
          The following test used the relative error criterion, though
          at certain points the relative error could be much higher than
          indicated.
             IEEE    -200, -4             10000     4.8e-16     1.3e-16
  
          Cephes Math Library Release 2.8:  June, 2000
          Copyright 1984, 1987, 1989, 1992, 2000 by Stephen L. Moshier
          Translated to AlgoPascal by Bochkanov Sergey (2005, 2006, 2007).
          *************************************************************************/
          public static double lngamma(double x,
              ref double sgngam)
          {
              double result = 0;
              double a = 0;
              double b = 0;
              double c = 0;
              double p = 0;
              double q = 0;
              double u = 0;
              double w = 0;
              double z = 0;
              int i = 0;
              double logpi = 0;
              double ls2pi = 0;
              double tmp = 0;
  
              sgngam = 0;
  
              sgngam = 1;
              logpi = 1.14472988584940017414;
              ls2pi = 0.91893853320467274178;
              if( (double)(x)<(double)(-34.0) )
              {
                  q = -x;
                  w = lngamma(q, ref tmp);
                  p = (int)Math.Floor(q);
                  i = (int)Math.Round(p);
                  if( i%2==0 )
                  {
                      sgngam = -1;
                  }
                  else
                  {
                      sgngam = 1;
                  }
                  z = q-p;
                  if( (double)(z)>(double)(0.5) )
                  {
                      p = p+1;
                      z = p-q;
                  }
                  z = q*Math.Sin(Math.PI*z);
                  result = logpi-Math.Log(z)-w;
                  return result;
              }
              if( (double)(x)<(double)(13) )
              {
                  z = 1;
                  p = 0;
                  u = x;
                  while( (double)(u)>=(double)(3) )
                  {
                      p = p-1;
                      u = x+p;
                      z = z*u;
                  }
                  while( (double)(u)<(double)(2) )
                  {
                      z = z/u;
                      p = p+1;
                      u = x+p;
                  }
                  if( (double)(z)<(double)(0) )
                  {
                      sgngam = -1;
                      z = -z;
                  }
                  else
                  {
                      sgngam = 1;
                  }
                  if( (double)(u)==(double)(2) )
                  {
                      result = Math.Log(z);
                      return result;
                  }
                  p = p-2;
                  x = x+p;
                  b = -1378.25152569120859100;
                  b = -38801.6315134637840924+x*b;
                  b = -331612.992738871184744+x*b;
                  b = -1162370.97492762307383+x*b;
                  b = -1721737.00820839662146+x*b;
                  b = -853555.664245765465627+x*b;
                  c = 1;
                  c = -351.815701436523470549+x*c;
                  c = -17064.2106651881159223+x*c;
                  c = -220528.590553854454839+x*c;
                  c = -1139334.44367982507207+x*c;
                  c = -2532523.07177582951285+x*c;
                  c = -2018891.41433532773231+x*c;
                  p = x*b/c;
                  result = Math.Log(z)+p;
                  return result;
              }
              q = (x-0.5)*Math.Log(x)-x+ls2pi;
              if( (double)(x)>(double)(100000000) )
              {
                  result = q;
                  return result;
              }
              p = 1/(x*x);
              if( (double)(x)>=(double)(1000.0) )
              {
                  q = q+((7.9365079365079365079365*0.0001*p-2.7777777777777777777778*0.001)*p+0.0833333333333333333333)/x;
              }
              else
              {
                  a = 8.11614167470508450300*0.0001;
                  a = -(5.95061904284301438324*0.0001)+p*a;
                  a = 7.93650340457716943945*0.0001+p*a;
                  a = -(2.77777777730099687205*0.001)+p*a;
                  a = 8.33333333333331927722*0.01+p*a;
                  q = q+a/x;
              }
              result = q;
              return result;
          }
  
  
          private static double gammastirf(double x)
          {
              double result = 0;
              double y = 0;
              double w = 0;
              double v = 0;
              double stir = 0;
  
              w = 1/x;
              stir = 7.87311395793093628397E-4;
              stir = -2.29549961613378126380E-4+w*stir;
              stir = -2.68132617805781232825E-3+w*stir;
              stir = 3.47222221605458667310E-3+w*stir;
              stir = 8.33333333333482257126E-2+w*stir;
              w = 1+w*stir;
              y = Math.Exp(x);
              if( (double)(x)>(double)(143.01608) )
              {
                  v = Math.Pow(x, 0.5*x-0.25);
                  y = v*(v/y);
              }
              else
              {
                  y = Math.Pow(x, x-0.5)/y;
              }
              result = 2.50662827463100050242*y*w;
              return result;
          }
  
  
      }
      public class normaldistr
      {
          /*************************************************************************
          Error function
  
          The integral is
  
                                    x
                                     -
                          2         | |          2
            erf(x)  =  --------     |    exp( - t  ) dt.
                       sqrt(pi)   | |
                                   -
                                    0
  
          For 0 <= |x| < 1, erf(x) = x * P4(x**2)/Q5(x**2); otherwise
          erf(x) = 1 - erfc(x).
  
  
          ACCURACY:
  
                               Relative error:
          arithmetic   domain     # trials      peak         rms
             IEEE      0,1         30000       3.7e-16     1.0e-16
  
          Cephes Math Library Release 2.8:  June, 2000
          Copyright 1984, 1987, 1988, 1992, 2000 by Stephen L. Moshier
          *************************************************************************/
          public static double errorfunction(double x)
          {
              double result = 0;
              double xsq = 0;
              double s = 0;
              double p = 0;
              double q = 0;
  
              s = Math.Sign(x);
              x = Math.Abs(x);
              if( (double)(x)<(double)(0.5) )
              {
                  xsq = x*x;
                  p = 0.007547728033418631287834;
                  p = -0.288805137207594084924010+xsq*p;
                  p = 14.3383842191748205576712+xsq*p;
                  p = 38.0140318123903008244444+xsq*p;
                  p = 3017.82788536507577809226+xsq*p;
                  p = 7404.07142710151470082064+xsq*p;
                  p = 80437.3630960840172832162+xsq*p;
                  q = 0.0;
                  q = 1.00000000000000000000000+xsq*q;
                  q = 38.0190713951939403753468+xsq*q;
                  q = 658.070155459240506326937+xsq*q;
                  q = 6379.60017324428279487120+xsq*q;
                  q = 34216.5257924628539769006+xsq*q;
                  q = 80437.3630960840172826266+xsq*q;
                  result = s*1.1283791670955125738961589031*x*p/q;
                  return result;
              }
              if( (double)(x)>=(double)(10) )
              {
                  result = s;
                  return result;
              }
              result = s*(1-errorfunctionc(x));
              return result;
          }
  
  
          /*************************************************************************
          Complementary error function
  
           1 - erf(x) =
  
                                    inf.
                                      -
                           2         | |          2
            erfc(x)  =  --------     |    exp( - t  ) dt
                        sqrt(pi)   | |
                                    -
                                     x
  
  
          For small x, erfc(x) = 1 - erf(x); otherwise rational
          approximations are computed.
  
  
          ACCURACY:
  
                               Relative error:
          arithmetic   domain     # trials      peak         rms
             IEEE      0,26.6417   30000       5.7e-14     1.5e-14
  
          Cephes Math Library Release 2.8:  June, 2000
          Copyright 1984, 1987, 1988, 1992, 2000 by Stephen L. Moshier
          *************************************************************************/
          public static double errorfunctionc(double x)
          {
              double result = 0;
              double p = 0;
              double q = 0;
  
              if( (double)(x)<(double)(0) )
              {
                  result = 2-errorfunctionc(-x);
                  return result;
              }
              if( (double)(x)<(double)(0.5) )
              {
                  result = 1.0-errorfunction(x);
                  return result;
              }
              if( (double)(x)>=(double)(10) )
              {
                  result = 0;
                  return result;
              }
              p = 0.0;
              p = 0.5641877825507397413087057563+x*p;
              p = 9.675807882987265400604202961+x*p;
              p = 77.08161730368428609781633646+x*p;
              p = 368.5196154710010637133875746+x*p;
              p = 1143.262070703886173606073338+x*p;
              p = 2320.439590251635247384768711+x*p;
              p = 2898.0293292167655611275846+x*p;
              p = 1826.3348842295112592168999+x*p;
              q = 1.0;
              q = 17.14980943627607849376131193+x*q;
              q = 137.1255960500622202878443578+x*q;
              q = 661.7361207107653469211984771+x*q;
              q = 2094.384367789539593790281779+x*q;
              q = 4429.612803883682726711528526+x*q;
              q = 6089.5424232724435504633068+x*q;
              q = 4958.82756472114071495438422+x*q;
              q = 1826.3348842295112595576438+x*q;
              result = Math.Exp(-math.sqr(x))*p/q;
              return result;
          }
  
  
          /*************************************************************************
          Normal distribution function
  
          Returns the area under the Gaussian probability density
          function, integrated from minus infinity to x:
  
                                     x
                                      -
                            1        | |          2
             ndtr(x)  = ---------    |    exp( - t /2 ) dt
                        sqrt(2pi)  | |
                                    -
                                   -inf.
  
                      =  ( 1 + erf(z) ) / 2
                      =  erfc(z) / 2
  
          where z = x/sqrt(2). Computation is via the functions
          erf and erfc.
  
  
          ACCURACY:
  
                               Relative error:
          arithmetic   domain     # trials      peak         rms
             IEEE     -13,0        30000       3.4e-14     6.7e-15
  
          Cephes Math Library Release 2.8:  June, 2000
          Copyright 1984, 1987, 1988, 1992, 2000 by Stephen L. Moshier
          *************************************************************************/
          public static double normaldistribution(double x)
          {
              double result = 0;
  
              result = 0.5*(errorfunction(x/1.41421356237309504880)+1);
              return result;
          }
  
  
          /*************************************************************************
          Inverse of the error function
  
          Cephes Math Library Release 2.8:  June, 2000
          Copyright 1984, 1987, 1988, 1992, 2000 by Stephen L. Moshier
          *************************************************************************/
          public static double inverf(double e)
          {
              double result = 0;
  
              result = invnormaldistribution(0.5*(e+1))/Math.Sqrt(2);
              return result;
          }
  
  
          /*************************************************************************
          Inverse of Normal distribution function
  
          Returns the argument, x, for which the area under the
          Gaussian probability density function (integrated from
          minus infinity to x) is equal to y.
  
  
          For small arguments 0 < y < exp(-2), the program computes
          z = sqrt( -2.0 * log(y) );  then the approximation is
          x = z - log(z)/z  - (1/z) P(1/z) / Q(1/z).
          There are two rational functions P/Q, one for 0 < y < exp(-32)
          and the other for y up to exp(-2).  For larger arguments,
          w = y - 0.5, and  x/sqrt(2pi) = w + w**3 R(w**2)/S(w**2)).
  
          ACCURACY:
  
                               Relative error:
          arithmetic   domain        # trials      peak         rms
             IEEE     0.125, 1        20000       7.2e-16     1.3e-16
             IEEE     3e-308, 0.135   50000       4.6e-16     9.8e-17
  
          Cephes Math Library Release 2.8:  June, 2000
          Copyright 1984, 1987, 1988, 1992, 2000 by Stephen L. Moshier
          *************************************************************************/
          public static double invnormaldistribution(double y0)
          {
              double result = 0;
              double expm2 = 0;
              double s2pi = 0;
              double x = 0;
              double y = 0;
              double z = 0;
              double y2 = 0;
              double x0 = 0;
              double x1 = 0;
              int code = 0;
              double p0 = 0;
              double q0 = 0;
              double p1 = 0;
              double q1 = 0;
              double p2 = 0;
              double q2 = 0;
  
              expm2 = 0.13533528323661269189;
              s2pi = 2.50662827463100050242;
              if( (double)(y0)<=(double)(0) )
              {
                  result = -math.maxrealnumber;
                  return result;
              }
              if( (double)(y0)>=(double)(1) )
              {
                  result = math.maxrealnumber;
                  return result;
              }
              code = 1;
              y = y0;
              if( (double)(y)>(double)(1.0-expm2) )
              {
                  y = 1.0-y;
                  code = 0;
              }
              if( (double)(y)>(double)(expm2) )
              {
                  y = y-0.5;
                  y2 = y*y;
                  p0 = -59.9633501014107895267;
                  p0 = 98.0010754185999661536+y2*p0;
                  p0 = -56.6762857469070293439+y2*p0;
                  p0 = 13.9312609387279679503+y2*p0;
                  p0 = -1.23916583867381258016+y2*p0;
                  q0 = 1;
                  q0 = 1.95448858338141759834+y2*q0;
                  q0 = 4.67627912898881538453+y2*q0;
                  q0 = 86.3602421390890590575+y2*q0;
                  q0 = -225.462687854119370527+y2*q0;
                  q0 = 200.260212380060660359+y2*q0;
                  q0 = -82.0372256168333339912+y2*q0;
                  q0 = 15.9056225126211695515+y2*q0;
                  q0 = -1.18331621121330003142+y2*q0;
                  x = y+y*y2*p0/q0;
                  x = x*s2pi;
                  result = x;
                  return result;
              }
              x = Math.Sqrt(-(2.0*Math.Log(y)));
              x0 = x-Math.Log(x)/x;
              z = 1.0/x;
              if( (double)(x)<(double)(8.0) )
              {
                  p1 = 4.05544892305962419923;
                  p1 = 31.5251094599893866154+z*p1;
                  p1 = 57.1628192246421288162+z*p1;
                  p1 = 44.0805073893200834700+z*p1;
                  p1 = 14.6849561928858024014+z*p1;
                  p1 = 2.18663306850790267539+z*p1;
                  p1 = -(1.40256079171354495875*0.1)+z*p1;
                  p1 = -(3.50424626827848203418*0.01)+z*p1;
                  p1 = -(8.57456785154685413611*0.0001)+z*p1;
                  q1 = 1;
                  q1 = 15.7799883256466749731+z*q1;
                  q1 = 45.3907635128879210584+z*q1;
                  q1 = 41.3172038254672030440+z*q1;
                  q1 = 15.0425385692907503408+z*q1;
                  q1 = 2.50464946208309415979+z*q1;
                  q1 = -(1.42182922854787788574*0.1)+z*q1;
                  q1 = -(3.80806407691578277194*0.01)+z*q1;
                  q1 = -(9.33259480895457427372*0.0001)+z*q1;
                  x1 = z*p1/q1;
              }
              else
              {
                  p2 = 3.23774891776946035970;
                  p2 = 6.91522889068984211695+z*p2;
                  p2 = 3.93881025292474443415+z*p2;
                  p2 = 1.33303460815807542389+z*p2;
                  p2 = 2.01485389549179081538*0.1+z*p2;
                  p2 = 1.23716634817820021358*0.01+z*p2;
                  p2 = 3.01581553508235416007*0.0001+z*p2;
                  p2 = 2.65806974686737550832*0.000001+z*p2;
                  p2 = 6.23974539184983293730*0.000000001+z*p2;
                  q2 = 1;
                  q2 = 6.02427039364742014255+z*q2;
                  q2 = 3.67983563856160859403+z*q2;
                  q2 = 1.37702099489081330271+z*q2;
                  q2 = 2.16236993594496635890*0.1+z*q2;
                  q2 = 1.34204006088543189037*0.01+z*q2;
                  q2 = 3.28014464682127739104*0.0001+z*q2;
                  q2 = 2.89247864745380683936*0.000001+z*q2;
                  q2 = 6.79019408009981274425*0.000000001+z*q2;
                  x1 = z*p2/q2;
              }
              x = x0-x1;
              if( code!=0 )
              {
                  x = -x;
              }
              result = x;
              return result;
          }
  
  
      }
      public class igammaf
      {
          /*************************************************************************
          Incomplete gamma integral
  
          The function is defined by
  
                                    x
                                     -
                            1       | |  -t  a-1
           igam(a,x)  =   -----     |   e   t   dt.
                           -      | |
                          | (a)    -
                                    0
  
  
          In this implementation both arguments must be positive.
          The integral is evaluated by either a power series or
          continued fraction expansion, depending on the relative
          values of a and x.
  
          ACCURACY:
  
                               Relative error:
          arithmetic   domain     # trials      peak         rms
             IEEE      0,30       200000       3.6e-14     2.9e-15
             IEEE      0,100      300000       9.9e-14     1.5e-14
  
          Cephes Math Library Release 2.8:  June, 2000
          Copyright 1985, 1987, 2000 by Stephen L. Moshier
          *************************************************************************/
          public static double incompletegamma(double a,
              double x)
          {
              double result = 0;
              double igammaepsilon = 0;
              double ans = 0;
              double ax = 0;
              double c = 0;
              double r = 0;
              double tmp = 0;
  
              igammaepsilon = 0.000000000000001;
              if( (double)(x)<=(double)(0) || (double)(a)<=(double)(0) )
              {
                  result = 0;
                  return result;
              }
              if( (double)(x)>(double)(1) && (double)(x)>(double)(a) )
              {
                  result = 1-incompletegammac(a, x);
                  return result;
              }
              ax = a*Math.Log(x)-x-gammafunc.lngamma(a, ref tmp);
              if( (double)(ax)<(double)(-709.78271289338399) )
              {
                  result = 0;
                  return result;
              }
              ax = Math.Exp(ax);
              r = a;
              c = 1;
              ans = 1;
              do
              {
                  r = r+1;
                  c = c*x/r;
                  ans = ans+c;
              }
              while( (double)(c/ans)>(double)(igammaepsilon) );
              result = ans*ax/a;
              return result;
          }
  
  
          /*************************************************************************
          Complemented incomplete gamma integral
  
          The function is defined by
  
  
           igamc(a,x)   =   1 - igam(a,x)
  
                                     inf.
                                       -
                              1       | |  -t  a-1
                        =   -----     |   e   t   dt.
                             -      | |
                            | (a)    -
                                      x
  
  
          In this implementation both arguments must be positive.
          The integral is evaluated by either a power series or
          continued fraction expansion, depending on the relative
          values of a and x.
  
          ACCURACY:
  
          Tested at random a, x.
                         a         x                      Relative error:
          arithmetic   domain   domain     # trials      peak         rms
             IEEE     0.5,100   0,100      200000       1.9e-14     1.7e-15
             IEEE     0.01,0.5  0,100      200000       1.4e-13     1.6e-15
  
          Cephes Math Library Release 2.8:  June, 2000
          Copyright 1985, 1987, 2000 by Stephen L. Moshier
          *************************************************************************/
          public static double incompletegammac(double a,
              double x)
          {
              double result = 0;
              double igammaepsilon = 0;
              double igammabignumber = 0;
              double igammabignumberinv = 0;
              double ans = 0;
              double ax = 0;
              double c = 0;
              double yc = 0;
              double r = 0;
              double t = 0;
              double y = 0;
              double z = 0;
              double pk = 0;
              double pkm1 = 0;
              double pkm2 = 0;
              double qk = 0;
              double qkm1 = 0;
              double qkm2 = 0;
              double tmp = 0;
  
              igammaepsilon = 0.000000000000001;
              igammabignumber = 4503599627370496.0;
              igammabignumberinv = 2.22044604925031308085*0.0000000000000001;
              if( (double)(x)<=(double)(0) || (double)(a)<=(double)(0) )
              {
                  result = 1;
                  return result;
              }
              if( (double)(x)<(double)(1) || (double)(x)<(double)(a) )
              {
                  result = 1-incompletegamma(a, x);
                  return result;
              }
              ax = a*Math.Log(x)-x-gammafunc.lngamma(a, ref tmp);
              if( (double)(ax)<(double)(-709.78271289338399) )
              {
                  result = 0;
                  return result;
              }
              ax = Math.Exp(ax);
              y = 1-a;
              z = x+y+1;
              c = 0;
              pkm2 = 1;
              qkm2 = x;
              pkm1 = x+1;
              qkm1 = z*x;
              ans = pkm1/qkm1;
              do
              {
                  c = c+1;
                  y = y+1;
                  z = z+2;
                  yc = y*c;
                  pk = pkm1*z-pkm2*yc;
                  qk = qkm1*z-qkm2*yc;
                  if( (double)(qk)!=(double)(0) )
                  {
                      r = pk/qk;
                      t = Math.Abs((ans-r)/r);
                      ans = r;
                  }
                  else
                  {
                      t = 1;
                  }
                  pkm2 = pkm1;
                  pkm1 = pk;
                  qkm2 = qkm1;
                  qkm1 = qk;
                  if( (double)(Math.Abs(pk))>(double)(igammabignumber) )
                  {
                      pkm2 = pkm2*igammabignumberinv;
                      pkm1 = pkm1*igammabignumberinv;
                      qkm2 = qkm2*igammabignumberinv;
                      qkm1 = qkm1*igammabignumberinv;
                  }
              }
              while( (double)(t)>(double)(igammaepsilon) );
              result = ans*ax;
              return result;
          }
  
  
          /*************************************************************************
          Inverse of complemented imcomplete gamma integral
  
          Given p, the function finds x such that
  
           igamc( a, x ) = p.
  
          Starting with the approximate value
  
                  3
           x = a t
  
           where
  
           t = 1 - d - ndtri(p) sqrt(d)
  
          and
  
           d = 1/9a,
  
          the routine performs up to 10 Newton iterations to find the
          root of igamc(a,x) - p = 0.
  
          ACCURACY:
  
          Tested at random a, p in the intervals indicated.
  
                         a        p                      Relative error:
          arithmetic   domain   domain     # trials      peak         rms
             IEEE     0.5,100   0,0.5       100000       1.0e-14     1.7e-15
             IEEE     0.01,0.5  0,0.5       100000       9.0e-14     3.4e-15
             IEEE    0.5,10000  0,0.5        20000       2.3e-13     3.8e-14
  
          Cephes Math Library Release 2.8:  June, 2000
          Copyright 1984, 1987, 1995, 2000 by Stephen L. Moshier
          *************************************************************************/
          public static double invincompletegammac(double a,
              double y0)
          {
              double result = 0;
              double igammaepsilon = 0;
              double iinvgammabignumber = 0;
              double x0 = 0;
              double x1 = 0;
              double x = 0;
              double yl = 0;
              double yh = 0;
              double y = 0;
              double d = 0;
              double lgm = 0;
              double dithresh = 0;
              int i = 0;
              int dir = 0;
              double tmp = 0;
  
              igammaepsilon = 0.000000000000001;
              iinvgammabignumber = 4503599627370496.0;
              x0 = iinvgammabignumber;
              yl = 0;
              x1 = 0;
              yh = 1;
              dithresh = 5*igammaepsilon;
              d = 1/(9*a);
              y = 1-d-normaldistr.invnormaldistribution(y0)*Math.Sqrt(d);
              x = a*y*y*y;
              lgm = gammafunc.lngamma(a, ref tmp);
              i = 0;
              while( i<10 )
              {
                  if( (double)(x)>(double)(x0) || (double)(x)<(double)(x1) )
                  {
                      d = 0.0625;
                      break;
                  }
                  y = incompletegammac(a, x);
                  if( (double)(y)<(double)(yl) || (double)(y)>(double)(yh) )
                  {
                      d = 0.0625;
                      break;
                  }
                  if( (double)(y)<(double)(y0) )
                  {
                      x0 = x;
                      yl = y;
                  }
                  else
                  {
                      x1 = x;
                      yh = y;
                  }
                  d = (a-1)*Math.Log(x)-x-lgm;
                  if( (double)(d)<(double)(-709.78271289338399) )
                  {
                      d = 0.0625;
                      break;
                  }
                  d = -Math.Exp(d);
                  d = (y-y0)/d;
                  if( (double)(Math.Abs(d/x))<(double)(igammaepsilon) )
                  {
                      result = x;
                      return result;
                  }
                  x = x-d;
                  i = i+1;
              }
              if( (double)(x0)==(double)(iinvgammabignumber) )
              {
                  if( (double)(x)<=(double)(0) )
                  {
                      x = 1;
                  }
                  while( (double)(x0)==(double)(iinvgammabignumber) )
                  {
                      x = (1+d)*x;
                      y = incompletegammac(a, x);
                      if( (double)(y)<(double)(y0) )
                      {
                          x0 = x;
                          yl = y;
                          break;
                      }
                      d = d+d;
                  }
              }
              d = 0.5;
              dir = 0;
              i = 0;
              while( i<400 )
              {
                  x = x1+d*(x0-x1);
                  y = incompletegammac(a, x);
                  lgm = (x0-x1)/(x1+x0);
                  if( (double)(Math.Abs(lgm))<(double)(dithresh) )
                  {
                      break;
                  }
                  lgm = (y-y0)/y0;
                  if( (double)(Math.Abs(lgm))<(double)(dithresh) )
                  {
                      break;
                  }
                  if( (double)(x)<=(double)(0.0) )
                  {
                      break;
                  }
                  if( (double)(y)>=(double)(y0) )
                  {
                      x1 = x;
                      yh = y;
                      if( dir<0 )
                      {
                          dir = 0;
                          d = 0.5;
                      }
                      else
                      {
                          if( dir>1 )
                          {
                              d = 0.5*d+0.5;
                          }
                          else
                          {
                              d = (y0-yl)/(yh-yl);
                          }
                      }
                      dir = dir+1;
                  }
                  else
                  {
                      x0 = x;
                      yl = y;
                      if( dir>0 )
                      {
                          dir = 0;
                          d = 0.5;
                      }
                      else
                      {
                          if( dir<-1 )
                          {
                              d = 0.5*d;
                          }
                          else
                          {
                              d = (y0-yl)/(yh-yl);
                          }
                      }
                      dir = dir-1;
                  }
                  i = i+1;
              }
              result = x;
              return result;
          }
  
  
      }
      public class elliptic
      {
          /*************************************************************************
          Complete elliptic integral of the first kind
  
          Approximates the integral
  
  
  
                     pi/2
                      -
                     | |
                     |           dt
          K(m)  =    |    ------------------
                     |                   2
                   | |    sqrt( 1 - m sin t )
                    -
                     0
  
          using the approximation
  
              P(x)  -  log x Q(x).
  
          ACCURACY:
  
                               Relative error:
          arithmetic   domain     # trials      peak         rms
             IEEE       0,1        30000       2.5e-16     6.8e-17
  
          Cephes Math Library, Release 2.8:  June, 2000
          Copyright 1984, 1987, 2000 by Stephen L. Moshier
          *************************************************************************/
          public static double ellipticintegralk(double m)
          {
              double result = 0;
  
              result = ellipticintegralkhighprecision(1.0-m);
              return result;
          }
  
  
          /*************************************************************************
          Complete elliptic integral of the first kind
  
          Approximates the integral
  
  
  
                     pi/2
                      -
                     | |
                     |           dt
          K(m)  =    |    ------------------
                     |                   2
                   | |    sqrt( 1 - m sin t )
                    -
                     0
  
          where m = 1 - m1, using the approximation
  
              P(x)  -  log x Q(x).
  
          The argument m1 is used rather than m so that the logarithmic
          singularity at m = 1 will be shifted to the origin; this
          preserves maximum accuracy.
  
          K(0) = pi/2.
  
          ACCURACY:
  
                               Relative error:
          arithmetic   domain     # trials      peak         rms
             IEEE       0,1        30000       2.5e-16     6.8e-17
  
          Cephes Math Library, Release 2.8:  June, 2000
          Copyright 1984, 1987, 2000 by Stephen L. Moshier
          *************************************************************************/
          public static double ellipticintegralkhighprecision(double m1)
          {
              double result = 0;
              double p = 0;
              double q = 0;
  
              if( (double)(m1)<=(double)(math.machineepsilon) )
              {
                  result = 1.3862943611198906188E0-0.5*Math.Log(m1);
              }
              else
              {
                  p = 1.37982864606273237150E-4;
                  p = p*m1+2.28025724005875567385E-3;
                  p = p*m1+7.97404013220415179367E-3;
                  p = p*m1+9.85821379021226008714E-3;
                  p = p*m1+6.87489687449949877925E-3;
                  p = p*m1+6.18901033637687613229E-3;
                  p = p*m1+8.79078273952743772254E-3;
                  p = p*m1+1.49380448916805252718E-2;
                  p = p*m1+3.08851465246711995998E-2;
                  p = p*m1+9.65735902811690126535E-2;
                  p = p*m1+1.38629436111989062502E0;
                  q = 2.94078955048598507511E-5;
                  q = q*m1+9.14184723865917226571E-4;
                  q = q*m1+5.94058303753167793257E-3;
                  q = q*m1+1.54850516649762399335E-2;
                  q = q*m1+2.39089602715924892727E-2;
                  q = q*m1+3.01204715227604046988E-2;
                  q = q*m1+3.73774314173823228969E-2;
                  q = q*m1+4.88280347570998239232E-2;
                  q = q*m1+7.03124996963957469739E-2;
                  q = q*m1+1.24999999999870820058E-1;
                  q = q*m1+4.99999999999999999821E-1;
                  result = p-q*Math.Log(m1);
              }
              return result;
          }
  
  
          /*************************************************************************
          Incomplete elliptic integral of the first kind F(phi|m)
  
          Approximates the integral
  
  
  
                         phi
                          -
                         | |
                         |           dt
          F(phi_\m)  =    |    ------------------
                         |                   2
                       | |    sqrt( 1 - m sin t )
                        -
                         0
  
          of amplitude phi and modulus m, using the arithmetic -
          geometric mean algorithm.
  
  
  
  
          ACCURACY:
  
          Tested at random points with m in [0, 1] and phi as indicated.
  
                               Relative error:
          arithmetic   domain     # trials      peak         rms
             IEEE     -10,10       200000      7.4e-16     1.0e-16
  
          Cephes Math Library Release 2.8:  June, 2000
          Copyright 1984, 1987, 2000 by Stephen L. Moshier
          *************************************************************************/
          public static double incompleteellipticintegralk(double phi,
              double m)
          {
              double result = 0;
              double a = 0;
              double b = 0;
              double c = 0;
              double e = 0;
              double temp = 0;
              double pio2 = 0;
              double t = 0;
              double k = 0;
              int d = 0;
              int md = 0;
              int s = 0;
              int npio2 = 0;
  
              pio2 = 1.57079632679489661923;
              if( (double)(m)==(double)(0) )
              {
                  result = phi;
                  return result;
              }
              a = 1-m;
              if( (double)(a)==(double)(0) )
              {
                  result = Math.Log(Math.Tan(0.5*(pio2+phi)));
                  return result;
              }
              npio2 = (int)Math.Floor(phi/pio2);
              if( npio2%2!=0 )
              {
                  npio2 = npio2+1;
              }
              if( npio2!=0 )
              {
                  k = ellipticintegralk(1-a);
                  phi = phi-npio2*pio2;
              }
              else
              {
                  k = 0;
              }
              if( (double)(phi)<(double)(0) )
              {
                  phi = -phi;
                  s = -1;
              }
              else
              {
                  s = 0;
              }
              b = Math.Sqrt(a);
              t = Math.Tan(phi);
              if( (double)(Math.Abs(t))>(double)(10) )
              {
                  e = 1.0/(b*t);
                  if( (double)(Math.Abs(e))<(double)(10) )
                  {
                      e = Math.Atan(e);
                      if( npio2==0 )
                      {
                          k = ellipticintegralk(1-a);
                      }
                      temp = k-incompleteellipticintegralk(e, m);
                      if( s<0 )
                      {
                          temp = -temp;
                      }
                      result = temp+npio2*k;
                      return result;
                  }
              }
              a = 1.0;
              c = Math.Sqrt(m);
              d = 1;
              md = 0;
              while( (double)(Math.Abs(c/a))>(double)(math.machineepsilon) )
              {
                  temp = b/a;
                  phi = phi+Math.Atan(t*temp)+md*Math.PI;
                  md = (int)((phi+pio2)/Math.PI);
                  t = t*(1.0+temp)/(1.0-temp*t*t);
                  c = 0.5*(a-b);
                  temp = Math.Sqrt(a*b);
                  a = 0.5*(a+b);
                  b = temp;
                  d = d+d;
              }
              temp = (Math.Atan(t)+md*Math.PI)/(d*a);
              if( s<0 )
              {
                  temp = -temp;
              }
              result = temp+npio2*k;
              return result;
          }
  
  
          /*************************************************************************
          Complete elliptic integral of the second kind
  
          Approximates the integral
  
  
                     pi/2
                      -
                     | |                 2
          E(m)  =    |    sqrt( 1 - m sin t ) dt
                   | |
                    -
                     0
  
          using the approximation
  
               P(x)  -  x log x Q(x).
  
          ACCURACY:
  
                               Relative error:
          arithmetic   domain     # trials      peak         rms
             IEEE       0, 1       10000       2.1e-16     7.3e-17
  
          Cephes Math Library, Release 2.8: June, 2000
          Copyright 1984, 1987, 1989, 2000 by Stephen L. Moshier
          *************************************************************************/
          public static double ellipticintegrale(double m)
          {
              double result = 0;
              double p = 0;
              double q = 0;
  
              alglib.ap.assert((double)(m)>=(double)(0) && (double)(m)<=(double)(1), "Domain error in EllipticIntegralE: m<0 or m>1");
              m = 1-m;
              if( (double)(m)==(double)(0) )
              {
                  result = 1;
                  return result;
              }
              p = 1.53552577301013293365E-4;
              p = p*m+2.50888492163602060990E-3;
              p = p*m+8.68786816565889628429E-3;
              p = p*m+1.07350949056076193403E-2;
              p = p*m+7.77395492516787092951E-3;
              p = p*m+7.58395289413514708519E-3;
              p = p*m+1.15688436810574127319E-2;
              p = p*m+2.18317996015557253103E-2;
              p = p*m+5.68051945617860553470E-2;
              p = p*m+4.43147180560990850618E-1;
              p = p*m+1.00000000000000000299E0;
              q = 3.27954898576485872656E-5;
              q = q*m+1.00962792679356715133E-3;
              q = q*m+6.50609489976927491433E-3;
              q = q*m+1.68862163993311317300E-2;
              q = q*m+2.61769742454493659583E-2;
              q = q*m+3.34833904888224918614E-2;
              q = q*m+4.27180926518931511717E-2;
              q = q*m+5.85936634471101055642E-2;
              q = q*m+9.37499997197644278445E-2;
              q = q*m+2.49999999999888314361E-1;
              result = p-q*m*Math.Log(m);
              return result;
          }
  
  
          /*************************************************************************
          Incomplete elliptic integral of the second kind
  
          Approximates the integral
  
  
                         phi
                          -
                         | |
                         |                   2
          E(phi_\m)  =    |    sqrt( 1 - m sin t ) dt
                         |
                       | |
                        -
                         0
  
          of amplitude phi and modulus m, using the arithmetic -
          geometric mean algorithm.
  
          ACCURACY:
  
          Tested at random arguments with phi in [-10, 10] and m in
          [0, 1].
                               Relative error:
          arithmetic   domain     # trials      peak         rms
             IEEE     -10,10      150000       3.3e-15     1.4e-16
  
          Cephes Math Library Release 2.8:  June, 2000
          Copyright 1984, 1987, 1993, 2000 by Stephen L. Moshier
          *************************************************************************/
          public static double incompleteellipticintegrale(double phi,
              double m)
          {
              double result = 0;
              double pio2 = 0;
              double a = 0;
              double b = 0;
              double c = 0;
              double e = 0;
              double temp = 0;
              double lphi = 0;
              double t = 0;
              double ebig = 0;
              int d = 0;
              int md = 0;
              int npio2 = 0;
              int s = 0;
  
              pio2 = 1.57079632679489661923;
              if( (double)(m)==(double)(0) )
              {
                  result = phi;
                  return result;
              }
              lphi = phi;
              npio2 = (int)Math.Floor(lphi/pio2);
              if( npio2%2!=0 )
              {
                  npio2 = npio2+1;
              }
              lphi = lphi-npio2*pio2;
              if( (double)(lphi)<(double)(0) )
              {
                  lphi = -lphi;
                  s = -1;
              }
              else
              {
                  s = 1;
              }
              a = 1.0-m;
              ebig = ellipticintegrale(m);
              if( (double)(a)==(double)(0) )
              {
                  temp = Math.Sin(lphi);
                  if( s<0 )
                  {
                      temp = -temp;
                  }
                  result = temp+npio2*ebig;
                  return result;
              }
              t = Math.Tan(lphi);
              b = Math.Sqrt(a);
              
              //
              // Thanks to Brian Fitzgerald <fitzgb@mml0.meche.rpi.edu>
              // for pointing out an instability near odd multiples of pi/2
              //
              if( (double)(Math.Abs(t))>(double)(10) )
              {
                  
                  //
                  // Transform the amplitude
                  //
                  e = 1.0/(b*t);
                  
                  //
                  // ... but avoid multiple recursions.
                  //
                  if( (double)(Math.Abs(e))<(double)(10) )
                  {
                      e = Math.Atan(e);
                      temp = ebig+m*Math.Sin(lphi)*Math.Sin(e)-incompleteellipticintegrale(e, m);
                      if( s<0 )
                      {
                          temp = -temp;
                      }
                      result = temp+npio2*ebig;
                      return result;
                  }
              }
              c = Math.Sqrt(m);
              a = 1.0;
              d = 1;
              e = 0.0;
              md = 0;
              while( (double)(Math.Abs(c/a))>(double)(math.machineepsilon) )
              {
                  temp = b/a;
                  lphi = lphi+Math.Atan(t*temp)+md*Math.PI;
                  md = (int)((lphi+pio2)/Math.PI);
                  t = t*(1.0+temp)/(1.0-temp*t*t);
                  c = 0.5*(a-b);
                  temp = Math.Sqrt(a*b);
                  a = 0.5*(a+b);
                  b = temp;
                  d = d+d;
                  e = e+c*Math.Sin(lphi);
              }
              temp = ebig/ellipticintegralk(m);
              temp = temp*((Math.Atan(t)+md*Math.PI)/(d*a));
              temp = temp+e;
              if( s<0 )
              {
                  temp = -temp;
              }
              result = temp+npio2*ebig;
              return result;
          }
  
  
      }
      public class hermite
      {
          /*************************************************************************
          Calculation of the value of the Hermite polynomial.
  
          Parameters:
              n   -   degree, n>=0
              x   -   argument
  
          Result:
              the value of the Hermite polynomial Hn at x
          *************************************************************************/
          public static double hermitecalculate(int n,
              double x)
          {
              double result = 0;
              int i = 0;
              double a = 0;
              double b = 0;
  
              result = 0;
              
              //
              // Prepare A and B
              //
              a = 1;
              b = 2*x;
              
              //
              // Special cases: N=0 or N=1
              //
              if( n==0 )
              {
                  result = a;
                  return result;
              }
              if( n==1 )
              {
                  result = b;
                  return result;
              }
              
              //
              // General case: N>=2
              //
              for(i=2; i<=n; i++)
              {
                  result = 2*x*b-2*(i-1)*a;
                  a = b;
                  b = result;
              }
              return result;
          }
  
  
          /*************************************************************************
          Summation of Hermite polynomials using Clenshaw's recurrence formula.
  
          This routine calculates
              c[0]*H0(x) + c[1]*H1(x) + ... + c[N]*HN(x)
  
          Parameters:
              n   -   degree, n>=0
              x   -   argument
  
          Result:
              the value of the Hermite polynomial at x
          *************************************************************************/
          public static double hermitesum(double[] c,
              int n,
              double x)
          {
              double result = 0;
              double b1 = 0;
              double b2 = 0;
              int i = 0;
  
              b1 = 0;
              b2 = 0;
              result = 0;
              for(i=n; i>=0; i--)
              {
                  result = 2*(x*b1-(i+1)*b2)+c[i];
                  b2 = b1;
                  b1 = result;
              }
              return result;
          }
  
  
          /*************************************************************************
          Representation of Hn as C[0] + C[1]*X + ... + C[N]*X^N
  
          Input parameters:
              N   -   polynomial degree, n>=0
  
          Output parameters:
              C   -   coefficients
          *************************************************************************/
          public static void hermitecoefficients(int n,
              ref double[] c)
          {
              int i = 0;
  
              c = new double[0];
  
              c = new double[n+1];
              for(i=0; i<=n; i++)
              {
                  c[i] = 0;
              }
              c[n] = Math.Exp(n*Math.Log(2));
              for(i=0; i<=n/2-1; i++)
              {
                  c[n-2*(i+1)] = -(c[n-2*i]*(n-2*i)*(n-2*i-1)/4/(i+1));
              }
          }
  
  
      }
      public class dawson
      {
          /*************************************************************************
          Dawson's Integral
  
          Approximates the integral
  
                                      x
                                      -
                               2     | |        2
           dawsn(x)  =  exp( -x  )   |    exp( t  ) dt
                                   | |
                                    -
                                    0
  
          Three different rational approximations are employed, for
          the intervals 0 to 3.25; 3.25 to 6.25; and 6.25 up.
  
          ACCURACY:
  
                               Relative error:
          arithmetic   domain     # trials      peak         rms
             IEEE      0,10        10000       6.9e-16     1.0e-16
  
          Cephes Math Library Release 2.8:  June, 2000
          Copyright 1984, 1987, 1989, 2000 by Stephen L. Moshier
          *************************************************************************/
          public static double dawsonintegral(double x)
          {
              double result = 0;
              double x2 = 0;
              double y = 0;
              int sg = 0;
              double an = 0;
              double ad = 0;
              double bn = 0;
              double bd = 0;
              double cn = 0;
              double cd = 0;
  
              sg = 1;
              if( (double)(x)<(double)(0) )
              {
                  sg = -1;
                  x = -x;
              }
              if( (double)(x)<(double)(3.25) )
              {
                  x2 = x*x;
                  an = 1.13681498971755972054E-11;
                  an = an*x2+8.49262267667473811108E-10;
                  an = an*x2+1.94434204175553054283E-8;
                  an = an*x2+9.53151741254484363489E-7;
                  an = an*x2+3.07828309874913200438E-6;
                  an = an*x2+3.52513368520288738649E-4;
                  an = an*x2+-8.50149846724410912031E-4;
                  an = an*x2+4.22618223005546594270E-2;
                  an = an*x2+-9.17480371773452345351E-2;
                  an = an*x2+9.99999999999999994612E-1;
                  ad = 2.40372073066762605484E-11;
                  ad = ad*x2+1.48864681368493396752E-9;
                  ad = ad*x2+5.21265281010541664570E-8;
                  ad = ad*x2+1.27258478273186970203E-6;
                  ad = ad*x2+2.32490249820789513991E-5;
                  ad = ad*x2+3.25524741826057911661E-4;
                  ad = ad*x2+3.48805814657162590916E-3;
                  ad = ad*x2+2.79448531198828973716E-2;
                  ad = ad*x2+1.58874241960120565368E-1;
                  ad = ad*x2+5.74918629489320327824E-1;
                  ad = ad*x2+1.00000000000000000539E0;
                  y = x*an/ad;
                  result = sg*y;
                  return result;
              }
              x2 = 1.0/(x*x);
              if( (double)(x)<(double)(6.25) )
              {
                  bn = 5.08955156417900903354E-1;
                  bn = bn*x2-2.44754418142697847934E-1;
                  bn = bn*x2+9.41512335303534411857E-2;
                  bn = bn*x2-2.18711255142039025206E-2;
                  bn = bn*x2+3.66207612329569181322E-3;
                  bn = bn*x2-4.23209114460388756528E-4;
                  bn = bn*x2+3.59641304793896631888E-5;
                  bn = bn*x2-2.14640351719968974225E-6;
                  bn = bn*x2+9.10010780076391431042E-8;
                  bn = bn*x2-2.40274520828250956942E-9;
                  bn = bn*x2+3.59233385440928410398E-11;
                  bd = 1.00000000000000000000E0;
                  bd = bd*x2-6.31839869873368190192E-1;
                  bd = bd*x2+2.36706788228248691528E-1;
                  bd = bd*x2-5.31806367003223277662E-2;
                  bd = bd*x2+8.48041718586295374409E-3;
                  bd = bd*x2-9.47996768486665330168E-4;
                  bd = bd*x2+7.81025592944552338085E-5;
                  bd = bd*x2-4.55875153252442634831E-6;
                  bd = bd*x2+1.89100358111421846170E-7;
                  bd = bd*x2-4.91324691331920606875E-9;
                  bd = bd*x2+7.18466403235734541950E-11;
                  y = 1.0/x+x2*bn/(bd*x);
                  result = sg*0.5*y;
                  return result;
              }
              if( (double)(x)>(double)(1.0E9) )
              {
                  result = sg*0.5/x;
                  return result;
              }
              cn = -5.90592860534773254987E-1;
              cn = cn*x2+6.29235242724368800674E-1;
              cn = cn*x2-1.72858975380388136411E-1;
              cn = cn*x2+1.64837047825189632310E-2;
              cn = cn*x2-4.86827613020462700845E-4;
              cd = 1.00000000000000000000E0;
              cd = cd*x2-2.69820057197544900361E0;
              cd = cd*x2+1.73270799045947845857E0;
              cd = cd*x2-3.93708582281939493482E-1;
              cd = cd*x2+3.44278924041233391079E-2;
              cd = cd*x2-9.73655226040941223894E-4;
              y = 1.0/x+x2*cn/(cd*x);
              result = sg*0.5*y;
              return result;
          }
  
  
      }
      public class trigintegrals
      {
          /*************************************************************************
          Sine and cosine integrals
  
          Evaluates the integrals
  
                                   x
                                   -
                                  |  cos t - 1
            Ci(x) = eul + ln x +  |  --------- dt,
                                  |      t
                                 -
                                  0
                      x
                      -
                     |  sin t
            Si(x) =  |  ----- dt
                     |    t
                    -
                     0
  
          where eul = 0.57721566490153286061 is Euler's constant.
          The integrals are approximated by rational functions.
          For x > 8 auxiliary functions f(x) and g(x) are employed
          such that
  
          Ci(x) = f(x) sin(x) - g(x) cos(x)
          Si(x) = pi/2 - f(x) cos(x) - g(x) sin(x)
  
  
          ACCURACY:
             Test interval = [0,50].
          Absolute error, except relative when > 1:
          arithmetic   function   # trials      peak         rms
             IEEE        Si        30000       4.4e-16     7.3e-17
             IEEE        Ci        30000       6.9e-16     5.1e-17
  
          Cephes Math Library Release 2.1:  January, 1989
          Copyright 1984, 1987, 1989 by Stephen L. Moshier
          *************************************************************************/
          public static void sinecosineintegrals(double x,
              ref double si,
              ref double ci)
          {
              double z = 0;
              double c = 0;
              double s = 0;
              double f = 0;
              double g = 0;
              int sg = 0;
              double sn = 0;
              double sd = 0;
              double cn = 0;
              double cd = 0;
              double fn = 0;
              double fd = 0;
              double gn = 0;
              double gd = 0;
  
              si = 0;
              ci = 0;
  
              if( (double)(x)<(double)(0) )
              {
                  sg = -1;
                  x = -x;
              }
              else
              {
                  sg = 0;
              }
              if( (double)(x)==(double)(0) )
              {
                  si = 0;
                  ci = -math.maxrealnumber;
                  return;
              }
              if( (double)(x)>(double)(1.0E9) )
              {
                  si = 1.570796326794896619-Math.Cos(x)/x;
                  ci = Math.Sin(x)/x;
                  return;
              }
              if( (double)(x)<=(double)(4) )
              {
                  z = x*x;
                  sn = -8.39167827910303881427E-11;
                  sn = sn*z+4.62591714427012837309E-8;
                  sn = sn*z-9.75759303843632795789E-6;
                  sn = sn*z+9.76945438170435310816E-4;
                  sn = sn*z-4.13470316229406538752E-2;
                  sn = sn*z+1.00000000000000000302E0;
                  sd = 2.03269266195951942049E-12;
                  sd = sd*z+1.27997891179943299903E-9;
                  sd = sd*z+4.41827842801218905784E-7;
                  sd = sd*z+9.96412122043875552487E-5;
                  sd = sd*z+1.42085239326149893930E-2;
                  sd = sd*z+9.99999999999999996984E-1;
                  s = x*sn/sd;
                  cn = 2.02524002389102268789E-11;
                  cn = cn*z-1.35249504915790756375E-8;
                  cn = cn*z+3.59325051419993077021E-6;
                  cn = cn*z-4.74007206873407909465E-4;
                  cn = cn*z+2.89159652607555242092E-2;
                  cn = cn*z-1.00000000000000000080E0;
                  cd = 4.07746040061880559506E-12;
                  cd = cd*z+3.06780997581887812692E-9;
                  cd = cd*z+1.23210355685883423679E-6;
                  cd = cd*z+3.17442024775032769882E-4;
                  cd = cd*z+5.10028056236446052392E-2;
                  cd = cd*z+4.00000000000000000080E0;
                  c = z*cn/cd;
                  if( sg!=0 )
                  {
                      s = -s;
                  }
                  si = s;
                  ci = 0.57721566490153286061+Math.Log(x)+c;
                  return;
              }
              s = Math.Sin(x);
              c = Math.Cos(x);
              z = 1.0/(x*x);
              if( (double)(x)<(double)(8) )
              {
                  fn = 4.23612862892216586994E0;
                  fn = fn*z+5.45937717161812843388E0;
                  fn = fn*z+1.62083287701538329132E0;
                  fn = fn*z+1.67006611831323023771E-1;
                  fn = fn*z+6.81020132472518137426E-3;
                  fn = fn*z+1.08936580650328664411E-4;
                  fn = fn*z+5.48900223421373614008E-7;
                  fd = 1.00000000000000000000E0;
                  fd = fd*z+8.16496634205391016773E0;
                  fd = fd*z+7.30828822505564552187E0;
                  fd = fd*z+1.86792257950184183883E0;
                  fd = fd*z+1.78792052963149907262E-1;
                  fd = fd*z+7.01710668322789753610E-3;
                  fd = fd*z+1.10034357153915731354E-4;
                  fd = fd*z+5.48900252756255700982E-7;
                  f = fn/(x*fd);
                  gn = 8.71001698973114191777E-2;
                  gn = gn*z+6.11379109952219284151E-1;
                  gn = gn*z+3.97180296392337498885E-1;
                  gn = gn*z+7.48527737628469092119E-2;
                  gn = gn*z+5.38868681462177273157E-3;
                  gn = gn*z+1.61999794598934024525E-4;
                  gn = gn*z+1.97963874140963632189E-6;
                  gn = gn*z+7.82579040744090311069E-9;
                  gd = 1.00000000000000000000E0;
                  gd = gd*z+1.64402202413355338886E0;
                  gd = gd*z+6.66296701268987968381E-1;
                  gd = gd*z+9.88771761277688796203E-2;
                  gd = gd*z+6.22396345441768420760E-3;
                  gd = gd*z+1.73221081474177119497E-4;
                  gd = gd*z+2.02659182086343991969E-6;
                  gd = gd*z+7.82579218933534490868E-9;
                  g = z*gn/gd;
              }
              else
              {
                  fn = 4.55880873470465315206E-1;
                  fn = fn*z+7.13715274100146711374E-1;
                  fn = fn*z+1.60300158222319456320E-1;
                  fn = fn*z+1.16064229408124407915E-2;
                  fn = fn*z+3.49556442447859055605E-4;
                  fn = fn*z+4.86215430826454749482E-6;
                  fn = fn*z+3.20092790091004902806E-8;
                  fn = fn*z+9.41779576128512936592E-11;
                  fn = fn*z+9.70507110881952024631E-14;
                  fd = 1.00000000000000000000E0;
                  fd = fd*z+9.17463611873684053703E-1;
                  fd = fd*z+1.78685545332074536321E-1;
                  fd = fd*z+1.22253594771971293032E-2;
                  fd = fd*z+3.58696481881851580297E-4;
                  fd = fd*z+4.92435064317881464393E-6;
                  fd = fd*z+3.21956939101046018377E-8;
                  fd = fd*z+9.43720590350276732376E-11;
                  fd = fd*z+9.70507110881952025725E-14;
                  f = fn/(x*fd);
                  gn = 6.97359953443276214934E-1;
                  gn = gn*z+3.30410979305632063225E-1;
                  gn = gn*z+3.84878767649974295920E-2;
                  gn = gn*z+1.71718239052347903558E-3;
                  gn = gn*z+3.48941165502279436777E-5;
                  gn = gn*z+3.47131167084116673800E-7;
                  gn = gn*z+1.70404452782044526189E-9;
                  gn = gn*z+3.85945925430276600453E-12;
                  gn = gn*z+3.14040098946363334640E-15;
                  gd = 1.00000000000000000000E0;
                  gd = gd*z+1.68548898811011640017E0;
                  gd = gd*z+4.87852258695304967486E-1;
                  gd = gd*z+4.67913194259625806320E-2;
                  gd = gd*z+1.90284426674399523638E-3;
                  gd = gd*z+3.68475504442561108162E-5;
                  gd = gd*z+3.57043223443740838771E-7;
                  gd = gd*z+1.72693748966316146736E-9;
                  gd = gd*z+3.87830166023954706752E-12;
                  gd = gd*z+3.14040098946363335242E-15;
                  g = z*gn/gd;
              }
              si = 1.570796326794896619-f*c-g*s;
              if( sg!=0 )
              {
                  si = -si;
              }
              ci = f*s-g*c;
          }
  
  
          /*************************************************************************
          Hyperbolic sine and cosine integrals
  
          Approximates the integrals
  
                                     x
                                     -
                                    | |   cosh t - 1
            Chi(x) = eul + ln x +   |    -----------  dt,
                                  | |          t
                                   -
                                   0
  
                        x
                        -
                       | |  sinh t
            Shi(x) =   |    ------  dt
                     | |       t
                      -
                      0
  
          where eul = 0.57721566490153286061 is Euler's constant.
          The integrals are evaluated by power series for x < 8
          and by Chebyshev expansions for x between 8 and 88.
          For large x, both functions approach exp(x)/2x.
          Arguments greater than 88 in magnitude return MAXNUM.
  
  
          ACCURACY:
  
          Test interval 0 to 88.
                               Relative error:
          arithmetic   function  # trials      peak         rms
             IEEE         Shi      30000       6.9e-16     1.6e-16
                 Absolute error, except relative when |Chi| > 1:
             IEEE         Chi      30000       8.4e-16     1.4e-16
  
          Cephes Math Library Release 2.8:  June, 2000
          Copyright 1984, 1987, 2000 by Stephen L. Moshier
          *************************************************************************/
          public static void hyperbolicsinecosineintegrals(double x,
              ref double shi,
              ref double chi)
          {
              double k = 0;
              double z = 0;
              double c = 0;
              double s = 0;
              double a = 0;
              int sg = 0;
              double b0 = 0;
              double b1 = 0;
              double b2 = 0;
  
              shi = 0;
              chi = 0;
  
              if( (double)(x)<(double)(0) )
              {
                  sg = -1;
                  x = -x;
              }
              else
              {
                  sg = 0;
              }
              if( (double)(x)==(double)(0) )
              {
                  shi = 0;
                  chi = -math.maxrealnumber;
                  return;
              }
              if( (double)(x)<(double)(8.0) )
              {
                  z = x*x;
                  a = 1.0;
                  s = 1.0;
                  c = 0.0;
                  k = 2.0;
                  do
                  {
                      a = a*z/k;
                      c = c+a/k;
                      k = k+1.0;
                      a = a/k;
                      s = s+a/k;
                      k = k+1.0;
                  }
                  while( (double)(Math.Abs(a/s))>=(double)(math.machineepsilon) );
                  s = s*x;
              }
              else
              {
                  if( (double)(x)<(double)(18.0) )
                  {
                      a = (576.0/x-52.0)/10.0;
                      k = Math.Exp(x)/x;
                      b0 = 1.83889230173399459482E-17;
                      b1 = 0.0;
                      chebiterationshichi(a, -9.55485532279655569575E-17, ref b0, ref b1, ref b2);
                      chebiterationshichi(a, 2.04326105980879882648E-16, ref b0, ref b1, ref b2);
                      chebiterationshichi(a, 1.09896949074905343022E-15, ref b0, ref b1, ref b2);
                      chebiterationshichi(a, -1.31313534344092599234E-14, ref b0, ref b1, ref b2);
                      chebiterationshichi(a, 5.93976226264314278932E-14, ref b0, ref b1, ref b2);
                      chebiterationshichi(a, -3.47197010497749154755E-14, ref b0, ref b1, ref b2);
                      chebiterationshichi(a, -1.40059764613117131000E-12, ref b0, ref b1, ref b2);
                      chebiterationshichi(a, 9.49044626224223543299E-12, ref b0, ref b1, ref b2);
                      chebiterationshichi(a, -1.61596181145435454033E-11, ref b0, ref b1, ref b2);
                      chebiterationshichi(a, -1.77899784436430310321E-10, ref b0, ref b1, ref b2);
                      chebiterationshichi(a, 1.35455469767246947469E-9, ref b0, ref b1, ref b2);
                      chebiterationshichi(a, -1.03257121792819495123E-9, ref b0, ref b1, ref b2);
                      chebiterationshichi(a, -3.56699611114982536845E-8, ref b0, ref b1, ref b2);
                      chebiterationshichi(a, 1.44818877384267342057E-7, ref b0, ref b1, ref b2);
                      chebiterationshichi(a, 7.82018215184051295296E-7, ref b0, ref b1, ref b2);
                      chebiterationshichi(a, -5.39919118403805073710E-6, ref b0, ref b1, ref b2);
                      chebiterationshichi(a, -3.12458202168959833422E-5, ref b0, ref b1, ref b2);
                      chebiterationshichi(a, 8.90136741950727517826E-5, ref b0, ref b1, ref b2);
                      chebiterationshichi(a, 2.02558474743846862168E-3, ref b0, ref b1, ref b2);
                      chebiterationshichi(a, 2.96064440855633256972E-2, ref b0, ref b1, ref b2);
                      chebiterationshichi(a, 1.11847751047257036625E0, ref b0, ref b1, ref b2);
                      s = k*0.5*(b0-b2);
                      b0 = -8.12435385225864036372E-18;
                      b1 = 0.0;
                      chebiterationshichi(a, 2.17586413290339214377E-17, ref b0, ref b1, ref b2);
                      chebiterationshichi(a, 5.22624394924072204667E-17, ref b0, ref b1, ref b2);
                      chebiterationshichi(a, -9.48812110591690559363E-16, ref b0, ref b1, ref b2);
                      chebiterationshichi(a, 5.35546311647465209166E-15, ref b0, ref b1, ref b2);
                      chebiterationshichi(a, -1.21009970113732918701E-14, ref b0, ref b1, ref b2);
                      chebiterationshichi(a, -6.00865178553447437951E-14, ref b0, ref b1, ref b2);
                      chebiterationshichi(a, 7.16339649156028587775E-13, ref b0, ref b1, ref b2);
                      chebiterationshichi(a, -2.93496072607599856104E-12, ref b0, ref b1, ref b2);
                      chebiterationshichi(a, -1.40359438136491256904E-12, ref b0, ref b1, ref b2);
                      chebiterationshichi(a, 8.76302288609054966081E-11, ref b0, ref b1, ref b2);
                      chebiterationshichi(a, -4.40092476213282340617E-10, ref b0, ref b1, ref b2);
                      chebiterationshichi(a, -1.87992075640569295479E-10, ref b0, ref b1, ref b2);
                      chebiterationshichi(a, 1.31458150989474594064E-8, ref b0, ref b1, ref b2);
                      chebiterationshichi(a, -4.75513930924765465590E-8, ref b0, ref b1, ref b2);
                      chebiterationshichi(a, -2.21775018801848880741E-7, ref b0, ref b1, ref b2);
                      chebiterationshichi(a, 1.94635531373272490962E-6, ref b0, ref b1, ref b2);
                      chebiterationshichi(a, 4.33505889257316408893E-6, ref b0, ref b1, ref b2);
                      chebiterationshichi(a, -6.13387001076494349496E-5, ref b0, ref b1, ref b2);
                      chebiterationshichi(a, -3.13085477492997465138E-4, ref b0, ref b1, ref b2);
                      chebiterationshichi(a, 4.97164789823116062801E-4, ref b0, ref b1, ref b2);
                      chebiterationshichi(a, 2.64347496031374526641E-2, ref b0, ref b1, ref b2);
                      chebiterationshichi(a, 1.11446150876699213025E0, ref b0, ref b1, ref b2);
                      c = k*0.5*(b0-b2);
                  }
                  else
                  {
                      if( (double)(x)<=(double)(88.0) )
                      {
                          a = (6336.0/x-212.0)/70.0;
                          k = Math.Exp(x)/x;
                          b0 = -1.05311574154850938805E-17;
                          b1 = 0.0;
                          chebiterationshichi(a, 2.62446095596355225821E-17, ref b0, ref b1, ref b2);
                          chebiterationshichi(a, 8.82090135625368160657E-17, ref b0, ref b1, ref b2);
                          chebiterationshichi(a, -3.38459811878103047136E-16, ref b0, ref b1, ref b2);
                          chebiterationshichi(a, -8.30608026366935789136E-16, ref b0, ref b1, ref b2);
                          chebiterationshichi(a, 3.93397875437050071776E-15, ref b0, ref b1, ref b2);
                          chebiterationshichi(a, 1.01765565969729044505E-14, ref b0, ref b1, ref b2);
                          chebiterationshichi(a, -4.21128170307640802703E-14, ref b0, ref b1, ref b2);
                          chebiterationshichi(a, -1.60818204519802480035E-13, ref b0, ref b1, ref b2);
                          chebiterationshichi(a, 3.34714954175994481761E-13, ref b0, ref b1, ref b2);
                          chebiterationshichi(a, 2.72600352129153073807E-12, ref b0, ref b1, ref b2);
                          chebiterationshichi(a, 1.66894954752839083608E-12, ref b0, ref b1, ref b2);
                          chebiterationshichi(a, -3.49278141024730899554E-11, ref b0, ref b1, ref b2);
                          chebiterationshichi(a, -1.58580661666482709598E-10, ref b0, ref b1, ref b2);
                          chebiterationshichi(a, -1.79289437183355633342E-10, ref b0, ref b1, ref b2);
                          chebiterationshichi(a, 1.76281629144264523277E-9, ref b0, ref b1, ref b2);
                          chebiterationshichi(a, 1.69050228879421288846E-8, ref b0, ref b1, ref b2);
                          chebiterationshichi(a, 1.25391771228487041649E-7, ref b0, ref b1, ref b2);
                          chebiterationshichi(a, 1.16229947068677338732E-6, ref b0, ref b1, ref b2);
                          chebiterationshichi(a, 1.61038260117376323993E-5, ref b0, ref b1, ref b2);
                          chebiterationshichi(a, 3.49810375601053973070E-4, ref b0, ref b1, ref b2);
                          chebiterationshichi(a, 1.28478065259647610779E-2, ref b0, ref b1, ref b2);
                          chebiterationshichi(a, 1.03665722588798326712E0, ref b0, ref b1, ref b2);
                          s = k*0.5*(b0-b2);
                          b0 = 8.06913408255155572081E-18;
                          b1 = 0.0;
                          chebiterationshichi(a, -2.08074168180148170312E-17, ref b0, ref b1, ref b2);
                          chebiterationshichi(a, -5.98111329658272336816E-17, ref b0, ref b1, ref b2);
                          chebiterationshichi(a, 2.68533951085945765591E-16, ref b0, ref b1, ref b2);
                          chebiterationshichi(a, 4.52313941698904694774E-16, ref b0, ref b1, ref b2);
                          chebiterationshichi(a, -3.10734917335299464535E-15, ref b0, ref b1, ref b2);
                          chebiterationshichi(a, -4.42823207332531972288E-15, ref b0, ref b1, ref b2);
                          chebiterationshichi(a, 3.49639695410806959872E-14, ref b0, ref b1, ref b2);
                          chebiterationshichi(a, 6.63406731718911586609E-14, ref b0, ref b1, ref b2);
                          chebiterationshichi(a, -3.71902448093119218395E-13, ref b0, ref b1, ref b2);
                          chebiterationshichi(a, -1.27135418132338309016E-12, ref b0, ref b1, ref b2);
                          chebiterationshichi(a, 2.74851141935315395333E-12, ref b0, ref b1, ref b2);
                          chebiterationshichi(a, 2.33781843985453438400E-11, ref b0, ref b1, ref b2);
                          chebiterationshichi(a, 2.71436006377612442764E-11, ref b0, ref b1, ref b2);
                          chebiterationshichi(a, -2.56600180000355990529E-10, ref b0, ref b1, ref b2);
                          chebiterationshichi(a, -1.61021375163803438552E-9, ref b0, ref b1, ref b2);
                          chebiterationshichi(a, -4.72543064876271773512E-9, ref b0, ref b1, ref b2);
                          chebiterationshichi(a, -3.00095178028681682282E-9, ref b0, ref b1, ref b2);
                          chebiterationshichi(a, 7.79387474390914922337E-8, ref b0, ref b1, ref b2);
                          chebiterationshichi(a, 1.06942765566401507066E-6, ref b0, ref b1, ref b2);
                          chebiterationshichi(a, 1.59503164802313196374E-5, ref b0, ref b1, ref b2);
                          chebiterationshichi(a, 3.49592575153777996871E-4, ref b0, ref b1, ref b2);
                          chebiterationshichi(a, 1.28475387530065247392E-2, ref b0, ref b1, ref b2);
                          chebiterationshichi(a, 1.03665693917934275131E0, ref b0, ref b1, ref b2);
                          c = k*0.5*(b0-b2);
                      }
                      else
                      {
                          if( sg!=0 )
                          {
                              shi = -math.maxrealnumber;
                          }
                          else
                          {
                              shi = math.maxrealnumber;
                          }
                          chi = math.maxrealnumber;
                          return;
                      }
                  }
              }
              if( sg!=0 )
              {
                  s = -s;
              }
              shi = s;
              chi = 0.57721566490153286061+Math.Log(x)+c;
          }
  
  
          private static void chebiterationshichi(double x,
              double c,
              ref double b0,
              ref double b1,
              ref double b2)
          {
              b2 = b1;
              b1 = b0;
              b0 = x*b1-b2+c;
          }
  
  
      }
      public class poissondistr
      {
          /*************************************************************************
          Poisson distribution
  
          Returns the sum of the first k+1 terms of the Poisson
          distribution:
  
            k         j
            --   -m  m
            >   e    --
            --       j!
           j=0
  
          The terms are not summed directly; instead the incomplete
          gamma integral is employed, according to the relation
  
          y = pdtr( k, m ) = igamc( k+1, m ).
  
          The arguments must both be positive.
          ACCURACY:
  
          See incomplete gamma function
  
          Cephes Math Library Release 2.8:  June, 2000
          Copyright 1984, 1987, 1995, 2000 by Stephen L. Moshier
          *************************************************************************/
          public static double poissondistribution(int k,
              double m)
          {
              double result = 0;
  
              alglib.ap.assert(k>=0 && (double)(m)>(double)(0), "Domain error in PoissonDistribution");
              result = igammaf.incompletegammac(k+1, m);
              return result;
          }
  
  
          /*************************************************************************
          Complemented Poisson distribution
  
          Returns the sum of the terms k+1 to infinity of the Poisson
          distribution:
  
           inf.       j
            --   -m  m
            >   e    --
            --       j!
           j=k+1
  
          The terms are not summed directly; instead the incomplete
          gamma integral is employed, according to the formula
  
          y = pdtrc( k, m ) = igam( k+1, m ).
  
          The arguments must both be positive.
  
          ACCURACY:
  
          See incomplete gamma function
  
          Cephes Math Library Release 2.8:  June, 2000
          Copyright 1984, 1987, 1995, 2000 by Stephen L. Moshier
          *************************************************************************/
          public static double poissoncdistribution(int k,
              double m)
          {
              double result = 0;
  
              alglib.ap.assert(k>=0 && (double)(m)>(double)(0), "Domain error in PoissonDistributionC");
              result = igammaf.incompletegamma(k+1, m);
              return result;
          }
  
  
          /*************************************************************************
          Inverse Poisson distribution
  
          Finds the Poisson variable x such that the integral
          from 0 to x of the Poisson density is equal to the
          given probability y.
  
          This is accomplished using the inverse gamma integral
          function and the relation
  
             m = igami( k+1, y ).
  
          ACCURACY:
  
          See inverse incomplete gamma function
  
          Cephes Math Library Release 2.8:  June, 2000
          Copyright 1984, 1987, 1995, 2000 by Stephen L. Moshier
          *************************************************************************/
          public static double invpoissondistribution(int k,
              double y)
          {
              double result = 0;
  
              alglib.ap.assert((k>=0 && (double)(y)>=(double)(0)) && (double)(y)<(double)(1), "Domain error in InvPoissonDistribution");
              result = igammaf.invincompletegammac(k+1, y);
              return result;
          }
  
  
      }
      public class bessel
      {
          /*************************************************************************
          Bessel function of order zero
  
          Returns Bessel function of order zero of the argument.
  
          The domain is divided into the intervals [0, 5] and
          (5, infinity). In the first interval the following rational
          approximation is used:
  
  
                 2         2
          (w - r  ) (w - r  ) P (w) / Q (w)
                1         2    3       8
  
                     2
          where w = x  and the two r's are zeros of the function.
  
          In the second interval, the Hankel asymptotic expansion
          is employed with two rational functions of degree 6/6
          and 7/7.
  
          ACCURACY:
  
                               Absolute error:
          arithmetic   domain     # trials      peak         rms
             IEEE      0, 30       60000       4.2e-16     1.1e-16
  
          Cephes Math Library Release 2.8:  June, 2000
          Copyright 1984, 1987, 1989, 2000 by Stephen L. Moshier
          *************************************************************************/
          public static double besselj0(double x)
          {
              double result = 0;
              double xsq = 0;
              double nn = 0;
              double pzero = 0;
              double qzero = 0;
              double p1 = 0;
              double q1 = 0;
  
              if( (double)(x)<(double)(0) )
              {
                  x = -x;
              }
              if( (double)(x)>(double)(8.0) )
              {
                  besselasympt0(x, ref pzero, ref qzero);
                  nn = x-Math.PI/4;
                  result = Math.Sqrt(2/Math.PI/x)*(pzero*Math.Cos(nn)-qzero*Math.Sin(nn));
                  return result;
              }
              xsq = math.sqr(x);
              p1 = 26857.86856980014981415848441;
              p1 = -40504123.71833132706360663322+xsq*p1;
              p1 = 25071582855.36881945555156435+xsq*p1;
              p1 = -8085222034853.793871199468171+xsq*p1;
              p1 = 1434354939140344.111664316553+xsq*p1;
              p1 = -136762035308817138.6865416609+xsq*p1;
              p1 = 6382059341072356562.289432465+xsq*p1;
              p1 = -117915762910761053603.8440800+xsq*p1;
              p1 = 493378725179413356181.6813446+xsq*p1;
              q1 = 1.0;
              q1 = 1363.063652328970604442810507+xsq*q1;
              q1 = 1114636.098462985378182402543+xsq*q1;
              q1 = 669998767.2982239671814028660+xsq*q1;
              q1 = 312304311494.1213172572469442+xsq*q1;
              q1 = 112775673967979.8507056031594+xsq*q1;
              q1 = 30246356167094626.98627330784+xsq*q1;
              q1 = 5428918384092285160.200195092+xsq*q1;
              q1 = 493378725179413356211.3278438+xsq*q1;
              result = p1/q1;
              return result;
          }
  
  
          /*************************************************************************
          Bessel function of order one
  
          Returns Bessel function of order one of the argument.
  
          The domain is divided into the intervals [0, 8] and
          (8, infinity). In the first interval a 24 term Chebyshev
          expansion is used. In the second, the asymptotic
          trigonometric representation is employed using two
          rational functions of degree 5/5.
  
          ACCURACY:
  
                               Absolute error:
          arithmetic   domain      # trials      peak         rms
             IEEE      0, 30       30000       2.6e-16     1.1e-16
  
          Cephes Math Library Release 2.8:  June, 2000
          Copyright 1984, 1987, 1989, 2000 by Stephen L. Moshier
          *************************************************************************/
          public static double besselj1(double x)
          {
              double result = 0;
              double s = 0;
              double xsq = 0;
              double nn = 0;
              double pzero = 0;
              double qzero = 0;
              double p1 = 0;
              double q1 = 0;
  
              s = Math.Sign(x);
              if( (double)(x)<(double)(0) )
              {
                  x = -x;
              }
              if( (double)(x)>(double)(8.0) )
              {
                  besselasympt1(x, ref pzero, ref qzero);
                  nn = x-3*Math.PI/4;
                  result = Math.Sqrt(2/Math.PI/x)*(pzero*Math.Cos(nn)-qzero*Math.Sin(nn));
                  if( (double)(s)<(double)(0) )
                  {
                      result = -result;
                  }
                  return result;
              }
              xsq = math.sqr(x);
              p1 = 2701.122710892323414856790990;
              p1 = -4695753.530642995859767162166+xsq*p1;
              p1 = 3413234182.301700539091292655+xsq*p1;
              p1 = -1322983480332.126453125473247+xsq*p1;
              p1 = 290879526383477.5409737601689+xsq*p1;
              p1 = -35888175699101060.50743641413+xsq*p1;
              p1 = 2316433580634002297.931815435+xsq*p1;
              p1 = -66721065689249162980.20941484+xsq*p1;
              p1 = 581199354001606143928.050809+xsq*p1;
              q1 = 1.0;
              q1 = 1606.931573481487801970916749+xsq*q1;
              q1 = 1501793.594998585505921097578+xsq*q1;
              q1 = 1013863514.358673989967045588+xsq*q1;
              q1 = 524371026216.7649715406728642+xsq*q1;
              q1 = 208166122130760.7351240184229+xsq*q1;
              q1 = 60920613989175217.46105196863+xsq*q1;
              q1 = 11857707121903209998.37113348+xsq*q1;
              q1 = 1162398708003212287858.529400+xsq*q1;
              result = s*x*p1/q1;
              return result;
          }
  
  
          /*************************************************************************
          Bessel function of integer order
  
          Returns Bessel function of order n, where n is a
          (possibly negative) integer.
  
          The ratio of jn(x) to j0(x) is computed by backward
          recurrence.  First the ratio jn/jn-1 is found by a
          continued fraction expansion.  Then the recurrence
          relating successive orders is applied until j0 or j1 is
          reached.
  
          If n = 0 or 1 the routine for j0 or j1 is called
          directly.
  
          ACCURACY:
  
                               Absolute error:
          arithmetic   range      # trials      peak         rms
             IEEE      0, 30        5000       4.4e-16     7.9e-17
  
  
          Not suitable for large n or x. Use jv() (fractional order) instead.
  
          Cephes Math Library Release 2.8:  June, 2000
          Copyright 1984, 1987, 2000 by Stephen L. Moshier
          *************************************************************************/
          public static double besseljn(int n,
              double x)
          {
              double result = 0;
              double pkm2 = 0;
              double pkm1 = 0;
              double pk = 0;
              double xk = 0;
              double r = 0;
              double ans = 0;
              int k = 0;
              int sg = 0;
  
              if( n<0 )
              {
                  n = -n;
                  if( n%2==0 )
                  {
                      sg = 1;
                  }
                  else
                  {
                      sg = -1;
                  }
              }
              else
              {
                  sg = 1;
              }
              if( (double)(x)<(double)(0) )
              {
                  if( n%2!=0 )
                  {
                      sg = -sg;
                  }
                  x = -x;
              }
              if( n==0 )
              {
                  result = sg*besselj0(x);
                  return result;
              }
              if( n==1 )
              {
                  result = sg*besselj1(x);
                  return result;
              }
              if( n==2 )
              {
                  if( (double)(x)==(double)(0) )
                  {
                      result = 0;
                  }
                  else
                  {
                      result = sg*(2.0*besselj1(x)/x-besselj0(x));
                  }
                  return result;
              }
              if( (double)(x)<(double)(math.machineepsilon) )
              {
                  result = 0;
                  return result;
              }
              k = 53;
              pk = 2*(n+k);
              ans = pk;
              xk = x*x;
              do
              {
                  pk = pk-2.0;
                  ans = pk-xk/ans;
                  k = k-1;
              }
              while( k!=0 );
              ans = x/ans;
              pk = 1.0;
              pkm1 = 1.0/ans;
              k = n-1;
              r = 2*k;
              do
              {
                  pkm2 = (pkm1*r-pk*x)/x;
                  pk = pkm1;
                  pkm1 = pkm2;
                  r = r-2.0;
                  k = k-1;
              }
              while( k!=0 );
              if( (double)(Math.Abs(pk))>(double)(Math.Abs(pkm1)) )
              {
                  ans = besselj1(x)/pk;
              }
              else
              {
                  ans = besselj0(x)/pkm1;
              }
              result = sg*ans;
              return result;
          }
  
  
          /*************************************************************************
          Bessel function of the second kind, order zero
  
          Returns Bessel function of the second kind, of order
          zero, of the argument.
  
          The domain is divided into the intervals [0, 5] and
          (5, infinity). In the first interval a rational approximation
          R(x) is employed to compute
            y0(x)  = R(x)  +   2 * log(x) * j0(x) / PI.
          Thus a call to j0() is required.
  
          In the second interval, the Hankel asymptotic expansion
          is employed with two rational functions of degree 6/6
          and 7/7.
  
  
  
          ACCURACY:
  
           Absolute error, when y0(x) < 1; else relative error:
  
          arithmetic   domain     # trials      peak         rms
             IEEE      0, 30       30000       1.3e-15     1.6e-16
  
          Cephes Math Library Release 2.8:  June, 2000
          Copyright 1984, 1987, 1989, 2000 by Stephen L. Moshier
          *************************************************************************/
          public static double bessely0(double x)
          {
              double result = 0;
              double nn = 0;
              double xsq = 0;
              double pzero = 0;
              double qzero = 0;
              double p4 = 0;
              double q4 = 0;
  
              if( (double)(x)>(double)(8.0) )
              {
                  besselasympt0(x, ref pzero, ref qzero);
                  nn = x-Math.PI/4;
                  result = Math.Sqrt(2/Math.PI/x)*(pzero*Math.Sin(nn)+qzero*Math.Cos(nn));
                  return result;
              }
              xsq = math.sqr(x);
              p4 = -41370.35497933148554125235152;
              p4 = 59152134.65686889654273830069+xsq*p4;
              p4 = -34363712229.79040378171030138+xsq*p4;
              p4 = 10255208596863.94284509167421+xsq*p4;
              p4 = -1648605817185729.473122082537+xsq*p4;
              p4 = 137562431639934407.8571335453+xsq*p4;
              p4 = -5247065581112764941.297350814+xsq*p4;
              p4 = 65874732757195549259.99402049+xsq*p4;
              p4 = -27502866786291095837.01933175+xsq*p4;
              q4 = 1.0;
              q4 = 1282.452772478993804176329391+xsq*q4;
              q4 = 1001702.641288906265666651753+xsq*q4;
              q4 = 579512264.0700729537480087915+xsq*q4;
              q4 = 261306575504.1081249568482092+xsq*q4;
              q4 = 91620380340751.85262489147968+xsq*q4;
              q4 = 23928830434997818.57439356652+xsq*q4;
              q4 = 4192417043410839973.904769661+xsq*q4;
              q4 = 372645883898616588198.9980+xsq*q4;
              result = p4/q4+2/Math.PI*besselj0(x)*Math.Log(x);
              return result;
          }
  
  
          /*************************************************************************
          Bessel function of second kind of order one
  
          Returns Bessel function of the second kind of order one
          of the argument.
  
          The domain is divided into the intervals [0, 8] and
          (8, infinity). In the first interval a 25 term Chebyshev
          expansion is used, and a call to j1() is required.
          In the second, the asymptotic trigonometric representation
          is employed using two rational functions of degree 5/5.
  
          ACCURACY:
  
                               Absolute error:
          arithmetic   domain      # trials      peak         rms
             IEEE      0, 30       30000       1.0e-15     1.3e-16
  
          Cephes Math Library Release 2.8:  June, 2000
          Copyright 1984, 1987, 1989, 2000 by Stephen L. Moshier
          *************************************************************************/
          public static double bessely1(double x)
          {
              double result = 0;
              double nn = 0;
              double xsq = 0;
              double pzero = 0;
              double qzero = 0;
              double p4 = 0;
              double q4 = 0;
  
              if( (double)(x)>(double)(8.0) )
              {
                  besselasympt1(x, ref pzero, ref qzero);
                  nn = x-3*Math.PI/4;
                  result = Math.Sqrt(2/Math.PI/x)*(pzero*Math.Sin(nn)+qzero*Math.Cos(nn));
                  return result;
              }
              xsq = math.sqr(x);
              p4 = -2108847.540133123652824139923;
              p4 = 3639488548.124002058278999428+xsq*p4;
              p4 = -2580681702194.450950541426399+xsq*p4;
              p4 = 956993023992168.3481121552788+xsq*p4;
              p4 = -196588746272214065.8820322248+xsq*p4;
              p4 = 21931073399177975921.11427556+xsq*p4;
              p4 = -1212297555414509577913.561535+xsq*p4;
              p4 = 26554738314348543268942.48968+xsq*p4;
              p4 = -99637534243069222259967.44354+xsq*p4;
              q4 = 1.0;
              q4 = 1612.361029677000859332072312+xsq*q4;
              q4 = 1563282.754899580604737366452+xsq*q4;
              q4 = 1128686837.169442121732366891+xsq*q4;
              q4 = 646534088126.5275571961681500+xsq*q4;
              q4 = 297663212564727.6729292742282+xsq*q4;
              q4 = 108225825940881955.2553850180+xsq*q4;
              q4 = 29549879358971486742.90758119+xsq*q4;
              q4 = 5435310377188854170800.653097+xsq*q4;
              q4 = 508206736694124324531442.4152+xsq*q4;
              result = x*p4/q4+2/Math.PI*(besselj1(x)*Math.Log(x)-1/x);
              return result;
          }
  
  
          /*************************************************************************
          Bessel function of second kind of integer order
  
          Returns Bessel function of order n, where n is a
          (possibly negative) integer.
  
          The function is evaluated by forward recurrence on
          n, starting with values computed by the routines
          y0() and y1().
  
          If n = 0 or 1 the routine for y0 or y1 is called
          directly.
  
          ACCURACY:
                               Absolute error, except relative
                               when y > 1:
          arithmetic   domain     # trials      peak         rms
             IEEE      0, 30       30000       3.4e-15     4.3e-16
  
          Cephes Math Library Release 2.8:  June, 2000
          Copyright 1984, 1987, 2000 by Stephen L. Moshier
          *************************************************************************/
          public static double besselyn(int n,
              double x)
          {
              double result = 0;
              int i = 0;
              double a = 0;
              double b = 0;
              double tmp = 0;
              double s = 0;
  
              s = 1;
              if( n<0 )
              {
                  n = -n;
                  if( n%2!=0 )
                  {
                      s = -1;
                  }
              }
              if( n==0 )
              {
                  result = bessely0(x);
                  return result;
              }
              if( n==1 )
              {
                  result = s*bessely1(x);
                  return result;
              }
              a = bessely0(x);
              b = bessely1(x);
              for(i=1; i<=n-1; i++)
              {
                  tmp = b;
                  b = 2*i/x*b-a;
                  a = tmp;
              }
  
gt:             result = s*b;