1 /*******************************************************************
2 ** m a t h 6 4 . c
3 ** Forth Inspired Command Language - 64 bit math support routines
4 ** Author: John Sadler (john_sadler@alum.mit.edu)
5 ** Created: 25 January 1998
6 ** Rev 2.03: Support for 128 bit DP math. This file really ouught to
7 ** be renamed!
8 ** $Id: math64.c,v 1.9 2001/12/05 07:21:34 jsadler Exp $
9 *******************************************************************/
10 /*
11 ** Copyright (c) 1997-2001 John Sadler (john_sadler@alum.mit.edu)
12 ** All rights reserved.
13 **
14 ** Get the latest Ficl release at http://ficl.sourceforge.net
15 **
16 ** I am interested in hearing from anyone who uses ficl. If you have
17 ** a problem, a success story, a defect, an enhancement request, or
18 ** if you would like to contribute to the ficl release, please
19 ** contact me by email at the address above.
20 **
21 ** L I C E N S E and D I S C L A I M E R
22 **
23 ** Redistribution and use in source and binary forms, with or without
24 ** modification, are permitted provided that the following conditions
25 ** are met:
26 ** 1. Redistributions of source code must retain the above copyright
27 ** notice, this list of conditions and the following disclaimer.
28 ** 2. Redistributions in binary form must reproduce the above copyright
29 ** notice, this list of conditions and the following disclaimer in the
30 ** documentation and/or other materials provided with the distribution.
31 **
32 ** THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND
33 ** ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
34 ** IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
35 ** ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE
36 ** FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
37 ** DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
38 ** OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
39 ** HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
40 ** LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
41 ** OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
42 ** SUCH DAMAGE.
43 */
44
45
46 #include "ficl.h"
47 #include "math64.h"
48
49
50 /**************************************************************************
51 m 6 4 A b s
52 ** Returns the absolute value of an DPINT
53 **************************************************************************/
m64Abs(DPINT x)54 DPINT m64Abs(DPINT x)
55 {
56 if (m64IsNegative(x))
57 x = m64Negate(x);
58
59 return x;
60 }
61
62
63 /**************************************************************************
64 m 6 4 F l o o r e d D i v I
65 **
66 ** FROM THE FORTH ANS...
67 ** Floored division is integer division in which the remainder carries
68 ** the sign of the divisor or is zero, and the quotient is rounded to
69 ** its arithmetic floor. Symmetric division is integer division in which
70 ** the remainder carries the sign of the dividend or is zero and the
71 ** quotient is the mathematical quotient rounded towards zero or
72 ** truncated. Examples of each are shown in tables 3.3 and 3.4.
73 **
74 ** Table 3.3 - Floored Division Example
75 ** Dividend Divisor Remainder Quotient
76 ** -------- ------- --------- --------
77 ** 10 7 3 1
78 ** -10 7 4 -2
79 ** 10 -7 -4 -2
80 ** -10 -7 -3 1
81 **
82 **
83 ** Table 3.4 - Symmetric Division Example
84 ** Dividend Divisor Remainder Quotient
85 ** -------- ------- --------- --------
86 ** 10 7 3 1
87 ** -10 7 -3 -1
88 ** 10 -7 3 -1
89 ** -10 -7 -3 1
90 **************************************************************************/
m64FlooredDivI(DPINT num,FICL_INT den)91 INTQR m64FlooredDivI(DPINT num, FICL_INT den)
92 {
93 INTQR qr;
94 UNSQR uqr;
95 int signRem = 1;
96 int signQuot = 1;
97
98 if (m64IsNegative(num))
99 {
100 num = m64Negate(num);
101 signQuot = -signQuot;
102 }
103
104 if (den < 0)
105 {
106 den = -den;
107 signRem = -signRem;
108 signQuot = -signQuot;
109 }
110
111 uqr = ficlLongDiv(m64CastIU(num), (FICL_UNS)den);
112 qr = m64CastQRUI(uqr);
113 if (signQuot < 0)
114 {
115 qr.quot = -qr.quot;
116 if (qr.rem != 0)
117 {
118 qr.quot--;
119 qr.rem = den - qr.rem;
120 }
121 }
122
123 if (signRem < 0)
124 qr.rem = -qr.rem;
125
126 return qr;
127 }
128
129
130 /**************************************************************************
131 m 6 4 I s N e g a t i v e
132 ** Returns TRUE if the specified DPINT has its sign bit set.
133 **************************************************************************/
m64IsNegative(DPINT x)134 int m64IsNegative(DPINT x)
135 {
136 return (x.hi < 0);
137 }
138
139
140 /**************************************************************************
141 m 6 4 M a c
142 ** Mixed precision multiply and accumulate primitive for number building.
143 ** Multiplies DPUNS u by FICL_UNS mul and adds FICL_UNS add. Mul is typically
144 ** the numeric base, and add represents a digit to be appended to the
145 ** growing number.
146 ** Returns the result of the operation
147 **************************************************************************/
m64Mac(DPUNS u,FICL_UNS mul,FICL_UNS add)148 DPUNS m64Mac(DPUNS u, FICL_UNS mul, FICL_UNS add)
149 {
150 DPUNS resultLo = ficlLongMul(u.lo, mul);
151 DPUNS resultHi = ficlLongMul(u.hi, mul);
152 resultLo.hi += resultHi.lo;
153 resultHi.lo = resultLo.lo + add;
154
155 if (resultHi.lo < resultLo.lo)
156 resultLo.hi++;
157
158 resultLo.lo = resultHi.lo;
159
160 return resultLo;
161 }
162
163
164 /**************************************************************************
165 m 6 4 M u l I
166 ** Multiplies a pair of FICL_INTs and returns an DPINT result.
167 **************************************************************************/
m64MulI(FICL_INT x,FICL_INT y)168 DPINT m64MulI(FICL_INT x, FICL_INT y)
169 {
170 DPUNS prod;
171 int sign = 1;
172
173 if (x < 0)
174 {
175 sign = -sign;
176 x = -x;
177 }
178
179 if (y < 0)
180 {
181 sign = -sign;
182 y = -y;
183 }
184
185 prod = ficlLongMul(x, y);
186 if (sign > 0)
187 return m64CastUI(prod);
188 else
189 return m64Negate(m64CastUI(prod));
190 }
191
192
193 /**************************************************************************
194 m 6 4 N e g a t e
195 ** Negates an DPINT by complementing and incrementing.
196 **************************************************************************/
m64Negate(DPINT x)197 DPINT m64Negate(DPINT x)
198 {
199 x.hi = ~x.hi;
200 x.lo = ~x.lo;
201 x.lo ++;
202 if (x.lo == 0)
203 x.hi++;
204
205 return x;
206 }
207
208
209 /**************************************************************************
210 m 6 4 P u s h
211 ** Push an DPINT onto the specified stack in the order required
212 ** by ANS Forth (most significant cell on top)
213 ** These should probably be macros...
214 **************************************************************************/
i64Push(FICL_STACK * pStack,DPINT i64)215 void i64Push(FICL_STACK *pStack, DPINT i64)
216 {
217 stackPushINT(pStack, i64.lo);
218 stackPushINT(pStack, i64.hi);
219 return;
220 }
221
u64Push(FICL_STACK * pStack,DPUNS u64)222 void u64Push(FICL_STACK *pStack, DPUNS u64)
223 {
224 stackPushINT(pStack, u64.lo);
225 stackPushINT(pStack, u64.hi);
226 return;
227 }
228
229
230 /**************************************************************************
231 m 6 4 P o p
232 ** Pops an DPINT off the stack in the order required by ANS Forth
233 ** (most significant cell on top)
234 ** These should probably be macros...
235 **************************************************************************/
i64Pop(FICL_STACK * pStack)236 DPINT i64Pop(FICL_STACK *pStack)
237 {
238 DPINT ret;
239 ret.hi = stackPopINT(pStack);
240 ret.lo = stackPopINT(pStack);
241 return ret;
242 }
243
u64Pop(FICL_STACK * pStack)244 DPUNS u64Pop(FICL_STACK *pStack)
245 {
246 DPUNS ret;
247 ret.hi = stackPopINT(pStack);
248 ret.lo = stackPopINT(pStack);
249 return ret;
250 }
251
252
253 /**************************************************************************
254 m 6 4 S y m m e t r i c D i v
255 ** Divide an DPINT by a FICL_INT and return a FICL_INT quotient and a
256 ** FICL_INT remainder. The absolute values of quotient and remainder are not
257 ** affected by the signs of the numerator and denominator (the operation
258 ** is symmetric on the number line)
259 **************************************************************************/
m64SymmetricDivI(DPINT num,FICL_INT den)260 INTQR m64SymmetricDivI(DPINT num, FICL_INT den)
261 {
262 INTQR qr;
263 UNSQR uqr;
264 int signRem = 1;
265 int signQuot = 1;
266
267 if (m64IsNegative(num))
268 {
269 num = m64Negate(num);
270 signRem = -signRem;
271 signQuot = -signQuot;
272 }
273
274 if (den < 0)
275 {
276 den = -den;
277 signQuot = -signQuot;
278 }
279
280 uqr = ficlLongDiv(m64CastIU(num), (FICL_UNS)den);
281 qr = m64CastQRUI(uqr);
282 if (signRem < 0)
283 qr.rem = -qr.rem;
284
285 if (signQuot < 0)
286 qr.quot = -qr.quot;
287
288 return qr;
289 }
290
291
292 /**************************************************************************
293 m 6 4 U M o d
294 ** Divides a DPUNS by base (an UNS16) and returns an UNS16 remainder.
295 ** Writes the quotient back to the original DPUNS as a side effect.
296 ** This operation is typically used to convert an DPUNS to a text string
297 ** in any base. See words.c:numberSignS, for example.
298 ** Mechanics: performs 4 ficlLongDivs, each of which produces 16 bits
299 ** of the quotient. C does not provide a way to divide an FICL_UNS by an
300 ** UNS16 and get an FICL_UNS quotient (ldiv is closest, but it's signed,
301 ** unfortunately), so I've used ficlLongDiv.
302 **************************************************************************/
303 #if (BITS_PER_CELL == 32)
304
305 #define UMOD_SHIFT 16
306 #define UMOD_MASK 0x0000ffff
307
308 #elif (BITS_PER_CELL == 64)
309
310 #define UMOD_SHIFT 32
311 #define UMOD_MASK 0x00000000ffffffff
312
313 #endif
314
m64UMod(DPUNS * pUD,UNS16 base)315 UNS16 m64UMod(DPUNS *pUD, UNS16 base)
316 {
317 DPUNS ud;
318 UNSQR qr;
319 DPUNS result;
320
321 result.hi = result.lo = 0;
322
323 ud.hi = 0;
324 ud.lo = pUD->hi >> UMOD_SHIFT;
325 qr = ficlLongDiv(ud, (FICL_UNS)base);
326 result.hi = qr.quot << UMOD_SHIFT;
327
328 ud.lo = (qr.rem << UMOD_SHIFT) | (pUD->hi & UMOD_MASK);
329 qr = ficlLongDiv(ud, (FICL_UNS)base);
330 result.hi |= qr.quot & UMOD_MASK;
331
332 ud.lo = (qr.rem << UMOD_SHIFT) | (pUD->lo >> UMOD_SHIFT);
333 qr = ficlLongDiv(ud, (FICL_UNS)base);
334 result.lo = qr.quot << UMOD_SHIFT;
335
336 ud.lo = (qr.rem << UMOD_SHIFT) | (pUD->lo & UMOD_MASK);
337 qr = ficlLongDiv(ud, (FICL_UNS)base);
338 result.lo |= qr.quot & UMOD_MASK;
339
340 *pUD = result;
341
342 return (UNS16)(qr.rem);
343 }
344
345
346 /**************************************************************************
347 ** Contributed by
348 ** Michael A. Gauland gaulandm@mdhost.cse.tek.com
349 **************************************************************************/
350 #if PORTABLE_LONGMULDIV != 0
351 /**************************************************************************
352 m 6 4 A d d
353 **
354 **************************************************************************/
m64Add(DPUNS x,DPUNS y)355 DPUNS m64Add(DPUNS x, DPUNS y)
356 {
357 DPUNS result;
358 int carry;
359
360 result.hi = x.hi + y.hi;
361 result.lo = x.lo + y.lo;
362
363
364 carry = ((x.lo | y.lo) & CELL_HI_BIT) && !(result.lo & CELL_HI_BIT);
365 carry |= ((x.lo & y.lo) & CELL_HI_BIT);
366
367 if (carry)
368 {
369 result.hi++;
370 }
371
372 return result;
373 }
374
375
376 /**************************************************************************
377 m 6 4 S u b
378 **
379 **************************************************************************/
m64Sub(DPUNS x,DPUNS y)380 DPUNS m64Sub(DPUNS x, DPUNS y)
381 {
382 DPUNS result;
383
384 result.hi = x.hi - y.hi;
385 result.lo = x.lo - y.lo;
386
387 if (x.lo < y.lo)
388 {
389 result.hi--;
390 }
391
392 return result;
393 }
394
395
396 /**************************************************************************
397 m 6 4 A S L
398 ** 64 bit left shift
399 **************************************************************************/
m64ASL(DPUNS x)400 DPUNS m64ASL( DPUNS x )
401 {
402 DPUNS result;
403
404 result.hi = x.hi << 1;
405 if (x.lo & CELL_HI_BIT)
406 {
407 result.hi++;
408 }
409
410 result.lo = x.lo << 1;
411
412 return result;
413 }
414
415
416 /**************************************************************************
417 m 6 4 A S R
418 ** 64 bit right shift (unsigned - no sign extend)
419 **************************************************************************/
m64ASR(DPUNS x)420 DPUNS m64ASR( DPUNS x )
421 {
422 DPUNS result;
423
424 result.lo = x.lo >> 1;
425 if (x.hi & 1)
426 {
427 result.lo |= CELL_HI_BIT;
428 }
429
430 result.hi = x.hi >> 1;
431 return result;
432 }
433
434
435 /**************************************************************************
436 m 6 4 O r
437 ** 64 bit bitwise OR
438 **************************************************************************/
m64Or(DPUNS x,DPUNS y)439 DPUNS m64Or( DPUNS x, DPUNS y )
440 {
441 DPUNS result;
442
443 result.hi = x.hi | y.hi;
444 result.lo = x.lo | y.lo;
445
446 return result;
447 }
448
449
450 /**************************************************************************
451 m 6 4 C o m p a r e
452 ** Return -1 if x < y; 0 if x==y, and 1 if x > y.
453 **************************************************************************/
m64Compare(DPUNS x,DPUNS y)454 int m64Compare(DPUNS x, DPUNS y)
455 {
456 int result;
457
458 if (x.hi > y.hi)
459 {
460 result = +1;
461 }
462 else if (x.hi < y.hi)
463 {
464 result = -1;
465 }
466 else
467 {
468 /* High parts are equal */
469 if (x.lo > y.lo)
470 {
471 result = +1;
472 }
473 else if (x.lo < y.lo)
474 {
475 result = -1;
476 }
477 else
478 {
479 result = 0;
480 }
481 }
482
483 return result;
484 }
485
486
487 /**************************************************************************
488 f i c l L o n g M u l
489 ** Portable versions of ficlLongMul and ficlLongDiv in C
490 ** Contributed by:
491 ** Michael A. Gauland gaulandm@mdhost.cse.tek.com
492 **************************************************************************/
ficlLongMul(FICL_UNS x,FICL_UNS y)493 DPUNS ficlLongMul(FICL_UNS x, FICL_UNS y)
494 {
495 DPUNS result = { 0, 0 };
496 DPUNS addend;
497
498 addend.lo = y;
499 addend.hi = 0; /* No sign extension--arguments are unsigned */
500
501 while (x != 0)
502 {
503 if ( x & 1)
504 {
505 result = m64Add(result, addend);
506 }
507 x >>= 1;
508 addend = m64ASL(addend);
509 }
510 return result;
511 }
512
513
514 /**************************************************************************
515 f i c l L o n g D i v
516 ** Portable versions of ficlLongMul and ficlLongDiv in C
517 ** Contributed by:
518 ** Michael A. Gauland gaulandm@mdhost.cse.tek.com
519 **************************************************************************/
ficlLongDiv(DPUNS q,FICL_UNS y)520 UNSQR ficlLongDiv(DPUNS q, FICL_UNS y)
521 {
522 UNSQR result;
523 DPUNS quotient;
524 DPUNS subtrahend;
525 DPUNS mask;
526
527 quotient.lo = 0;
528 quotient.hi = 0;
529
530 subtrahend.lo = y;
531 subtrahend.hi = 0;
532
533 mask.lo = 1;
534 mask.hi = 0;
535
536 while ((m64Compare(subtrahend, q) < 0) &&
537 (subtrahend.hi & CELL_HI_BIT) == 0)
538 {
539 mask = m64ASL(mask);
540 subtrahend = m64ASL(subtrahend);
541 }
542
543 while (mask.lo != 0 || mask.hi != 0)
544 {
545 if (m64Compare(subtrahend, q) <= 0)
546 {
547 q = m64Sub( q, subtrahend);
548 quotient = m64Or(quotient, mask);
549 }
550 mask = m64ASR(mask);
551 subtrahend = m64ASR(subtrahend);
552 }
553
554 result.quot = quotient.lo;
555 result.rem = q.lo;
556 return result;
557 }
558
559 #endif
560
561