*This class contains a Java version of the FORTRAN routine *rnor which generates standard normal random numbers. *The FORTRAN version of rnor was written by David Kahaner and *George Marsaglia. *

*This class was translated by a statistician from the FORTRAN routine. *It is NOT an official translation. If an official translation *becomes available, it should be used rather than this version. *

*If you have suggestions for improving this *code, please contact Steve Verrill at steve@swst.org * *@author Steve Verrill *@version .5 --- March 9, 2004 * */ public class Normal extends Object { /* As written the routine assumes that the the machine has words that are at least 32 bits in length. If the machine has words that are 64 bits in length, then mdig (see below) could be altered. In this case however the values of m1 and m2 would have to be altered to be in accord with the FORTRAN statements: MDIG = 64 M1 = 2**(MDIG-2) + (2**(MDIG-2)-1) M2 = 2**(MDIG/2) Also m1 and m2 would have to be declared to be long. */ // 2^31 - 1: static final int m1 = 2147483647; // 2^16: static final int m2 = 65536; static final double aa = 12.37586; static final double b = .4878992; static final double c = 12.67706; static final double c1 = .9689279; static final double c2 = 1.301198; static final double pc = .1958303E-1; static final double xn = 2.776994; static final double v[] = {0.0,.3409450, .4573146, .5397793, .6062427, .6631691, .7136975, .7596125, .8020356, .8417227, .8792102, .9148948, .9490791, .9820005, 1.0138492, 1.0447810, 1.0749254, 1.1043917, 1.1332738, 1.1616530, 1.1896010, 1.2171815, 1.2444516, 1.2714635, 1.2982650, 1.3249008, 1.3514125, 1.3778399, 1.4042211, 1.4305929, 1.4569915, 1.4834526, 1.5100121, 1.5367061, 1.5635712, 1.5906454, 1.6179680, 1.6455802, 1.6735255, 1.7018503, 1.7306045, 1.7598422, 1.7896223, 1.8200099, 1.8510770, 1.8829044, 1.9155830, 1.9492166, 1.9839239, 2.0198430, 2.0571356, 2.0959930, 2.1366450, 2.1793713, 2.2245175, 2.2725185, 2.3239338, 2.3795007, 2.4402218, 2.5075117, 2.5834658, 2.6713916, 2.7769943, 2.7769943, 2.7769943, 2.7769943}; int i1; int j1; int m[] = new int[18]; double rmax; double w[] = new double[66]; Uniform uniform; // constructor public Normal(int iseed) { int jseed,k0,k1,j0,i; double dummy; jseed = Math.min(Math.abs(iseed),m1); if (jseed%2 == 0) { jseed--; } k0 = 9069%m2; k1 = 9069/m2; j0 = jseed%m2; j1 = jseed/m2; for (i = 1; i <= 17; i++) { jseed = j0*k0; j1 = (jseed/m2 + j0*k1 + j1*k0)%(m2/2); j0 = jseed%m2; m[i] = j0 + m2*j1; } i1 = 5; j1 = 17; rmax = 1.0/m1; uniform = new Uniform(iseed); // To make results match those of the FORTRAN version // we include the following dummy call: dummy = uniform.uni(); for (i = 1; i <= 65; i++) { w[i] = rmax*v[i]; } } public double rnor() { int i,j,brk; double rnorm; double x,y,s; double tttttt; i = m[i1] - m[j1]; if (i < 0) { i += m1; } m[j1] = i; i1--; if (i1 == 0) { i1 = 17; } j1--; if (j1 == 0) { j1 = 17; } j = i%64 + 1; rnorm = i*w[j+1]; if ( ((i/m2)/2)*2 == (i/m2) ) { rnorm = -rnorm; } if (Math.abs(rnorm) <= v[j]) { return rnorm; } x = (Math.abs(rnorm) - v[j])/(v[j+1] - v[j]); y = uniform.uni(); s = x + y; if (s > c2) { if (rnorm < 0.0) { rnorm = b*x - b; } else { rnorm = b - b*x; } return rnorm; } else { if (s <= c1) { return rnorm; } if (y > c - aa*Math.exp(-.5*(b-b*x)*(b-b*x))) { if (rnorm < 0.0) { rnorm = b*x - b; } else { rnorm = b - b*x; } return rnorm; } if (Math.exp(-.5*v[j+1]*v[j+1]) + y*pc/v[j+1] <= Math.exp(-.5*rnorm*rnorm)) { return rnorm; } brk = 0; while (brk == 0) { // A version of the following while loop was added on 10/23/89 // by Steve Verrill to the FORTRAN code to ensure that we are not // attempting to take the log of 0.0 (something that the original code // was occasionally doing) tttttt = 0.0; while (tttttt == 0.0) { tttttt = uniform.uni(); } s = xn - Math.log(tttttt)/xn; // A version of the following while loop was added on 10/23/89 // by Steve Verrill to the FORTRAN code to ensure that we are not // attempting to take the log of 0.0 (something that the original code // was occasionally doing) tttttt = 0.0; while (tttttt == 0.0) { tttttt = uniform.uni(); } if (3.855849 + Math.log(tttttt) - xn*s <= -.5*s*s) { brk = 1; } } if (rnorm < 0.0) { return -s; } else { return s; } } } }