diff rDiff/src/locfit/Source/libtube.c @ 0:0f80a5141704

version 0.3 uploaded

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/rDiff/src/locfit/Source/libtube.c	Thu Feb 14 23:38:36 2013 -0500
@@ -0,0 +1,708 @@
+/*
+ * Copyright 1996-2006 Catherine Loader.
+ */
+
+#include "mex.h"
+/*
+ * Copyright 1996-2006 Catherine Loader.
+ */
+/*
+ *   Copyright (c) 1998-2006 Catherine Loader.
+ *   This file contains functions to compute the constants
+ *   appearing in the tube formula.
+ */
+
+#include <stdio.h>
+#include <math.h>
+#include "tube.h"
+
+static double *fd, *ft;
+static int globm, (*wdf)(), use_covar, kap_terms;
+
+int k0_reqd(d,n,uc)
+int d, n, uc;
+{ int m;
+  m = d*(d+1)+1;
+  if (uc) return(2*m*m);
+     else return(2*n*m);
+}
+
+void assignk0(z,d,n) /* z should be n*(2*d*d+2*d+2); */
+double *z;
+int d, n;
+{ ft = z; z += n*(d*(d+1)+1);
+  fd = z; z += n*(d*(d+1)+1);
+}
+
+/* Residual projection of y to the columns of A,
+ * (I - A(R^TR)^{-1}A^T)y
+ * R should be from the QR-decomp. of A.
+ */
+void rproject(y,A,R,n,p)
+double *y, *A, *R;
+int n, p;
+{ double v[1+TUBE_MXDIM];
+  int i, j;
+
+  for (i=0; i<p; i++) v[i] = innerprod(&A[i*n],y,n);
+  qrsolv(R,v,n,p);
+  for (i=0; i<n; i++)
+    for (j=0; j<p; j++)
+      y[i] -= A[j*n+i]*v[j];
+}
+
+double k2c(lij,A,m,dd,d)
+double *lij, *A;
+int m, d, dd;
+{ int i, j, k, l;
+  double sum, *bk, v[TUBE_MXDIM];
+
+  for (i=0; i<dd*d; i++)
+    chol_hsolve(fd,&lij[i*m],m,dd+1);
+  for (i=0; i<dd*d; i++)
+    for (j=0; j<dd*d; j++)
+      lij[i*m+j+d+1] -= innerprod(&lij[i*m],&lij[j*m],dd+1);
+
+  sum = 0;
+  for (i=0; i<dd; i++)
+    for (j=0; j<i; j++)
+    { bk = &lij[i*d*m + j*d + d+1];
+      for (k=0; k<dd; k++)
+      { v[0] = 0;
+        for (l=0; l<dd; l++) v[l+1] = bk[k*m+l];
+        chol_solve(fd,v,m,dd+1);
+        for (l=0; l<dd; l++) bk[k*m+l] = v[l+1];
+      }
+      for (k=0; k<dd; k++)
+      { v[0] = 0;
+        for (l=0; l<dd; l++) v[l+1] = bk[l*m+k];
+        chol_solve(fd,v,m,dd+1);
+        for (l=0; l<dd; l++) bk[l*m+k] = v[l+1];
+      }
+      sum += bk[i*m+j] - bk[j*m+i];
+    }
+  return(sum*fd[0]*fd[0]);
+}
+
+double k2x(lij,A,m,d,dd)
+double *lij, *A;
+int m, d, dd;
+{ int i, j, k;
+  double s, v[1+TUBE_MXDIM], *ll;
+
+/* residual projections of lij onto A = [l,l1,...,ld] */
+  for (i=0; i<d; i++)
+    for (j=i; j<d; j++)
+    { ll = &lij[(i*dd+j)*m];
+      rproject(ll,A,fd,m,d+1);
+      if (i!=j) memcpy(&lij[(j*dd+i)*m],ll,m*sizeof(double));
+    }
+
+/* compute lij[j][i] = e_i^T (A^T A)^{-1} B_j^T */
+  for (k=0; k<m; k++)
+    for (j=0; j<d; j++)
+    { v[0] = 0;
+      for (i=0; i<d; i++) v[i+1] = lij[(j*dd+i)*m+k];
+      qrsolv(fd,v,m,d+1);
+      for (i=0; i<d; i++) lij[(j*dd+i)*m+k] = v[i+1];
+    } 
+
+/* finally, add up to get the kappa2 term */
+  s = 0;
+  for (j=0; j<d; j++)
+    for (k=0; k<j; k++)
+      s += innerprod(&lij[(j*dd+j)*m],&lij[(k*dd+k)*m],m)
+         - innerprod(&lij[(j*dd+k)*m],&lij[(k*dd+j)*m],m);
+
+  return(s*fd[0]*fd[0]);
+}
+
+void d2c(ll,nn,li,ni,lij,nij,M,m,dd,d)
+double *ll, *nn, *li, *ni, *lij, *nij, *M;
+int m, dd, d;
+{ int i, j, k, l, t, u, v, w;
+  double z;
+
+  for (i=0; i<dd; i++)
+    for (j=0; j<dd; j++)
+    { for (k=0; k<d; k++)
+      { for (l=0; l<d; l++)
+        { z = M[i*d+k]*M[j*d+l];
+          if (z != 0.0)
+          { nij[(i*d+j)*m] += z*lij[(k*d+l)*m];
+            for (t=0; t<d; t++) /* need d, not dd here */
+              for (u=0; u<d; u++)
+                nij[(i*d+j)*m+t+1] += z*M[t*d+u]*lij[(k*d+l)*m+u+1];
+            for (t=0; t<dd; t++)
+              for (u=0; u<dd; u++)
+              { for (v=0; v<d; v++)
+                  for (w=0; w<d; w++)
+                    nij[(i*d+j)*m+(t*d+u)+d+1] +=
+                      z*M[t*d+v]*M[u*d+w]*lij[(k*d+l)*m+(v*d+w)+d+1];
+                for (v=0; v<d; v++)
+                  nij[(i*d+j)*m+(t*d+u)+d+1] += z*M[(v+1)*d*d+t*d+u]*lij[(k*d+l)*m+v+1];
+              }
+          }
+        }
+
+        z = M[(k+1)*d*d+i*d+j];
+        if (z!=0.0)
+        { nij[(i*d+j)*m] += z*li[k*m];
+          for (t=0; t<d; t++)
+            for (u=0; u<d; u++)
+              nij[(i*d+j)*m+t+1] += z*M[t*d+u]*li[k*m+u+1];
+          for (t=0; t<dd; t++)
+            for (u=0; u<dd; u++)
+            { for (v=0; v<d; v++)
+                for (w=0; w<d; w++)
+                  nij[(i*d+j)*m+(t*d+u)+d+1] += z*M[t*d+v]*M[u*d+w]*lij[(v*d+w)*m+k+1];
+              for (v=0; v<d; v++)
+                nij[(i*d+j)*m+(t*d+u)+d+1] += z*M[(v+1)*d*d+t*d+u]*li[k*m+v+1];
+            }
+        }
+      }
+    }
+}
+
+void d2x(li,lij,nij,M,m,dd,d)
+double *li, *lij, *nij, *M;
+int m, dd, d;
+{ int i, j, k, l, z;
+  double t;
+  for (i=0; i<dd; i++)
+    for (j=0; j<dd; j++)
+    { for (k=0; k<d; k++)
+      { for (l=0; l<d; l++)
+        { t = M[i*d+k] * M[j*d+l];
+          if (t != 0.0)
+          { for (z=0; z<m; z++)
+              nij[(i*d+j)*m+z] += t*lij[(k*d+l)*m+z];
+          }
+        }
+        t = M[(k+1)*d*d+i*d+j];
+        if (t!=0.0)
+          for (z=0; z<m; z++)
+            nij[(i*d+j)*m+z] += t*li[k*m+z];
+      }
+    }
+}
+
+int k0x(x,d,kap,M)
+double *x, *kap, *M;
+int d;
+{ double det, *lij, *nij, z;
+  int j, m, r;
+
+  r = 1 + ((d>=2) & (kap_terms >= 3));
+  m = globm = wdf(x,ft,r);
+
+  memcpy(fd,ft,m*(d+1)*sizeof(double));
+  if (use_covar) chol_dec(fd,m,d+1);
+            else qr(fd,m,d+1,NULL);
+
+  det = 1;
+  for (j=1; j<=d; j++)
+    det *= fd[j*(m+1)]/fd[0];
+  kap[0] = det;
+  if (kap_terms == 1) return(1);
+  kap[1] = 0.0;
+  if ((kap_terms == 2) | (d<=1)) return(2);
+
+  lij = &ft[(d+1)*m];
+  nij = &fd[(d+1)*m];
+  memcpy(nij,lij,m*d*d*sizeof(double));
+  z = (use_covar) ?  k2c(nij,ft,m,d,d) : k2x(nij,ft,m,d,d);
+  kap[2] = z*det;
+  if ((kap_terms == 3) | (d==2)) return(3);
+
+  kap[3] = 0;
+  return(4);
+}
+
+void d1c(li,ni,m,d,M)
+double *li, *ni, *M;
+int m, d;
+{ int i, j, k, l;
+  double t;
+
+  fd[0] = ft[0];
+  for (i=0; i<d; i++)
+  { t = 0;
+    for (j=0; j<d; j++) t += M[i*d+j]*li[j*m];
+    fd[i+1] = ni[i*m] = t;
+
+    for (j=0; j<d; j++)
+    { t = 0;
+      for (k=0; k<d; k++)
+        for (l=0; l<d; l++)
+          t += li[k*m+l+1] * M[i*d+k] * M[j*d+l];
+      ni[i*m+j+1] = t;
+    }
+  }
+}
+
+void d1x(li,ni,m,d,M)
+double *li, *ni, *M;
+int m, d;
+{ int i, j, k;
+  memcpy(fd,ft,m*sizeof(double));
+  setzero(ni,m*d);
+  for (j=0; j<d; j++)
+    for (k=0; k<d; k++)
+      if (M[j*d+k]!=0)
+        for (i=0; i<m; i++) ni[j*m+i] += M[j*d+k]*li[k*m+i];
+}
+
+int l1x(x,d,lap,M)
+double *x, *lap, *M;
+int d;
+{ double det, sumcj, *u, v[TUBE_MXDIM];
+  double *ll, *li, *lij, *ni, *nij;
+  int i, j, m;
+  if (kap_terms<=1) return(0);
+  m = globm;
+  li  = &ft[m]; lij = &ft[(d+1)*m];
+  ni  = &fd[m]; nij = &fd[(d+1)*m];
+  setzero(ni,m*d);
+  setzero(nij,m*d*d);
+
+  if (use_covar) d1c(li,ni,m,d,M);
+            else d1x(li,ni,m,d,M);
+
+/* the last (d+1) columns of nij are free, use for an extra copy of ni */
+  ll = &fd[d*d*m];
+  u = &ll[d*m];
+  if (use_covar)
+    memcpy(u,&ni[(d-1)*m],d*sizeof(double)); /* cov(ld, (l,l1,...ld-1)) */
+  else
+    memcpy(ll,fd,(d+1)*m*sizeof(double));
+
+  if (use_covar) chol_dec(fd,m,d+1);
+            else qr(fd,m,d+1,NULL);
+  det = 1;
+  for (j=1; j<d; j++)
+    det *= fd[(m+1)*j]/fd[0];
+  lap[0] = det;
+  if ((kap_terms==2) | (d<=1)) return(1);
+
+  sumcj = 0.0;
+  if (use_covar)
+  { d2c(ft,fd,li,ni,lij,nij,M,m,d-1,d);
+    chol_solve(fd,u,m,d);
+    for (i=0; i<d-1; i++)
+    { v[0] = 0;
+      for (j=0; j<d-1; j++)
+        v[j+1] = nij[(i*d+j)*m+d] - innerprod(u,&nij[(i*d+j)*m],d);
+      chol_solve(fd,v,m,d);
+      sumcj -= v[i+1];
+    }
+  }
+  else
+  { d2x(li,lij,nij,M,m,d-1,d);
+    rproject(u,ll,fd,m,d);
+    for (i=0; i<d-1; i++)
+    { v[0] = 0;
+      for (j=0; j<d-1; j++) v[j+1] = innerprod(&nij[(i*d+j)*m],u,m);
+      qrsolv(fd,v,m,d);
+      sumcj -= v[i+1];
+    }
+  }
+
+  lap[1] = sumcj*det*fd[0]/fd[(m+1)*d];
+  if ((kap_terms==3) | (d==2)) return(2);
+
+  if (use_covar) lap[2] = k2c(nij,ll,m,d-1,d)*det;
+            else lap[2] = k2x(nij,ll,m,d-1,d)*det;
+  return(3);
+}
+
+int m0x(x,d,m0,M)
+double *x, *m0, *M;
+int d;
+{ double det, *li, *ni, *lij, *nij, *ll, *u1, *u2;
+  double om, so, co, sumcj, v[TUBE_MXDIM];
+  int m, i, j;
+  
+  if ((kap_terms<=2) | (d<=1)) return(0);
+
+  m = globm;
+  li  = &ft[m]; lij = &ft[(d+1)*m];
+  ni  = &fd[m]; nij = &fd[(d+1)*m];
+  setzero(ni,m*d); 
+  setzero(nij,m*d*d);
+
+  if (use_covar) d1c(li,ni,m,d,M);
+            else d1x(li,ni,m,d,M);
+
+/* the last (d+1) columns of nij are free, use for an extra copy of ni */
+  ll = &fd[d*d*m];
+  u1 = &ll[d*m];
+  u2 = &ll[(d-1)*m];
+  if (use_covar)
+  { memcpy(u1,&ni[(d-1)*m],d*sizeof(double));
+    memcpy(u2,&ni[(d-2)*m],d*sizeof(double));
+  }
+  else
+    memcpy(ll,fd,(d+1)*m*sizeof(double));
+
+  if (use_covar) chol_dec(fd,m,d+1);
+            else qr(fd,m,d+1,NULL);
+  det = 1;
+  for (j=1; j<d-1; j++)
+    det *= fd[j*(m+1)]/fd[0];
+  om = atan2(fd[d*(m+1)],-fd[d*(m+1)-1]);
+  m0[0] = det*om;
+  if ((kap_terms==3) | (d==2)) return(1);
+
+  so = sin(om)/fd[d*(m+1)];
+  co = (1-cos(om))/fd[(d-1)*(m+1)];
+  sumcj = 0.0;
+  if (use_covar)
+  { d2c(ft,fd,li,ni,lij,nij,M,m,d-2,d);
+    chol_solve(fd,u1,m,d);
+    chol_solve(fd,u2,m,d-1);
+    for (i=0; i<d-2; i++)
+    { v[0] = 0;
+      for (j=0; j<d-2; j++)
+        v[j+1] =
+          so*(nij[(i*d+j)*m+d]-innerprod(u1,&nij[(i*d+j)*m],d))
+         +co*(nij[(i*d+j)*m+d-1]-innerprod(u2,&nij[(i*d+j)*m],d-1));
+      qrsolv(fd,v,m,d-1);
+      sumcj -= v[i+1];
+    }
+  }
+  else
+  { d2x(li,lij,nij,M,m,d-2,d);
+    rproject(u1,ll,fd,m,d);
+    rproject(u2,ll,fd,m,d-1); /* now, u1, u2 are unnormalized n1*, n2* */
+    for (i=0; i<m; i++)
+      u1[i] = so*u1[i] + co*u2[i];  /* for n1*, n2* */
+    for (i=0; i<d-2; i++)
+    { v[0] = 0;
+      for (j=0; j<d-2; j++)
+        v[j+1] = innerprod(&nij[(i*d+j)*m],u1,m);
+      qrsolv(fd,v,m,d-1);
+      sumcj -= v[i+1];
+    }
+  }
+
+  m0[1] = sumcj*det*fd[0];
+  return(2);
+}
+
+int n0x(x,d,n0,M)
+double *x, *n0, *M;
+int d;
+{ double det, *li, *ni, *a0, *a1, *a2;
+  int j, m;
+
+  if ((kap_terms <= 3) | (d <= 2)) return(0);
+
+  m = globm;
+  li  = &ft[m];
+  ni  = &fd[m];
+
+  if (use_covar) d1c(li,ni,m,d,M);
+            else d1x(li,ni,m,d,M);
+
+  det = 1;
+  if (use_covar) chol_dec(fd,m,d+1);
+            else qr(fd,m,d+1,NULL);
+  for (j=1; j<d-2; j++)
+    det *= fd[j*(m+1)]/fd[0];
+
+  a0 = &ni[(d-3)*m+d-2];
+  a1 = &ni[(d-2)*m+d-2];
+  a2 = &ni[(d-1)*m+d-2];
+
+  a0[0] = a1[1]*a2[2];
+  a0[1] =-a1[0]*a2[2];
+  a0[2] = a1[0]*a2[1]-a1[1]*a2[0];
+  a1[0] = 0;
+  a1[1] = a2[2];
+  a1[2] =-a2[1];
+  a2[0] = a2[1] = 0.0; a2[2] = 1.0;
+  rn3(a0); rn3(a1);
+  n0[0] = det*sptarea(a0,a1,a2);
+  return(1);
+}
+
+int kodf(ll,ur,mg,kap,lap)
+double *ll, *ur, *kap, *lap;
+int *mg;
+{ double x[1], *l0, *l1, t, sum;
+  int i, j, n;
+
+  sum = 0.0;
+  for (i=0; i<=mg[0]; i++)
+  { if (i&1) { l1 = fd; l0 = ft; }
+        else { l1 = ft; l0 = fd; }
+    x[0] = ll[0] + (ur[0]-ll[0])*i/mg[0];
+    n = wdf(x,l0,0);
+
+    t = sqrt(innerprod(l0,l0,n));
+    for (j=0; j<n; j++) l0[j] /= t;
+
+    if (i>0)
+    { t = 0.0;
+      for (j=0; j<n; j++) t += (l1[j]-l0[j])*(l1[j]-l0[j]);
+      sum += 2*asin(sqrt(t)/2);
+    }
+  }
+  kap[0] = sum;
+  if (kap_terms<=1) return(1);
+  kap[1] = 0.0;
+  lap[0] = 2.0;
+  return(2);
+}
+
+int tube_constants(f,d,m,ev,mg,fl,kap,wk,terms,uc)
+double *fl, *kap, *wk;
+int d, m, ev, *mg, (*f)(), terms, uc;
+{ int aw, deb=0;
+  double k0[4], l0[3], m0[2], n0[1], z[TUBE_MXDIM];
+
+  wdf = f;
+
+  aw = (wk==NULL);
+  if (aw) {
+    wk = (double *)calloc(k0_reqd(d,m,uc),sizeof(double));
+    if ( wk == NULL ) {
+      printf("Problem allocating memory for wk\n");fflush(stdout);
+    }
+  }
+  assignk0(wk,d,m);
+
+  k0[0] = k0[1] = k0[2] = k0[3] = 0.0;
+  l0[0] = l0[1] = l0[2] = 0.0;
+  m0[0] = m0[1] = 0.0;
+  n0[0] = 0.0;
+
+  use_covar = uc;
+  kap_terms = terms;
+  if ((kap_terms <=0) | (kap_terms >= 5))
+    mut_printf("Warning: terms = %2d\n",kap_terms);
+
+  switch(ev)
+  {
+    case IMONTE:
+      monte(k0x,fl,&fl[d],d,k0,mg[0]);
+      break;
+    case ISPHERIC:
+      if (d==2) integ_disc(k0x,l1x,fl,k0,l0,mg);
+      if (d==3) integ_sphere(k0x,l1x,fl,k0,l0,mg);
+      break;
+    case ISIMPSON:
+      if (use_covar) simpson4(k0x,l1x,m0x,n0x,fl,&fl[d],d,k0,l0,m0,n0,mg,z);
+                else simpson4(k0x,l1x,m0x,n0x,fl,&fl[d],d,k0,l0,m0,n0,mg,z);
+      break;
+    case IDERFREE:
+      kodf(fl,&fl[d],mg,k0,l0);
+      break;
+    default:
+      mut_printf("Unknown integration type in tube_constants().\n");
+  }
+
+  if (deb>0)
+  { mut_printf("constants:\n");
+    mut_printf("  k0: %8.5f %8.5f %8.5f %8.5f\n",k0[0],k0[1],k0[2],k0[3]);
+    mut_printf("  l0: %8.5f %8.5f %8.5f\n",l0[0],l0[1],l0[2]);
+    mut_printf("  m0: %8.5f %8.5f\n",m0[0],m0[1]);
+    mut_printf("  n0: %8.5f\n",n0[0]);
+    if (d==2) mut_printf("  check: %8.5f\n",(k0[0]+k0[2]+l0[1]+m0[0])/(2*PI));
+    if (d==3) mut_printf("  check: %8.5f\n",(l0[0]+l0[2]+m0[1]+n0[0])/(4*PI));
+  }
+
+  if (aw) free(wk);
+
+  kap[0] = k0[0];
+  if (kap_terms==1) return(1);
+  kap[1] = l0[0]/2;
+  if ((kap_terms==2) | (d==1)) return(2);
+  kap[2] = (k0[2]+l0[1]+m0[0])/(2*PI);
+  if ((kap_terms==3) | (d==2)) return(3);
+  kap[3] = (l0[2]+m0[1]+n0[0])/(4*PI);
+  return(4);
+}
+/*
+ * Copyright 1996-2006 Catherine Loader.
+ */
+/*
+ *   Copyright (c) 1998-2006 Catherine Loader.
+ *
+ *   Computes the critical values from constants kappa0 etc
+ *   and significance level.
+ */
+
+#include <math.h>
+#include "tube.h"
+
+#define LOGPI 1.144729885849400174143427
+
+/* area(d) = 2 pi^(d/2) / Gamma(d/2)
+ *         = surface area of unit sphere in R^d
+ */
+static double A[10] =
+  { 1, /* d=0, whatever */
+    2,
+    6.2831853071795864770, /* 2*pi */
+    12.566370614359172954, /* 4*pi */
+    19.739208802178717238, /* 2*pi^2 */
+    26.318945069571622985, /* 8/3*pi^2 */
+    31.006276680299820177, /* pi^3 */
+    33.073361792319808190, /* 16/15*pi^3 */
+    32.469697011334145747, /* 1/3*pi^4 */
+    29.686580124648361825  /* 32/105*pi^4 */
+  };
+
+double area(d)
+int d;
+{ if (d<10) return(A[d]);
+  return(2*exp(d*LOGPI/2.0-mut_lgammai(d)));
+}
+
+double tailp_uniform(c,k0,m,d,s,n)
+double c, *k0, n;
+int m, d, s;
+{ int i;
+  double p;
+  p = 0.0;
+  for (i=0; i<m; i++) if (k0[i] != 0.0)
+    p += k0[i] * ibeta(1-c*c,(n-d+i-1)/2.0,(d+1-i)/2.0) / area(d+1-i);
+  return( (s==TWO_SIDED) ? 2*p : p );
+}
+
+double tailp_gaussian(c,k0,m,d,s,n)
+double c, *k0, n;
+int m, d, s;
+{ int i;
+  double p;
+  p = 0.0;
+  for (i=0; i<m; i++) if (k0[i] != 0.0)
+    p += k0[i] * (1-pchisq(c*c,(double) d+1-i)) / area(d+1-i);
+  return( (s==TWO_SIDED) ? 2*p : p );
+}
+
+double tailp_tprocess(c,k0,m,d,s,n)
+double c, *k0, n;
+int m, d, s;
+{ int i;
+  double p;
+  p = 0.0;
+  for (i=0; i<m; i++) if (k0[i] != 0.0)
+    p += k0[i] * (1-pf(c*c/(d+1-i),(double) d+1-i, n)) / area(d+1-i);
+  return( (s==TWO_SIDED) ? 2*p : p );
+}
+
+double taild_uniform(c,k0,m,d,s,n)
+double c, *k0, n;
+int m, d, s;
+{ int i;
+  double p;
+  p = 0.0;
+  for (i=0; i<m; i++) if (k0[i] != 0.0)
+    p += k0[i] * 2*c*dbeta(1-c*c,(n-d+i-1)/2.0,(d+1-i)/2.0,0) / area(d+1-i);
+  return( (s==TWO_SIDED) ? 2*p : p );
+}
+
+double taild_gaussian(c,k0,m,d,s,n)
+double c, *k0, n;
+int m, d, s;
+{ int i;
+  double p;
+  p = 0.0;
+  for (i=0; i<m; i++) if (k0[i] != 0.0)
+    p += k0[i] * 2*c*dchisq(c*c,(double) d+1-i,0) / area(d+1-i);
+  return( (s==TWO_SIDED) ? 2*p : p );
+}
+
+double taild_tprocess(c,k0,m,d,s,n)
+double c, *k0, n;
+int m, d, s;
+{ int i;
+  double p;
+  p = 0.0;
+  for (i=0; i<m; i++) if (k0[i] != 0.0)
+    p += k0[i] * 2*c*df(c*c/(d+1-i),(double) d+1-i, n,0) / ((d+1-i)*area(d+1-i));
+  return( (s==TWO_SIDED) ? 2*p : p );
+}
+
+double tailp(c,k0,m,d,s,nu, process)
+double c, *k0, nu;
+int m, d, s, process;
+{ switch(process)
+  { case UNIF:  return(tailp_uniform(c,k0,m,d,s,nu));
+    case GAUSS: return(tailp_gaussian(c,k0,m,d,s,nu));
+    case TPROC: return(tailp_tprocess(c,k0,m,d,s,nu));
+  }
+  mut_printf("taild: unknown process.\n");
+  return(0.0);
+}
+
+double taild(c,k0,m,d,s,nu, process)
+double c, *k0, nu;
+int m, d, s, process;
+{ switch(process)
+  { case UNIF:  return(taild_uniform(c,k0,m,d,s,nu));
+    case GAUSS: return(taild_gaussian(c,k0,m,d,s,nu));
+    case TPROC: return(taild_tprocess(c,k0,m,d,s,nu));
+  }
+  mut_printf("taild: unknown process.\n");
+  return(0.0);
+}
+
+double critval(alpha,k0,m,d,s,nu,process)
+double alpha, *k0, nu;
+int m, d, s, process;
+{ double c, cn, c0, c1, tp, td;
+  int j, maxit;
+  double (*tpf)(), (*tdf)();
+
+  maxit = 20;
+  if (m<0)
+  { mut_printf("critval: no terms?\n");
+    return(2.0);
+  }
+  if (m>d+1) m = d+1;
+  if ((alpha<=0) | (alpha>=1))
+  { mut_printf("critval: invalid alpha %8.5f\n",alpha);
+    return(2.0);
+  }
+  if (alpha>0.5)
+    mut_printf("critval: A mighty large tail probability alpha=%8.5f\n",alpha);
+  if (m==0) { d = 0; k0[0] = 1; m = 1; }
+
+  switch(process)
+  { case UNIF:
+      c = 0.5; c0 = 0.0; c1 = 1.0;
+      tpf = tailp_uniform;
+      tdf = taild_uniform;
+      break;
+    case GAUSS:
+      c = 2.0; c0 = 0.0; c1 = 0.0;
+      tpf = tailp_gaussian;
+      tdf = taild_gaussian;
+      break;
+    case TPROC:
+      c = 2.0; c0 = 0.0; c1 = 0.0;
+      tpf = tailp_tprocess;
+      tdf = taild_tprocess;
+      break;
+    default:
+      mut_printf("critval: unknown process.\n");
+      return(0.0);
+  }
+
+  for (j=0; j<maxit; j++)
+  { tp = tpf(c,k0,m,d,s,nu)-alpha;
+    td = tdf(c,k0,m,d,s,nu);
+    if (tp>0) c0 = c;
+    if (tp<0) c1 = c;
+    cn = c + tp/td;
+    if (cn<c0) cn = (c+c0)/2;
+    if ((c1>0.0) && (cn>c1)) cn = (c+c1)/2;
+    c = cn;
+    if (fabs(tp/alpha)<1.0e-10) return(c);
+  }
+  return(c);
+}