#include <iostream>
#include <stdlib.h>
#include <string.h>
#include <math.h>

int myatoi(char *s, int *i) {
  float wert;
  int   rc=0;
  char  *endp;

  wert=strtod(s, &endp);
  if (s==endp) {
    rc=0;
  }
  else {
    *i=(int)floor(wert);
    rc=1;
  }
  return(rc);
}

int myatof(char *s, float *f) {
  double wert;
  int    rc=0;
  char   *endp;

  wert=strtod(s, &endp);
  if (s==endp) {
    rc=0;
  }
  else {
    *f=(float)wert;
    rc=1;
  }
  return(rc);
}

char *getsig(char *s) {
  char *a, *b, *rc;
  /* allocates memory for a string in */
  /* doublequotes and returns a pointer to it */
  if (strchr(s, '"')) {
    a=strchr(s, '"');
    b=strrchr(s, '"');
    if (a<b) {
      rc=new char[b-a];
      strncpy(rc, (a+1), b-a-1);
      rc[b-a-1]='\0';
    }
    else {
      rc=new char[1];
      strcpy(rc,"\0");
    }
  }
  else {
    rc=new char[1];
    strcpy(rc,"\0");
  }
  return rc;
}

float prho(int n, float is, int *ifault) {
/*
    Algorithm AS 89   Appl. Statist. (1975) Vol.24, No. 3, P377.
       
    To evaluate the probability of obtaining a value greater than or
    equal to is, where is=(n**3-n)*(1-r)/6, r=Spearman's rho and n
    must be greater than 1

    n is the number of data pairs in the correlation
    The calculated probability is for one-sided tests.
    For two-sided tests multiply the probability value with 2.

    Auxiliary function required: ALNORM = algorithm AS66
*/

  float  c1=0.2274, c2=0.2531, c3=0.1745, c4=0.0758, c5=0.1033, c6=0.3932;
  float  c7=0.0879, c8=0.0151, c9=0.0072, c10=0.0831, c11=0.0131, c12=0.00046;
  float  rc;
  double b, x, y, u;
  long   i, js, ifr, nfac, m, n1, mt, nn, ise;
  long   l[7];

  extern float alnorm(float x, int upper);

  /* Test admissibility of arguments and initialize */
  rc=1.0;
  *ifault=1;
  if (n<=1) return(rc);
  *ifault=0;
  if (is<=0) return(rc);
  rc=0.0;
  if (is > n*(n*n-1)/3) return(rc);
  js=(long)is;
  if (js!= 2*(js/2)) js=js+1;
  /*  if (n>6) goto mark6; */
  goto mark6;

  /* Exact evaluation of probability */
  nfac=1;
  for (i=1; i<=n; i++) {
    nfac*=i;
    l[i]=i;
  }
mark1:
  rc=1.0/(float)(nfac);
  if (js != n*(n*n-1)/3) return(rc);
  ifr=0;
  for (m=1; m<=nfac; m++) {
    ise=0;
    for (i=1; i<=n; i++) {
      ise+=(i-l[i])*(i-l[i]);
    }
  mark2:
    if (js <= ise) ifr++;
    n1=n;
  mark3:
    mt=l[1];
    nn=n1-1;
    for (i=1; i<=nn; i++) {
      l[i]=l[i+1];
    }
  mark4:
    l[n1]=mt;
    if ((l[n1]!=n1) || (n1==2)) goto mark5;
    n1--;
    if (m!=nfac) goto mark3;
  }
mark5:
  rc=(float)(ifr)/(float)(nfac);
  return(rc);
  /* Evaluation by Edgeworth series expansion */
mark6:
  b=1.0/(double)(n);
  x=(6.0*((double)(js)-1.0)*b/(1.0/(b*b)-1.0)-1.0)*sqrt(1.0/b-1.0);
  y=x*x;
  u=x*b*(c1+b*(c2+c3*b)+y*(-c4+b*(c5+c6*b)-y*b*(c7+c8*b-y*(c9-c10*b+y*b*(c11-c12*y)))));
  /* Call to algorithm AS 66 */
  rc=u/exp(y/2.0)+alnorm(x,1);
  if (rc < 0.0) rc=0.0;
  if (rc > 1.0) rc=1.0;
  return(rc);
}


/*
      Algorithm AS66 Applied Statistics (1973) vol22 no.3
      Evaluates the tail area of the standardised normal curve
      from x to infinity if upper is .true. or
      from minus infinity to x if upper is .false.
*/
float stdnorm(float x) {
  return (1.0/sqrt(2*3.141592653))*exp(-0.5*x*x);
}
float alnorm(float x, int upper) {
  float q;

  float stdnorm(float x);
  extern float qtrap(float (*func)(float x), float a, float b);
  
  if (x==0.0) q=0;
  else q=qtrap(&stdnorm, 0, x);

  if (upper) q=1.0-q-0.5;
  else q+=0.5;
 
  return(q);
}


#define EPS 1.0e-5
#define JMAX 20
float qtrap(float (*func)(float), float a, float b)
{
  float trapzd(float (*func)(float), float a, float b, int n);
  void nrerror(char error_text[]);
  int j;
  float s,olds;

  olds = -1.0e30;
  for (j=1;j<=JMAX;j++) {
    s=trapzd(func,a,b,j);
    if (fabs(s-olds) < EPS*fabs(olds)) return s;
    olds=s;
  }
  nrerror("Too many steps in routine qtrap");
  return 0.0;
}
#undef EPS
#undef JMAX

#define FUNC(x) ((*func)(x))
float trapzd(float (*func)(float), float a, float b, int n)
{
  float x,tnm,sum,del;
  static float s;
  int it,j;
  
  if (n == 1) {
    return (s=0.5*(b-a)*(FUNC(a)+FUNC(b)));
  } else {
    for (it=1,j=1;j<n-1;j++) it <<= 1;
    tnm=it;
    del=(b-a)/tnm;
    x=a+0.5*del;
    for (sum=0.0,j=1;j<=it;j++,x+=del) sum += FUNC(x);
    s=0.5*(s+(b-a)*sum/tnm);
    return s;
  }
}
#undef FUNC