rts-int128.c
1 //! @file rts-int128.c
2 //! @author J. Marcel van der Veer
3
4 //! @section Copyright
5 //!
6 //! This file is part of Algol68G - an Algol 68 compiler-interpreter.
7 //! Copyright 2001-2024 J. Marcel van der Veer [algol68g@xs4all.nl].
8
9 //! @section License
10 //!
11 //! This program is free software; you can redistribute it and/or modify it
12 //! under the terms of the GNU General Public License as published by the
13 //! Free Software Foundation; either version 3 of the License, or
14 //! (at your option) any later version.
15 //!
16 //! This program is distributed in the hope that it will be useful, but
17 //! WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
18 //! or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for
19 //! more details. You should have received a copy of the GNU General Public
20 //! License along with this program. If not, see [http://www.gnu.org/licenses/].
21
22 //! @section Synopsis
23 //!
24 //! 128-bit INT support.
25
26 // A68G keeps this code since some (for instance 32-bit) platforms do not
27 // natively support a 128-bit int type.
28
29 #include "a68g.h"
30
31 #if (A68_LEVEL >= 3)
32
33 #include "a68g-prelude.h"
34 #include "a68g-genie.h"
35 #include "a68g-double.h"
36
37 void m64to128 (DOUBLE_NUM_T * w, UNSIGNED_T u, UNSIGNED_T v)
38 {
39 // Knuth Algorithm M, multiprecision multiplication of natural numbers.
40 #define M (0xffffffff)
41 #define N 32
42 UNSIGNED_T hu = u >> N, lu = u & M, hv = v >> N, lv = v & M;
43 UNSIGNED_T t = lu * lv;
44 UNSIGNED_T w3 = t & M, k = t >> N;
45 t = hu * lv + k;
46 UNSIGNED_T w2 = t & M, w1 = t >> N;
47 t = lu * hv + w2;
48 k = t >> N;
49 HW (*w) = hu * hv + w1 + k;
50 LW (*w) = (t << N) + w3;
51 #undef M
52 #undef N
53 }
54
55 void m128to128 (NODE_T * p, MOID_T * m, DOUBLE_NUM_T * w, DOUBLE_NUM_T u, DOUBLE_NUM_T v)
56 {
57 // Knuth Algorithm M, multiprecision multiplication of natural numbers.
58 UNSIGNED_T hu = HW (u), lu = LW (u), hv = HW (v), lv = LW (v);
59 DOUBLE_NUM_T k, t, h, w1, w2, w3;
60 if (lu == 0 || lv == 0) {
61 set_lw (t, 0);
62 } else {
63 m64to128 (&t, lu, lv);
64 }
65 set_lw (w3, LW (t));
66 set_lw (k, HW (t));
67 if (hu == 0 || lv == 0) {
68 set_lw (t, 0);
69 } else {
70 m64to128 (&t, hu, lv);
71 }
72 add_double (p, m, t, t, k);
73 set_lw (w2, LW (t));
74 set_lw (w1, HW (t));
75 if (lu == 0 || hv == 0) {
76 set_lw (t, 0);
77 } else {
78 m64to128 (&t, lu, hv);
79 }
80 add_double (p, m, t, t, w2);
81 set_lw (k, HW (t));
82 if (hu == 0 || hv == 0) {
83 set_lw (h, 0);
84 } else {
85 m64to128 (&h, hu, hv);
86 }
87 add_double (p, m, h, h, w1);
88 add_double (p, m, h, h, k);
89 set_hw (*w, LW (t));
90 add_double (p, m, *w, *w, w3);
91 PRELUDE_ERROR (MODCHK (p, m, HW (h) != 0 || LW (h) != 0), p, ERROR_MATH, M_LONG_INT)
92 }
93
94 DOUBLE_NUM_T double_udiv (NODE_T * p, MOID_T * m, DOUBLE_NUM_T n, DOUBLE_NUM_T d, int mode)
95 {
96 // A bit naive long division.
97 DOUBLE_NUM_T q, r;
98 // Special cases.
99 PRELUDE_ERROR (IS_ZERO (d), p, ERROR_DIVISION_BY_ZERO, M_LONG_INT);
100 if (IS_ZERO (n)) {
101 if (mode == 0) {
102 set_lw (q, 0);
103 return q;
104 } else {
105 set_lw (r, 0);
106 return r;
107 }
108 }
109 // Assuming random n and d, circa half of the divisions is trivial:
110 // n / d = 1 when n = d, n / d = 0 when n < d.
111 if (EQ (n, d)) {
112 if (mode == 0) {
113 set_lw (q, 1);
114 return q;
115 } else {
116 set_lw (r, 0);
117 return r;
118 }
119 } else if (GT (d, n)) {
120 if (mode == 0) {
121 set_lw (q, 0);
122 return q;
123 } else {
124 return n;
125 }
126 }
127 // Halfword divide.
128 if (HW (n) == 0 && HW (d) == 0) {
129 if (mode == 0) {
130 set_lw (q, LW (n) / LW (d));
131 return q;
132 } else {
133 set_lw (r, LW (n) % LW (d));
134 return r;
135 }
136 }
137 // We now know that n and d both have > 64 bits.
138 // Full divide.
139 set_lw (q, 0);
140 set_lw (r, 0);
141 for (int k = 128; k > 0; k--) {
142 UNSIGNED_T carry = (LW (q) & D_SIGN) ? 0x1 : 0x0;
143 LW (q) <<= 1;
144 HW (q) = (HW (q) << 1) | carry;
145 // Left-shift r
146 carry = (LW (r) & D_SIGN) ? 0x1 : 0x0;
147 LW (r) <<= 1;
148 HW (r) = (HW (r) << 1) | carry;
149 // r[0] = n[k]
150 if (HW (n) & D_SIGN) {
151 LW (r) |= 0x1;
152 }
153 carry = (LW (n) & D_SIGN) ? 0x1 : 0x0;
154 LW (n) <<= 1;
155 HW (n) = (HW (n) << 1) | carry;
156 // if r >= d
157 if (GE (r, d)) {
158 // r = r - d
159 sub_double (p, m, r, r, d);
160 // q[k] = 1
161 LW (q) |= 0x1;
162 }
163 }
164 if (mode == 0) {
165 return q;
166 } else {
167 return r;
168 }
169 }
170
171 DOUBLE_NUM_T double_uadd (NODE_T * p, MOID_T * m, DOUBLE_NUM_T u, DOUBLE_NUM_T v)
172 {
173 DOUBLE_NUM_T w;
174 (void) p;
175 add_double (p, m, w, u, v);
176 return w;
177 }
178
179 DOUBLE_NUM_T double_usub (NODE_T * p, MOID_T * m, DOUBLE_NUM_T u, DOUBLE_NUM_T v)
180 {
181 DOUBLE_NUM_T w;
182 (void) p;
183 sub_double (p, m, w, u, v);
184 return w;
185 }
186
187 DOUBLE_NUM_T double_umul (NODE_T * p, MOID_T * m, DOUBLE_NUM_T u, DOUBLE_NUM_T v)
188 {
189 DOUBLE_NUM_T w;
190 m128to128 (p, m, &w, u, v);
191 return w;
192 }
193
194 // Signed integer.
195
196 DOUBLE_NUM_T double_sadd (NODE_T * p, DOUBLE_NUM_T u, DOUBLE_NUM_T v)
197 {
198 int neg_u = D_NEG (u), neg_v = D_NEG (v);
199 if (neg_u) {
200 u = neg_double_int (u);
201 }
202 if (neg_v) {
203 v = neg_double_int (v);
204 }
205 DOUBLE_NUM_T w;
206 set_lw (w, 0);
207 if (!neg_u && !neg_v) {
208 w = double_uadd (p, M_LONG_INT, u, v);
209 PRELUDE_ERROR (D_NEG (w), p, ERROR_MATH, M_LONG_INT);
210 } else if (neg_u && neg_v) {
211 w = neg_double_int (double_sadd (p, u, v));
212 } else if (neg_u) {
213 w = double_ssub (p, v, u);
214 } else if (neg_v) {
215 w = double_ssub (p, u, v);
216 }
217 return w;
218 }
219
220 DOUBLE_NUM_T double_ssub (NODE_T * p, DOUBLE_NUM_T u, DOUBLE_NUM_T v)
221 {
222 int neg_u = D_NEG (u), neg_v = D_NEG (v);
223 if (neg_u) {
224 u = neg_double_int (u);
225 }
226 if (neg_v) {
227 v = neg_double_int (v);
228 }
229 DOUBLE_NUM_T w;
230 set_lw (w, 0);
231 if (!neg_u && !neg_v) {
232 if (D_LT (u, v)) {
233 w = neg_double_int (double_usub (p, M_LONG_INT, v, u));
234 } else {
235 w = double_usub (p, M_LONG_INT, u, v);
236 }
237 } else if (neg_u && neg_v) {
238 w = double_ssub (p, v, u);
239 } else if (neg_u) {
240 w = neg_double_int (double_sadd (p, u, v));
241 } else if (neg_v) {
242 w = double_sadd (p, u, v);
243 }
244 return w;
245 }
246
247 DOUBLE_NUM_T double_smul (NODE_T * p, DOUBLE_NUM_T u, DOUBLE_NUM_T v)
248 {
249 int neg_u = D_NEG (u), neg_v = D_NEG (v);
250 DOUBLE_NUM_T w;
251 if (neg_u) {
252 u = neg_double_int (u);
253 }
254 if (neg_v) {
255 v = neg_double_int (v);
256 }
257 w = double_umul (p, M_LONG_INT, u, v);
258 if (neg_u != neg_v) {
259 w = neg_double_int (w);
260 }
261 return w;
262 }
263
264 DOUBLE_NUM_T double_sdiv (NODE_T * p, DOUBLE_NUM_T u, DOUBLE_NUM_T v, int mode)
265 {
266 int neg_u = D_NEG (u), neg_v = D_NEG (v);
267 if (neg_u) {
268 u = neg_double_int (u);
269 }
270 if (neg_v) {
271 v = neg_double_int (v);
272 }
273 DOUBLE_NUM_T w = double_udiv (p, M_LONG_INT, u, v, mode);
274 if (mode == 0 && neg_u != neg_v) {
275 w = neg_double_int (w);
276 } else if (mode == 1 && D_NEG (w)) {
277 w = double_sadd (p, w, v);
278 }
279 return w;
280 }
281
282 #endif
