1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
|
/*
* mf_arith.c -- optimized arithmetic routines for UNIX METAFONT 84
*
* These functions account for a significant percentage of the processing
* time used by METAFONT, and have been recoded in assembly language for
* some classes of machines.
*
* As yet I don't have code for SPARC, so I'm just using a C version
* by richards, in hopes that cc does better than pc on the brute-force
* algorithm found in mf.web.
*
* function makefraction(p, q: integer): fraction;
* { Calculate the function floor( (p * 2^28) / q + 0.5 ) }
* { if both p and q are positive. If not, then the value }
* { of makefraction is calculated as though both *were* }
* { positive, then the result sign adjusted. }
* { (e.g. makefraction ALWAYS rounds away from zero) }
* { In case of an overflow, return the largest possible }
* { value (2^32-1) with the correct sign, and set global }
* { variable "aritherror" to 1. Note that -2^32 is }
* { considered to be an illegal product for this type of }
* { arithmetic! }
*
* function takefraction(f: fraction; q: integer): integer;
* { Calculate the function floor( (p * q) / 2^28 + 0.5 ) }
* { takefraction() rounds in a similar manner as }
* { makefraction() above. }
*
*/
#include "h00vars.h"
extern bool aritherror; /* to be set on overflow */
#define EL_GORDO 0x7fffffff /* 2^31-1 */
#define FRACTION_ONE 0x10000000 /* 2^28 */
#define FRACTION_HALF 0x08000000 /* 2^27 */
#define FRACTION_FOUR 0x40000000 /* 2^30 */
int makefraction(p, q)
register int p, q;
{
int negative;
register int be_careful;
register int f;
int n;
if (p < 0) {
negative = 1;
p = -p;
} else negative = 0;
if (q < 0) {
negative = !negative;
q = -q;
}
n = p / q;
p = p % q;
if (n >= 8) {
aritherror = TRUE;
return (negative? -EL_GORDO : EL_GORDO);
}
n = (n-1)*FRACTION_ONE;
f = 1;
do {
be_careful = p - q;
p = be_careful + p;
if (p >= 0) {
f = f + f + 1;
} else {
f <<= 1;
p += q;
}
} while (f < FRACTION_ONE);
be_careful = p - q;
if ((be_careful + p) >= 0)
f += 1;
return (negative? -(f+n) : (f+n));
}
int takefraction(q, f)
register int q;
register int f;
{
int n, negative;
register int p, be_careful;
if (f < 0) {
negative = 1;
f = -f;
} else negative = 0;
if (q < 0) {
negative = !negative;
q = -q;
}
if (f < FRACTION_ONE)
n = 0;
else {
n = f / FRACTION_ONE;
f = f % FRACTION_ONE;
if (q < (EL_GORDO / n))
n = n * q;
else {
aritherror = TRUE;
n = EL_GORDO;
}
}
f += FRACTION_ONE;
p = FRACTION_HALF;
if (q < FRACTION_FOUR)
do {
if (f & 0x1)
p = (p+q) >> 1;
else
p >>= 1;
f >>= 1;
} while (f > 1);
else
do {
if (f & 0x1)
p = p + ((q-p) >> 1);
else
p >>= 1;
f >>= 1;
} while (f > 1);
be_careful = n - EL_GORDO;
if ((be_careful + p) > 0) {
aritherror = TRUE;
n = EL_GORDO - p;
}
return(negative? -(n+p) : (n+p));
}
|
