1 /*
2 * m a t h 6 4 . c
3 * Forth Inspired Command Language - 64 bit math support routines
4 * Authors: Michael A. Gauland (gaulandm@mdhost.cse.tek.com)
5 * Larry Hastings (larry@hastings.org)
6 * John Sadler (john_sadler@alum.mit.edu)
7 * Created: 25 January 1998
8 * Rev 2.03: Support for 128 bit DP math. This file really ouught to
9 * be renamed!
10 * $Id: double.c,v 1.2 2010/09/12 15:18:07 asau Exp $
11 */
12 /*
13 * Copyright (c) 1997-2001 John Sadler (john_sadler@alum.mit.edu)
14 * All rights reserved.
15 *
16 * Get the latest Ficl release at http://ficl.sourceforge.net
17 *
18 * I am interested in hearing from anyone who uses Ficl. If you have
19 * a problem, a success story, a defect, an enhancement request, or
20 * if you would like to contribute to the Ficl release, please
21 * contact me by email at the address above.
22 *
23 * L I C E N S E and D I S C L A I M E R
24 *
25 * Redistribution and use in source and binary forms, with or without
26 * modification, are permitted provided that the following conditions
27 * are met:
28 * 1. Redistributions of source code must retain the above copyright
29 * notice, this list of conditions and the following disclaimer.
30 * 2. Redistributions in binary form must reproduce the above copyright
31 * notice, this list of conditions and the following disclaimer in the
32 * documentation and/or other materials provided with the distribution.
33 *
34 * THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND
35 * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
36 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
37 * ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE
38 * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
39 * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
40 * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
41 * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
42 * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
43 * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
44 * SUCH DAMAGE.
45 */
46
47 #include "ficl.h"
48
49 #if FICL_PLATFORM_HAS_2INTEGER
50 ficl2UnsignedQR
ficl2UnsignedDivide(ficl2Unsigned q,ficlUnsigned y)51 ficl2UnsignedDivide(ficl2Unsigned q, ficlUnsigned y)
52 {
53 ficl2UnsignedQR result;
54
55 result.quotient = q / y;
56 /*
57 * Once we have the quotient, it's cheaper to calculate the
58 * remainder this way than with % (mod). --lch
59 */
60 result.remainder = (ficlInteger)(q - (result.quotient * y));
61
62 return (result);
63 }
64
65 #else /* FICL_PLATFORM_HAS_2INTEGER */
66
67 #define FICL_CELL_HIGH_BIT ((uintmax_t)1 << (FICL_BITS_PER_CELL-1))
68 #define UMOD_SHIFT (FICL_BITS_PER_CELL / 2)
69 #define UMOD_MASK ((1L << (FICL_BITS_PER_CELL / 2)) - 1)
70
71 /*
72 * ficl2IntegerIsNegative
73 * Returns TRUE if the specified ficl2Unsigned has its sign bit set.
74 */
75 int
ficl2IntegerIsNegative(ficl2Integer x)76 ficl2IntegerIsNegative(ficl2Integer x)
77 {
78 return (x.high < 0);
79 }
80
81 /*
82 * ficl2IntegerNegate
83 * Negates an ficl2Unsigned by complementing and incrementing.
84 */
85 ficl2Integer
ficl2IntegerNegate(ficl2Integer x)86 ficl2IntegerNegate(ficl2Integer x)
87 {
88 x.high = ~x.high;
89 x.low = ~x.low;
90 x.low ++;
91 if (x.low == 0)
92 x.high++;
93
94 return (x);
95 }
96
97 /*
98 * ficl2UnsignedMultiplyAccumulate
99 * Mixed precision multiply and accumulate primitive for number building.
100 * Multiplies ficl2Unsigned u by ficlUnsigned mul and adds ficlUnsigned add.
101 * Mul is typically the numeric base, and add represents a digit to be
102 * appended to the growing number.
103 * Returns the result of the operation
104 */
105 ficl2Unsigned
ficl2UnsignedMultiplyAccumulate(ficl2Unsigned u,ficlUnsigned mul,ficlUnsigned add)106 ficl2UnsignedMultiplyAccumulate(ficl2Unsigned u, ficlUnsigned mul,
107 ficlUnsigned add)
108 {
109 ficl2Unsigned resultLo = ficl2UnsignedMultiply(u.low, mul);
110 ficl2Unsigned resultHi = ficl2UnsignedMultiply(u.high, mul);
111 resultLo.high += resultHi.low;
112 resultHi.low = resultLo.low + add;
113
114 if (resultHi.low < resultLo.low)
115 resultLo.high++;
116
117 resultLo.low = resultHi.low;
118
119 return (resultLo);
120 }
121
122 /*
123 * ficl2IntegerMultiply
124 * Multiplies a pair of ficlIntegers and returns an ficl2Integer result.
125 */
126 ficl2Integer
ficl2IntegerMultiply(ficlInteger x,ficlInteger y)127 ficl2IntegerMultiply(ficlInteger x, ficlInteger y)
128 {
129 ficl2Unsigned prod;
130 ficl2Integer result;
131 int sign = 1;
132
133 if (x < 0) {
134 sign = -sign;
135 x = -x;
136 }
137
138 if (y < 0) {
139 sign = -sign;
140 y = -y;
141 }
142
143 prod = ficl2UnsignedMultiply(x, y);
144 FICL_2INTEGER_SET(FICL_2UNSIGNED_GET_HIGH(prod),
145 FICL_2UNSIGNED_GET_LOW(prod), result);
146 if (sign > 0)
147 return (result);
148 else
149 return (ficl2IntegerNegate(result));
150 }
151
152 ficl2Integer
ficl2IntegerDecrement(ficl2Integer x)153 ficl2IntegerDecrement(ficl2Integer x)
154 {
155 if (x.low == INTMAX_MIN)
156 x.high--;
157 x.low--;
158
159 return (x);
160 }
161
162 ficl2Unsigned
ficl2UnsignedAdd(ficl2Unsigned x,ficl2Unsigned y)163 ficl2UnsignedAdd(ficl2Unsigned x, ficl2Unsigned y)
164 {
165 ficl2Unsigned result;
166
167 result.high = x.high + y.high;
168 result.low = x.low + y.low;
169
170 if (result.low < y.low)
171 result.high++;
172
173 return (result);
174 }
175
176 /*
177 * ficl2UnsignedMultiply
178 * Contributed by:
179 * Michael A. Gauland gaulandm@mdhost.cse.tek.com
180 */
181 ficl2Unsigned
ficl2UnsignedMultiply(ficlUnsigned x,ficlUnsigned y)182 ficl2UnsignedMultiply(ficlUnsigned x, ficlUnsigned y)
183 {
184 ficl2Unsigned result = { 0, 0 };
185 ficl2Unsigned addend;
186
187 addend.low = y;
188 addend.high = 0; /* No sign extension--arguments are unsigned */
189
190 while (x != 0) {
191 if (x & 1) {
192 result = ficl2UnsignedAdd(result, addend);
193 }
194 x >>= 1;
195 addend = ficl2UnsignedArithmeticShiftLeft(addend);
196 }
197 return (result);
198 }
199
200 /*
201 * ficl2UnsignedSubtract
202 */
203 ficl2Unsigned
ficl2UnsignedSubtract(ficl2Unsigned x,ficl2Unsigned y)204 ficl2UnsignedSubtract(ficl2Unsigned x, ficl2Unsigned y)
205 {
206 ficl2Unsigned result;
207
208 result.high = x.high - y.high;
209 result.low = x.low - y.low;
210
211 if (x.low < y.low) {
212 result.high--;
213 }
214
215 return (result);
216 }
217
218 /*
219 * ficl2UnsignedArithmeticShiftLeft
220 * 64 bit left shift
221 */
222 ficl2Unsigned
ficl2UnsignedArithmeticShiftLeft(ficl2Unsigned x)223 ficl2UnsignedArithmeticShiftLeft(ficl2Unsigned x)
224 {
225 ficl2Unsigned result;
226
227 result.high = x.high << 1;
228 if (x.low & FICL_CELL_HIGH_BIT) {
229 result.high++;
230 }
231
232 result.low = x.low << 1;
233
234 return (result);
235 }
236
237 /*
238 * ficl2UnsignedArithmeticShiftRight
239 * 64 bit right shift (unsigned - no sign extend)
240 */
241 ficl2Unsigned
ficl2UnsignedArithmeticShiftRight(ficl2Unsigned x)242 ficl2UnsignedArithmeticShiftRight(ficl2Unsigned x)
243 {
244 ficl2Unsigned result;
245
246 result.low = x.low >> 1;
247 if (x.high & 1) {
248 result.low |= FICL_CELL_HIGH_BIT;
249 }
250
251 result.high = x.high >> 1;
252 return (result);
253 }
254
255 /*
256 * ficl2UnsignedOr
257 * 64 bit bitwise OR
258 */
259 ficl2Unsigned
ficl2UnsignedOr(ficl2Unsigned x,ficl2Unsigned y)260 ficl2UnsignedOr(ficl2Unsigned x, ficl2Unsigned y)
261 {
262 ficl2Unsigned result;
263
264 result.high = x.high | y.high;
265 result.low = x.low | y.low;
266
267 return (result);
268 }
269
270 /*
271 * ficl2UnsignedCompare
272 * Return -1 if x < y; 0 if x==y, and 1 if x > y.
273 */
274 int
ficl2UnsignedCompare(ficl2Unsigned x,ficl2Unsigned y)275 ficl2UnsignedCompare(ficl2Unsigned x, ficl2Unsigned y)
276 {
277 if (x.high > y.high)
278 return (1);
279 if (x.high < y.high)
280 return (-1);
281
282 /* High parts are equal */
283
284 if (x.low > y.low)
285 return (1);
286 else if (x.low < y.low)
287 return (-1);
288
289 return (0);
290 }
291
292 /*
293 * ficl2UnsignedDivide
294 * Portable versions of ficl2Multiply and ficl2Divide in C
295 * Contributed by:
296 * Michael A. Gauland gaulandm@mdhost.cse.tek.com
297 */
298 ficl2UnsignedQR
ficl2UnsignedDivide(ficl2Unsigned q,ficlUnsigned y)299 ficl2UnsignedDivide(ficl2Unsigned q, ficlUnsigned y)
300 {
301 ficl2UnsignedQR result;
302 ficl2Unsigned quotient;
303 ficl2Unsigned subtrahend;
304 ficl2Unsigned mask;
305
306 quotient.low = 0;
307 quotient.high = 0;
308
309 subtrahend.low = y;
310 subtrahend.high = 0;
311
312 mask.low = 1;
313 mask.high = 0;
314
315 while ((ficl2UnsignedCompare(subtrahend, q) < 0) &&
316 (subtrahend.high & FICL_CELL_HIGH_BIT) == 0) {
317 mask = ficl2UnsignedArithmeticShiftLeft(mask);
318 subtrahend = ficl2UnsignedArithmeticShiftLeft(subtrahend);
319 }
320
321 while (mask.low != 0 || mask.high != 0) {
322 if (ficl2UnsignedCompare(subtrahend, q) <= 0) {
323 q = ficl2UnsignedSubtract(q, subtrahend);
324 quotient = ficl2UnsignedOr(quotient, mask);
325 }
326 mask = ficl2UnsignedArithmeticShiftRight(mask);
327 subtrahend = ficl2UnsignedArithmeticShiftRight(subtrahend);
328 }
329
330 result.quotient = quotient;
331 result.remainder = q.low;
332 return (result);
333 }
334 #endif /* !FICL_PLATFORM_HAS_2INTEGER */
335
336 /*
337 * ficl2IntegerDivideFloored
338 *
339 * FROM THE FORTH ANS...
340 * Floored division is integer division in which the remainder carries
341 * the sign of the divisor or is zero, and the quotient is rounded to
342 * its arithmetic floor. Symmetric division is integer division in which
343 * the remainder carries the sign of the dividend or is zero and the
344 * quotient is the mathematical quotient rounded towards zero or
345 * truncated. Examples of each are shown in tables 3.3 and 3.4.
346 *
347 * Table 3.3 - Floored Division Example
348 * Dividend Divisor Remainder Quotient
349 * -------- ------- --------- --------
350 * 10 7 3 1
351 * -10 7 4 -2
352 * 10 -7 -4 -2
353 * -10 -7 -3 1
354 *
355 *
356 * Table 3.4 - Symmetric Division Example
357 * Dividend Divisor Remainder Quotient
358 * -------- ------- --------- --------
359 * 10 7 3 1
360 * -10 7 -3 -1
361 * 10 -7 3 -1
362 * -10 -7 -3 1
363 */
364 ficl2IntegerQR
ficl2IntegerDivideFloored(ficl2Integer num,ficlInteger den)365 ficl2IntegerDivideFloored(ficl2Integer num, ficlInteger den)
366 {
367 ficl2IntegerQR qr;
368 ficl2UnsignedQR uqr;
369 ficl2Unsigned u;
370 int signRem = 1;
371 int signQuot = 1;
372
373 if (ficl2IntegerIsNegative(num)) {
374 num = ficl2IntegerNegate(num);
375 signQuot = -signQuot;
376 }
377
378 if (den < 0) {
379 den = -den;
380 signRem = -signRem;
381 signQuot = -signQuot;
382 }
383
384 FICL_2UNSIGNED_SET(FICL_2UNSIGNED_GET_HIGH(num),
385 FICL_2UNSIGNED_GET_LOW(num), u);
386 uqr = ficl2UnsignedDivide(u, (ficlUnsigned)den);
387 qr = FICL_2UNSIGNEDQR_TO_2INTEGERQR(uqr);
388 if (signQuot < 0) {
389 qr.quotient = ficl2IntegerNegate(qr.quotient);
390 if (qr.remainder != 0) {
391 qr.quotient = ficl2IntegerDecrement(qr.quotient);
392 qr.remainder = den - qr.remainder;
393 }
394 }
395
396 if (signRem < 0)
397 qr.remainder = -qr.remainder;
398
399 return (qr);
400 }
401
402 /*
403 * ficl2IntegerDivideSymmetric
404 * Divide an ficl2Unsigned by a ficlInteger and return a ficlInteger quotient
405 * and a ficlInteger remainder. The absolute values of quotient and remainder
406 * are not affected by the signs of the numerator and denominator
407 * (the operation is symmetric on the number line)
408 */
409 ficl2IntegerQR
ficl2IntegerDivideSymmetric(ficl2Integer num,ficlInteger den)410 ficl2IntegerDivideSymmetric(ficl2Integer num, ficlInteger den)
411 {
412 ficl2IntegerQR qr;
413 ficl2UnsignedQR uqr;
414 ficl2Unsigned u;
415 int signRem = 1;
416 int signQuot = 1;
417
418 if (ficl2IntegerIsNegative(num)) {
419 num = ficl2IntegerNegate(num);
420 signRem = -signRem;
421 signQuot = -signQuot;
422 }
423
424 if (den < 0) {
425 den = -den;
426 signQuot = -signQuot;
427 }
428
429 FICL_2UNSIGNED_SET(FICL_2UNSIGNED_GET_HIGH(num),
430 FICL_2UNSIGNED_GET_LOW(num), u);
431 uqr = ficl2UnsignedDivide(u, (ficlUnsigned)den);
432 qr = FICL_2UNSIGNEDQR_TO_2INTEGERQR(uqr);
433 if (signRem < 0)
434 qr.remainder = -qr.remainder;
435
436 if (signQuot < 0)
437 qr.quotient = ficl2IntegerNegate(qr.quotient);
438
439 return (qr);
440 }
441