question: 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 ) ; }
target:  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.