/*

 *  Mathlib : A C Library of Special Functions

 *  Copyright (C) 1998-2015 The R Core Team

 *  based on AS243 (C) 1989 Royal Statistical Society

 *

 *  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; 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.

 *

 *  You should have received a copy of the GNU General Public License

 *  along with this program; if not, a copy is available at

 *  https://www.R-project.org/Licenses/

 */

/*  Algorithm AS 243  Lenth,R.V. (1989). Appl. Statist., Vol.38, 185-189.

 *  ----------------

 *  Cumulative probability at t of the non-central t-distribution

 *  with df degrees of freedom (may be fractional) and non-centrality

 *  parameter delta.

 *

 *  NOTE

 *

 *    Requires the following auxiliary routines:

 *

 *  lgammafn(x) - log gamma function

 *  pbeta(x, a, b)  - incomplete beta function

 *  pnorm(x)  - normal distribution function

 *

 *  CONSTANTS

 *

 *    M_SQRT_2dPI  = 1/ {gamma(1.5) * sqrt(2)} = sqrt(2 / pi)

 *    M_LN_SQRT_PI = ln(sqrt(pi)) = ln(pi)/2

 */

#include "nmath.h"

#include "dpq.h"

/*----------- DEBUGGING -------------

 *

 *  make CFLAGS='-DDEBUG_pnt -g'

 * -- Feb.3, 1999; M.Maechler:

  - For 't > ncp > 20' (or so)  the result is completely WRONG!  <== no longer true

  - but for ncp > 100

 */

double pnt(double t, double df, double ncp, int lower_tail, int log_p)

{

    double albeta, a, b, del, errbd, lambda, rxb, tt, x;

    LDOUBLE geven, godd, p, q, s, tnc, xeven, xodd;

    int it, negdel;

    /* note - itrmax and errmax may be changed to suit one's needs. */

    const int itrmax = 1000;

    const static double errmax = 1.e-12;

    if (df <= 0.0) ML_WARN_return_NAN;

    if(ncp == 0.0) return pt(t, df, lower_tail, log_p);

    if(!R_FINITE(t))

  return (t < 0) ? R_DT_0 : R_DT_1;

    if (t >= 0.) {

  negdel = FALSE; tt = t;  del = ncp;

    }

    else {

  /* We deal quickly with left tail if extreme,

     since pt(q, df, ncp) <= pt(0, df, ncp) = \Phi(-ncp) */

  if (ncp > 40 && (!log_p || !lower_tail)) return R_DT_0;

  negdel = TRUE;  tt = -t; del = -ncp;

    }

    if (df > 4e5 || del*del > 2*M_LN2*(-(DBL_MIN_EXP))) {

  /*-- 2nd part: if del > 37.62, then p=0 below

    FIXME: test should depend on `df', `tt' AND `del' ! */

  /* Approx. from  Abramowitz & Stegun 26.7.10 (p.949) */

  s = 1./(4.*df);

  return pnorm((double)(tt*(1. - s)), del,

         sqrt((double) (1. + tt*tt*2.*s)),

         lower_tail != negdel, log_p);

    }

    /* initialize twin series */

    /* Guenther, J. (1978). Statist. Computn. Simuln. vol.6, 199. */

    x = t * t;

    rxb = df/(x + df);/* := (1 - x) {x below} -- but more accurately */

    x = x / (x + df);/* in [0,1) */

#ifdef DEBUG_pnt

    REprintf("pnt(t=%7g, df=%7g, ncp=%7g) ==> x= %10g:",t,df,ncp, x);

#endif

    if (x > 0.) {/* <==>  t != 0 */

  lambda = del * del;

  p = .5 * exp(-.5 * lambda);

#ifdef DEBUG_pnt

  REprintf("\t p=%10Lg\n",p);

#endif

  if(p == 0.) { /* underflow! */

      /*========== really use an other algorithm for this case !!! */

      ML_WARNING(ME_UNDERFLOW, "pnt");

      ML_WARNING(ME_RANGE, "pnt"); /* |ncp| too large */

      return R_DT_0;

  }

#ifdef DEBUG_pnt

  REprintf("it  1e5*(godd,   geven)|          p           q           s"

         /* 1.3 1..4..7.9 1..4..7.9|1..4..7.901 1..4..7.901 1..4..7.901 */

     "        pnt(*)     errbd\n");

         /* 1..4..7..0..34 1..4..7.9*/

#endif

  q = M_SQRT_2dPI * p * del;

  s = .5 - p;

  /* s = 0.5 - p = 0.5*(1 - exp(-.5 L)) =  -0.5*expm1(-.5 L)) */

  if(s < 1e-7)

      s = -0.5 * expm1(-0.5 * lambda);

  a = .5;

  b = .5 * df;

  /* rxb = (1 - x) ^ b   [ ~= 1 - b*x for tiny x --> see 'xeven' below]

   *       where '(1 - x)' =: rxb {accurately!} above */

  rxb = pow(rxb, b);

  albeta = M_LN_SQRT_PI + lgammafn(b) - lgammafn(.5 + b);

  xodd = pbeta(x, a, b, /*lower*/TRUE, /*log_p*/FALSE);

  godd = 2. * rxb * exp(a * log(x) - albeta);

  tnc = b * x;

  xeven = (tnc < DBL_EPSILON) ? tnc : 1. - rxb;

  geven = tnc * rxb;

  tnc = p * xodd + q * xeven;

  /* repeat until convergence or iteration limit */

  for(it = 1; it <= itrmax; it++) {

      a += 1.;

      xodd  -= godd;

      xeven -= geven;

      godd  *= x * (a + b - 1.) / a;

      geven *= x * (a + b - .5) / (a + .5);

      p *= lambda / (2 * it);

      q *= lambda / (2 * it + 1);

      tnc += p * xodd + q * xeven;

      s -= p;

      /* R 2.4.0 added test for rounding error here. */

      if(s < -1.e-10) { /* happens e.g. for (t,df,ncp)=(40,10,38.5), after 799 it.*/

    ML_WARNING(ME_PRECISION, "pnt");

#ifdef DEBUG_pnt

    REprintf("s = %#14.7Lg < 0 !!! ---> non-convergence!!\n", s);

#endif

    goto finis;

      }

      if(s <= 0 && it > 1) goto finis;

      errbd = (double)(2. * s * (xodd - godd));

#ifdef DEBUG_pnt

      REprintf("%3d %#9.4g %#9.4g|%#11.4Lg %#11.4Lg %#11.4Lg %#14.10Lg %#9.4g\n",

         it, 1e5*(double)godd, 1e5*(double)geven, p, q, s, tnc, errbd);

#endif

      if(fabs(errbd) < errmax) goto finis;/*convergence*/

  }

  /* non-convergence:*/

  ML_WARNING(ME_NOCONV, "pnt");

    }

    else { /* x = t = 0 */

  tnc = 0.;

    }

 finis:

    tnc += pnorm(- del, 0., 1., /*lower*/TRUE, /*log_p*/FALSE);

    lower_tail = lower_tail != negdel; /* xor */

    if(tnc > 1 - 1e-10 && lower_tail)

  ML_WARNING(ME_PRECISION, "pnt{final}");

    return R_DT_val(fmin2((double)tnc, 1.) /* Precaution */);

}