1% arithmetic.w 2% 3% Copyright 2009-2010 Taco Hoekwater <taco@@luatex.org> 4% 5% This file is part of LuaTeX. 6% 7% LuaTeX is free software; you can redistribute it and/or modify it under 8% the terms of the GNU General Public License as published by the Free 9% Software Foundation; either version 2 of the License, or (at your 10% option) any later version. 11% 12% LuaTeX is distributed in the hope that it will be useful, but WITHOUT 13% ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or 14% FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public 15% License for more details. 16% 17% You should have received a copy of the GNU General Public License along 18% with LuaTeX; if not, see <http://www.gnu.org/licenses/>. 19 20\def\MP{MetaPost} 21 22@ @c 23 24 25#include "ptexlib.h" 26 27@ The principal computations performed by \TeX\ are done entirely in terms of 28integers less than $2^{31}$ in magnitude; and divisions are done only when both 29dividend and divisor are nonnegative. Thus, the arithmetic specified in this 30program can be carried out in exactly the same way on a wide variety of 31computers, including some small ones. Why? Because the arithmetic 32calculations need to be spelled out precisely in order to guarantee that 33\TeX\ will produce identical output on different machines. If some 34quantities were rounded differently in different implementations, we would 35find that line breaks and even page breaks might occur in different places. 36Hence the arithmetic of \TeX\ has been designed with care, and systems that 37claim to be implementations of \TeX82 should follow precisely the 38@:TeX82}{\TeX82@> 39calculations as they appear in the present program. 40 41(Actually there are three places where \TeX\ uses |div| with a possibly negative 42numerator. These are harmless; see |div| in the index. Also if the user 43sets the \.{\\time} or the \.{\\year} to a negative value, some diagnostic 44information will involve negative-numerator division. The same remarks 45apply for |mod| as well as for |div|.) 46 47Here is a routine that calculates half of an integer, using an 48unambiguous convention with respect to signed odd numbers. 49 50@c 51int half(int x) 52{ 53 if (odd(x)) 54 return ((x + 1) / 2); 55 else 56 return (x / 2); 57} 58 59 60@ The following function is used to create a scaled integer from a given decimal 61fraction $(.d_0d_1\ldots d_{k-1})$, where |0<=k<=17|. The digit $d_i$ is 62given in |dig[i]|, and the calculation produces a correctly rounded result. 63 64@c 65scaled round_decimals(int k) 66{ /* converts a decimal fraction */ 67 int a; /* the accumulator */ 68 a = 0; 69 while (k-- > 0) { 70 a = (a + dig[k] * two) / 10; 71 } 72 return ((a + 1) / 2); 73} 74 75 76@ Conversely, here is a procedure analogous to |print_int|. If the output 77of this procedure is subsequently read by \TeX\ and converted by the 78|round_decimals| routine above, it turns out that the original value will 79be reproduced exactly; the ``simplest'' such decimal number is output, 80but there is always at least one digit following the decimal point. 81 82The invariant relation in the \&{repeat} loop is that a sequence of 83decimal digits yet to be printed will yield the original number if and only if 84they form a fraction~$f$ in the range $s-\delta\L10\cdot2^{16}f<s$. 85We can stop if and only if $f=0$ satisfies this condition; the loop will 86terminate before $s$ can possibly become zero. 87 88@c 89void print_scaled(scaled s) 90{ /* prints scaled real, rounded to five digits */ 91 scaled delta; /* amount of allowable inaccuracy */ 92 char buffer[20]; 93 int i = 0; 94 if (s < 0) { 95 print_char('-'); 96 negate(s); /* print the sign, if negative */ 97 } 98 print_int(s / unity); /* print the integer part */ 99 buffer[i++] = '.'; 100 s = 10 * (s % unity) + 5; 101 delta = 10; 102 do { 103 if (delta > unity) 104 s = s + 0100000 - 50000; /* round the last digit */ 105 buffer[i++] = '0' + (s / unity); 106 s = 10 * (s % unity); 107 delta = delta * 10; 108 } while (s > delta); 109 buffer[i++] = '\0'; 110 tprint(buffer); 111} 112 113@ Physical sizes that a \TeX\ user specifies for portions of documents are 114represented internally as scaled points. Thus, if we define an `sp' (scaled 115@^sp@> 116point) as a unit equal to $2^{-16}$ printer's points, every dimension 117inside of \TeX\ is an integer number of sp. There are exactly 1184,736,286.72 sp per inch. Users are not allowed to specify dimensions 119larger than $2^{30}-1$ sp, which is a distance of about 18.892 feet (5.7583 120meters); two such quantities can be added without overflow on a 32-bit 121computer. 122 123The present implementation of \TeX\ does not check for overflow when 124@^overflow in arithmetic@> 125dimensions are added or subtracted. This could be done by inserting a 126few dozen tests of the form `\ignorespaces|if x>=010000000000 then 127@t\\{report\_overflow}@>|', but the chance of overflow is so remote that 128such tests do not seem worthwhile. 129 130\TeX\ needs to do only a few arithmetic operations on scaled quantities, 131other than addition and subtraction, and the following subroutines do most of 132the work. A single computation might use several subroutine calls, and it is 133desirable to avoid producing multiple error messages in case of arithmetic 134overflow; so the routines set the global variable |arith_error| to |true| 135instead of reporting errors directly to the user. Another global variable, 136|tex_remainder|, holds the remainder after a division. 137 138@c 139boolean arith_error; /* has arithmetic overflow occurred recently? */ 140scaled tex_remainder; /* amount subtracted to get an exact division */ 141 142 143@ The first arithmetical subroutine we need computes $nx+y$, where |x| 144and~|y| are |scaled| and |n| is an integer. We will also use it to 145multiply integers. 146 147@c 148scaled mult_and_add(int n, scaled x, scaled y, scaled max_answer) 149{ 150 if (n == 0) 151 return y; 152 if (n < 0) { 153 negate(x); 154 negate(n); 155 } 156 if (((x <= (max_answer - y) / n) && (-x <= (max_answer + y) / n))) { 157 return (n * x + y); 158 } else { 159 arith_error = true; 160 return 0; 161 } 162} 163 164@ We also need to divide scaled dimensions by integers. 165@c 166scaled x_over_n(scaled x, int n) 167{ 168 boolean negative; /* should |tex_remainder| be negated? */ 169 negative = false; 170 if (n == 0) { 171 arith_error = true; 172 tex_remainder = x; 173 return 0; 174 } else { 175 if (n < 0) { 176 negate(x); 177 negate(n); 178 negative = true; 179 } 180 if (x >= 0) { 181 tex_remainder = x % n; 182 if (negative) 183 negate(tex_remainder); 184 return (x / n); 185 } else { 186 tex_remainder = -((-x) % n); 187 if (negative) 188 negate(tex_remainder); 189 return (-((-x) / n)); 190 } 191 } 192} 193 194 195@ Then comes the multiplication of a scaled number by a fraction |n/d|, 196where |n| and |d| are nonnegative integers |<=@t$2^{16}$@>| and |d| is 197positive. It would be too dangerous to multiply by~|n| and then divide 198by~|d|, in separate operations, since overflow might well occur; and it 199would be too inaccurate to divide by |d| and then multiply by |n|. Hence 200this subroutine simulates 1.5-precision arithmetic. 201 202@c 203scaled xn_over_d(scaled x, int n, int d) 204{ 205 nonnegative_integer t, u, v, xx, dd; /* intermediate quantities */ 206 boolean positive = true; /* was |x>=0|? */ 207 if (x < 0) { 208 negate(x); 209 positive = false; 210 } 211 xx = (nonnegative_integer) x; 212 dd = (nonnegative_integer) d; 213 t = ((xx % 0100000) * (nonnegative_integer) n); 214 u = ((xx / 0100000) * (nonnegative_integer) n + (t / 0100000)); 215 v = (u % dd) * 0100000 + (t % 0100000); 216 if (u / dd >= 0100000) 217 arith_error = true; 218 else 219 u = 0100000 * (u / dd) + (v / dd); 220 if (positive) { 221 tex_remainder = (int) (v % dd); 222 return (scaled) u; 223 } else { 224 /* casts are for ms cl */ 225 tex_remainder = -(int) (v % dd); 226 return -(scaled) (u); 227 } 228} 229 230 231@ The next subroutine is used to compute the ``badness'' of glue, when a 232total~|t| is supposed to be made from amounts that sum to~|s|. According 233to {\sl The \TeX book}, the badness of this situation is $100(t/s)^3$; 234however, badness is simply a heuristic, so we need not squeeze out the 235last drop of accuracy when computing it. All we really want is an 236approximation that has similar properties. 237@:TeXbook}{\sl The \TeX book@> 238 239The actual method used to compute the badness is easier to read from the 240program than to describe in words. It produces an integer value that is a 241reasonably close approximation to $100(t/s)^3$, and all implementations 242of \TeX\ should use precisely this method. Any badness of $2^{13}$ or more is 243treated as infinitely bad, and represented by 10000. 244 245It is not difficult to prove that $$\hbox{|badness(t+1,s)>=badness(t,s) 246>=badness(t,s+1)|}.$$ The badness function defined here is capable of 247computing at most 1095 distinct values, but that is plenty. 248 249@c 250halfword badness(scaled t, scaled s) 251{ /* compute badness, given |t>=0| */ 252 int r; /* approximation to $\alpha t/s$, where $\alpha^3\approx 253 100\cdot2^{18}$ */ 254 if (t == 0) { 255 return 0; 256 } else if (s <= 0) { 257 return inf_bad; 258 } else { 259 if (t <= 7230584) 260 r = (t * 297) / s; /* $297^3=99.94\times2^{18}$ */ 261 else if (s >= 1663497) 262 r = t / (s / 297); 263 else 264 r = t; 265 if (r > 1290) 266 return inf_bad; /* $1290^3<2^{31}<1291^3$ */ 267 else 268 return ((r * r * r + 0400000) / 01000000); 269 /* that was $r^3/2^{18}$, rounded to the nearest integer */ 270 } 271} 272 273 274@ When \TeX\ ``packages'' a list into a box, it needs to calculate the 275proportionality ratio by which the glue inside the box should stretch 276or shrink. This calculation does not affect \TeX's decision making, 277so the precise details of rounding, etc., in the glue calculation are not 278of critical importance for the consistency of results on different computers. 279 280We shall use the type |glue_ratio| for such proportionality ratios. 281A glue ratio should take the same amount of memory as an 282|integer| (usually 32 bits) if it is to blend smoothly with \TeX's 283other data structures. Thus |glue_ratio| should be equivalent to 284|short_real| in some implementations of PASCAL. Alternatively, 285it is possible to deal with glue ratios using nothing but fixed-point 286arithmetic; see {\sl TUGboat \bf3},1 (March 1982), 10--27. (But the 287routines cited there must be modified to allow negative glue ratios.) 288@^system dependencies@> 289 290 291@ This section is (almost) straight from MetaPost. I had to change 292the types (use |integer| instead of |fraction|), but that should 293not have any influence on the actual calculations (the original 294comments refer to quantities like |fraction_four| ($2^{30}$), and 295that is the same as the numeric representation of |max_dimen|). 296 297I've copied the low-level variables and routines that are needed, but 298only those (e.g. |m_log|), not the accompanying ones like |m_exp|. Most 299of the following low-level numeric routines are only needed within the 300calculation of |norm_rand|. I've been forced to rename |make_fraction| 301to |make_frac| because TeX already has a routine by that name with 302a wholly different function (it creates a |fraction_noad| for math 303typesetting) -- Taco 304 305And now let's complete our collection of numeric utility routines 306by considering random number generation. 307\MP{} generates pseudo-random numbers with the additive scheme recommended 308in Section 3.6 of {\sl The Art of Computer Programming}; however, the 309results are random fractions between 0 and |fraction_one-1|, inclusive. 310 311There's an auxiliary array |randoms| that contains 55 pseudo-random 312fractions. Using the recurrence $x_n=(x_{n-55}-x_{n-31})\bmod 2^{28}$, 313we generate batches of 55 new $x_n$'s at a time by calling |new_randoms|. 314The global variable |j_random| tells which element has most recently 315been consumed. 316 317@c 318static int randoms[55]; /* the last 55 random values generated */ 319static int j_random; /* the number of unused |randoms| */ 320scaled random_seed; /* the default random seed */ 321 322@ A small bit of metafont is needed. 323 324@c 325#define fraction_half 01000000000 /* $2^{27}$, represents 0.50000000 */ 326#define fraction_one 02000000000 /* $2^{28}$, represents 1.00000000 */ 327#define fraction_four 010000000000 /* $2^{30}$, represents 4.00000000 */ 328#define el_gordo 017777777777 /* $2^{31}-1$, the largest value that \MP\ likes */ 329 330@ The |make_frac| routine produces the |fraction| equivalent of 331|p/q|, given integers |p| and~|q|; it computes the integer 332$f=\lfloor2^{28}p/q+{1\over2}\rfloor$, when $p$ and $q$ are 333positive. If |p| and |q| are both of the same scaled type |t|, 334the ``type relation'' |make_frac(t,t)=fraction| is valid; 335and it's also possible to use the subroutine ``backwards,'' using 336the relation |make_frac(t,fraction)=t| between scaled types. 337 338If the result would have magnitude $2^{31}$ or more, |make_frac| 339sets |arith_error:=true|. Most of \MP's internal computations have 340been designed to avoid this sort of error. 341 342If this subroutine were programmed in assembly language on a typical 343machine, we could simply compute |(@t$2^{28}$@>*p)div q|, since a 344double-precision product can often be input to a fixed-point division 345instruction. But when we are restricted to PASCAL arithmetic it 346is necessary either to resort to multiple-precision maneuvering 347or to use a simple but slow iteration. The multiple-precision technique 348would be about three times faster than the code adopted here, but it 349would be comparatively long and tricky, involving about sixteen 350additional multiplications and divisions. 351 352This operation is part of \MP's ``inner loop''; indeed, it will 353consume nearly 10\%! of the running time (exclusive of input and output) 354if the code below is left unchanged. A machine-dependent recoding 355will therefore make \MP\ run faster. The present implementation 356is highly portable, but slow; it avoids multiplication and division 357except in the initial stage. System wizards should be careful to 358replace it with a routine that is guaranteed to produce identical 359results in all cases. 360@^system dependencies@> 361 362As noted below, a few more routines should also be replaced by machine-dependent 363code, for efficiency. But when a procedure is not part of the ``inner loop,'' 364such changes aren't advisable; simplicity and robustness are 365preferable to trickery, unless the cost is too high. 366 367@c 368static int make_frac(int p, int q) 369{ 370 int f; /* the fraction bits, with a leading 1 bit */ 371 int n; /* the integer part of $\vert p/q\vert$ */ 372 register int be_careful; /* disables certain compiler optimizations */ 373 boolean negative = false; /* should the result be negated? */ 374 if (p < 0) { 375 negate(p); 376 negative = true; 377 } 378 if (q <= 0) { 379#ifdef DEBUG 380 if (q == 0) 381 confusion("/"); 382#endif 383 negate(q); 384 negative = !negative; 385 } 386 n = p / q; 387 p = p % q; 388 if (n >= 8) { 389 arith_error = true; 390 if (negative) 391 return (-el_gordo); 392 else 393 return el_gordo; 394 } else { 395 n = (n - 1) * fraction_one; 396 /* Compute $f=\lfloor 2^{28}(1+p/q)+{1\over2}\rfloor$ */ 397 /* The |repeat| loop here preserves the following invariant relations 398 between |f|, |p|, and~|q|: 399 (i)~|0<=p<q|; (ii)~$fq+p=2^k(q+p_0)$, where $k$ is an integer and 400 $p_0$ is the original value of~$p$. 401 402 Notice that the computation specifies 403 |(p-q)+p| instead of |(p+p)-q|, because the latter could overflow. 404 Let us hope that optimizing compilers do not miss this point; a 405 special variable |be_careful| is used to emphasize the necessary 406 order of computation. Optimizing compilers should keep |be_careful| 407 in a register, not store it in memory. 408 */ 409 f = 1; 410 do { 411 be_careful = p - q; 412 p = be_careful + p; 413 if (p >= 0) 414 f = f + f + 1; 415 else { 416 f += f; 417 p = p + q; 418 } 419 } while (f < fraction_one); 420 be_careful = p - q; 421 if (be_careful + p >= 0) 422 incr(f); 423 424 if (negative) 425 return (-(f + n)); 426 else 427 return (f + n); 428 } 429} 430 431@ @c 432static int take_frac(int q, int f) 433{ 434 int p; /* the fraction so far */ 435 int n; /* additional multiple of $q$ */ 436 register int be_careful; /* disables certain compiler optimizations */ 437 boolean negative = false; /* should the result be negated? */ 438 /* Reduce to the case that |f>=0| and |q>0| */ 439 if (f < 0) { 440 negate(f); 441 negative = true; 442 } 443 if (q < 0) { 444 negate(q); 445 negative = !negative; 446 } 447 448 if (f < fraction_one) { 449 n = 0; 450 } else { 451 n = f / fraction_one; 452 f = f % fraction_one; 453 if (q <= el_gordo / n) { 454 n = n * q; 455 } else { 456 arith_error = true; 457 n = el_gordo; 458 } 459 } 460 f = f + fraction_one; 461 /* Compute $p=\lfloor qf/2^{28}+{1\over2}\rfloor-q$ */ 462 /* The invariant relations in this case are (i)~$\lfloor(qf+p)/2^k\rfloor 463 =\lfloor qf_0/2^{28}+{1\over2}\rfloor$, where $k$ is an integer and 464 $f_0$ is the original value of~$f$; (ii)~$2^k\L f<2^{k+1}$. 465 */ 466 p = fraction_half; /* that's $2^{27}$; the invariants hold now with $k=28$ */ 467 if (q < fraction_four) { 468 do { 469 if (odd(f)) 470 p = halfp(p + q); 471 else 472 p = halfp(p); 473 f = halfp(f); 474 } while (f != 1); 475 } else { 476 do { 477 if (odd(f)) 478 p = p + halfp(q - p); 479 else 480 p = halfp(p); 481 f = halfp(f); 482 } while (f != 1); 483 } 484 485 be_careful = n - el_gordo; 486 if (be_careful + p > 0) { 487 arith_error = true; 488 n = el_gordo - p; 489 } 490 if (negative) 491 return (-(n + p)); 492 else 493 return (n + p); 494} 495 496 497 498@ The subroutines for logarithm and exponential involve two tables. 499The first is simple: |two_to_the[k]| equals $2^k$. The second involves 500a bit more calculation, which the author claims to have done correctly: 501|spec_log[k]| is $2^{27}$ times $\ln\bigl(1/(1-2^{-k})\bigr)= 5022^{-k}+{1\over2}2^{-2k}+{1\over3}2^{-3k}+\cdots\,$, rounded to the 503nearest integer. 504 505@c 506static int two_to_the[31]; /* powers of two */ 507static int spec_log[29]; /* special logarithms */ 508 509@ @c 510void initialize_arithmetic(void) 511{ 512 int k; 513 two_to_the[0] = 1; 514 for (k = 1; k <= 30; k++) 515 two_to_the[k] = 2 * two_to_the[k - 1]; 516 spec_log[1] = 93032640; 517 spec_log[2] = 38612034; 518 spec_log[3] = 17922280; 519 spec_log[4] = 8662214; 520 spec_log[5] = 4261238; 521 spec_log[6] = 2113709; 522 spec_log[7] = 1052693; 523 spec_log[8] = 525315; 524 spec_log[9] = 262400; 525 spec_log[10] = 131136; 526 spec_log[11] = 65552; 527 spec_log[12] = 32772; 528 spec_log[13] = 16385; 529 for (k = 14; k <= 27; k++) 530 spec_log[k] = two_to_the[27 - k]; 531 spec_log[28] = 1; 532} 533 534@ @c 535static int m_log(int x) 536{ 537 int y, z; /* auxiliary registers */ 538 int k; /* iteration counter */ 539 if (x <= 0) { 540 /* Handle non-positive logarithm */ 541 print_err("Logarithm of "); 542 print_scaled(x); 543 tprint(" has been replaced by 0"); 544 help2("Since I don't take logs of non-positive numbers,", 545 "I'm zeroing this one. Proceed, with fingers crossed."); 546 error(); 547 return 0; 548 } else { 549 y = 1302456956 + 4 - 100; /* $14\times2^{27}\ln2\approx1302456956.421063$ */ 550 z = 27595 + 6553600; /* and $2^{16}\times .421063\approx 27595$ */ 551 while (x < fraction_four) { 552 x += x; 553 y = y - 93032639; 554 z = z - 48782; 555 } /* $2^{27}\ln2\approx 93032639.74436163$ 556 and $2^{16}\times.74436163\approx 48782$ */ 557 y = y + (z / unity); 558 k = 2; 559 while (x > fraction_four + 4) { 560 /* Increase |k| until |x| can be multiplied by a 561 factor of $2^{-k}$, and adjust $y$ accordingly */ 562 z = ((x - 1) / two_to_the[k]) + 1; /* $z=\lceil x/2^k\rceil$ */ 563 while (x < fraction_four + z) { 564 z = halfp(z + 1); 565 k = k + 1; 566 } 567 y = y + spec_log[k]; 568 x = x - z; 569 } 570 return (y / 8); 571 } 572} 573 574 575 576@ The following somewhat different subroutine tests rigorously if $ab$ is 577greater than, equal to, or less than~$cd$, 578given integers $(a,b,c,d)$. In most cases a quick decision is reached. 579The result is $+1$, 0, or~$-1$ in the three respective cases. 580 581@c 582static int ab_vs_cd(int a, int b, int c, int d) 583{ 584 int q, r; /* temporary registers */ 585 /* Reduce to the case that |a,c>=0|, |b,d>0| */ 586 if (a < 0) { 587 negate(a); 588 negate(b); 589 } 590 if (c < 0) { 591 negate(c); 592 negate(d); 593 } 594 if (d <= 0) { 595 if (b >= 0) 596 return (((a == 0 || b == 0) && (c == 0 || d == 0)) ? 0 : 1); 597 if (d == 0) 598 return (a == 0 ? 0 : -1); 599 q = a; 600 a = c; 601 c = q; 602 q = -b; 603 b = -d; 604 d = q; 605 } else if (b <= 0) { 606 if (b < 0 && a > 0) 607 return -1; 608 return (c == 0 ? 0 : -1); 609 } 610 611 while (1) { 612 q = a / d; 613 r = c / b; 614 if (q != r) 615 return (q > r ? 1 : -1); 616 q = a % d; 617 r = c % b; 618 if (r == 0) 619 return (q == 0 ? 0 : 1); 620 if (q == 0) 621 return -1; 622 a = b; 623 b = q; 624 c = d; 625 d = r; /* now |a>d>0| and |c>b>0| */ 626 } 627} 628 629 630 631@ To consume a random integer, the program below will say `|next_random|' 632and then it will fetch |randoms[j_random]|. 633 634@c 635#define next_random() do { \ 636 if (j_random==0) new_randoms(); else decr(j_random); \ 637 } while (0) 638 639static void new_randoms(void) 640{ 641 int k; /* index into |randoms| */ 642 int x; /* accumulator */ 643 for (k = 0; k <= 23; k++) { 644 x = randoms[k] - randoms[k + 31]; 645 if (x < 0) 646 x = x + fraction_one; 647 randoms[k] = x; 648 } 649 for (k = 24; k <= 54; k++) { 650 x = randoms[k] - randoms[k - 24]; 651 if (x < 0) 652 x = x + fraction_one; 653 randoms[k] = x; 654 } 655 j_random = 54; 656} 657 658 659@ To initialize the |randoms| table, we call the following routine. 660 661@c 662void init_randoms(int seed) 663{ 664 int j, jj, k; /* more or less random integers */ 665 int i; /* index into |randoms| */ 666 j = abs(seed); 667 while (j >= fraction_one) 668 j = halfp(j); 669 k = 1; 670 for (i = 0; i <= 54; i++) { 671 jj = k; 672 k = j - k; 673 j = jj; 674 if (k < 0) 675 k = k + fraction_one; 676 randoms[(i * 21) % 55] = j; 677 } 678 new_randoms(); 679 new_randoms(); 680 new_randoms(); /* ``warm up'' the array */ 681} 682 683 684@ To produce a uniform random number in the range |0<=u<x| or |0>=u>x| 685or |0=u=x|, given a |scaled| value~|x|, we proceed as shown here. 686 687Note that the call of |take_frac| will produce the values 0 and~|x| 688with about half the probability that it will produce any other particular 689values between 0 and~|x|, because it rounds its answers. 690 691@c 692int unif_rand(int x) 693{ 694 int y; /* trial value */ 695 next_random(); 696 y = take_frac(abs(x), randoms[j_random]); 697 if (y == abs(x)) 698 return 0; 699 else if (x > 0) 700 return y; 701 else 702 return -y; 703} 704 705 706@ Finally, a normal deviate with mean zero and unit standard deviation 707can readily be obtained with the ratio method (Algorithm 3.4.1R in 708{\sl The Art of Computer Programming\/}). 709 710@c 711int norm_rand(void) 712{ 713 int x, u, l; /* what the book would call $2^{16}X$, $2^{28}U$, and $-2^{24}\ln U$ */ 714 do { 715 do { 716 next_random(); 717 x = take_frac(112429, randoms[j_random] - fraction_half); 718 /* $2^{16}\sqrt{8/e}\approx 112428.82793$ */ 719 next_random(); 720 u = randoms[j_random]; 721 } while (abs(x) >= u); 722 x = make_frac(x, u); 723 l = 139548960 - m_log(u); /* $2^{24}\cdot12\ln2\approx139548959.6165$ */ 724 } while (ab_vs_cd(1024, l, x, x) < 0); 725 return x; 726} 727 728@ This function could also be expressed as a macro, but it is a useful 729 breakpoint for debugging. 730 731@c 732int fix_int(int val, int min, int max) 733{ 734 return (val < min ? min : (val > max ? max : val)); 735} 736