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
;