/****************************************************************************
*                   chi2.c
*
*  This module contains the function for the chi square distribution.
*
*  All functions in this module are taken from the Cephes math library.
*
*    Cephes Math Library Release 2.0:  April, 1987
*    Copyright 1984, 1987 by Stephen L. Moshier
*    Direct inquiries to 30 Frost Street, Cambridge, MA 02140
*
*  from Persistence of Vision(tm) Ray Tracer
*  Copyright 1996 Persistence of Vision Team
*---------------------------------------------------------------------------
*  NOTICE: This source code file is provided so that users may experiment
*  with enhancements to POV-Ray and to port the software to platforms other
*  than those supported by the POV-Ray Team.  There are strict rules under
*  which you are permitted to use this file.  The rules are in the file
*  named POVLEGAL.DOC which should be distributed with this file. If
*  POVLEGAL.DOC is not available or for more info please contact the POV-Ray
*  Team Coordinator by leaving a message in CompuServe's Graphics Developer's
*  Forum.  The latest version of POV-Ray may be found there as well.
*
* This program is based on the popular DKB raytracer version 2.12.
* DKBTrace was originally written by David K. Buck.
* DKBTrace Ver 2.0-2.12 were written by David K. Buck & Aaron A. Collins.
*
*****************************************************************************/

#include "frame.h"
#include "povproto.h"
#include "chi2.h"

/*
Cephes Math Library Release 2.0:  April, 1987
Copyright 1984, 1987 by Stephen L. Moshier
Direct inquiries to 30 Frost Street, Cambridge, MA 02140
*/

/*****************************************************************************
* Local preprocessor defines
******************************************************************************/

#define MAXLGM 2.556348e305
#define BIG    1.44115188075855872E+17
#define MACHEP 1.38777878078144567553E-17 /* 2**-56 */
#define MAXLOG 8.8029691931113054295988E1 /* log(2**127) */
#define MAXNUM 1.701411834604692317316873e38  /* 2**127 */
#define LOGPI 1.14472988584940017414



/*****************************************************************************
* Local typedefs
******************************************************************************/



/*****************************************************************************
* Local variables
******************************************************************************/

static int sgngam = 0;

/*
 * A[]: Stirling's formula expansion of log gamma
 * B[], C[]: log gamma function between 2 and 3
 */

static DBL A[] =
{
   8.11614167470508450300E-4,
  -5.95061904284301438324E-4,
   7.93650340457716943945E-4,
  -2.77777777730099687205E-3,
   8.33333333333331927722E-2
};

static DBL B[] =
{
  -1.37825152569120859100E3,
  -3.88016315134637840924E4,
  -3.31612992738871184744E5,
  -1.16237097492762307383E6,
  -1.72173700820839662146E6,
  -8.53555664245765465627E5
};

static DBL C[] =
{
   1.00000000000000000000E0,
  -3.51815701436523470549E2,
  -1.70642106651881159223E4,
  -2.20528590553854454839E5,
  -1.13933444367982507207E6,
  -2.53252307177582951285E6,
  -2.01889141433532773231E6
};

/* log(sqrt(2pi)) */

static DBL LS2PI = 0.91893853320467274178;

/* sqrt(2pi) */

static DBL s2pi = 2.50662827463100050242E0;

/* approximation for 0 <= |y - 0.5| <= 3/8 */

static DBL P0[5] =
{
  -5.99633501014107895267E1,
   9.80010754185999661536E1,
  -5.66762857469070293439E1,
   1.39312609387279679503E1,
  -1.23916583867381258016E0,
};

static DBL Q0[8] =
{
/* 1.00000000000000000000E0,*/
   1.95448858338141759834E0,
   4.67627912898881538453E0,
   8.63602421390890590575E1,
  -2.25462687854119370527E2,
   2.00260212380060660359E2,
  -8.20372256168333339912E1,
   1.59056225126211695515E1,
  -1.18331621121330003142E0,
};

/*
 * Approximation for interval z = sqrt(-2 log y ) between 2 and 8
 * i.e., y between exp(-2) = .135 and exp(-32) = 1.27e-14.
 */

