diff options
Diffstat (limited to 'Build/source/texk/web2c/luatexdir/tex/arithmetic.c')
-rw-r--r-- | Build/source/texk/web2c/luatexdir/tex/arithmetic.c | 815 |
1 files changed, 815 insertions, 0 deletions
diff --git a/Build/source/texk/web2c/luatexdir/tex/arithmetic.c b/Build/source/texk/web2c/luatexdir/tex/arithmetic.c new file mode 100644 index 00000000000..2782c253400 --- /dev/null +++ b/Build/source/texk/web2c/luatexdir/tex/arithmetic.c @@ -0,0 +1,815 @@ +/* + +arithmetic.w + +Copyright 2009-2010 Taco Hoekwater <taco@@luatex.org> + +This file is part of LuaTeX. + +LuaTeX is free software; you can redistribute it and/or modify it under +the terms of the GNU General Public License as published by the Free +Software Foundation; either version 2 of the License, or (at your +option) any later version. + +LuaTeX is distributed in the hope that it will be useful, but WITHOUT +ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or +FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public +License for more details. + +You should have received a copy of the GNU General Public License along +with LuaTeX; if not, see <http://www.gnu.org/licenses/>. + +*/ + +#include "ptexlib.h" + +/*tex + +The principal computations performed by \TeX\ are done entirely in terms of +integers less than $2^{31}$ in magnitude; and divisions are done only when both +dividend and divisor are nonnegative. Thus, the arithmetic specified in this +program can be carried out in exactly the same way on a wide variety of +computers, including some small ones. Why? Because the arithmetic calculations +need to be spelled out precisely in order to guarantee that \TeX\ will produce +identical output on different machines. If some quantities were rounded +differently in different implementations, we would find that line breaks and even +page breaks might occur in different places. Hence the arithmetic of \TeX\ has +been designed with care, and systems that claim to be implementations of \TeX82 +should follow precisely the @:TeX82}{\TeX82@> calculations as they appear in the +present program. + +Actually there are three places where \TeX\ uses |div| with a possibly negative +numerator. These are harmless; see |div| in the index. Also if the user sets the +\.{\\time} or the \.{\\year} to a negative value, some diagnostic information +will involve negative-numerator division. The same remarks apply for |mod| as +well as for |div|. + +Here is a routine that calculates half of an integer, using an unambiguous +convention with respect to signed odd numbers. + +*/ + +int half(int x) +{ + if (odd(x)) + return ((x + 1) / 2); + else + return (x / 2); +} + +/*tex + +The following function is used to create a scaled integer from a given decimal +fraction $(.d_0d_1\ldots d_{k-1})$, where |0<=k<=17|. The digit $d_i$ is +given in |dig[i]|, and the calculation produces a correctly rounded result. + +*/ + +scaled round_decimals(int k) +{ + int a = 0; + while (k-- > 0) { + a = (a + dig[k] * two) / 10; + } + return ((a + 1) / 2); +} + +/*tex + +Conversely, here is a procedure analogous to |print_int|. If the output of this +procedure is subsequently read by \TeX\ and converted by the |round_decimals| +routine above, it turns out that the original value will be reproduced exactly; +the ``simplest'' such decimal number is output, but there is always at least one +digit following the decimal point. + +The invariant relation in the \&{repeat} loop is that a sequence of decimal +digits yet to be printed will yield the original number if and only if they form +a fraction~$f$ in the range $s-\delta\L10\cdot2^{16}f<s$. We can stop if and only +if $f=0$ satisfies this condition; the loop will terminate before $s$ can +possibly become zero. + +The next one prints a scaled real, rounded to five digits. + +*/ + +void print_scaled(scaled s) +{ + /*tex The amount of allowable inaccuracy: */ + scaled delta; + char buffer[20]; + int i = 0; + if (s < 0) { + /*tex Print the sign, if negative. */ + print_char('-'); + negate(s); + } + /*tex Print the integer part. */ + print_int(s / unity); + buffer[i++] = '.'; + s = 10 * (s % unity) + 5; + delta = 10; + do { + if (delta > unity) { + /*tex Round the last digit. */ + s = s + 0100000 - 50000; + } + buffer[i++] = '0' + (s / unity); + s = 10 * (s % unity); + delta = delta * 10; + } while (s > delta); + buffer[i++] = '\0'; + tprint(buffer); +} + +/*tex + +Physical sizes that a \TeX\ user specifies for portions of documents are +represented internally as scaled points. Thus, if we define an `sp' (scaled +@^sp@> point) as a unit equal to $2^{-16}$ printer's points, every dimension +inside of \TeX\ is an integer number of sp. There are exactly 4,736,286.72 sp per +inch. Users are not allowed to specify dimensions larger than $2^{30}-1$ sp, +which is a distance of about 18.892 feet (5.7583 meters); two such quantities can +be added without overflow on a 32-bit computer. + +The present implementation of \TeX\ does not check for overflow when @^overflow +in arithmetic@> dimensions are added or subtracted. This could be done by +inserting a few dozen tests of the form `\ignorespaces|if x>=010000000000 then +@t\\{report\_overflow}@>|', but the chance of overflow is so remote that such +tests do not seem worthwhile. + +\TeX\ needs to do only a few arithmetic operations on scaled quantities, other +than addition and subtraction, and the following subroutines do most of the work. +A single computation might use several subroutine calls, and it is desirable to +avoid producing multiple error messages in case of arithmetic overflow; so the +routines set the global variable |arith_error| to |true| instead of reporting +errors directly to the user. Another global variable, |tex_remainder|, holds the +remainder after a division. + +*/ + +/*tex Has arithmetic overflow occurred recently? */ + +boolean arith_error; + +/*tex The amount subtracted to get an exact division. */ + +scaled tex_remainder; + +/*tex + + The first arithmetical subroutine we need computes $nx+y$, where |x| +and~|y| are |scaled| and |n| is an integer. We will also use it to +multiply integers. + +*/ + +scaled mult_and_add(int n, scaled x, scaled y, scaled max_answer) +{ + if (n == 0) + return y; + if (n < 0) { + negate(x); + negate(n); + } + if (((x <= (max_answer - y) / n) && (-x <= (max_answer + y) / n))) { + return (n * x + y); + } else { + arith_error = true; + return 0; + } +} + +/*tex + +We also need to divide scaled dimensions by integers. + +*/ + +scaled x_over_n(scaled x, int n) +{ + /*tex Should |tex_remainder| be negated? */ + boolean negative = false; + if (n == 0) { + arith_error = true; + tex_remainder = x; + return 0; + } else { + if (n < 0) { + negate(x); + negate(n); + negative = true; + } + if (x >= 0) { + tex_remainder = x % n; + if (negative) + negate(tex_remainder); + return (x / n); + } else { + tex_remainder = -((-x) % n); + if (negative) + negate(tex_remainder); + return (-((-x) / n)); + } + } +} + +/*tex + +Then comes the multiplication of a scaled number by a fraction |n/d|, where |n| +and |d| are nonnegative integers |<=@t$2^{16}$@>| and |d| is positive. It would +be too dangerous to multiply by~|n| and then divide by~|d|, in separate +operations, since overflow might well occur; and it would be too inaccurate to +divide by |d| and then multiply by |n|. Hence this subroutine simulates +1.5-precision arithmetic. + +*/ + +scaled xn_over_d(scaled x, int n, int d) +{ + nonnegative_integer t, u, v, xx, dd; + boolean positive = true; + if (x < 0) { + negate(x); + positive = false; + } + xx = (nonnegative_integer) x; + dd = (nonnegative_integer) d; + t = ((xx % 0100000) * (nonnegative_integer) n); + u = ((xx / 0100000) * (nonnegative_integer) n + (t / 0100000)); + v = (u % dd) * 0100000 + (t % 0100000); + if (u / dd >= 0100000) + arith_error = true; + else + u = 0100000 * (u / dd) + (v / dd); + if (positive) { + tex_remainder = (int) (v % dd); + return (scaled) u; + } else { + /*tex The casts are for ms cl. */ + tex_remainder = -(int) (v % dd); + return -(scaled) (u); + } +} + +/*tex + +The next subroutine is used to compute the ``badness'' of glue, when a total~|t| +is supposed to be made from amounts that sum to~|s|. According to {\sl The \TeX +book}, the badness of this situation is $100(t/s)^3$; however, badness is simply +a heuristic, so we need not squeeze out the last drop of accuracy when computing +it. All we really want is an approximation that has similar properties. +@:TeXbook}{\sl The \TeX book@> + +The actual method used to compute the badness is easier to read from the program +than to describe in words. It produces an integer value that is a reasonably +close approximation to $100(t/s)^3$, and all implementations of \TeX\ should use +precisely this method. Any badness of $2^{13}$ or more is treated as infinitely +bad, and represented by 10000. + +It is not difficult to prove that $$\hbox{|badness(t+1,s)>=badness(t,s) +>= badness(t,s+1)|}.$$ The badness function defined here is capable of computing +at most 1095 distinct values, but that is plenty. + +*/ + +halfword badness(scaled t, scaled s) +{ + /*tex Approximation to $\alpha t/s$, where $\alpha^3\approx 100\cdot2^{18}$ */ + int r; + if (t == 0) { + return 0; + } else if (s <= 0) { + return inf_bad; + } else { + /*tex $297^3=99.94\times2^{18}$ */ + if (t <= 7230584) { + r = (t * 297) / s; + } else if (s >= 1663497) { + r = t / (s / 297); + } else { + r = t; + } + if (r > 1290) { + /*tex $1290^3<2^{31}<1291^3$ */ + return inf_bad; + } else { + /*tex This is $r^3/2^{18}$, rounded to the nearest integer. */ + return ((r * r * r + 0400000) / 01000000); + } + } +} + +/*tex + +When \TeX\ ``packages'' a list into a box, it needs to calculate the +proportionality ratio by which the glue inside the box should stretch or shrink. +This calculation does not affect \TeX's decision making, so the precise details +of rounding, etc., in the glue calculation are not of critical importance for the +consistency of results on different computers. + +We shall use the type |glue_ratio| for such proportionality ratios. A glue ratio +should take the same amount of memory as an |integer| (usually 32 bits) if it is +to blend smoothly with \TeX's other data structures. Thus |glue_ratio| should be +equivalent to |short_real| in some implementations of PASCAL. Alternatively, it +is possible to deal with glue ratios using nothing but fixed-point arithmetic; +see {\sl TUGboat \bf3},1 (March 1982), 10--27. (But the routines cited there must +be modified to allow negative glue ratios.) @^system dependencies@> + +*/ + +/* + +This section is (almost) straight from MetaPost. I (Taco) had to change the types +(use |integer| instead of |fraction|), but that should not have any influence on +the actual calculations (the original comments refer to quantities like +|fraction_four| ($2^{30}$), and that is the same as the numeric representation of +|max_dimen|). + +I've copied the low-level variables and routines that are needed, but only those +(e.g. |m_log|), not the accompanying ones like |m_exp|. Most of the following +low-level numeric routines are only needed within the calculation of |norm_rand|. +I've been forced to rename |make_fraction| to |make_frac| because TeX already has +a routine by that name with a wholly different function (it creates a +|fraction_noad| for math typesetting) + +And now let's complete our collection of numeric utility routines by considering +random number generation. \MP{} generates pseudo-random numbers with the additive +scheme recommended in Section 3.6 of {\sl The Art of Computer Programming}; +however, the results are random fractions between 0 and |fraction_one-1|, +inclusive. + +There's an auxiliary array |randoms| that contains 55 pseudo-random fractions. +Using the recurrence $x_n=(x_{n-55}-x_{n-31})\bmod 2^{28}$, we generate batches +of 55 new $x_n$'s at a time by calling |new_randoms|. The global variable +|j_random| tells which element has most recently been consumed. + +*/ + +/*tex The last 55 random values generated: */ + +static int randoms[55]; + +/*tex The number of unused |randoms|: */ + +static int j_random; + +/*tex The default random seed: */ + +scaled random_seed; + +/*tex A small bit of \METAPOST\ is needed. */ + +#define fraction_half 01000000000 /* $2^{27} $, represents 0.50000000 */ +#define fraction_one 02000000000 /* $2^{28} $, represents 1.00000000 */ +#define fraction_four 010000000000 /* $2^{30} $, represents 4.00000000 */ +#define el_gordo 017777777777 /* $2^{31}-1$, the largest value that \MP\ likes */ + +/*tex + +The |make_frac| routine produces the |fraction| equivalent of |p/q|, given +integers |p| and~|q|; it computes the integer +$f=\lfloor2^{28}p/q+{1\over2}\rfloor$, when $p$ and $q$ are positive. If |p| and +|q| are both of the same scaled type |t|, the ``type relation'' +|make_frac(t,t)=fraction| is valid; and it's also possible to use the subroutine +``backwards,'' using the relation |make_frac(t,fraction)=t| between scaled types. + +If the result would have magnitude $2^{31}$ or more, |make_frac| sets +|arith_error:=true|. Most of \MP's internal computations have been designed to +avoid this sort of error. + +If this subroutine were programmed in assembly language on a typical machine, we +could simply compute |(@t$2^{28}$@>*p)div q|, since a double-precision product +can often be input to a fixed-point division instruction. But when we are +restricted to PASCAL arithmetic it is necessary either to resort to +multiple-precision maneuvering or to use a simple but slow iteration. The +multiple-precision technique would be about three times faster than the code +adopted here, but it would be comparatively long and tricky, involving about +sixteen additional multiplications and divisions. + +This operation is part of \MP's ``inner loop''; indeed, it will consume nearly +10\%! of the running time (exclusive of input and output) if the code below is +left unchanged. A machine-dependent recoding will therefore make \MP\ run faster. +The present implementation is highly portable, but slow; it avoids multiplication +and division except in the initial stage. System wizards should be careful to +replace it with a routine that is guaranteed to produce identical results in all +cases. @^system dependencies@> + +As noted below, a few more routines should also be replaced by machine-dependent +code, for efficiency. But when a procedure is not part of the ``inner loop,'' +such changes aren't advisable; simplicity and robustness are preferable to +trickery, unless the cost is too high. + +*/ + +static int make_frac(int p, int q) +{ + /*tex The fraction bits, with a leading 1 bit: */ + int f; + /*tex The integer part of $\vert p/q\vert$: */ + int n; + /*tex Disables certain compiler optimizations: */ + register int be_careful; + /*tex Should the result be negated? */ + boolean negative = false; + if (p < 0) { + negate(p); + negative = true; + } + if (q <= 0) { + negate(q); + negative = !negative; + } + n = p / q; + p = p % q; + if (n >= 8) { + arith_error = true; + if (negative) + return (-el_gordo); + else + return el_gordo; + } else { + n = (n - 1) * fraction_one; + /*tex_remainder + + Compute $f=\lfloor 2^{28}(1+p/q)+{1\over2}\rfloor$. The |repeat| loop + here preserves the following invariant relations between |f|, |p|, + and~|q|: (i)~|0<=p<q|; (ii)~$fq+p=2^k(q+p_0)$, where $k$ is an + integer and $p_0$ is the original value of~$p$. + + Notice that the computation specifies |(p-q)+p| instead of |(p+p)-q|, + because the latter could overflow. Let us hope that optimizing + compilers do not miss this point; a special variable |be_careful| is + used to emphasize the necessary order of computation. Optimizing + compilers should keep |be_careful| in a register, not store it in + memory. + + */ + f = 1; + do { + be_careful = p - q; + p = be_careful + p; + if (p >= 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); + } +} + +static int take_frac(int q, int f) +{ + /*tex The fraction so far: */ + int p; + /*tex Additional multiple of $q$: */ + int n; + /*tex Disables certain compiler optimizations. */ + register int be_careful; + /*tex Should the result be negated? */ + boolean negative = false; + /*tex 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; + /*tex + + 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}$. + + Here |p| becomes $2^{27}$; the invariants hold now with $k=28$: + + */ + p = fraction_half; + 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); +} + +/*tex + +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. + +*/ + +/*tex The powers of two: */ + +static int two_to_the[31]; + +/*tex Special logarithms: */ + +static int spec_log[29]; + +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; +} + +static int m_log(int x) +{ + /*tex Auxiliary registers: */ + int y, z; + /*tex Iteration counter: */ + int k; + if (x <= 0) { + /*tex 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 { + /*tex $14\times2^{27}\ln2\approx1302456956.421063$ */ + y = 1302456956 + 4 - 100; + /*tex $2^{16}\times .421063\approx 27595$ */ + z = 27595 + 6553600; + while (x < fraction_four) { + x += x; + /*tex $2^{27}\ln2\approx 93032639.74436163$ */ + y = y - 93032639; + /*tex $2^{16}\times.74436163\approx 48782$ */ + z = z - 48782; + } + + y = y + (z / unity); + k = 2; + while (x > fraction_four + 4) { + /*tex + Increase |k| until |x| can be multiplied by a factor of $2^{-k}$, + and adjust $y$ accordingly. Here $z=\lceil x/2^k\rceil$. + */ + z = ((x - 1) / two_to_the[k]) + 1; + while (x < fraction_four + z) { + z = halfp(z + 1); + k = k + 1; + } + y = y + spec_log[k]; + x = x - z; + } + return (y / 8); + } +} + +/*tex + +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. + +*/ + +static int ab_vs_cd(int a, int b, int c, int d) +{ + int q, r; + /*tex Reduce to the case that |a,c>=0| and |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; + /*tex Now |a>d>0| and |c>b>0|. */ + } +} + +/*tex + +To consume a random integer, the program below will say `|next_random|' and then +it will fetch |randoms[j_random]|. + +*/ + +#define next_random() do { \ + if (j_random==0) \ + new_randoms(); \ + else \ + decr(j_random); \ +} while (0) + +static void new_randoms(void) +{ + /*tex The index into |randoms|. */ + int k; + /*tex The accumulator. */ + int x; + 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; +} + +/*tex + +To initialize the |randoms| table, we call the following routine. + +*/ + +void init_randoms(int seed) +{ + /*tex Three more or less random integers. */ + int j, jj, k; + /*tex The index into |randoms|. */ + int i; + 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; + } + /*tex We ``warm up'' the array. */ + new_randoms(); + new_randoms(); + new_randoms(); +} + +/*tex + +To produce a uniform random number in the range |0<=u<x| or |0>=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. + +*/ + +int unif_rand(int x) +{ + int y; + 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; +} + +/*tex + +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\/}. + +*/ + +int norm_rand(void) +{ + /*tex What the book would call $2^{16}X$, $2^{28}U$, and $-2^{24}\ln U$. */ + int x, u, l; + do { + do { + next_random(); + x = take_frac(112429, randoms[j_random] - fraction_half); + /*tex Which is $2^{16}\sqrt{8/e}\approx 112428.82793$. */ + next_random(); + u = randoms[j_random]; + } while (abs(x) >= u); + x = make_frac(x, u); + /*tex More fuzzyness: $2^{24}\cdot12\ln2\approx139548959.6165$. */ + l = 139548960 - m_log(u); + } while (ab_vs_cd(1024, l, x, x) < 0); + return x; +} + +/*tex + +This function could also be expressed as a macro, but it is a useful breakpoint +for debugging. + +*/ + +int fix_int(int val, int min, int max) +{ + return (val < min ? min : (val > max ? max : val)); +} |