# include <math.h>
# include <stdio.h>
# include <stdlib.h>

# include "asa005.h"

/******************************************************************************/

double prncst ( double st, int idf, double d, int *ifault )

/******************************************************************************/
/*
  Purpose:

    prncst() computes the lower tail of noncentral T distribution.

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    23 October 2010

  Author:

    Original FORTRAN77 version by BE Cooper.
    C version by John Burkardt.

  Reference:

    BE Cooper,
    Algorithm AS 5:
    The Integral of the Non-Central T-Distribution,
    Applied Statistics,
    Volume 17, Number 2, 1968, page 193.

  Input:

    double ST, the argument.

    int IDF, the number of degrees of freedom.

    double D, the noncentrality parameter.

  Output:

    int *IFAULT, error flag.
    0, no error occurred.
    nonzero, an error occurred.

    double PRNCST, the value of the lower tail of
    the noncentral T distribution.

  Local:

    double G1, 1.0 / sqrt(2.0 * pi)

    double G2, 1.0 / (2.0 * pi)

    double G3, sqrt(2.0 * pi)
*/
{
  double a;
  double ak;
  double b;
  double da;
  double drb;
  double emin = 12.5;
  double f;
  double fk;
  double fkm1;
  double fmkm1;
  double fmkm2;
  double g1 = 0.3989422804;
  double g2 = 0.1591549431;
  double g3 = 2.5066282746;
  int ioe;
  int k;
  double rb;
  double sum;
  double value;

  f = ( double ) ( idf );
/*
  For very large IDF, use the normal approximation.
*/
  if ( 100 < idf )
  {
    *ifault = 1;

    a = sqrt ( 0.5 * f ) * exp ( lgamma ( 0.5 * ( f - 1.0 ) ) 
      - lgamma ( 0.5 * f ) ) * d;

    value = alnorm ( ( st - a ) / sqrt ( f * ( 1.0 + d * d ) 
      / ( f - 2.0 ) - a * a ), 0 );
    return value;
  }

  *ifault = 0;
  ioe = ( idf % 2 );
  a = st / sqrt ( f );
  b = f / ( f + st * st );
  rb = sqrt ( b );
  da = d * a;
  drb = d * rb;

  if ( idf == 1 )
  {
    value = alnorm ( drb, 1 ) + 2.0 * tfn ( drb, a );
    return value;
  }

  sum = 0.0;

  if ( fabs ( drb ) < emin )
  {
    fmkm2 = a * rb * exp ( - 0.5 * drb * drb ) 
    * alnorm ( a * drb, 0 ) * g1;
  }
  else
  {
    fmkm2 = 0.0;
  }

  fmkm1 = b * da * fmkm2;

  if ( fabs ( d ) < emin )
  {
    fmkm1 = fmkm1 + b * a * g2 * exp ( - 0.5 * d * d );
  }

  if ( ioe == 0 )
  {
    sum = fmkm2;
  }
  else
  {
    sum = fmkm1;
  }

  ak = 1.0;
  fk = 2.0;

  for ( k = 2; k <= idf - 2; k = k + 2 )
  {
    fkm1 = fk - 1.0;
    fmkm2 = b * ( da * ak * fmkm1 + fmkm2 ) * fkm1 / fk;
    ak = 1.0 / ( ak * fkm1 );
    fmkm1 = b * ( da * ak * fmkm2 + fmkm1 ) * fk / ( fk + 1.0 );

    if ( ioe == 0 )
    {
      sum = sum + fmkm2;
    }
    else
    {
      sum = sum + fmkm1;
    }
    ak = 1.0 / ( ak * fk );
    fk = fk + 2.0;
  }

  if ( ioe == 0 )
  {
    value = alnorm ( d, 1 ) + sum * g3;
  }
  else
  {
    value = alnorm ( drb, 1 ) + 2.0 * ( sum + tfn ( drb, a ) );
  }

  return value;
}