static DBL P1[9] =
{
   4.05544892305962419923E0,
   3.15251094599893866154E1,
   5.71628192246421288162E1,
   4.40805073893200834700E1,
   1.46849561928858024014E1,
   2.18663306850790267539E0,
  -1.40256079171354495875E-1,
  -3.50424626827848203418E-2,
  -8.57456785154685413611E-4,
};

static DBL Q1[8] =
{
/*  1.00000000000000000000E0,*/
   1.57799883256466749731E1,
   4.53907635128879210584E1,
   4.13172038254672030440E1,
   1.50425385692907503408E1,
   2.50464946208309415979E0,
  -1.42182922854787788574E-1,
  -3.80806407691578277194E-2,
  -9.33259480895457427372E-4,
};

/*
 * Approximation for interval z = sqrt(-2 log y ) between 8 and 64
 * i.e., y between exp(-32) = 1.27e-14 and exp(-2048) = 3.67e-890.
 */

static DBL P2[9] =
{
  3.23774891776946035970E0,
  6.91522889068984211695E0,
  3.93881025292474443415E0,
  1.33303460815807542389E0,
  2.01485389549179081538E-1,
  1.23716634817820021358E-2,
  3.01581553508235416007E-4,
  2.65806974686737550832E-6,
  6.23974539184983293730E-9,
};

static DBL Q2[8] =
{
/*  1.00000000000000000000E0,*/
  6.02427039364742014255E0,
  3.67983563856160859403E0,
  1.37702099489081330271E0,
  2.16236993594496635890E-1,
  1.34204006088543189037E-2,
  3.28014464682127739104E-4,
  2.89247864745380683936E-6,
  6.79019408009981274425E-9,
};


/*****************************************************************************
* Static functions
******************************************************************************/

static DBL igami PARAMS((DBL a, DBL y0));
static DBL lgam PARAMS((DBL x));
static DBL polevl PARAMS((DBL x, DBL * coef, int N));
static DBL p1evl PARAMS((DBL x, DBL * coef, int N));
static DBL igamc PARAMS((DBL a, DBL x));
static DBL igam PARAMS((DBL a, DBL x));
static DBL ndtri PARAMS((DBL y0));



/*              chdtri()
 *
 *  Inverse of complemented Chi-square distribution
 *
 *
 *
 * SYNOPSIS:
 *
 * DBL df, x, y, chdtri();
 *
 * x = chdtri( df, y );
 *
 *
 *
 *
 * DESCRIPTION:
 *
 * 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 igami.c.
 *
 * ERROR MESSAGES:
 *
 *   message         condition      value returned
 * chdtri domain   y < 0 or y > 1        0.0
 *                     v < 1
 *
 */

DBL chdtri(df, y)
DBL df, y;
{
  DBL x;

  if ((y < 0.0) || (y > 1.0) || (df < 1.0))
  {
    Error("Illegal values fd=%f and y=%f in chdtri().\n", df, y);

    return (0.0);
  }

  x = igami(0.5 * df, y);

  return (2.0 * x);
}



/*              lgam()
 *
 *  Natural logarithm of gamma function
 *
 *
 *
 * SYNOPSIS:
 *
 * DBL x, y, lgam();
 * extern int sgngam;
 *
 * y = lgam( x );
 *
 *
 *
 * DESCRIPTION:
 *
 * Returns the base e (2.718...) logarithm of the absolute
 * value of the gamma function of the argument.
 * The sign (+1 or -1) of the gamma function is returned in a
 * global (extern) variable named sgngam.
 *
 * For arguments greater than 13, the logarithm of the gamma
 * function is approximated by the logarithmic version of
 * Stirling's formula using a polynomial approximation of
 * degree 4. Arguments between -33 and +33 are reduced by
 * recurrence to the interval [2,3] of a rational approximation.
 * The cosecant reflection formula is employed for arguments
 * less than -33.
 *
 * Arguments greater than MAXLGM return MAXNUM and an error
 * message.  MAXLGM = 2.035093e36 for DEC
 * arithmetic or 2.556348e305 for IEEE arithmetic.
 *
 *
 *
 * ACCURACY:
 *
 *
 * arithmetic      domain        # trials     peak         rms
 *    DEC     0, 3                  7000     5.2e-17     1.3e-17
 *    DEC     2.718, 2.035e36       5000     3.9e-17     9.9e-18
 *    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
 *
 */

