1 /**************************************************************************/
2 /*                                                                        */
3 /*                                 OCaml                                  */
4 /*                                                                        */
5 /*             Xavier Leroy, projet Cristal, INRIA Rocquencourt           */
6 /*                                                                        */
7 /*   Copyright 2003 Institut National de Recherche en Informatique et     */
8 /*     en Automatique.                                                    */
9 /*                                                                        */
10 /*   All rights reserved.  This file is distributed under the terms of    */
11 /*   the GNU Lesser General Public License version 2.1, with the          */
12 /*   special exception on linking described in the file LICENSE.          */
13 /*                                                                        */
14 /**************************************************************************/
15 
16 /**** Generic operations on digits ****/
17 
18 /* These macros can be defined in the machine-specific include file.
19    Below are the default definitions (in plain C).
20    Except for BngMult, all macros are guaranteed to evaluate their
21    arguments exactly once. */
22 
23 #ifndef BngAdd2
24 /* res = arg1 + arg2.  carryout = carry out. */
25 #define BngAdd2(res,carryout,arg1,arg2) {                                   \
26   bngdigit tmp1, tmp2;                                                      \
27   tmp1 = arg1;                                                              \
28   tmp2 = tmp1 + (arg2);                                                     \
29   carryout = (tmp2 < tmp1);                                                 \
30   res = tmp2;                                                               \
31 }
32 #endif
33 
34 #ifndef BngAdd2Carry
35 /* res = arg1 + arg2 + carryin.  carryout = carry out. */
36 #define BngAdd2Carry(res,carryout,arg1,arg2,carryin) {                      \
37   bngdigit tmp1, tmp2, tmp3;                                                \
38   tmp1 = arg1;                                                              \
39   tmp2 = tmp1 + (arg2);                                                     \
40   tmp3 = tmp2 + (carryin);                                                  \
41   carryout = (tmp2 < tmp1) + (tmp3 < tmp2);                                 \
42   res = tmp3;                                                               \
43 }
44 #endif
45 
46 #ifndef BngAdd3
47 /* res = arg1 + arg2 + arg3.  Each carry increments carryaccu. */
48 #define BngAdd3(res,carryaccu,arg1,arg2,arg3) {                             \
49   bngdigit tmp1, tmp2, tmp3;                                                \
50   tmp1 = arg1;                                                              \
51   tmp2 = tmp1 + (arg2);                                                     \
52   carryaccu += (tmp2 < tmp1);                                               \
53   tmp3 = tmp2 + (arg3);                                                     \
54   carryaccu += (tmp3 < tmp2);                                               \
55   res = tmp3;                                                               \
56 }
57 #endif
58 
59 #ifndef BngSub2
60 /* res = arg1 - arg2.  carryout = carry out. */
61 #define BngSub2(res,carryout,arg1,arg2) {                                   \
62   bngdigit tmp1, tmp2;                                                      \
63   tmp1 = arg1;                                                              \
64   tmp2 = arg2;                                                              \
65   res = tmp1 - tmp2;                                                        \
66   carryout = (tmp1 < tmp2);                                                 \
67 }
68 #endif
69 
70 #ifndef BngSub2Carry
71 /* res = arg1 - arg2 - carryin.  carryout = carry out. */
72 #define BngSub2Carry(res,carryout,arg1,arg2,carryin) {                      \
73   bngdigit tmp1, tmp2, tmp3;                                                \
74   tmp1 = arg1;                                                              \
75   tmp2 = arg2;                                                              \
76   tmp3 = tmp1 - tmp2;                                                       \
77   res = tmp3 - (carryin);                                                   \
78   carryout = (tmp1 < tmp2) + (tmp3 < carryin);                              \
79 }
80 #endif
81 
82 #ifndef BngSub3
83 /* res = arg1 - arg2 - arg3.  Each carry increments carryaccu. */
84 #define BngSub3(res,carryaccu,arg1,arg2,arg3) {                             \
85   bngdigit tmp1, tmp2, tmp3, tmp4;                                          \
86   tmp1 = arg1;                                                              \
87   tmp2 = arg2;                                                              \
88   tmp3 = arg3;                                                              \
89   tmp4 = tmp1 - tmp2;                                                       \
90   res = tmp4 - tmp3;                                                        \
91   carryaccu += (tmp1 < tmp2) + (tmp4 < tmp3);                               \
92 }
93 #endif
94 
95 #define BngLowHalf(d) ((d) & (((bngdigit)1 << BNG_BITS_PER_HALF_DIGIT) - 1))
96 #define BngHighHalf(d) ((d) >> BNG_BITS_PER_HALF_DIGIT)
97 
98 #ifndef BngMult
99 /* resl = low  digit of product arg1 * arg2
100    resh = high digit of product arg1 * arg2. */
101 #if SIZEOF_PTR == 4 && defined(ARCH_UINT64_TYPE)
102 #define BngMult(resh,resl,arg1,arg2) {                                      \
103   ARCH_UINT64_TYPE p = (ARCH_UINT64_TYPE)(arg1) * (ARCH_UINT64_TYPE)(arg2); \
104   resh = p >> 32;                                                           \
105   resl = p;                                                                 \
106 }
107 #else
108 #define BngMult(resh,resl,arg1,arg2) {                                      \
109   bngdigit p11 = BngLowHalf(arg1) * BngLowHalf(arg2);                       \
110   bngdigit p12 = BngLowHalf(arg1) * BngHighHalf(arg2);                      \
111   bngdigit p21 = BngHighHalf(arg1) * BngLowHalf(arg2);                      \
112   bngdigit p22 = BngHighHalf(arg1) * BngHighHalf(arg2);                     \
113   resh = p22 + (p12 >> BNG_BITS_PER_HALF_DIGIT)                             \
114              + (p21 >> BNG_BITS_PER_HALF_DIGIT);                            \
115   BngAdd3(resl, resh,                                                       \
116      p11, p12 << BNG_BITS_PER_HALF_DIGIT, p21 << BNG_BITS_PER_HALF_DIGIT);  \
117 }
118 #endif
119 #endif
120 
121 #ifndef BngDiv
122 /* Divide the double-width number nh:nl by d.
123    Require d != 0 and nh < d.
124    Store quotient in quo, remainder in rem.
125    Can be slow if d is not normalized. */
126 #define BngDiv(quo,rem,nh,nl,d) bng_div_aux(&(quo),&(rem),nh,nl,d)
127 #define BngDivNeedsNormalization
128 
bng_div_aux(bngdigit * quo,bngdigit * rem,bngdigit nh,bngdigit nl,bngdigit d)129 static void bng_div_aux(bngdigit * quo, bngdigit * rem,
130                         bngdigit nh, bngdigit nl, bngdigit d)
131 {
132   bngdigit dl, dh, ql, qh, pl, ph, nsaved;
133 
134   dl = BngLowHalf(d);
135   dh = BngHighHalf(d);
136   /* Under-estimate the top half of the quotient (qh) */
137   qh = nh / (dh + 1);
138   /* Shift nh:nl right by BNG_BITS_PER_HALF_DIGIT bits,
139      so that we focus on the top 1.5 digits of the numerator.
140      Then, subtract (qh * d) from nh:nl. */
141   nsaved = BngLowHalf(nl);
142   ph = qh * dh;
143   pl = qh * dl;
144   nh -= ph; /* Subtract before shifting so that carry propagates for free */
145   nl = (nl >> BNG_BITS_PER_HALF_DIGIT) | (nh << BNG_BITS_PER_HALF_DIGIT);
146   nh = (nh >> BNG_BITS_PER_HALF_DIGIT);
147   nh -= (nl < pl);  /* Borrow */
148   nl -= pl;
149   /* Adjust estimate qh until nh:nl < 0:d */
150   while (nh != 0 || nl >= d) {
151     nh -= (nl < d); /* Borrow */
152     nl -= d;
153     qh++;
154   }
155   /* Under-estimate the bottom half of the quotient (ql) */
156   ql = nl / (dh + 1);
157   /* Shift nh:nl left by BNG_BITS_PER_HALF_DIGIT bits, restoring the
158      low bits we saved earlier, so that we focus on the bottom 1.5 digit
159      of the numerator.  Then, subtract (ql * d) from nh:nl. */
160   ph = ql * dh;
161   pl = ql * dl;
162   nl -= ph; /* Subtract before shifting so that carry propagates for free */
163   nh = (nl >> BNG_BITS_PER_HALF_DIGIT);
164   nl = (nl << BNG_BITS_PER_HALF_DIGIT) | nsaved;
165   nh -= (nl < pl);  /* Borrow */
166   nl -= pl;
167   /* Adjust estimate ql until nh:nl < 0:d */
168   while (nh != 0 || nl >= d) {
169     nh -= (nl < d); /* Borrow */
170     nl -= d;
171     ql++;
172   }
173   /* We're done */
174   *quo = (qh << BNG_BITS_PER_HALF_DIGIT) | ql;
175   *rem = nl;
176 }
177 
178 #endif
179