3 % Copyright
2009-2010 Taco Hoekwater
<taco@@luatex.org
>
5 % This file is part of LuaTeX.
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.
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.
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
/>.
27 @ The principal computations performed by \TeX\ are done entirely in terms of
28 integers less than $
2^
{31}$ in magnitude
; and divisions are done only when both
29 dividend and divisor are nonnegative. Thus
, the arithmetic specified in this
30 program can be carried out in exactly the same way on a wide variety of
31 computers
, including some small ones. Why? Because the arithmetic
32 calculations need to be spelled out precisely in order to guarantee that
33 \TeX\ will produce identical output on different machines. If some
34 quantities were rounded differently in different implementations
, we would
35 find that line breaks and even page breaks might occur in different places.
36 Hence the arithmetic of \TeX\ has been designed with care
, and systems that
37 claim to be implementations of \TeX82 should follow precisely the
39 calculations as they appear in the present program.
41 (Actually there are three places where \TeX\ uses |div| with a possibly negative
42 numerator. These are harmless
; see |div| in the index. Also if the user
43 sets the \.
{\\time
} or the \.
{\\year
} to a negative value
, some diagnostic
44 information will involve negative-numerator division. The same remarks
45 apply for |mod| as well as for |div|.
)
47 Here is a routine that calculates half of an integer
, using an
48 unambiguous convention with respect to signed odd numbers.
60 @ The following function is used to create a scaled integer from a given decimal
61 fraction $
(.d_0d_1\ldots d_
{k-1
})$
, where |
0<=k
<=17|. The digit $d_i$ is
62 given in |dig
[i
]|
, and the calculation produces a correctly rounded result.
65 scaled round_decimals
(int k
)
66 { /* converts a decimal fraction
*/
67 int a
; /* the accumulator
*/
70 a
= (a
+ dig
[k
] * two
) / 10;
76 @ Conversely
, here is a procedure analogous to |print_int|. If the output
77 of this procedure is subsequently read by \TeX\ and converted by the
78 |round_decimals| routine above
, it turns out that the original value will
79 be reproduced exactly
; the ``simplest'' such decimal number is output
,
80 but there is always at least one digit following the decimal point.
82 The invariant relation in the \
&{repeat} loop is that a sequence of
83 decimal digits yet to be printed will yield the original number if and only if
84 they form a fraction~$f$ in the range $s-\delta\L10\cdot2^
{16}f
<s$.
85 We can stop if and only if $f
=0$ satisfies this condition
; the loop will
86 terminate before $s$ can possibly become zero.
89 void print_scaled
(scaled s
)
90 { /* prints scaled real
, rounded to five digits
*/
91 scaled delta
; /* amount of allowable inaccuracy
*/
96 negate
(s
); /* print the sign
, if negative
*/
98 print_int
(s
/ unity
); /* print the integer part
*/
100 s
= 10 * (s
% unity
) + 5;
104 s
= s
+ 0100000 - 50000; /* round the last digit
*/
105 buffer
[i
++] = '
0'
+ (s
/ unity
);
106 s
= 10 * (s
% unity
);
113 @ Physical sizes that a \TeX\ user specifies for portions of documents are
114 represented internally as scaled points. Thus
, if we define an `sp'
(scaled
116 point
) as a unit equal to $
2^
{-16}$ printer's points
, every dimension
117 inside of \TeX\ is an integer number of sp. There are exactly
118 4,736,286.72 sp per inch. Users are not allowed to specify dimensions
119 larger than $
2^
{30}-1$ sp
, which is a distance of about
18.892 feet
(5.7583
120 meters
); two such quantities can be added without overflow on a
32-bit
123 The present implementation of \TeX\ does not check for overflow when
124 @^overflow in arithmetic@
>
125 dimensions are added or subtracted. This could be done by inserting a
126 few dozen tests of the form `\ignorespaces|if x
>=010000000000 then
127 @t\\
{report\_overflow
}@
>|'
, but the chance of overflow is so remote that
128 such tests do not seem worthwhile.
130 \TeX\ needs to do only a few arithmetic operations on scaled quantities
,
131 other than addition and subtraction
, and the following subroutines do most of
132 the work. A single computation might use several subroutine calls
, and it is
133 desirable to avoid producing multiple error messages in case of arithmetic
134 overflow
; so the routines set the global variable |arith_error| to |true|
135 instead of reporting errors directly to the user. Another global variable
,
136 |tex_remainder|
, holds the remainder after a division.
139 boolean arith_error
; /* has arithmetic overflow occurred recently?
*/
140 scaled tex_remainder
; /* amount subtracted to get an exact division
*/
143 @ The first arithmetical subroutine we need computes $nx
+y$
, where |x|
144 and~|y| are |scaled| and |n| is an integer. We will also use it to
148 scaled mult_and_add
(int n
, scaled x
, scaled y
, scaled max_answer
)
156 if
(((x
<= (max_answer
- y
) / n
) && (-x <= (max_answer + y) / n))) {
164 @ We also need to divide scaled dimensions by integers.
166 scaled x_over_n
(scaled x
, int n
)
168 boolean negative
; /* should |tex_remainder| be negated?
*/
181 tex_remainder
= x
% n
;
183 negate
(tex_remainder
);
186 tex_remainder
= -((-x
) % n
);
188 negate
(tex_remainder
);
189 return
(-((-x
) / n
));
195 @ Then comes the multiplication of a scaled number by a fraction |n
/d|
,
196 where |n| and |d| are nonnegative integers |
<=@t$
2^
{16}$@
>| and |d| is
197 positive. It would be too dangerous to multiply by~|n| and then divide
198 by~|d|
, in separate operations
, since overflow might well occur
; and it
199 would be too inaccurate to divide by |d| and then multiply by |n|. Hence
200 this subroutine simulates
1.5-precision arithmetic.
203 scaled xn_over_d
(scaled x
, int n
, int d
)
205 nonnegative_integer t
, u
, v
, xx
, dd
; /* intermediate quantities
*/
206 boolean positive
= true
; /* was |x
>=0|?
*/
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)
219 u
= 0100000 * (u
/ dd
) + (v
/ dd
);
221 tex_remainder
= (int
) (v
% dd
);
224 /* casts are for ms cl
*/
225 tex_remainder
= -(int
) (v
% dd
);
226 return
-(scaled
) (u
);
231 @ The next subroutine is used to compute the ``badness'' of glue
, when a
232 total~|t| is supposed to be made from amounts that sum to~|s|. According
233 to
{\sl The \TeX book
}, the badness of this situation is $
100(t
/s
)^
3$
;
234 however
, badness is simply a heuristic
, so we need not squeeze out the
235 last drop of accuracy when computing it. All we really want is an
236 approximation that has similar properties.
237 @
:TeXbook
}{\sl The \TeX book@
>
239 The actual method used to compute the badness is easier to read from the
240 program than to describe in words. It produces an integer value that is a
241 reasonably close approximation to $
100(t
/s
)^
3$
, and all implementations
242 of \TeX\ should use precisely this method. Any badness of $
2^
{13}$ or more is
243 treated as infinitely bad
, and represented by
10000.
245 It 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
247 computing at most
1095 distinct values
, but that is plenty.
250 halfword badness
(scaled t
, scaled s
)
251 { /* compute badness
, given |t
>=0|
*/
252 int r
; /* approximation to $\alpha t
/s$
, where $\alpha^
3\approx
260 r
= (t
* 297) / s
; /* $
297^
3=99.94\times2^
{18}$
*/
261 else if
(s
>= 1663497)
266 return inf_bad
; /* $
1290^
3<2^
{31}<1291^
3$
*/
268 return
((r
* r
* r
+ 0400000) / 01000000);
269 /* that was $r^
3/2^
{18}$
, rounded to the nearest integer
*/
274 @ When \TeX\ ``packages'' a list into a box
, it needs to calculate the
275 proportionality ratio by which the glue inside the box should stretch
276 or shrink. This calculation does not affect \TeX's decision making
,
277 so the precise details of rounding
, etc.
, in the glue calculation are not
278 of critical importance for the consistency of results on different computers.
280 We shall use the type |glue_ratio| for such proportionality ratios.
281 A 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
283 other data structures. Thus |glue_ratio| should be equivalent to
284 |short_real| in some implementations of
PASCAL. Alternatively
,
285 it is possible to deal with glue ratios using nothing but fixed-point
286 arithmetic
; see
{\sl TUGboat \bf3
},1 (March
1982), 10--27.
(But the
287 routines cited there must be modified to allow negative glue ratios.
)
288 @^system dependencies@
>
291 @ This section is
(almost
) straight from MetaPost. I had to change
292 the types
(use |integer| instead of |fraction|
), but that should
293 not have any influence on the actual calculations
(the original
294 comments refer to quantities like |fraction_four|
($
2^
{30}$
), and
295 that is the same as the numeric representation of |max_dimen|
).
297 I've copied the low-level variables and routines that are needed
, but
298 only those
(e.g. |m_log|
), not the accompanying ones like |m_exp|. Most
299 of the following low-level numeric routines are only needed within the
300 calculation of |norm_rand|. I've been forced to rename |make_fraction|
301 to |make_frac| because TeX already has a routine by that name with
302 a wholly different function
(it creates a |fraction_noad| for math
305 And now let's complete our collection of numeric utility routines
306 by considering random number generation.
307 \MP
{} generates pseudo-random numbers with the additive scheme recommended
308 in Section
3.6 of
{\sl The Art of Computer Programming
}; however
, the
309 results are random fractions between
0 and |fraction_one-1|
, inclusive.
311 There's an auxiliary array |randoms| that contains
55 pseudo-random
312 fractions. Using the recurrence $x_n
=(x_
{n-55
}-x_
{n-31
})\bmod
2^
{28}$
,
313 we generate batches of
55 new $x_n$'s at a time by calling |new_randoms|.
314 The global variable |j_random| tells which element has most recently
318 static int randoms
[55]; /* the last
55 random values generated
*/
319 static int j_random
; /* the number of unused |randoms|
*/
320 scaled random_seed
; /* the default random seed
*/
322 @ A small bit of metafont is needed.
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
*/
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
333 positive. If |p| and |q| are both of the same scaled type |t|
,
334 the ``type relation'' |make_frac
(t
,t
)=fraction| is valid
;
335 and it's also possible to use the subroutine ``backwards
,'' using
336 the relation |make_frac
(t
,fraction
)=t| between scaled types.
338 If the result would have magnitude $
2^
{31}$ or more
, |make_frac|
339 sets |arith_error
:=true|. Most of \MP's internal computations have
340 been designed to avoid this sort of error.
342 If this subroutine were programmed in assembly language on a typical
343 machine
, we could simply compute |
(@t$
2^
{28}$@
>*p
)div q|
, since a
344 double-precision product can often be input to a fixed-point division
345 instruction. But when we are restricted to
PASCAL arithmetic it
346 is necessary either to resort to multiple-precision maneuvering
347 or to use a simple but slow iteration. The multiple-precision technique
348 would be about three times faster than the code adopted here
, but it
349 would be comparatively long and tricky
, involving about sixteen
350 additional multiplications and divisions.
352 This operation is part of \MP's ``inner loop''
; indeed
, it will
353 consume nearly
10\
%! of the running time
(exclusive of input and output
)
354 if the code below is left unchanged. A machine-dependent recoding
355 will therefore make \MP\ run faster. The present implementation
356 is highly portable
, but slow
; it avoids multiplication and division
357 except in the initial stage. System wizards should be careful to
358 replace it with a routine that is guaranteed to produce identical
359 results in all cases.
360 @^system dependencies@
>
362 As noted below
, a few more routines should also be replaced by machine-dependent
363 code
, for efficiency. But when a procedure is not part of the ``inner loop
,''
364 such changes aren't advisable
; simplicity and robustness are
365 preferable to trickery
, unless the cost is too high.
368 static int make_frac
(int p
, int q
)
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?
*/
384 negative
= !negative
;
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$.
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.
419 } while
(f
< fraction_one
);
421 if
(be_careful
+ p
>= 0)
432 static int take_frac
(int q
, int f
)
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|
*/
445 negative
= !negative
;
448 if
(f
< fraction_one
) {
451 n
= f
/ fraction_one
;
452 f
= f
% fraction_one
;
453 if
(q
<= el_gordo
/ n
) {
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}$.
466 p
= fraction_half
; /* that's $
2^
{27}$
; the invariants hold now with $k
=28$
*/
467 if
(q
< fraction_four
) {
478 p
= p
+ halfp
(q
- p
);
485 be_careful
= n
- el_gordo
;
486 if
(be_careful
+ p
> 0) {
498 @ The subroutines for logarithm and exponential involve two tables.
499 The first is simple
: |two_to_the
[k
]| equals $
2^k$. The second involves
500 a 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
)=
502 2^
{-k
}+{1\over2
}2^
{-2k
}+{1\over3
}2^
{-3k
}+\cdots\
,$
, rounded to the
506 static int two_to_the
[31]; /* powers of two
*/
507 static int spec_log
[29]; /* special logarithms
*/
510 void initialize_arithmetic
(void
)
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
];
535 static int m_log
(int x
)
537 int y
, z
; /* auxiliary registers
*/
538 int k
; /* iteration counter
*/
540 /* Handle non-positive logarithm
*/
541 print_err
("Logarithm of ");
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.");
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
) {
555 } /* $
2^
{27}\ln2\approx
93032639.74436163$
556 and $
2^
{16}\times
.74436163\approx
48782$
*/
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
) {
576 @ The following somewhat different subroutine tests rigorously if $ab$ is
577 greater than
, equal to
, or less than~$cd$
,
578 given integers $
(a
,b
,c
,d
)$. In most cases a quick decision is reached.
579 The result is $
+1$
, 0, or~$
-1$ in the three respective cases.
582 static int ab_vs_cd
(int a
, int b
, int c
, int d
)
584 int q
, r
; /* temporary registers
*/
585 /* Reduce to the case that |a
,c
>=0|
, |b
,d
>0|
*/
596 return
(((a
== 0 || b
== 0) && (c == 0 || d == 0)) ? 0 : 1);
598 return
(a
== 0 ?
0 : -1);
608 return
(c
== 0 ?
0 : -1);
615 return
(q
> r ?
1 : -1);
619 return
(q
== 0 ?
0 : 1);
625 d
= r
; /* now |a
>d
>0| and |c
>b
>0|
*/
631 @ To consume a random integer
, the program below will say `|next_random|'
632 and then it will fetch |randoms
[j_random
]|.
635 #define next_random
() do
{ \
636 if
(j_random
==0) new_randoms
(); else decr
(j_random
); \
639 static void new_randoms
(void
)
641 int k
; /* index into |randoms|
*/
642 int x
; /* accumulator
*/
643 for
(k
= 0; k
<= 23; k
++) {
644 x
= randoms
[k
] - randoms
[k
+ 31];
646 x
= x
+ fraction_one
;
649 for
(k
= 24; k
<= 54; k
++) {
650 x
= randoms
[k
] - randoms
[k
- 24];
652 x
= x
+ fraction_one
;
659 @ To initialize the |randoms| table
, we call the following routine.
662 void init_randoms
(int seed
)
664 int j
, jj
, k
; /* more or less random integers
*/
665 int i
; /* index into |randoms|
*/
667 while
(j
>= fraction_one
)
670 for
(i
= 0; i
<= 54; i
++) {
675 k
= k
+ fraction_one
;
676 randoms
[(i
* 21) % 55] = j
;
680 new_randoms
(); /* ``warm up'' the array
*/
684 @ To produce a uniform random number in the range |
0<=u
<x| or |
0>=u
>x|
685 or |
0=u
=x|
, given a |scaled| value~|x|
, we proceed as shown here.
687 Note that the call of |take_frac| will produce the values
0 and~|x|
688 with about half the probability that it will produce any other particular
689 values between
0 and~|x|
, because it rounds its answers.
694 int y
; /* trial value
*/
696 y
= take_frac
(abs
(x
), randoms
[j_random
]);
706 @ Finally
, a normal deviate with mean zero and unit standard deviation
707 can readily be obtained with the ratio method
(Algorithm
3.4.1R in
708 {\sl The Art of Computer Programming\
/}).
713 int x
, u
, l
; /* what the book would call $
2^
{16}X$
, $
2^
{28}U$
, and $
-2^
{24}\ln U$
*/
717 x
= take_frac
(112429, randoms
[j_random
] - fraction_half
);
718 /* $
2^
{16}\sqrt
{8/e
}\approx
112428.82793$
*/
720 u
= randoms
[j_random
];
721 } while
(abs
(x
) >= u
);
723 l
= 139548960 - m_log
(u
); /* $
2^
{24}\cdot12\ln2\approx139548959.6165$
*/
724 } while
(ab_vs_cd
(1024, l
, x
, x
) < 0);
728 @ This function could also be expressed as a macro
, but it is a useful
729 breakpoint for debugging.
732 int fix_int
(int val
, int min
, int max
)
734 return
(val
< min ? min
: (val
> max ? max
: val
));