static DBL lgam(x)
DBL x;
{
  DBL p, q, w, z;
  int i;

  sgngam = 1;

  if (x < -34.0)
  {
    q = -x;
    w = lgam(q);  /* note this modifies sgngam! */
    p = floor(q);

    if (p == q)
    {
      goto loverf;
    }

    i = p;

    if ((i & 1) == 0)
    {
      sgngam = -1;
    }
    else
    {
      sgngam = 1;
    }

    z = q - p;

    if (z > 0.5)
    {
      p += 1.0;

      z = p - q;
    }

    z = q * sin(M_PI * z);

    if (z == 0.0)
    {
      goto loverf;
    }

/*  z = log(M_PI) - log( z ) - w;*/

    z = LOGPI - log(z) - w;

    return (z);
  }

  if (x < 13.0)
  {
    z = 1.0;

    while (x >= 3.0)
    {
      x -= 1.0;

      z *= x;
    }

    while (x < 2.0)
    {
      if (x == 0.0)
      {
        goto loverf;
      }

      z /= x;

      x += 1.0;
    }

    if (z < 0.0)
    {
      sgngam = -1;

      z = -z;
    }
    else
    {
      sgngam = 1;
    }

    if (x == 2.0)
    {
      return (log(z));
    }

    x -= 2.0;

    p = x * polevl(x, B, 5) / p1evl(x, C, 6);

    return (log(z) + p);
  }

  if (x > MAXLGM)
  {
    loverf:
/*
    mtherr("lgam", OVERFLOW);
*/
    return (sgngam * MAXNUM);
  }

  q = (x - 0.5) * log(x) - x + LS2PI;

  if (x > 1.0e8)
  {
    return (q);
  }

  p = 1.0 / (x * x);

  if (x >= 1000.0)
  {
    q += ((7.9365079365079365079365e-4 * p -
           2.7777777777777777777778e-3) * p +
           0.0833333333333333333333) / x;
  }
  else
  {
    q += polevl(p, A, 4) / x;
  }

  return (q);
}



/*              igamc()
 *
 *  Complemented incomplete gamma integral
 *
 *
 *
 * SYNOPSIS:
 *
 * DBL a, x, y, igamc();
 *
 * y = igamc( a, x );
 *
 *
 *
 * DESCRIPTION:
 *
 * 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:
 *
 *                      Relative error:
 * arithmetic   domain     # trials      peak         rms
 *    DEC       0,30         2000       2.7e-15     4.0e-16
 *    IEEE      0,30        60000       1.4e-12     6.3e-15
 *
 */

static DBL igamc(a, x)
DBL a, x;
{
  DBL ans, c, yc, ax, y, z;
  DBL pk, pkm1, pkm2, qk, qkm1, qkm2;
  DBL r, t;

  if ((x <= 0) || (a <= 0))
  {
    return (1.0);
  }

  if ((x < 1.0) || (x < a))
  {
    return (1.0 - igam(a, x));
  }

  ax = a * log(x) - x - lgam(a);

  if (ax < -MAXLOG)
  {
/*
    mtherr("igamc", UNDERFLOW);
*/
    return (0.0);
  }

  ax = exp(ax);

/* continued fraction */

  y = 1.0 - a;
  z = x + y + 1.0;
  c = 0.0;

  pkm2 = 1.0;
  qkm2 = x;
  pkm1 = x + 1.0;
  qkm1 = z * x;

  ans = pkm1 / qkm1;

  do
  {
    c += 1.0;
    y += 1.0;
    z += 2.0;

    yc = y * c;

    pk = pkm1 * z - pkm2 * yc;
    qk = qkm1 * z - qkm2 * yc;

    if (qk != 0)
    {
      r = pk / qk;
      t = fabs((ans - r) / r);
      ans = r;
    }
    else
    {
      t = 1.0;
    }

    pkm2 = pkm1;
    pkm1 = pk;
    qkm2 = qkm1;
    qkm1 = qk;

    if (fabs(pk) > BIG)
    {
      pkm2 /= BIG;
      pkm1 /= BIG;
      qkm2 /= BIG;
      qkm1 /= BIG;
    }
  }
  while (t > MACHEP);

  return (ans * ax);
}



