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
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
|
/***********************************************************************/
/* */
/* OCaml */
/* */
/* Xavier Leroy, projet Cristal, INRIA Rocquencourt */
/* */
/* Copyright 2003 Institut National de Recherche en Informatique et */
/* en Automatique. All rights reserved. This file is distributed */
/* under the terms of the GNU Library General Public License, with */
/* the special exception on linking described in file ../../LICENSE. */
/* */
/***********************************************************************/
/**** Generic operations on digits ****/
/* These macros can be defined in the machine-specific include file.
Below are the default definitions (in plain C).
Except for BngMult, all macros are guaranteed to evaluate their
arguments exactly once. */
#ifndef BngAdd2
/* res = arg1 + arg2. carryout = carry out. */
#define BngAdd2(res,carryout,arg1,arg2) { \
bngdigit tmp1, tmp2; \
tmp1 = arg1; \
tmp2 = tmp1 + (arg2); \
carryout = (tmp2 < tmp1); \
res = tmp2; \
}
#endif
#ifndef BngAdd2Carry
/* res = arg1 + arg2 + carryin. carryout = carry out. */
#define BngAdd2Carry(res,carryout,arg1,arg2,carryin) { \
bngdigit tmp1, tmp2, tmp3; \
tmp1 = arg1; \
tmp2 = tmp1 + (arg2); \
tmp3 = tmp2 + (carryin); \
carryout = (tmp2 < tmp1) + (tmp3 < tmp2); \
res = tmp3; \
}
#endif
#ifndef BngAdd3
/* res = arg1 + arg2 + arg3. Each carry increments carryaccu. */
#define BngAdd3(res,carryaccu,arg1,arg2,arg3) { \
bngdigit tmp1, tmp2, tmp3; \
tmp1 = arg1; \
tmp2 = tmp1 + (arg2); \
carryaccu += (tmp2 < tmp1); \
tmp3 = tmp2 + (arg3); \
carryaccu += (tmp3 < tmp2); \
res = tmp3; \
}
#endif
#ifndef BngSub2
/* res = arg1 - arg2. carryout = carry out. */
#define BngSub2(res,carryout,arg1,arg2) { \
bngdigit tmp1, tmp2; \
tmp1 = arg1; \
tmp2 = arg2; \
res = tmp1 - tmp2; \
carryout = (tmp1 < tmp2); \
}
#endif
#ifndef BngSub2Carry
/* res = arg1 - arg2 - carryin. carryout = carry out. */
#define BngSub2Carry(res,carryout,arg1,arg2,carryin) { \
bngdigit tmp1, tmp2, tmp3; \
tmp1 = arg1; \
tmp2 = arg2; \
tmp3 = tmp1 - tmp2; \
res = tmp3 - (carryin); \
carryout = (tmp1 < tmp2) + (tmp3 < carryin); \
}
#endif
#ifndef BngSub3
/* res = arg1 - arg2 - arg3. Each carry increments carryaccu. */
#define BngSub3(res,carryaccu,arg1,arg2,arg3) { \
bngdigit tmp1, tmp2, tmp3, tmp4; \
tmp1 = arg1; \
tmp2 = arg2; \
tmp3 = arg3; \
tmp4 = tmp1 - tmp2; \
res = tmp4 - tmp3; \
carryaccu += (tmp1 < tmp2) + (tmp4 < tmp3); \
}
#endif
#define BngLowHalf(d) ((d) & (((bngdigit)1 << BNG_BITS_PER_HALF_DIGIT) - 1))
#define BngHighHalf(d) ((d) >> BNG_BITS_PER_HALF_DIGIT)
#ifndef BngMult
/* resl = low digit of product arg1 * arg2
resh = high digit of product arg1 * arg2. */
#if SIZEOF_PTR == 4 && defined(ARCH_UINT64_TYPE)
#define BngMult(resh,resl,arg1,arg2) { \
ARCH_UINT64_TYPE p = (ARCH_UINT64_TYPE)(arg1) * (ARCH_UINT64_TYPE)(arg2); \
resh = p >> 32; \
resl = p; \
}
#else
#define BngMult(resh,resl,arg1,arg2) { \
bngdigit p11 = BngLowHalf(arg1) * BngLowHalf(arg2); \
bngdigit p12 = BngLowHalf(arg1) * BngHighHalf(arg2); \
bngdigit p21 = BngHighHalf(arg1) * BngLowHalf(arg2); \
bngdigit p22 = BngHighHalf(arg1) * BngHighHalf(arg2); \
resh = p22 + (p12 >> BNG_BITS_PER_HALF_DIGIT) \
+ (p21 >> BNG_BITS_PER_HALF_DIGIT); \
BngAdd3(resl, resh, \
p11, p12 << BNG_BITS_PER_HALF_DIGIT, p21 << BNG_BITS_PER_HALF_DIGIT); \
}
#endif
#endif
#ifndef BngDiv
/* Divide the double-width number nh:nl by d.
Require d != 0 and nh < d.
Store quotient in quo, remainder in rem.
Can be slow if d is not normalized. */
#define BngDiv(quo,rem,nh,nl,d) bng_div_aux(&(quo),&(rem),nh,nl,d)
#define BngDivNeedsNormalization
static void bng_div_aux(bngdigit * quo, bngdigit * rem,
bngdigit nh, bngdigit nl, bngdigit d)
{
bngdigit dl, dh, ql, qh, pl, ph, nsaved;
dl = BngLowHalf(d);
dh = BngHighHalf(d);
/* Under-estimate the top half of the quotient (qh) */
qh = nh / (dh + 1);
/* Shift nh:nl right by BNG_BITS_PER_HALF_DIGIT bits,
so that we focus on the top 1.5 digits of the numerator.
Then, subtract (qh * d) from nh:nl. */
nsaved = BngLowHalf(nl);
ph = qh * dh;
pl = qh * dl;
nh -= ph; /* Subtract before shifting so that carry propagates for free */
nl = (nl >> BNG_BITS_PER_HALF_DIGIT) | (nh << BNG_BITS_PER_HALF_DIGIT);
nh = (nh >> BNG_BITS_PER_HALF_DIGIT);
nh -= (nl < pl); /* Borrow */
nl -= pl;
/* Adjust estimate qh until nh:nl < 0:d */
while (nh != 0 || nl >= d) {
nh -= (nl < d); /* Borrow */
nl -= d;
qh++;
}
/* Under-estimate the bottom half of the quotient (ql) */
ql = nl / (dh + 1);
/* Shift nh:nl left by BNG_BITS_PER_HALF_DIGIT bits, restoring the
low bits we saved earlier, so that we focus on the bottom 1.5 digit
of the numerator. Then, subtract (ql * d) from nh:nl. */
ph = ql * dh;
pl = ql * dl;
nl -= ph; /* Subtract before shifting so that carry propagates for free */
nh = (nl >> BNG_BITS_PER_HALF_DIGIT);
nl = (nl << BNG_BITS_PER_HALF_DIGIT) | nsaved;
nh -= (nl < pl); /* Borrow */
nl -= pl;
/* Adjust estimate ql until nh:nl < 0:d */
while (nh != 0 || nl >= d) {
nh -= (nl < d); /* Borrow */
nl -= d;
ql++;
}
/* We're done */
*quo = (qh << BNG_BITS_PER_HALF_DIGIT) | ql;
*rem = nl;
}
#endif
|
