`123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113` ``````/* ADinf.c \$Revision: 1.1 \$ \$Date: 2014/06/09 10:18:10 \$ Original C code by G. and J. Marsaglia R interface by Adrian Baddeley */ #include double ADinf(double z); /* A procedure for evaluating the limiting distribution of the Anderson-Darling statistic A_n=-n-(1/n)[ln(x_1(1-x_n)+3ln(x_2(1-x_{n-1})+5ln(x_3(1-x_{n-2})+... +(2n-1)ln(x_n(1-x_1))] where x_1infty} Pr[A_n150.) return 0.; a=2.22144146907918*exp(-t)/sqrt(t); b=3.93740248643060*2.*cPhi(sqrt(2*t));/* initialization requires cPhi */ /*if you have erfc(), replace 2*cPhi(sqrt(2*t)) with erfc(sqrt(t))*/ r=z*.125; f=a+b*r; for(i=1;i<200;i++) { c=((i-.5-t)*b+t*a)/i; a=b; b=c; r*=z/(8*i+8); if(fabs(r)<1e-40 || fabs(c)<1.e-40) return f; fnew=f+c*r; if(f==fnew) return f; f=fnew; } return f; } double ADinf(double z){ int j; double ad,adnew,r; if(z<.01) return 0.; /* avoids exponent limits; ADinf(.01)=.528e-52 */ r=1./z; ad=r*ADf(z,0); for(j=1;j<100;j++){ r*=(.5-j)/j; adnew=ad+(4*j+1)*r*ADf(z,j); if(ad==adnew) {return ad;} ad=adnew; } return ad; } /* Complementary normal distribution function cPhi(x) = integral from x to infinity of phi(x)=exp(-.5*t^2)/sqrt(2*pi) 13-15 digit accuracy for abs(x)<16. Stores R(0),R(2),R(4),...,R(16), with cPhi(x)=R(x)*phi(x), phi normal density, then uses Taylor series for R(z+h)=R(z)+hR'(z)+(1/2)h^2R''(z)+... with -10) ? s: 1-s); } } /* end i loop */ /* If not converged, return last estimate */ return ((x>0) ? s: 1-s); } /* R interface */ void ADprobExactInf(double *a, int *na, double *prob) { int i, m; m = *na; for(i = 0; i < m; i++) prob[i] = ADinf(a[i]); } ``````