/*              igam.c
 *
 *  Incomplete gamma integral
 *
 *
 *
 * SYNOPSIS:
 *
 * DBL a, x, y, igam();
 *
 * y = igam( a, x );
 *
 *
 *
 * DESCRIPTION:
 *
 * 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
 *    DEC       0,30         4000       4.4e-15     6.3e-16
 *    IEEE      0,30        10000       3.6e-14     5.1e-15
 *
 */

/* left tail of incomplete gamma function:
 *
 *          inf.      k
 *   a  -x   -       x
 *  x  e     >   ----------
 *           -     -
 *          k=0   | (a+k+1)
 *
 */

static DBL igam(a, x)
DBL a, x;
{
  DBL ans, ax, c, r;

  if ((x <= 0) || (a <= 0))
  {
    return (0.0);
  }

  if ((x > 1.0) && (x > a))
  {
    return (1.0 - igamc(a, x));
  }

/* Compute  x**a * exp(-x) / gamma(a)  */
  ax = a * log(x) - x - lgam(a);

  if (ax < -MAXLOG)
  {
/*
    mtherr("igam", UNDERFLOW);
*/
    return (0.0);
  }

  ax = exp(ax);

/* power series */
  r = a;
  c = 1.0;
  ans = 1.0;

  do
  {
    r += 1.0;
    c *= x / r;
    ans += c;
  }
  while (c / ans > MACHEP);

  return (ans * ax / a);
}



/*              igami()
 *
 *      Inverse of complemented imcomplete gamma integral
 *
 *
 *
 * SYNOPSIS:
 *
 * DBL a, x, y, igami();
 *
 * x = igami( a, y );
 *
 *
 *
 * DESCRIPTION:
 *
 * Given y, the function finds x such that
 *
 *  igamc( a, x ) = y.
 *
 * Starting with the approximate value
 *
 *         3
 *  x = a t
 *
 *  where
 *
 *  t = 1 - d - ndtri(y) sqrt(d)
 *
 * and
 *
 *  d = 1/9a,
 *
 * the routine performs up to 10 Newton iterations to find the
 * root of igamc(a,x) - y = 0.
 *
 *
 * ACCURACY:
 *
 * Tested for a ranging from 0.5 to 30 and x from 0 to 0.5.
 *
 *                      Relative error:
 * arithmetic   domain     # trials      peak         rms
 *    DEC       0,0.5         3400       8.8e-16     1.3e-16
 *    IEEE      0,0.5        10000       1.1e-14     1.0e-15
 *
 */

static DBL igami(a, y0)
DBL a, y0;
{
  DBL d, y, x0, lgm;
  int i;

/* approximation to inverse function */
  d = 1.0 / (9.0 * a);
  y = (1.0 - d - ndtri(y0) * sqrt(d));

  x0 = a * y * y * y;

  lgm = lgam(a);

  for (i = 0; i < 10; i++)
  {
    if (x0 <= 0.0)
    {
/*
      mtherr("igami", UNDERFLOW);
*/
      return (0.0);
    }

    y = igamc(a, x0);

/* compute the derivative of the function at this point */
    d = (a - 1.0) * log(x0) - x0 - lgm;

    if (d < -MAXLOG)
    {
/*
      mtherr("igami", UNDERFLOW);
*/
      goto done;
    }

    d = -exp(d);

/* compute the step to the next approximation of x */
    if (d == 0.0)
    {
      goto done;
    }

    d = (y - y0) / d;

    x0 = x0 - d;

    if (i < 3)
    {
      continue;
    }

    if (fabs(d / x0) < 2.0 * MACHEP)
    {
      goto done;
    }
  }

done:

  return (x0);
}



