= 0) f = f + f + 1; else { f += f; p = p + q; } } while (f < fraction_one); be_careful = p - q; if (be_careful + p >= 0) incr(f); if (negative) return (-(f + n)); else return (f + n); } } @ @c static int take_frac(int q, int f) { int p; /* the fraction so far */ int n; /* additional multiple of $q$ */ register int be_careful; /* disables certain compiler optimizations */ boolean negative = false; /* should the result be negated? */ /* Reduce to the case that |f>=0| and |q>0| */ if (f < 0) { negate(f); negative = true; } if (q < 0) { negate(q); negative = !negative; } if (f < fraction_one) { n = 0; } else { n = f / fraction_one; f = f % fraction_one; if (q <= el_gordo / n) { n = n * q; } else { arith_error = true; n = el_gordo; } } f = f + fraction_one; /* Compute $p=\lfloor qf/2^{28}+{1\over2}\rfloor-q$ */ /* The invariant relations in this case are (i)~$\lfloor(qf+p)/2^k\rfloor =\lfloor qf_0/2^{28}+{1\over2}\rfloor$, where $k$ is an integer and $f_0$ is the original value of~$f$; (ii)~$2^k\L f<2^{k+1}$. */ p = fraction_half; /* that's $2^{27}$; the invariants hold now with $k=28$ */ if (q < fraction_four) { do { if (odd(f)) p = halfp(p + q); else p = halfp(p); f = halfp(f); } while (f != 1); } else { do { if (odd(f)) p = p + halfp(q - p); else p = halfp(p); f = halfp(f); } while (f != 1); } be_careful = n - el_gordo; if (be_careful + p > 0) { arith_error = true; n = el_gordo - p; } if (negative) return (-(n + p)); else return (n + p); } @ The subroutines for logarithm and exponential involve two tables. The first is simple: |two_to_the[k]| equals $2^k$. The second involves a bit more calculation, which the author claims to have done correctly: |spec_log[k]| is $2^{27}$ times $\ln\bigl(1/(1-2^{-k})\bigr)= 2^{-k}+{1\over2}2^{-2k}+{1\over3}2^{-3k}+\cdots\,$, rounded to the nearest integer. @c static int two_to_the[31]; /* powers of two */ static int spec_log[29]; /* special logarithms */ @ @c void initialize_arithmetic(void) { int k; two_to_the[0] = 1; for (k = 1; k <= 30; k++) two_to_the[k] = 2 * two_to_the[k - 1]; spec_log[1] = 93032640; spec_log[2] = 38612034; spec_log[3] = 17922280; spec_log[4] = 8662214; spec_log[5] = 4261238; spec_log[6] = 2113709; spec_log[7] = 1052693; spec_log[8] = 525315; spec_log[9] = 262400; spec_log[10] = 131136; spec_log[11] = 65552; spec_log[12] = 32772; spec_log[13] = 16385; for (k = 14; k <= 27; k++) spec_log[k] = two_to_the[27 - k]; spec_log[28] = 1; } @ @c static int m_log(int x) { int y, z; /* auxiliary registers */ int k; /* iteration counter */ if (x <= 0) { /* Handle non-positive logarithm */ print_err("Logarithm of "); print_scaled(x); tprint(" has been replaced by 0"); help2("Since I don't take logs of non-positive numbers,", "I'm zeroing this one. Proceed, with fingers crossed."); error(); return 0; } else { y = 1302456956 + 4 - 100; /* $14\times2^{27}\ln2\approx1302456956.421063$ */ z = 27595 + 6553600; /* and $2^{16}\times .421063\approx 27595$ */ while (x < fraction_four) { x += x; y = y - 93032639; z = z - 48782; } /* $2^{27}\ln2\approx 93032639.74436163$ and $2^{16}\times.74436163\approx 48782$ */ y = y + (z / unity); k = 2; while (x > fraction_four + 4) { /* Increase |k| until |x| can be multiplied by a factor of $2^{-k}$, and adjust $y$ accordingly */ z = ((x - 1) / two_to_the[k]) + 1; /* $z=\lceil x/2^k\rceil$ */ while (x < fraction_four + z) { z = halfp(z + 1); k = k + 1; } y = y + spec_log[k]; x = x - z; } return (y / 8); } } @ The following somewhat different subroutine tests rigorously if $ab$ is greater than, equal to, or less than~$cd$, given integers $(a,b,c,d)$. In most cases a quick decision is reached. The result is $+1$, 0, or~$-1$ in the three respective cases. @c static int ab_vs_cd(int a, int b, int c, int d) { int q, r; /* temporary registers */ /* Reduce to the case that |a,c>=0|, |b,d>0| */ if (a < 0) { negate(a); negate(b); } if (c < 0) { negate(c); negate(d); } if (d <= 0) { if (b >= 0) return (((a == 0 || b == 0) && (c == 0 || d == 0)) ? 0 : 1); if (d == 0) return (a == 0 ? 0 : -1); q = a; a = c; c = q; q = -b; b = -d; d = q; } else if (b <= 0) { if (b < 0 && a > 0) return -1; return (c == 0 ? 0 : -1); } while (1) { q = a / d; r = c / b; if (q != r) return (q > r ? 1 : -1); q = a % d; r = c % b; if (r == 0) return (q == 0 ? 0 : 1); if (q == 0) return -1; a = b; b = q; c = d; d = r; /* now |a>d>0| and |c>b>0| */ } } @ To consume a random integer, the program below will say `|next_random|' and then it will fetch |randoms[j_random]|. @c #define next_random() do { \ if (j_random==0) new_randoms(); else decr(j_random); \ } while (0) static void new_randoms(void) { int k; /* index into |randoms| */ int x; /* accumulator */ for (k = 0; k <= 23; k++) { x = randoms[k] - randoms[k + 31]; if (x < 0) x = x + fraction_one; randoms[k] = x; } for (k = 24; k <= 54; k++) { x = randoms[k] - randoms[k - 24]; if (x < 0) x = x + fraction_one; randoms[k] = x; } j_random = 54; } @ To initialize the |randoms| table, we call the following routine. @c void init_randoms(int seed) { int j, jj, k; /* more or less random integers */ int i; /* index into |randoms| */ j = abs(seed); while (j >= fraction_one) j = halfp(j); k = 1; for (i = 0; i <= 54; i++) { jj = k; k = j - k; j = jj; if (k < 0) k = k + fraction_one; randoms[(i * 21) % 55] = j; } new_randoms(); new_randoms(); new_randoms(); /* ``warm up'' the array */ } @ To produce a uniform random number in the range |0<=u=u>x| or |0=u=x|, given a |scaled| value~|x|, we proceed as shown here. Note that the call of |take_frac| will produce the values 0 and~|x| with about half the probability that it will produce any other particular values between 0 and~|x|, because it rounds its answers. @c int unif_rand(int x) { int y; /* trial value */ next_random(); y = take_frac(abs(x), randoms[j_random]); if (y == abs(x)) return 0; else if (x > 0) return y; else return -y; } @ Finally, a normal deviate with mean zero and unit standard deviation can readily be obtained with the ratio method (Algorithm 3.4.1R in {\sl The Art of Computer Programming\/}). @c int norm_rand(void) { int x, u, l; /* what the book would call $2^{16}X$, $2^{28}U$, and $-2^{24}\ln U$ */ do { do { next_random(); x = take_frac(112429, randoms[j_random] - fraction_half); /* $2^{16}\sqrt{8/e}\approx 112428.82793$ */ next_random(); u = randoms[j_random]; } while (abs(x) >= u); x = make_frac(x, u); l = 139548960 - m_log(u); /* $2^{24}\cdot12\ln2\approx139548959.6165$ */ } while (ab_vs_cd(1024, l, x, x) < 0); return x; } @ This function could also be expressed as a macro, but it is a useful breakpoint for debugging. @c int fix_int(int val, int min, int max) { return (val < min ? min : (val > max ? max : val)); }