/*-Algorithm AS 177
 * Expected Normal Order Statistics (Exact and Approximate),
 * by J.P. Royston, 1982.
 * Applied Statistics, 31(2):161-165.
 *
 * Translation to C by James Darrell McCauley, mccauley@ecn.purdue.edu.
 *
 * The functions nscor1() and nscor2() calculate the expected values of
 * normal order statistics in exact or approximate form, respectively.
 *
 */

#define NSTEP 721
#define H 0.025
 
#include<math.h>
#include<stdio.h>

void nscor1  (double s[], int n, int n2, double work[], int *ifault)




/* exact calculation of normal scores */

{
  double ani, c, c1, d, scor, alnfac();
  int i,j;

  *ifault = 3;
  if (n2 != n / 2)
    return;
  *ifault = 1;
  if (n <= 1)
    return;
  *ifault = 0;
  if (n > 2000)
    *ifault = 2;
  /* calculate the natural log of factorial(n) */
  c1 = alnfac (n);
  d = c1 - log ((double) n);
  /* accumulate ordinates for calculation of integral for rankits */
  for (i = 0; i < n2; ++i)
  {
    ani = (double) n - i - 1;
    c = c1 - d;
    for (scor=0.0, j = 0; j < NSTEP; ++j)
      scor += work[0*NSTEP+j] *
	exp (work[1*NSTEP+j] + work[2*NSTEP+j] * i
	     + work[3*NSTEP+j] * ani + c);
    s[i] = scor * H;
    d += log ((double)(i + 1.0) / ani);
  }
  return;
}

void init  (double work[])

{
  double xstart = -9.0, pi2 = -0.918938533, xx, alnorm();
  int i;

  xx = xstart;
  /* set up arrays for calculation of integral */
  for (i = 0; i < NSTEP; ++i)
  {
    work[0*NSTEP+i] = xx;
    work[1*NSTEP+i] = pi2 - xx * xx * 0.5;
    work[2*NSTEP+i] = log (alnorm (xx,1));
    work[3*NSTEP+i] = log (alnorm (xx,0));
    xx = xstart + H * (i + 1.0);
  }
  return;
}

double alnfac  (int j)



/*-Algorithm AS 177.2 Appl. Statist. (1982) Vol.31, No.2
 * Natural logarithm of factorial for non-negative argument
 */

{
  static double r[7] = {0.0, 0.0, 0.69314718056, 1.79175946923,
  3.17805383035, 4.78749174278, 6.57925121101};
  double w, z;

  if (j == 1)
    return (double) 1.0;
  else if (j <= 7)
    return r[j];
  w = (double) j + 1;
  z = 1.0 / (w * w);
    return (w - 0.5) * log (w) - w + 0.918938522305 +
    (((4.0 - 3.0 * z) * z - 14.0) * z + 420.0) / (5040.0 * w);
}

void nscor2  (double s[], int n, int n2, int *ifault)




/*-Algorithm AS 177.3 Appl. Statist. (1982) Vol.31, No.2
 * Approximation for Rankits
 */

{
  static double eps[4] = {0.419885, 0.450536, 0.456936, 0.468488};
  static double dl1[4] = {0.112063, 0.121770, 0.239299, 0.215159};
  static double dl2[4] = {0.080122, 0.111348, -0.211867, -0.115049};
  static double gam[4] = {0.474798, 0.469051, 0.208597, 0.259784};
  static double lam[4] = {0.282765, 0.304856, 0.407708, 0.414093};
  static double bb = -0.283833, d = -0.106136, b1 = 0.5641896;
  double e1, e2, l1, correc(), ppnd16();
  int i, k;

  *ifault = 3;
  if (n2 != n / 2)
    return;
  *ifault = 1;
  if (n <= 1)
    return;
  *ifault = 0;
  if (n > 2000)
    *ifault = 2;
  s[0] = b1;
  if (n == 2)
    return;
  /* calculate normal areas for 3 largest rankits */
  k=(n2<3) ? n2 : 3;
  for (i = 0; i < k; ++i)
  {
    e1 = (1.0 + i - eps[i]) / (n + gam[i]);
    e2 = pow (e1, lam[i]);
    s[i] = e1 + e2 * (dl1[i] + e2 * dl2[i]) / n - correc (1 + i, n);
  }
  if (n2 != k)
  {
    /* calculate normal areas for remaining rankits */
    for (i = 3; i < n2; ++i)
    {
      l1 = lam[3] + bb / (1.0 + i + d);
      e1 = (1.0 + i - eps[3]) / (n + gam[3]);
      e2 = pow (e1, l1);
      s[i] = e1 + e2 * (dl1[3] + e2 * dl2[3]) / n - correc (1+i, n);
    }
  }
  /* convert normal tail areas to normal deviates */
  for (i = 0; i < n2; ++i)
    s[i] = -ppnd16 (s[i]);
  return;
}

double correc  (int i, int n)



/*-Algorithm AS 177.4 Appl. Statist. (1982) Vol.31, No.2
 * Calculates correction for tail area of noraml distribution
 * corresponding to ith largest rankit in sample size n.
 */

{
  static double c1[7] = {9.5, 28.7, 1.9, 0.0, -7.0, -6.2, -1.6};
  static double c2[7] = {-6.195e3, -9.569e3, -6.728e3, -17.614e3,
  -8.278e3, -3.570e3, 1.075e3};
  static double c3[7] = {9.338e4, 1.7516e5, 4.1040e5, 2.157e6,
  2.376e6, 2.065e6, 2.065e6};
  static double mic = 1.0e-6, c14 = 1.9e-5;
  double an, ret_val;

  ret_val = c14;
  if (i * n == 4)
    return ret_val;
  ret_val = 0.0;
  if (i < 1 || i > 7)
    return ret_val;
  else if (i != 4 && n > 20)
    return ret_val;
  else if (i == 4 && n > 40)
    return ret_val;
  else
  {
    an = 1.0 / (double) (n * n);
    ret_val = (c1[i-1] + an * (c2[i-1] + an * c3[i-1])) * mic;
    return ret_val;
  }
}