Section 99: Arithmetic with scaled dimensions
The principal computations performed by are done entirely in terms of integers less than 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 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 has been designed with care, and systems that claim to be implementations of should follow precisely the calculations as they appear in the present program.
(Actually there are three places where 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.)
Section 100
Here is a routine that calculates half of an integer, using an unambiguous convention with respect to signed odd numbers.
// << Start file |arithmetic.c|, 1382 >>
int half(int x) {
int res;
if (odd(x)) {
res = (x + 1) / 2;
}
else {
res = x / 2;
}
return res;
}
Section 101
Fixed-point arithmetic is done on scaled integers that are multiples of . In other words, a binary point is assumed to be sixteen bit positions from the right end of a binary computer word.
#define UNITY 65536 // 2^{16}, represents 1.00000
#define TWO 131072 // 2^{17}, represents 2.00000
⟨ Types in the outer block 18 ⟩+≡
typedef int scaled; // this is used for scaled integer
typedef unsigned char small_number; // this type is self-explanatory
Section 102
The following function is used to create a scaled integer from a given decimal fraction , where 0 k 17. The digit is given in dig[i], and the calculation produces a correctly rounded result.
// converts a decimal fraction
scaled round_decimals(small_number k) {
int a = 0; // the accumulator
while (k > 0) {
decr(k);
a = (a + dig[k] * TWO) / 10;
}
return (a + 1) / 2;
}
Section 103
Conversely, here is a procedure analogous to print_int. If the output of this procedure is subsequently read by 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 in the range . We can stop if and only if satisfies this condition; the loop will terminate before can possibly become zero.
// prints scaled real, rounded to five digits
void print_scaled(scaled s) {
scaled delta; // amount of allowable inaccuracy
if (s < 0) {
// print the sign, if negative
print_char('-');
negate(s);
}
print_int(s / UNITY); // print the integer part
print_char('.');
s = (s % UNITY)*10 + 5;
delta = 10;
do {
if (delta > UNITY) {
s += 32768 - 50000; // round the last digit
}
print_char('0' + (s / UNITY));
s = 10 * (s % UNITY);
delta *= 10;
} while (s > delta);
}
Section 104
Physical sizes that a user specifies for portions of documents are represented internally as scaled points. Thus, if we define an ‘sp’ (scaled point) as a unit equal to printer’s points, every dimension inside of 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 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 does not check for overflow when dimensions are added or subtracted. This could be done by inserting a few dozen tests of the form ‘if x 0x40000000 then report_overflow’, but the chance of overflow is so remote that such tests do not seem worthwhile.
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, remainder, holds the remainder after a division.
The global variable remainder has been renamed remainder_ to avoid conflict with a function of the same name in math library.
⟨ Global variables 13 ⟩+≡
bool arith_error; // has arithmetic overflow occurred recently?
scaled remainder_; // amount subtracted to get an exact division
Section 105
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.
// << Start file |arithmetic.h|, 1381 >>
#define nx_plus_y(A, B, C) mult_and_add((A), (B), (C), 0x3fffffff)
#define mult_integers(A, B) mult_and_add((A), (B), 0, 0x7fffffff)
scaled mult_and_add(int n, scaled x, scaled y, scaled max_answer) {
scaled res;
if (n < 0) {
negate(x);
negate(n);
}
if (n == 0) {
res = y;
}
else if ((x <= (max_answer - y) / n)
&& (-x <= (max_answer + y) / n))
{
res = n * x + y;
}
else {
arith_error = true;
res = 0;
}
return res;
}
Section 106
We also need to divide scaled dimensions by integers.
scaled x_over_n(scaled x, int n) {
bool negative = false; // should |remainder_| be negated?
scaled quotient;
if (n == 0) {
arith_error = true;
quotient = 0;
remainder_ = x;
}
else {
if (n < 0) {
negate(x);
negate(n);
negative = true;
}
if (x >= 0) {
quotient = x / n;
remainder_ = x % n;
}
else {
quotient = -((-x) / n);
remainder_ = -((-x) % n);
}
}
if (negative) {
negate(remainder_);
}
return quotient;
}
Section 107
Then comes the multiplication of a scaled number by a fraction nd, where n and d are nonnegative integers 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) {
bool positive = true; // was |x >= 0|?
int t, u, v; // intermediate quantities
scaled quotient;
if (x < 0) {
negate(x);
positive = false;
}
t = (x % 32768) * n;
u = (x / 32768) * n + (t / 32768);
v = (u % d) * 32768 + (t % 32768);
if (u / d >= 32768) {
arith_error = true;
}
else {
u = 32768 * (u / d) + (v / d);
}
if (positive) {
quotient = u;
remainder_ = v % d;
}
else {
quotient = -u;
remainder_ = -(v % d);
}
return quotient;
}
Section 108
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 The TeXbook, the badness of this situation is ; 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.
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 , and all implementations of should use precisely this method. Any badness of or more is treated as infinitely bad, and represented by 10000.
It is not difficult to prove that
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.
#define INF_BAD 10000 // infinitely bad value
// compute badness, given |t >= 0|
halfword badness(scaled t, scaled s) {
int r; // approximation to \alpha*t/s, where alpha^3\approx 100*2^{18}
halfword badness;
if (t == 0) {
badness = 0;
}
else if (s <= 0) {
badness = INF_BAD;
}
else {
if (t <= 7230584) {
r = (t * 297) / s; // 297^3 = 99.94 * 2^{18}
}
else if (s >= 1663497) {
r = t / (s / 297);
}
else {
r = t;
}
if (r > 1290) {
badness = INF_BAD; // 1290^3 < 2^{31} < 1291^3
}
else {
badness = (r*r*r + 0x20000) / 0x40000;
}
} // that was r^3/2^{18}, rounded to the nearest integer
return badness;
}
Section 109
When “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 ’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 int (usually 32 bits) if it is to blend smoothly with ’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 TUGboat 3,1 (March 1982), 10–27. (But the routines cited there must be modified to allow negative glue ratios.)
float and unfloat are unused. float_constant is directly used by appending
.0
to a value. Finally, glue_ration is made into a double.
#define set_glue_ratio_zero(X) (X) = 0.0 // store the representation of zero ratio
#define set_glue_ratio_one(X) (X) = 1.0 // store the representation of unit ratio
⟨ Types in the outer block 18 ⟩+≡
typedef double glue_ratio; // one-word representation of a glue expansion factor