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
46 * $FreeBSD: src/sys/boot/ficl/math64.c,v 1.4 2002/04/09 17:45:11 dcs Exp $
47 * $DragonFly: src/sys/boot/ficl/math64.c,v 1.3 2003/11/10 06:08:33 dillon Exp $
54 /**************************************************************************
56 ** Returns the absolute value of an DPINT
57 **************************************************************************/
67 /**************************************************************************
68 m 6 4 F l o o r e d D i v I
70 ** FROM THE FORTH ANS...
71 ** Floored division is integer division in which the remainder carries
72 ** the sign of the divisor or is zero, and the quotient is rounded to
73 ** its arithmetic floor. Symmetric division is integer division in which
74 ** the remainder carries the sign of the dividend or is zero and the
75 ** quotient is the mathematical quotient rounded towards zero or
76 ** truncated. Examples of each are shown in tables 3.3 and 3.4.
78 ** Table 3.3 - Floored Division Example
79 ** Dividend Divisor Remainder Quotient
80 ** -------- ------- --------- --------
87 ** Table 3.4 - Symmetric Division Example
88 ** Dividend Divisor Remainder Quotient
89 ** -------- ------- --------- --------
94 **************************************************************************/
95 INTQR m64FlooredDivI(DPINT num, FICL_INT den)
102 if (m64IsNegative(num))
104 num = m64Negate(num);
105 signQuot = -signQuot;
112 signQuot = -signQuot;
115 uqr = ficlLongDiv(m64CastIU(num), (FICL_UNS)den);
116 qr = m64CastQRUI(uqr);
123 qr.rem = den - qr.rem;
134 /**************************************************************************
135 m 6 4 I s N e g a t i v e
136 ** Returns TRUE if the specified DPINT has its sign bit set.
137 **************************************************************************/
138 int m64IsNegative(DPINT x)
144 /**************************************************************************
146 ** Mixed precision multiply and accumulate primitive for number building.
147 ** Multiplies DPUNS u by FICL_UNS mul and adds FICL_UNS add. Mul is typically
148 ** the numeric base, and add represents a digit to be appended to the
150 ** Returns the result of the operation
151 **************************************************************************/
152 DPUNS m64Mac(DPUNS u, FICL_UNS mul, FICL_UNS add)
154 DPUNS resultLo = ficlLongMul(u.lo, mul);
155 DPUNS resultHi = ficlLongMul(u.hi, mul);
156 resultLo.hi += resultHi.lo;
157 resultHi.lo = resultLo.lo + add;
159 if (resultHi.lo < resultLo.lo)
162 resultLo.lo = resultHi.lo;
168 /**************************************************************************
170 ** Multiplies a pair of FICL_INTs and returns an DPINT result.
171 **************************************************************************/
172 DPINT m64MulI(FICL_INT x, FICL_INT y)
189 prod = ficlLongMul(x, y);
191 return m64CastUI(prod);
193 return m64Negate(m64CastUI(prod));
197 /**************************************************************************
199 ** Negates an DPINT by complementing and incrementing.
200 **************************************************************************/
201 DPINT m64Negate(DPINT x)
213 /**************************************************************************
215 ** Push an DPINT onto the specified stack in the order required
216 ** by ANS Forth (most significant cell on top)
217 ** These should probably be macros...
218 **************************************************************************/
219 void i64Push(FICL_STACK *pStack, DPINT i64)
221 stackPushINT(pStack, i64.lo);
222 stackPushINT(pStack, i64.hi);
226 void u64Push(FICL_STACK *pStack, DPUNS u64)
228 stackPushINT(pStack, u64.lo);
229 stackPushINT(pStack, u64.hi);
234 /**************************************************************************
236 ** Pops an DPINT off the stack in the order required by ANS Forth
237 ** (most significant cell on top)
238 ** These should probably be macros...
239 **************************************************************************/
240 DPINT i64Pop(FICL_STACK *pStack)
243 ret.hi = stackPopINT(pStack);
244 ret.lo = stackPopINT(pStack);
248 DPUNS u64Pop(FICL_STACK *pStack)
251 ret.hi = stackPopINT(pStack);
252 ret.lo = stackPopINT(pStack);
257 /**************************************************************************
258 m 6 4 S y m m e t r i c D i v
259 ** Divide an DPINT by a FICL_INT and return a FICL_INT quotient and a
260 ** FICL_INT remainder. The absolute values of quotient and remainder are not
261 ** affected by the signs of the numerator and denominator (the operation
262 ** is symmetric on the number line)
263 **************************************************************************/
264 INTQR m64SymmetricDivI(DPINT num, FICL_INT den)
271 if (m64IsNegative(num))
273 num = m64Negate(num);
275 signQuot = -signQuot;
281 signQuot = -signQuot;
284 uqr = ficlLongDiv(m64CastIU(num), (FICL_UNS)den);
285 qr = m64CastQRUI(uqr);
296 /**************************************************************************
298 ** Divides a DPUNS by base (an UNS16) and returns an UNS16 remainder.
299 ** Writes the quotient back to the original DPUNS as a side effect.
300 ** This operation is typically used to convert an DPUNS to a text string
301 ** in any base. See words.c:numberSignS, for example.
302 ** Mechanics: performs 4 ficlLongDivs, each of which produces 16 bits
303 ** of the quotient. C does not provide a way to divide an FICL_UNS by an
304 ** UNS16 and get an FICL_UNS quotient (ldiv is closest, but it's signed,
305 ** unfortunately), so I've used ficlLongDiv.
306 **************************************************************************/
307 #if (BITS_PER_CELL == 32)
309 #define UMOD_SHIFT 16
310 #define UMOD_MASK 0x0000ffff
312 #elif (BITS_PER_CELL == 64)
314 #define UMOD_SHIFT 32
315 #define UMOD_MASK 0x00000000ffffffff
319 UNS16 m64UMod(DPUNS *pUD, UNS16 base)
325 result.hi = result.lo = 0;
328 ud.lo = pUD->hi >> UMOD_SHIFT;
329 qr = ficlLongDiv(ud, (FICL_UNS)base);
330 result.hi = qr.quot << UMOD_SHIFT;
332 ud.lo = (qr.rem << UMOD_SHIFT) | (pUD->hi & UMOD_MASK);
333 qr = ficlLongDiv(ud, (FICL_UNS)base);
334 result.hi |= qr.quot & UMOD_MASK;
336 ud.lo = (qr.rem << UMOD_SHIFT) | (pUD->lo >> UMOD_SHIFT);
337 qr = ficlLongDiv(ud, (FICL_UNS)base);
338 result.lo = qr.quot << UMOD_SHIFT;
340 ud.lo = (qr.rem << UMOD_SHIFT) | (pUD->lo & UMOD_MASK);
341 qr = ficlLongDiv(ud, (FICL_UNS)base);
342 result.lo |= qr.quot & UMOD_MASK;
346 return (UNS16)(qr.rem);
350 /**************************************************************************
352 ** Michael A. Gauland gaulandm@mdhost.cse.tek.com
353 **************************************************************************/
354 #if PORTABLE_LONGMULDIV != 0
355 /**************************************************************************
358 **************************************************************************/
359 DPUNS m64Add(DPUNS x, DPUNS y)
364 result.hi = x.hi + y.hi;
365 result.lo = x.lo + y.lo;
368 carry = ((x.lo | y.lo) & CELL_HI_BIT) && !(result.lo & CELL_HI_BIT);
369 carry |= ((x.lo & y.lo) & CELL_HI_BIT);
380 /**************************************************************************
383 **************************************************************************/
384 DPUNS m64Sub(DPUNS x, DPUNS y)
388 result.hi = x.hi - y.hi;
389 result.lo = x.lo - y.lo;
400 /**************************************************************************
403 **************************************************************************/
404 DPUNS m64ASL( DPUNS x )
408 result.hi = x.hi << 1;
409 if (x.lo & CELL_HI_BIT)
414 result.lo = x.lo << 1;
420 /**************************************************************************
422 ** 64 bit right shift (unsigned - no sign extend)
423 **************************************************************************/
424 DPUNS m64ASR( DPUNS x )
428 result.lo = x.lo >> 1;
431 result.lo |= CELL_HI_BIT;
434 result.hi = x.hi >> 1;
439 /**************************************************************************
442 **************************************************************************/
443 DPUNS m64Or( DPUNS x, DPUNS y )
447 result.hi = x.hi | y.hi;
448 result.lo = x.lo | y.lo;
454 /**************************************************************************
456 ** Return -1 if x < y; 0 if x==y, and 1 if x > y.
457 **************************************************************************/
458 int m64Compare(DPUNS x, DPUNS y)
466 else if (x.hi < y.hi)
472 /* High parts are equal */
477 else if (x.lo < y.lo)
491 /**************************************************************************
492 f i c l L o n g M u l
493 ** Portable versions of ficlLongMul and ficlLongDiv in C
495 ** Michael A. Gauland gaulandm@mdhost.cse.tek.com
496 **************************************************************************/
497 DPUNS ficlLongMul(FICL_UNS x, FICL_UNS y)
499 DPUNS result = { 0, 0 };
503 addend.hi = 0; /* No sign extension--arguments are unsigned */
509 result = m64Add(result, addend);
512 addend = m64ASL(addend);
518 /**************************************************************************
519 f i c l L o n g D i v
520 ** Portable versions of ficlLongMul and ficlLongDiv in C
522 ** Michael A. Gauland gaulandm@mdhost.cse.tek.com
523 **************************************************************************/
524 UNSQR ficlLongDiv(DPUNS q, FICL_UNS y)
540 while ((m64Compare(subtrahend, q) < 0) &&
541 (subtrahend.hi & CELL_HI_BIT) == 0)
544 subtrahend = m64ASL(subtrahend);
547 while (mask.lo != 0 || mask.hi != 0)
549 if (m64Compare(subtrahend, q) <= 0)
551 q = m64Sub( q, subtrahend);
552 quotient = m64Or(quotient, mask);
555 subtrahend = m64ASR(subtrahend);
558 result.quot = quotient.lo;