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 <stdint.h>
48
49 #include "ficl.h"
50
51
52 #if FICL_PLATFORM_HAS_2INTEGER
53
54
55
ficl2UnsignedDivide(ficl2Unsigned q,ficlUnsigned y)56 ficl2UnsignedQR ficl2UnsignedDivide(ficl2Unsigned q, ficlUnsigned y)
57 {
58 ficl2UnsignedQR result;
59
60 result.quotient = q / y;
61 /*
62 ** Once we have the quotient, it's cheaper to calculate the
63 ** remainder this way than with % (mod). --lch
64 */
65 result.remainder = (ficlInteger)(q - (result.quotient * y));
66
67 return result;
68 }
69
70
71 #else /* FICL_PLATFORM_HAS_2INTEGER */
72
73
74 #define FICL_CELL_HIGH_BIT ((uintmax_t)1 << (FICL_BITS_PER_CELL-1))
75 #define UMOD_SHIFT (FICL_BITS_PER_CELL / 2)
76 #define UMOD_MASK ((1L << (FICL_BITS_PER_CELL / 2)) - 1)
77
78
79 /**************************************************************************
80 ficl2IntegerIsNegative
81 ** Returns TRUE if the specified ficl2Unsigned has its sign bit set.
82 **************************************************************************/
ficl2IntegerIsNegative(ficl2Integer x)83 int ficl2IntegerIsNegative(ficl2Integer x)
84 {
85 return (x.high < 0);
86 }
87
88
89 /**************************************************************************
90 ficl2IntegerNegate
91 ** Negates an ficl2Unsigned by complementing and incrementing.
92 **************************************************************************/
ficl2IntegerNegate(ficl2Integer x)93 ficl2Integer ficl2IntegerNegate(ficl2Integer x)
94 {
95 x.high = ~x.high;
96 x.low = ~x.low;
97 x.low ++;
98 if (x.low == 0)
99 x.high++;
100
101 return x;
102 }
103
104 /**************************************************************************
105 ficl2UnsignedMultiplyAccumulate
106 ** Mixed precision multiply and accumulate primitive for number building.
107 ** Multiplies ficl2Unsigned u by ficlUnsigned mul and adds ficlUnsigned add. Mul is typically
108 ** the numeric base, and add represents a digit to be appended to the
109 ** growing number.
110 ** Returns the result of the operation
111 **************************************************************************/
ficl2UnsignedMultiplyAccumulate(ficl2Unsigned u,ficlUnsigned mul,ficlUnsigned add)112 ficl2Unsigned ficl2UnsignedMultiplyAccumulate(ficl2Unsigned u, ficlUnsigned mul, ficlUnsigned add)
113 {
114 ficl2Unsigned resultLo = ficl2UnsignedMultiply(u.low, mul);
115 ficl2Unsigned resultHi = ficl2UnsignedMultiply(u.high, mul);
116 resultLo.high += resultHi.low;
117 resultHi.low = resultLo.low + add;
118
119 if (resultHi.low < resultLo.low)
120 resultLo.high++;
121
122 resultLo.low = resultHi.low;
123
124 return resultLo;
125 }
126
127
128 /**************************************************************************
129 ficl2IntegerMultiply
130 ** Multiplies a pair of ficlIntegers and returns an ficl2Integer result.
131 **************************************************************************/
ficl2IntegerMultiply(ficlInteger x,ficlInteger y)132 ficl2Integer ficl2IntegerMultiply(ficlInteger x, ficlInteger y)
133 {
134 ficl2Unsigned prod;
135 int sign = 1;
136
137 if (x < 0)
138 {
139 sign = -sign;
140 x = -x;
141 }
142
143 if (y < 0)
144 {
145 sign = -sign;
146 y = -y;
147 }
148
149 prod = ficl2UnsignedMultiply(x, y);
150 if (sign > 0)
151 return FICL_2UNSIGNED_TO_2INTEGER(prod);
152 else
153 return ficl2IntegerNegate(FICL_2UNSIGNED_TO_2INTEGER(prod));
154 }
155
156
157
ficl2IntegerDecrement(ficl2Integer x)158 ficl2Integer ficl2IntegerDecrement(ficl2Integer x)
159 {
160 if (x.low == INTMAX_MIN)
161 x.high--;
162 x.low--;
163
164 return x;
165 }
166
167
ficl2UnsignedAdd(ficl2Unsigned x,ficl2Unsigned y)168 ficl2Unsigned ficl2UnsignedAdd(ficl2Unsigned x, ficl2Unsigned y)
169 {
170 ficl2Unsigned result;
171
172 result.high = x.high + y.high;
173 result.low = x.low + y.low;
174
175 if (result.low < y.low)
176 {
177 result.high++;
178 }
179
180 return result;
181 }
182
183 /**************************************************************************
184 ficl2UnsignedMultiply
185 ** Contributed by:
186 ** Michael A. Gauland gaulandm@mdhost.cse.tek.com
187 **************************************************************************/
ficl2UnsignedMultiply(ficlUnsigned x,ficlUnsigned y)188 ficl2Unsigned ficl2UnsignedMultiply(ficlUnsigned x, ficlUnsigned y)
189 {
190 ficl2Unsigned result = { 0, 0 };
191 ficl2Unsigned addend;
192
193 addend.low = y;
194 addend.high = 0; /* No sign extension--arguments are unsigned */
195
196 while (x != 0)
197 {
198 if ( x & 1)
199 {
200 result = ficl2UnsignedAdd(result, addend);
201 }
202 x >>= 1;
203 addend = ficl2UnsignedArithmeticShiftLeft(addend);
204 }
205 return result;
206 }
207
208
209
210 /**************************************************************************
211 ficl2UnsignedSubtract
212 **
213 **************************************************************************/
ficl2UnsignedSubtract(ficl2Unsigned x,ficl2Unsigned y)214 ficl2Unsigned ficl2UnsignedSubtract(ficl2Unsigned x, ficl2Unsigned y)
215 {
216 ficl2Unsigned result;
217
218 result.high = x.high - y.high;
219 result.low = x.low - y.low;
220
221 if (x.low < y.low)
222 {
223 result.high--;
224 }
225
226 return result;
227 }
228
229
230 /**************************************************************************
231 ficl2UnsignedArithmeticShiftLeft
232 ** 64 bit left shift
233 **************************************************************************/
ficl2UnsignedArithmeticShiftLeft(ficl2Unsigned x)234 ficl2Unsigned ficl2UnsignedArithmeticShiftLeft( ficl2Unsigned x )
235 {
236 ficl2Unsigned result;
237
238 result.high = x.high << 1;
239 if (x.low & FICL_CELL_HIGH_BIT)
240 {
241 result.high++;
242 }
243
244 result.low = x.low << 1;
245
246 return result;
247 }
248
249
250 /**************************************************************************
251 ficl2UnsignedArithmeticShiftRight
252 ** 64 bit right shift (unsigned - no sign extend)
253 **************************************************************************/
ficl2UnsignedArithmeticShiftRight(ficl2Unsigned x)254 ficl2Unsigned ficl2UnsignedArithmeticShiftRight( ficl2Unsigned x )
255 {
256 ficl2Unsigned result;
257
258 result.low = x.low >> 1;
259 if (x.high & 1)
260 {
261 result.low |= FICL_CELL_HIGH_BIT;
262 }
263
264 result.high = x.high >> 1;
265 return result;
266 }
267
268
269 /**************************************************************************
270 ficl2UnsignedOr
271 ** 64 bit bitwise OR
272 **************************************************************************/
ficl2UnsignedOr(ficl2Unsigned x,ficl2Unsigned y)273 ficl2Unsigned ficl2UnsignedOr( ficl2Unsigned x, ficl2Unsigned y )
274 {
275 ficl2Unsigned result;
276
277 result.high = x.high | y.high;
278 result.low = x.low | y.low;
279
280 return result;
281 }
282
283
284 /**************************************************************************
285 ficl2UnsignedCompare
286 ** Return -1 if x < y; 0 if x==y, and 1 if x > y.
287 **************************************************************************/
ficl2UnsignedCompare(ficl2Unsigned x,ficl2Unsigned y)288 int ficl2UnsignedCompare(ficl2Unsigned x, ficl2Unsigned y)
289 {
290 if (x.high > y.high)
291 return 1;
292 if (x.high < y.high)
293 return -1;
294
295 /* High parts are equal */
296
297 if (x.low > y.low)
298 return 1;
299 else if (x.low < y.low)
300 return -1;
301
302 return 0;
303 }
304
305
306
307 /**************************************************************************
308 ficl2UnsignedDivide
309 ** Portable versions of ficl2Multiply and ficl2Divide in C
310 ** Contributed by:
311 ** Michael A. Gauland gaulandm@mdhost.cse.tek.com
312 **************************************************************************/
ficl2UnsignedDivide(ficl2Unsigned q,ficlUnsigned y)313 ficl2UnsignedQR ficl2UnsignedDivide(ficl2Unsigned q, ficlUnsigned y)
314 {
315 ficl2UnsignedQR result;
316 ficl2Unsigned quotient;
317 ficl2Unsigned subtrahend;
318 ficl2Unsigned mask;
319
320 quotient.low = 0;
321 quotient.high = 0;
322
323 subtrahend.low = y;
324 subtrahend.high = 0;
325
326 mask.low = 1;
327 mask.high = 0;
328
329 while ((ficl2UnsignedCompare(subtrahend, q) < 0) &&
330 (subtrahend.high & FICL_CELL_HIGH_BIT) == 0)
331 {
332 mask = ficl2UnsignedArithmeticShiftLeft(mask);
333 subtrahend = ficl2UnsignedArithmeticShiftLeft(subtrahend);
334 }
335
336 while (mask.low != 0 || mask.high != 0)
337 {
338 if (ficl2UnsignedCompare(subtrahend, q) <= 0)
339 {
340 q = ficl2UnsignedSubtract( q, subtrahend);
341 quotient = ficl2UnsignedOr(quotient, mask);
342 }
343 mask = ficl2UnsignedArithmeticShiftRight(mask);
344 subtrahend = ficl2UnsignedArithmeticShiftRight(subtrahend);
345 }
346
347 result.quotient = quotient;
348 result.remainder = q.low;
349 return result;
350 }
351
352 #endif /* !FICL_PLATFORM_HAS_2INTEGER */
353
354
355
356 /**************************************************************************
357 ficl2IntegerAbsoluteValue
358 ** Returns the absolute value of an ficl2Unsigned
359 **************************************************************************/
ficl2IntegerAbsoluteValue(ficl2Integer x)360 ficl2Integer ficl2IntegerAbsoluteValue(ficl2Integer x)
361 {
362 if (ficl2IntegerIsNegative(x))
363 return ficl2IntegerNegate(x);
364 return x;
365 }
366
367
368 /**************************************************************************
369 ficl2IntegerDivideFloored
370 **
371 ** FROM THE FORTH ANS...
372 ** Floored division is integer division in which the remainder carries
373 ** the sign of the divisor or is zero, and the quotient is rounded to
374 ** its arithmetic floor. Symmetric division is integer division in which
375 ** the remainder carries the sign of the dividend or is zero and the
376 ** quotient is the mathematical quotient rounded towards zero or
377 ** truncated. Examples of each are shown in tables 3.3 and 3.4.
378 **
379 ** Table 3.3 - Floored Division Example
380 ** Dividend Divisor Remainder Quotient
381 ** -------- ------- --------- --------
382 ** 10 7 3 1
383 ** -10 7 4 -2
384 ** 10 -7 -4 -2
385 ** -10 -7 -3 1
386 **
387 **
388 ** Table 3.4 - Symmetric Division Example
389 ** Dividend Divisor Remainder Quotient
390 ** -------- ------- --------- --------
391 ** 10 7 3 1
392 ** -10 7 -3 -1
393 ** 10 -7 3 -1
394 ** -10 -7 -3 1
395 **************************************************************************/
ficl2IntegerDivideFloored(ficl2Integer num,ficlInteger den)396 ficl2IntegerQR ficl2IntegerDivideFloored(ficl2Integer num, ficlInteger den)
397 {
398 ficl2IntegerQR qr;
399 ficl2UnsignedQR uqr;
400 int signRem = 1;
401 int signQuot = 1;
402
403 if (ficl2IntegerIsNegative(num))
404 {
405 num = ficl2IntegerNegate(num);
406 signQuot = -signQuot;
407 }
408
409 if (den < 0)
410 {
411 den = -den;
412 signRem = -signRem;
413 signQuot = -signQuot;
414 }
415
416 uqr = ficl2UnsignedDivide(FICL_2INTEGER_TO_2UNSIGNED(num), (ficlUnsigned)den);
417 qr = FICL_2UNSIGNEDQR_TO_2INTEGERQR(uqr);
418 if (signQuot < 0)
419 {
420 qr.quotient = ficl2IntegerNegate(qr.quotient);
421 if (qr.remainder != 0)
422 {
423 qr.quotient = ficl2IntegerDecrement(qr.quotient);
424 qr.remainder = den - qr.remainder;
425 }
426 }
427
428 if (signRem < 0)
429 qr.remainder = -qr.remainder;
430
431 return qr;
432 }
433
434
435
436 /**************************************************************************
437 ficl2IntegerDivideSymmetric
438 ** Divide an ficl2Unsigned by a ficlInteger and return a ficlInteger quotient and a
439 ** ficlInteger remainder. The absolute values of quotient and remainder are not
440 ** affected by the signs of the numerator and denominator (the operation
441 ** is symmetric on the number line)
442 **************************************************************************/
ficl2IntegerDivideSymmetric(ficl2Integer num,ficlInteger den)443 ficl2IntegerQR ficl2IntegerDivideSymmetric(ficl2Integer num, ficlInteger den)
444 {
445 ficl2IntegerQR qr;
446 ficl2UnsignedQR uqr;
447 int signRem = 1;
448 int signQuot = 1;
449
450 if (ficl2IntegerIsNegative(num))
451 {
452 num = ficl2IntegerNegate(num);
453 signRem = -signRem;
454 signQuot = -signQuot;
455 }
456
457 if (den < 0)
458 {
459 den = -den;
460 signQuot = -signQuot;
461 }
462
463 uqr = ficl2UnsignedDivide(FICL_2INTEGER_TO_2UNSIGNED(num), (ficlUnsigned)den);
464 qr = FICL_2UNSIGNEDQR_TO_2INTEGERQR(uqr);
465 if (signRem < 0)
466 qr.remainder = -qr.remainder;
467
468 if (signQuot < 0)
469 qr.quotient = ficl2IntegerNegate(qr.quotient);
470
471 return qr;
472 }
473
474
475