math64.c revision 40843 
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**
7*******************************************************************/
8
9#include "ficl.h"
10#include "math64.h"
11
12
13/**************************************************************************
14                        m 6 4 A b s
15** Returns the absolute value of an INT64
16**************************************************************************/
17INT64 m64Abs(INT64 x)
18{
19    if (m64IsNegative(x))
20        x = m64Negate(x);
21
22    return x;
23}
24
25
26/**************************************************************************
27                        m 6 4 F l o o r e d D i v I
28**
29** FROM THE FORTH ANS...
30** Floored division is integer division in which the remainder carries
31** the sign of the divisor or is zero, and the quotient is rounded to
32** its arithmetic floor. Symmetric division is integer division in which
33** the remainder carries the sign of the dividend or is zero and the
34** quotient is the mathematical quotient rounded towards zero or
35** truncated. Examples of each are shown in tables 3.3 and 3.4.
36**
37** Table 3.3 - Floored Division Example
38** Dividend        Divisor Remainder       Quotient
39** --------        ------- ---------       --------
40**  10                7       3                1
41** -10                7       4               -2
42**  10               -7      -4               -2
43** -10               -7      -3                1
44**
45**
46** Table 3.4 - Symmetric Division Example
47** Dividend        Divisor Remainder       Quotient
48** --------        ------- ---------       --------
49**  10                7       3                1
50** -10                7      -3               -1
51**  10               -7       3               -1
52** -10               -7      -3                1
53**************************************************************************/
54INTQR m64FlooredDivI(INT64 num, INT32 den)
55{
56    INTQR qr;
57    UNSQR uqr;
58    int signRem = 1;
59    int signQuot = 1;
60
61    if (m64IsNegative(num))
62    {
63        num = m64Negate(num);
64        signQuot = -signQuot;
65    }
66
67    if (den < 0)
68    {
69        den      = -den;
70        signRem  = -signRem;
71        signQuot = -signQuot;
72    }
73
74    uqr = ficlLongDiv(m64CastIU(num), (UNS32)den);
75    qr = m64CastQRUI(uqr);
76    if (signQuot < 0)
77    {
78        qr.quot = -qr.quot;
79        if (qr.rem != 0)
80        {
81            qr.quot--;
82            qr.rem = den - qr.rem;
83        }
84    }
85
86    if (signRem < 0)
87        qr.rem = -qr.rem;
88
89    return qr;
90}
91
92
93/**************************************************************************
94                        m 6 4 I s N e g a t i v e
95** Returns TRUE if the specified INT64 has its sign bit set.
96**************************************************************************/
97int m64IsNegative(INT64 x)
98{
99    return (x.hi < 0);
100}
101
102
103/**************************************************************************
104                        m 6 4 M a c
105** Mixed precision multiply and accumulate primitive for number building.
106** Multiplies UNS64 u by UNS32 mul and adds UNS32 add. Mul is typically
107** the numeric base, and add represents a digit to be appended to the
108** growing number.
109** Returns the result of the operation
110**************************************************************************/
111UNS64 m64Mac(UNS64 u, UNS32 mul, UNS32 add)
112{
113    UNS64 resultLo = ficlLongMul(u.lo, mul);
114    UNS64 resultHi = ficlLongMul(u.hi, mul);
115    resultLo.hi += resultHi.lo;
116    resultHi.lo = resultLo.lo + add;
117
118    if (resultHi.lo < resultLo.lo)
119        resultLo.hi++;
120
121    resultLo.lo = resultHi.lo;
122
123    return resultLo;
124}
125
126
127/**************************************************************************
128                        m 6 4 M u l I
129** Multiplies a pair of INT32s and returns an INT64 result.
130**************************************************************************/
131INT64 m64MulI(INT32 x, INT32 y)
132{
133    UNS64 prod;
134    int sign = 1;
135
136    if (x < 0)
137    {
138        sign = -sign;
139        x = -x;
140    }
141
142    if (y < 0)
143    {
144        sign = -sign;
145        y = -y;
146    }
147
148    prod = ficlLongMul(x, y);
149    if (sign > 0)
150        return m64CastUI(prod);
151    else
152        return m64Negate(m64CastUI(prod));
153}
154
155
156/**************************************************************************
157                        m 6 4 N e g a t e
158** Negates an INT64 by complementing and incrementing.
159**************************************************************************/
160INT64 m64Negate(INT64 x)
161{
162    x.hi = ~x.hi;
163    x.lo = ~x.lo;
164    x.lo ++;
165    if (x.lo == 0)
166        x.hi++;
167
168    return x;
169}
170
171
172/**************************************************************************
173                        m 6 4 P u s h
174** Push an INT64 onto the specified stack in the order required
175** by ANS Forth (most significant cell on top)
176** These should probably be macros...
177**************************************************************************/
178void  i64Push(FICL_STACK *pStack, INT64 i64)
179{
180    stackPushINT32(pStack, i64.lo);
181    stackPushINT32(pStack, i64.hi);
182    return;
183}
184
185void  u64Push(FICL_STACK *pStack, UNS64 u64)
186{
187    stackPushINT32(pStack, u64.lo);
188    stackPushINT32(pStack, u64.hi);
189    return;
190}
191
192
193/**************************************************************************
194                        m 6 4 P o p
195** Pops an INT64 off the stack in the order required by ANS Forth
196** (most significant cell on top)
197** These should probably be macros...
198**************************************************************************/
199INT64 i64Pop(FICL_STACK *pStack)
200{
201    INT64 ret;
202    ret.hi = stackPopINT32(pStack);
203    ret.lo = stackPopINT32(pStack);
204    return ret;
205}
206
207UNS64 u64Pop(FICL_STACK *pStack)
208{
209    UNS64 ret;
210    ret.hi = stackPopINT32(pStack);
211    ret.lo = stackPopINT32(pStack);
212    return ret;
213}
214
215
216/**************************************************************************
217                        m 6 4 S y m m e t r i c D i v
218** Divide an INT64 by an INT32 and return an INT32 quotient and an INT32
219** remainder. The absolute values of quotient and remainder are not
220** affected by the signs of the numerator and denominator (the operation
221** is symmetric on the number line)
222**************************************************************************/
223INTQR m64SymmetricDivI(INT64 num, INT32 den)
224{
225    INTQR qr;
226    UNSQR uqr;
227    int signRem = 1;
228    int signQuot = 1;
229
230    if (m64IsNegative(num))
231    {
232        num = m64Negate(num);
233        signRem  = -signRem;
234        signQuot = -signQuot;
235    }
236
237    if (den < 0)
238    {
239        den      = -den;
240        signQuot = -signQuot;
241    }
242
243    uqr = ficlLongDiv(m64CastIU(num), (UNS32)den);
244    qr = m64CastQRUI(uqr);
245    if (signRem < 0)
246        qr.rem = -qr.rem;
247
248    if (signQuot < 0)
249        qr.quot = -qr.quot;
250
251    return qr;
252}
253
254
255/**************************************************************************
256                        m 6 4 U M o d
257** Divides an UNS64 by base (an UNS16) and returns an UNS16 remainder.
258** Writes the quotient back to the original UNS64 as a side effect.
259** This operation is typically used to convert an UNS64 to a text string
260** in any base. See words.c:numberSignS, for example.
261** Mechanics: performs 4 ficlLongDivs, each of which produces 16 bits
262** of the quotient. C does not provide a way to divide an UNS32 by an
263** UNS16 and get an UNS32 quotient (ldiv is closest, but it's signed,
264** unfortunately), so I've used ficlLongDiv.
265**************************************************************************/
266UNS16 m64UMod(UNS64 *pUD, UNS16 base)
267{
268    UNS64 ud;
269    UNSQR qr;
270    UNS64 result;
271
272    result.hi = result.lo = 0;
273
274    ud.hi = 0;
275    ud.lo = pUD->hi >> 16;
276    qr = ficlLongDiv(ud, (UNS32)base);
277    result.hi = qr.quot << 16;
278
279    ud.lo = (qr.rem << 16) | (pUD->hi & 0x0000ffff);
280    qr = ficlLongDiv(ud, (UNS32)base);
281    result.hi |= qr.quot & 0x0000ffff;
282
283    ud.lo = (qr.rem << 16) | (pUD->lo >> 16);
284    qr = ficlLongDiv(ud, (UNS32)base);
285    result.lo = qr.quot << 16;
286
287    ud.lo = (qr.rem << 16) | (pUD->lo & 0x0000ffff);
288    qr = ficlLongDiv(ud, (UNS32)base);
289    result.lo |= qr.quot & 0x0000ffff;
290
291    *pUD = result;
292
293    return (UNS16)(qr.rem);
294}
295
296
297