private static void swilk(boolean[] init,double[] x,int n,int n1,int n2,double[] a,double[] w,double[] pw,int[] ifault){
  pw[0]=1.0;
  if (w[0] >= 0.0) {
    w[0]=1.0;
  }
  double an=n;
  ifault[0]=3;
  int nn2=n / 2;
  if (n2 < nn2) {
    return;
  }
  ifault[0]=1;
  if (n < 3) {
    return;
  }
  if (!init[0]) {
    if (n == 3) {
      a[1]=SQRTH;
    }
 else {
      double an25=an + 0.25;
      double summ2=0.0;
      for (int i=1; i <= n2; ++i) {
        a[i]=ppnd((i - TH) / an25);
        summ2+=a[i] * a[i];
      }
      summ2*=2.0;
      double ssumm2=Math.sqrt(summ2);
      double rsn=1.0 / Math.sqrt(an);
      double a1=poly(C1,6,rsn) - a[1] / ssumm2;
      int i1;
      double fac;
      if (n > 5) {
        i1=3;
        double a2=-a[2] / ssumm2 + poly(C2,6,rsn);
        fac=Math.sqrt((summ2 - 2.0 * a[1] * a[1] - 2.0 * a[2] * a[2]) / (1.0 - 2.0 * a1 * a1 - 2.0 * a2 * a2));
        a[1]=a1;
        a[2]=a2;
      }
 else {
        i1=2;
        fac=Math.sqrt((summ2 - 2.0 * a[1] * a[1]) / (1.0 - 2.0 * a1 * a1));
        a[1]=a1;
      }
      for (int i=i1; i <= nn2; ++i) {
        a[i]=-a[i] / fac;
      }
    }
    init[0]=true;
  }
  if (n1 < 3) {
    return;
  }
  int ncens=n - n1;
  ifault[0]=4;
  if (ncens < 0 || (ncens > 0 && n < 20)) {
    return;
  }
  ifault[0]=5;
  double delta=ncens / an;
  if (delta > 0.8) {
    return;
  }
  double w1, xx;
  if (w[0] < 0.0) {
    w1=1.0 + w[0];
    ifault[0]=0;
  }
 else {
    ifault[0]=6;
    double range=x[n1] - x[1];
    if (range < SMALL) {
      return;
    }
    ifault[0]=7;
    xx=x[1] / range;
    double sx=xx;
    double sa=-a[1];
    int j=n - 1;
    for (int i=2; i <= n1; ++i) {
      double xi=x[i] / range;
      sx+=xi;
      if (i != j) {
        sa+=sign(1,i - j) * a[Math.min(i,j)];
      }
      xx=xi;
      --j;
    }
    ifault[0]=0;
    if (n > 5000) {
      ifault[0]=2;
    }
    sa/=n1;
    sx/=n1;
    double ssa=0.0;
    double ssx=0.0;
    double sax=0.0;
    j=n;
    double asa;
    for (int i=1; i <= n1; ++i) {
      if (i != j) {
        asa=sign(1,i - j) * a[Math.min(i,j)] - sa;
      }
 else {
        asa=-sa;
      }
      double xsx=x[i] / range - sx;
      ssa+=asa * asa;
      ssx+=xsx * xsx;
      sax+=asa * xsx;
      --j;
    }
    double ssassx=Math.sqrt(ssa * ssx);
    w1=(ssassx - sax) * (ssassx + sax) / (ssa * ssx);
  }
  w[0]=1.0 - w1;
  if (n == 3) {
    pw[0]=PI6 * (Math.asin(Math.sqrt(w[0])) - STQR);
    return;
  }
  double y=Math.log(w1);
  xx=Math.log(an);
  double m=0.0;
  double s=1.0;
  if (n <= 11) {
    double gamma=poly(G,2,an);
    if (y >= gamma) {
      pw[0]=SMALL;
      return;
    }
    y=-Math.log(gamma - y);
    m=poly(C3,4,an);
    s=Math.exp(poly(C4,4,an));
  }
 else {
    m=poly(C5,4,xx);
    s=Math.exp(poly(C6,3,xx));
  }
  if (ncens > 0) {
    double ld=-Math.log(delta);
    double bf=1.0 + xx * BF1;
    double z90f=Z90 + bf * Math.pow(poly(C7,2,Math.pow(XX90,xx)),ld);
    double z95f=Z95 + bf * Math.pow(poly(C8,2,Math.pow(XX95,xx)),ld);
    double z99f=Z99 + bf * Math.pow(poly(C9,2,xx),ld);
    double zfm=(z90f + z95f + z99f) / 3.0;
    double zsd=(Z90 * (z90f - zfm) + Z95 * (z95f - zfm) + Z99 * (z99f - zfm)) / ZSS;
    double zbar=zfm - zsd * ZM;
    m+=zbar * s;
    s*=zsd;
  }
  pw[0]=alnorm((y - m) / s,UPPER);
}

