1 /*******************************************************************
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
8 ** $Id: math64.c,v 1.9 2001/12/05 07:21:34 jsadler Exp $
9 *******************************************************************/
11 ** Copyright (c) 1997-2001 John Sadler (john_sadler@alum.mit.edu)
12 ** All rights reserved.
14 ** Get the latest Ficl release at http://ficl.sourceforge.net
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.
21 ** L I C E N S E and D I S C L A I M E R
23 ** Redistribution and use in source and binary forms, with or without
24 ** modification, are permitted provided that the following conditions
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.
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
51 /**************************************************************************
53 ** Returns the absolute value of an DPINT
54 **************************************************************************/
64 /**************************************************************************
65 m 6 4 F l o o r e d D i v I
67 ** FROM THE FORTH ANS...
68 ** Floored division is integer division in which the remainder carries
69 ** the sign of the divisor or is zero, and the quotient is rounded to
70 ** its arithmetic floor. Symmetric division is integer division in which
71 ** the remainder carries the sign of the dividend or is zero and the
72 ** quotient is the mathematical quotient rounded towards zero or
73 ** truncated. Examples of each are shown in tables 3.3 and 3.4.
75 ** Table 3.3 - Floored Division Example
76 ** Dividend Divisor Remainder Quotient
77 ** -------- ------- --------- --------
84 ** Table 3.4 - Symmetric Division Example
85 ** Dividend Divisor Remainder Quotient
86 ** -------- ------- --------- --------
91 **************************************************************************/
92 INTQR m64FlooredDivI(DPINT num, FICL_INT den)
99 if (m64IsNegative(num))
101 num = m64Negate(num);
102 signQuot = -signQuot;
109 signQuot = -signQuot;
112 uqr = ficlLongDiv(m64CastIU(num), (FICL_UNS)den);
113 qr = m64CastQRUI(uqr);
120 qr.rem = den - qr.rem;
131 /**************************************************************************
132 m 6 4 I s N e g a t i v e
133 ** Returns TRUE if the specified DPINT has its sign bit set.
134 **************************************************************************/
135 int m64IsNegative(DPINT x)
141 /**************************************************************************
143 ** Mixed precision multiply and accumulate primitive for number building.
144 ** Multiplies DPUNS u by FICL_UNS mul and adds FICL_UNS add. Mul is typically
145 ** the numeric base, and add represents a digit to be appended to the
147 ** Returns the result of the operation
148 **************************************************************************/
149 DPUNS m64Mac(DPUNS u, FICL_UNS mul, FICL_UNS add)
151 DPUNS resultLo = ficlLongMul(u.lo, mul);
152 DPUNS resultHi = ficlLongMul(u.hi, mul);
153 resultLo.hi += resultHi.lo;
154 resultHi.lo = resultLo.lo + add;
156 if (resultHi.lo < resultLo.lo)
159 resultLo.lo = resultHi.lo;
165 /**************************************************************************
167 ** Multiplies a pair of FICL_INTs and returns an DPINT result.
168 **************************************************************************/
169 DPINT m64MulI(FICL_INT x, FICL_INT y)
186 prod = ficlLongMul(x, y);
188 return m64CastUI(prod);
190 return m64Negate(m64CastUI(prod));
194 /**************************************************************************
196 ** Negates an DPINT by complementing and incrementing.
197 **************************************************************************/
198 DPINT m64Negate(DPINT x)
210 /**************************************************************************
212 ** Push an DPINT onto the specified stack in the order required
213 ** by ANS Forth (most significant cell on top)
214 ** These should probably be macros...
215 **************************************************************************/
216 void i64Push(FICL_STACK *pStack, DPINT i64)
218 stackPushINT(pStack, i64.lo);
219 stackPushINT(pStack, i64.hi);
223 void u64Push(FICL_STACK *pStack, DPUNS u64)
225 stackPushINT(pStack, u64.lo);
226 stackPushINT(pStack, u64.hi);
231 /**************************************************************************
233 ** Pops an DPINT off the stack in the order required by ANS Forth
234 ** (most significant cell on top)
235 ** These should probably be macros...
236 **************************************************************************/
237 DPINT i64Pop(FICL_STACK *pStack)
240 ret.hi = stackPopINT(pStack);
241 ret.lo = stackPopINT(pStack);
245 DPUNS u64Pop(FICL_STACK *pStack)
248 ret.hi = stackPopINT(pStack);
249 ret.lo = stackPopINT(pStack);
254 /**************************************************************************
255 m 6 4 S y m m e t r i c D i v
256 ** Divide an DPINT by a FICL_INT and return a FICL_INT quotient and a
257 ** FICL_INT remainder. The absolute values of quotient and remainder are not
258 ** affected by the signs of the numerator and denominator (the operation
259 ** is symmetric on the number line)
260 **************************************************************************/
261 INTQR m64SymmetricDivI(DPINT num, FICL_INT den)
268 if (m64IsNegative(num))
270 num = m64Negate(num);
272 signQuot = -signQuot;
278 signQuot = -signQuot;
281 uqr = ficlLongDiv(m64CastIU(num), (FICL_UNS)den);
282 qr = m64CastQRUI(uqr);
293 /**************************************************************************
295 ** Divides a DPUNS by base (an UNS16) and returns an UNS16 remainder.
296 ** Writes the quotient back to the original DPUNS as a side effect.
297 ** This operation is typically used to convert an DPUNS to a text string
298 ** in any base. See words.c:numberSignS, for example.
299 ** Mechanics: performs 4 ficlLongDivs, each of which produces 16 bits
300 ** of the quotient. C does not provide a way to divide an FICL_UNS by an
301 ** UNS16 and get an FICL_UNS quotient (ldiv is closest, but it's signed,
302 ** unfortunately), so I've used ficlLongDiv.
303 **************************************************************************/
304 #if (BITS_PER_CELL == 32)
306 #define UMOD_SHIFT 16
307 #define UMOD_MASK 0x0000ffff
309 #elif (BITS_PER_CELL == 64)
311 #define UMOD_SHIFT 32
312 #define UMOD_MASK 0x00000000ffffffff
316 UNS16 m64UMod(DPUNS *pUD, UNS16 base)
322 result.hi = result.lo = 0;
325 ud.lo = pUD->hi >> UMOD_SHIFT;
326 qr = ficlLongDiv(ud, (FICL_UNS)base);
327 result.hi = qr.quot << UMOD_SHIFT;
329 ud.lo = (qr.rem << UMOD_SHIFT) | (pUD->hi & UMOD_MASK);
330 qr = ficlLongDiv(ud, (FICL_UNS)base);
331 result.hi |= qr.quot & UMOD_MASK;
333 ud.lo = (qr.rem << UMOD_SHIFT) | (pUD->lo >> UMOD_SHIFT);
334 qr = ficlLongDiv(ud, (FICL_UNS)base);
335 result.lo = qr.quot << UMOD_SHIFT;
337 ud.lo = (qr.rem << UMOD_SHIFT) | (pUD->lo & UMOD_MASK);
338 qr = ficlLongDiv(ud, (FICL_UNS)base);
339 result.lo |= qr.quot & UMOD_MASK;
343 return (UNS16)(qr.rem);
347 /**************************************************************************
349 ** Michael A. Gauland gaulandm@mdhost.cse.tek.com
350 **************************************************************************/
351 #if PORTABLE_LONGMULDIV != 0
352 /**************************************************************************
355 **************************************************************************/
356 DPUNS m64Add(DPUNS x, DPUNS y)
361 result.hi = x.hi + y.hi;
362 result.lo = x.lo + y.lo;
365 carry = ((x.lo | y.lo) & CELL_HI_BIT) && !(result.lo & CELL_HI_BIT);
366 carry |= ((x.lo & y.lo) & CELL_HI_BIT);
377 /**************************************************************************
380 **************************************************************************/
381 DPUNS m64Sub(DPUNS x, DPUNS y)
385 result.hi = x.hi - y.hi;
386 result.lo = x.lo - y.lo;
397 /**************************************************************************
400 **************************************************************************/
401 DPUNS m64ASL( DPUNS x )
405 result.hi = x.hi << 1;
406 if (x.lo & CELL_HI_BIT)
411 result.lo = x.lo << 1;
417 /**************************************************************************
419 ** 64 bit right shift (unsigned - no sign extend)
420 **************************************************************************/
421 DPUNS m64ASR( DPUNS x )
425 result.lo = x.lo >> 1;
428 result.lo |= CELL_HI_BIT;
431 result.hi = x.hi >> 1;
436 /**************************************************************************
439 **************************************************************************/
440 DPUNS m64Or( DPUNS x, DPUNS y )
444 result.hi = x.hi | y.hi;
445 result.lo = x.lo | y.lo;
451 /**************************************************************************
453 ** Return -1 if x < y; 0 if x==y, and 1 if x > y.
454 **************************************************************************/
455 int m64Compare(DPUNS x, DPUNS y)
463 else if (x.hi < y.hi)
469 /* High parts are equal */
474 else if (x.lo < y.lo)
488 /**************************************************************************
489 f i c l L o n g M u l
490 ** Portable versions of ficlLongMul and ficlLongDiv in C
492 ** Michael A. Gauland gaulandm@mdhost.cse.tek.com
493 **************************************************************************/
494 DPUNS ficlLongMul(FICL_UNS x, FICL_UNS y)
496 DPUNS result = { 0, 0 };
500 addend.hi = 0; /* No sign extension--arguments are unsigned */
506 result = m64Add(result, addend);
509 addend = m64ASL(addend);
515 /**************************************************************************
516 f i c l L o n g D i v
517 ** Portable versions of ficlLongMul and ficlLongDiv in C
519 ** Michael A. Gauland gaulandm@mdhost.cse.tek.com
520 **************************************************************************/
521 UNSQR ficlLongDiv(DPUNS q, FICL_UNS y)
537 while ((m64Compare(subtrahend, q) < 0) &&
538 (subtrahend.hi & CELL_HI_BIT) == 0)
541 subtrahend = m64ASL(subtrahend);
544 while (mask.lo != 0 || mask.hi != 0)
546 if (m64Compare(subtrahend, q) <= 0)
548 q = m64Sub( q, subtrahend);
549 quotient = m64Or(quotient, mask);
552 subtrahend = m64ASR(subtrahend);
555 result.quot = quotient.lo;