/*              ndtri.c
 *
 *  Inverse of Normal distribution function
 *
 *
 *
 * SYNOPSIS:
 *
 * DBL x, y, ndtri();
 *
 * x = ndtri( y );
 *
 *
 *
 * DESCRIPTION:
 *
 * 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
 *    DEC      0.125, 1         5500       9.5e-17     2.1e-17
 *    DEC      6e-39, 0.135     3500       5.7e-17     1.3e-17
 *    IEEE     0.125, 1        20000       7.2e-16     1.3e-16
 *    IEEE     3e-308, 0.135   50000       4.6e-16     9.8e-17
 *
 *
 * ERROR MESSAGES:
 *
 *   message         condition    value returned
 * ndtri domain       x <= 0        -MAXNUM
 * ndtri domain       x >= 1         MAXNUM
 *
 */

static DBL ndtri(y0)
DBL y0;
{
  DBL x, y, z, y2, x0, x1;
  int code;

  if (y0 <= 0.0)
  {
/*
    mtherr("ndtri", DOMAIN);
*/
    return (-MAXNUM);
  }

  if (y0 >= 1.0)
  {
/*
    mtherr("ndtri", DOMAIN);
*/
    return (MAXNUM);
  }

  code = 1;

  y = y0;

  if (y > (1.0 - 0.13533528323661269189)) /* 0.135... = exp(-2) */
  {
    y = 1.0 - y;
    code = 0;
  }

  if (y > 0.13533528323661269189)
  {
    y = y - 0.5;
    y2 = y * y;
    x = y + y * (y2 * polevl(y2, P0, 4) / p1evl(y2, Q0, 8));
    x = x * s2pi;

    return (x);
  }

  x = sqrt(-2.0 * log(y));
  x0 = x - log(x) / x;

  z = 1.0 / x;

  if (x < 8.0)  /* y > exp(-32) = 1.2664165549e-14 */
  {
    x1 = z * polevl(z, P1, 8) / p1evl(z, Q1, 8);
  }
  else
  {
    x1 = z * polevl(z, P2, 8) / p1evl(z, Q2, 8);
  }

  x = x0 - x1;

  if (code != 0)
  {
    x = -x;
  }

  return (x);
}



/*              polevl.c
 *              p1evl.c
 *
 *  Evaluate polynomial
 *
 *
 *
 * SYNOPSIS:
 *
 * int N;
 * DBL x, y, coef[N+1], polevl[];
 *
 * y = polevl( x, coef, N );
 *
 *
 *
 * DESCRIPTION:
 *
 * Evaluates polynomial of degree N:
 *
 *                     2          N
 * y  =  C  + C x + C x  +...+ C x
 *        0    1     2          N
 *
 * Coefficients are stored in reverse order:
 *
 * coef[0] = C  , ..., coef[N] = C  .
 *            N                   0
 *
 *  The function p1evl() assumes that coef[N] = 1.0 and is
 * omitted from the array.  Its calling arguments are
 * otherwise the same as polevl().
 *
 *
 * SPEED:
 *
 * In the interest of speed, there are no checks for out
 * of bounds arithmetic.  This routine is used by most of
 * the functions in the library.  Depending on available
 * equipment features, the user may wish to rewrite the
 * program in microcode or assembly language.
 *
 */

static DBL polevl(x, coef, N)
DBL x;
DBL coef[];
int N;
{
  DBL ans;
  int i;
  DBL *p;

  p = coef;
  ans = *p++;
  i = N;

  do
  {
    ans = ans * x + *p++;
  }
  while (--i);

  return (ans);
}

/*              p1evl() */
/*                                          N
 * Evaluate polynomial when coefficient of x  is 1.0.
 * Otherwise same as polevl.
 */

static DBL p1evl(x, coef, N)
DBL x;
DBL coef[];
int N;
{
  DBL ans;
  DBL *p;
  int i;

  p = coef;
  ans = x + *p++;
  i = N - 1;

  do
  {
    ans = ans * x + *p++;
  }
  while (--i);

  return (ans);
}