private static void swilk ( boolean [ ] init , double [ ] x , int n , int n1 , int n2 , double [ ] a , double [ ] w , double [ ] pw , int [ ] ifault ) { pw [ 0 ] = 1.0 ; if ( w [ 0 ] >= 0.0 ) { w [ 0 ] = 1.0 ; } double an = n ; ifault [ 0 ] = 3 ; int nn2 = n / 2 ; if ( n2 < nn2 ) { return ; } ifault [ 0 ] = 1 ; if ( n < 3 ) { return ; } if ( ! init [ 0 ] ) { if ( n == 3 ) { a [ 1 ] = SQRTH ; } else { double an25 = an + 0.25 ; double summ2 = 0.0 ; for ( int i = 1 ; i <= n2 ; ++ i ) { a [ i ] = ppnd ( ( i - TH ) / an25 ) ; summ2 += a [ i ] * a [ i ] ; } summ2 *= 2.0 ; double ssumm2 = Math . sqrt ( summ2 ) ; double rsn = 1.0 / Math . sqrt ( an ) ; double a1 = poly ( C1 , 6 , rsn ) - a [ 1 ] / ssumm2 ; int i1 ; double fac ; if ( n > 5 ) { i1 = 3 ; double a2 = - a [ 2 ] / ssumm2 + poly ( C2 , 6 , rsn ) ; fac = Math . sqrt ( ( summ2 - 2.0 * a [ 1 ] * a [ 1 ] - 2.0 * a [ 2 ] * a [ 2 ] ) / ( 1.0 - 2.0 * a1 * a1 - 2.0 * a2 * a2 ) ) ; a [ 1 ] = a1 ; a [ 2 ] = a2 ; } else { i1 = 2 ; fac = Math . sqrt ( ( summ2 - 2.0 * a [ 1 ] * a [ 1 ] ) / ( 1.0 - 2.0 * a1 * a1 ) ) ; a [ 1 ] = a1 ; } for ( int i = i1 ; i <= nn2 ; ++ i ) { a [ i ] = - a [ i ] / fac ; } } init [ 0 ] = true ; } if ( n1 < 3 ) { return ; } int ncens = n - n1 ; ifault [ 0 ] = 4 ; if ( ncens < 0 || ( ncens > 0 && n < 20 ) ) { return ; } ifault [ 0 ] = 5 ; double delta = ncens / an ; if ( delta > 0.8 ) { return ; } double w1 , xx ; if ( w [ 0 ] < 0.0 ) { w1 = 1.0 + w [ 0 ] ; ifault [ 0 ] = 0 ; } else { ifault [ 0 ] = 6 ; double range = x [ n1 ] - x [ 1 ] ; if ( range < SMALL ) { return ; } ifault [ 0 ] = 7 ; xx = x [ 1 ] / range ; double sx = xx ; double sa = - a [ 1 ] ; int j = n - 1 ; for ( int i = 2 ; i <= n1 ; ++ i ) { double xi = x [ i ] / range ; sx += xi ; if ( i != j ) { sa += sign ( 1 , i - j ) * a [ Math . min ( i , j ) ] ; } xx = xi ; -- j ; } ifault [ 0 ] = 0 ; if ( n > 5000 ) { ifault [ 0 ] = 2 ; } sa /= n1 ; sx /= n1 ; double ssa = 0.0 ; double ssx = 0.0 ; double sax = 0.0 ; j = n ; double asa ; for ( int i = 1 ; i <= n1 ; ++ i ) { if ( i != j ) { asa = sign ( 1 , i - j ) * a [ Math . min ( i , j ) ] - sa ; } else { asa = - sa ; } double xsx = x [ i ] / range - sx ; ssa += asa * asa ; ssx += xsx * xsx ; sax += asa * xsx ; -- j ; } double ssassx = Math . sqrt ( ssa * ssx ) ; w1 = ( ssassx - sax ) * ( ssassx + sax ) / ( ssa * ssx ) ; } w [ 0 ] = 1.0 - w1 ; if ( n == 3 ) { pw [ 0 ] = PI6 * ( Math . asin ( Math . sqrt ( w [ 0 ] ) ) - STQR ) ; return ; } double y = Math . log ( w1 ) ; xx = Math . log ( an ) ; double m = 0.0 ; double s = 1.0 ; if ( n <= 11 ) { double gamma = poly ( G , 2 , an ) ; if ( y >= gamma ) { pw [ 0 ] = SMALL ; return ; } y = - Math . log ( gamma - y ) ; m = poly ( C3 , 4 , an ) ; s = Math . exp ( poly ( C4 , 4 , an ) ) ; } else { m = poly ( C5 , 4 , xx ) ; s = Math . exp ( poly ( C6 , 3 , xx ) ) ; } if ( ncens > 0 ) { double ld = - Math . log ( delta ) ; double bf = 1.0 + xx * BF1 ; double z90f = Z90 + bf * Math . pow ( poly ( C7 , 2 , Math . pow ( XX90 , xx ) ) , ld ) ; double z95f = Z95 + bf * Math . pow ( poly ( C8 , 2 , Math . pow ( XX95 , xx ) ) , ld ) ; double z99f = Z99 + bf * Math . pow ( poly ( C9 , 2 , xx ) , ld ) ; double zfm = ( z90f + z95f + z99f ) / 3.0 ; double zsd = ( Z90 * ( z90f - zfm ) + Z95 * ( z95f - zfm ) + Z99 * ( z99f - zfm ) ) / ZSS ; double zbar = zfm - zsd * ZM ; m += zbar * s ; s *= zsd ; } pw [ 0 ] = alnorm ( ( y - m ) / s , UPPER ) ; }
algorithm as r94 appl . statist . ( 1995 ) vol . 44 , no . 4 < p > calculates shapiro - wilk normality test and p - value for sample sizes 3 < = n < = 5000 . handles censored or uncensored data . corrects as 181 , which was found to be inaccurate for n > 50 . < p > note : semi - strange porting kludge alert . fortran allows subroutine arguments to be modified by the called routine ( passed by reference , not value ) , and the original code for this routine makes use of that feature to return multiple results . to avoid changing the code any more than necessary , i ' ve used java arrays to simulate this pass - by - reference feature . specifically , in the original code w , pw , and ifault are output results , not input parameters . pass in double [ 1 ] arrays for w and pw , and extract the computed values from the [ 0 ] element on return . the argument init is both input and output ; use a boolean [ 1 ] array and initialize [ 0 ] to false before the first call . the routine will update the value to true to record that initialization has been performed , to speed up subsequent calls on the same data set . note that although the contents of a [ ] will be computed by the routine on the first call , the caller must still allocate the array space and pass the unfilled array in to the subroutine . the routine will set the contents but not allocate the space . < p > as described above with the constants , the data arrays x [ ] and a [ ] are referenced with a base element of 1 ( like fortran ) instead of 0 ( like java ) to avoid screwing up the algorithm . to pass in 100 data points , declare x [ 101 ] and fill elements x [ 1 ] through x [ 100 ] with data . x [ 0 ] will be ignored .
