summaryrefslogtreecommitdiffstats
path: root/stand/ficl/math64.c
diff options
context:
space:
mode:
Diffstat (limited to 'stand/ficl/math64.c')
-rw-r--r--stand/ficl/math64.c561
1 files changed, 561 insertions, 0 deletions
diff --git a/stand/ficl/math64.c b/stand/ficl/math64.c
new file mode 100644
index 0000000..6e50458
--- /dev/null
+++ b/stand/ficl/math64.c
@@ -0,0 +1,561 @@
+/*******************************************************************
+** m a t h 6 4 . c
+** Forth Inspired Command Language - 64 bit math support routines
+** Author: John Sadler (john_sadler@alum.mit.edu)
+** Created: 25 January 1998
+** Rev 2.03: Support for 128 bit DP math. This file really ouught to
+** be renamed!
+** $Id: math64.c,v 1.9 2001/12/05 07:21:34 jsadler Exp $
+*******************************************************************/
+/*
+** Copyright (c) 1997-2001 John Sadler (john_sadler@alum.mit.edu)
+** All rights reserved.
+**
+** Get the latest Ficl release at http://ficl.sourceforge.net
+**
+** I am interested in hearing from anyone who uses ficl. If you have
+** a problem, a success story, a defect, an enhancement request, or
+** if you would like to contribute to the ficl release, please
+** contact me by email at the address above.
+**
+** L I C E N S E and D I S C L A I M E R
+**
+** Redistribution and use in source and binary forms, with or without
+** modification, are permitted provided that the following conditions
+** are met:
+** 1. Redistributions of source code must retain the above copyright
+** notice, this list of conditions and the following disclaimer.
+** 2. Redistributions in binary form must reproduce the above copyright
+** notice, this list of conditions and the following disclaimer in the
+** documentation and/or other materials provided with the distribution.
+**
+** THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND
+** ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+** IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
+** ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE
+** FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
+** DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
+** OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
+** HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
+** LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
+** OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
+** SUCH DAMAGE.
+*/
+
+/* $FreeBSD$ */
+
+#include "ficl.h"
+#include "math64.h"
+
+
+/**************************************************************************
+ m 6 4 A b s
+** Returns the absolute value of an DPINT
+**************************************************************************/
+DPINT m64Abs(DPINT x)
+{
+ if (m64IsNegative(x))
+ x = m64Negate(x);
+
+ return x;
+}
+
+
+/**************************************************************************
+ m 6 4 F l o o r e d D i v I
+**
+** FROM THE FORTH ANS...
+** Floored division is integer division in which the remainder carries
+** the sign of the divisor or is zero, and the quotient is rounded to
+** its arithmetic floor. Symmetric division is integer division in which
+** the remainder carries the sign of the dividend or is zero and the
+** quotient is the mathematical quotient rounded towards zero or
+** truncated. Examples of each are shown in tables 3.3 and 3.4.
+**
+** Table 3.3 - Floored Division Example
+** Dividend Divisor Remainder Quotient
+** -------- ------- --------- --------
+** 10 7 3 1
+** -10 7 4 -2
+** 10 -7 -4 -2
+** -10 -7 -3 1
+**
+**
+** Table 3.4 - Symmetric Division Example
+** Dividend Divisor Remainder Quotient
+** -------- ------- --------- --------
+** 10 7 3 1
+** -10 7 -3 -1
+** 10 -7 3 -1
+** -10 -7 -3 1
+**************************************************************************/
+INTQR m64FlooredDivI(DPINT num, FICL_INT den)
+{
+ INTQR qr;
+ UNSQR uqr;
+ int signRem = 1;
+ int signQuot = 1;
+
+ if (m64IsNegative(num))
+ {
+ num = m64Negate(num);
+ signQuot = -signQuot;
+ }
+
+ if (den < 0)
+ {
+ den = -den;
+ signRem = -signRem;
+ signQuot = -signQuot;
+ }
+
+ uqr = ficlLongDiv(m64CastIU(num), (FICL_UNS)den);
+ qr = m64CastQRUI(uqr);
+ if (signQuot < 0)
+ {
+ qr.quot = -qr.quot;
+ if (qr.rem != 0)
+ {
+ qr.quot--;
+ qr.rem = den - qr.rem;
+ }
+ }
+
+ if (signRem < 0)
+ qr.rem = -qr.rem;
+
+ return qr;
+}
+
+
+/**************************************************************************
+ m 6 4 I s N e g a t i v e
+** Returns TRUE if the specified DPINT has its sign bit set.
+**************************************************************************/
+int m64IsNegative(DPINT x)
+{
+ return (x.hi < 0);
+}
+
+
+/**************************************************************************
+ m 6 4 M a c
+** Mixed precision multiply and accumulate primitive for number building.
+** Multiplies DPUNS u by FICL_UNS mul and adds FICL_UNS add. Mul is typically
+** the numeric base, and add represents a digit to be appended to the
+** growing number.
+** Returns the result of the operation
+**************************************************************************/
+DPUNS m64Mac(DPUNS u, FICL_UNS mul, FICL_UNS add)
+{
+ DPUNS resultLo = ficlLongMul(u.lo, mul);
+ DPUNS resultHi = ficlLongMul(u.hi, mul);
+ resultLo.hi += resultHi.lo;
+ resultHi.lo = resultLo.lo + add;
+
+ if (resultHi.lo < resultLo.lo)
+ resultLo.hi++;
+
+ resultLo.lo = resultHi.lo;
+
+ return resultLo;
+}
+
+
+/**************************************************************************
+ m 6 4 M u l I
+** Multiplies a pair of FICL_INTs and returns an DPINT result.
+**************************************************************************/
+DPINT m64MulI(FICL_INT x, FICL_INT y)
+{
+ DPUNS prod;
+ int sign = 1;
+
+ if (x < 0)
+ {
+ sign = -sign;
+ x = -x;
+ }
+
+ if (y < 0)
+ {
+ sign = -sign;
+ y = -y;
+ }
+
+ prod = ficlLongMul(x, y);
+ if (sign > 0)
+ return m64CastUI(prod);
+ else
+ return m64Negate(m64CastUI(prod));
+}
+
+
+/**************************************************************************
+ m 6 4 N e g a t e
+** Negates an DPINT by complementing and incrementing.
+**************************************************************************/
+DPINT m64Negate(DPINT x)
+{
+ x.hi = ~x.hi;
+ x.lo = ~x.lo;
+ x.lo ++;
+ if (x.lo == 0)
+ x.hi++;
+
+ return x;
+}
+
+
+/**************************************************************************
+ m 6 4 P u s h
+** Push an DPINT onto the specified stack in the order required
+** by ANS Forth (most significant cell on top)
+** These should probably be macros...
+**************************************************************************/
+void i64Push(FICL_STACK *pStack, DPINT i64)
+{
+ stackPushINT(pStack, i64.lo);
+ stackPushINT(pStack, i64.hi);
+ return;
+}
+
+void u64Push(FICL_STACK *pStack, DPUNS u64)
+{
+ stackPushINT(pStack, u64.lo);
+ stackPushINT(pStack, u64.hi);
+ return;
+}
+
+
+/**************************************************************************
+ m 6 4 P o p
+** Pops an DPINT off the stack in the order required by ANS Forth
+** (most significant cell on top)
+** These should probably be macros...
+**************************************************************************/
+DPINT i64Pop(FICL_STACK *pStack)
+{
+ DPINT ret;
+ ret.hi = stackPopINT(pStack);
+ ret.lo = stackPopINT(pStack);
+ return ret;
+}
+
+DPUNS u64Pop(FICL_STACK *pStack)
+{
+ DPUNS ret;
+ ret.hi = stackPopINT(pStack);
+ ret.lo = stackPopINT(pStack);
+ return ret;
+}
+
+
+/**************************************************************************
+ m 6 4 S y m m e t r i c D i v
+** Divide an DPINT by a FICL_INT and return a FICL_INT quotient and a
+** FICL_INT remainder. The absolute values of quotient and remainder are not
+** affected by the signs of the numerator and denominator (the operation
+** is symmetric on the number line)
+**************************************************************************/
+INTQR m64SymmetricDivI(DPINT num, FICL_INT den)
+{
+ INTQR qr;
+ UNSQR uqr;
+ int signRem = 1;
+ int signQuot = 1;
+
+ if (m64IsNegative(num))
+ {
+ num = m64Negate(num);
+ signRem = -signRem;
+ signQuot = -signQuot;
+ }
+
+ if (den < 0)
+ {
+ den = -den;
+ signQuot = -signQuot;
+ }
+
+ uqr = ficlLongDiv(m64CastIU(num), (FICL_UNS)den);
+ qr = m64CastQRUI(uqr);
+ if (signRem < 0)
+ qr.rem = -qr.rem;
+
+ if (signQuot < 0)
+ qr.quot = -qr.quot;
+
+ return qr;
+}
+
+
+/**************************************************************************
+ m 6 4 U M o d
+** Divides a DPUNS by base (an UNS16) and returns an UNS16 remainder.
+** Writes the quotient back to the original DPUNS as a side effect.
+** This operation is typically used to convert an DPUNS to a text string
+** in any base. See words.c:numberSignS, for example.
+** Mechanics: performs 4 ficlLongDivs, each of which produces 16 bits
+** of the quotient. C does not provide a way to divide an FICL_UNS by an
+** UNS16 and get an FICL_UNS quotient (ldiv is closest, but it's signed,
+** unfortunately), so I've used ficlLongDiv.
+**************************************************************************/
+#if (BITS_PER_CELL == 32)
+
+#define UMOD_SHIFT 16
+#define UMOD_MASK 0x0000ffff
+
+#elif (BITS_PER_CELL == 64)
+
+#define UMOD_SHIFT 32
+#define UMOD_MASK 0x00000000ffffffff
+
+#endif
+
+UNS16 m64UMod(DPUNS *pUD, UNS16 base)
+{
+ DPUNS ud;
+ UNSQR qr;
+ DPUNS result;
+
+ result.hi = result.lo = 0;
+
+ ud.hi = 0;
+ ud.lo = pUD->hi >> UMOD_SHIFT;
+ qr = ficlLongDiv(ud, (FICL_UNS)base);
+ result.hi = qr.quot << UMOD_SHIFT;
+
+ ud.lo = (qr.rem << UMOD_SHIFT) | (pUD->hi & UMOD_MASK);
+ qr = ficlLongDiv(ud, (FICL_UNS)base);
+ result.hi |= qr.quot & UMOD_MASK;
+
+ ud.lo = (qr.rem << UMOD_SHIFT) | (pUD->lo >> UMOD_SHIFT);
+ qr = ficlLongDiv(ud, (FICL_UNS)base);
+ result.lo = qr.quot << UMOD_SHIFT;
+
+ ud.lo = (qr.rem << UMOD_SHIFT) | (pUD->lo & UMOD_MASK);
+ qr = ficlLongDiv(ud, (FICL_UNS)base);
+ result.lo |= qr.quot & UMOD_MASK;
+
+ *pUD = result;
+
+ return (UNS16)(qr.rem);
+}
+
+
+/**************************************************************************
+** Contributed by
+** Michael A. Gauland gaulandm@mdhost.cse.tek.com
+**************************************************************************/
+#if PORTABLE_LONGMULDIV != 0
+/**************************************************************************
+ m 6 4 A d d
+**
+**************************************************************************/
+DPUNS m64Add(DPUNS x, DPUNS y)
+{
+ DPUNS result;
+ int carry;
+
+ result.hi = x.hi + y.hi;
+ result.lo = x.lo + y.lo;
+
+
+ carry = ((x.lo | y.lo) & CELL_HI_BIT) && !(result.lo & CELL_HI_BIT);
+ carry |= ((x.lo & y.lo) & CELL_HI_BIT);
+
+ if (carry)
+ {
+ result.hi++;
+ }
+
+ return result;
+}
+
+
+/**************************************************************************
+ m 6 4 S u b
+**
+**************************************************************************/
+DPUNS m64Sub(DPUNS x, DPUNS y)
+{
+ DPUNS result;
+
+ result.hi = x.hi - y.hi;
+ result.lo = x.lo - y.lo;
+
+ if (x.lo < y.lo)
+ {
+ result.hi--;
+ }
+
+ return result;
+}
+
+
+/**************************************************************************
+ m 6 4 A S L
+** 64 bit left shift
+**************************************************************************/
+DPUNS m64ASL( DPUNS x )
+{
+ DPUNS result;
+
+ result.hi = x.hi << 1;
+ if (x.lo & CELL_HI_BIT)
+ {
+ result.hi++;
+ }
+
+ result.lo = x.lo << 1;
+
+ return result;
+}
+
+
+/**************************************************************************
+ m 6 4 A S R
+** 64 bit right shift (unsigned - no sign extend)
+**************************************************************************/
+DPUNS m64ASR( DPUNS x )
+{
+ DPUNS result;
+
+ result.lo = x.lo >> 1;
+ if (x.hi & 1)
+ {
+ result.lo |= CELL_HI_BIT;
+ }
+
+ result.hi = x.hi >> 1;
+ return result;
+}
+
+
+/**************************************************************************
+ m 6 4 O r
+** 64 bit bitwise OR
+**************************************************************************/
+DPUNS m64Or( DPUNS x, DPUNS y )
+{
+ DPUNS result;
+
+ result.hi = x.hi | y.hi;
+ result.lo = x.lo | y.lo;
+
+ return result;
+}
+
+
+/**************************************************************************
+ m 6 4 C o m p a r e
+** Return -1 if x < y; 0 if x==y, and 1 if x > y.
+**************************************************************************/
+int m64Compare(DPUNS x, DPUNS y)
+{
+ int result;
+
+ if (x.hi > y.hi)
+ {
+ result = +1;
+ }
+ else if (x.hi < y.hi)
+ {
+ result = -1;
+ }
+ else
+ {
+ /* High parts are equal */
+ if (x.lo > y.lo)
+ {
+ result = +1;
+ }
+ else if (x.lo < y.lo)
+ {
+ result = -1;
+ }
+ else
+ {
+ result = 0;
+ }
+ }
+
+ return result;
+}
+
+
+/**************************************************************************
+ f i c l L o n g M u l
+** Portable versions of ficlLongMul and ficlLongDiv in C
+** Contributed by:
+** Michael A. Gauland gaulandm@mdhost.cse.tek.com
+**************************************************************************/
+DPUNS ficlLongMul(FICL_UNS x, FICL_UNS y)
+{
+ DPUNS result = { 0, 0 };
+ DPUNS addend;
+
+ addend.lo = y;
+ addend.hi = 0; /* No sign extension--arguments are unsigned */
+
+ while (x != 0)
+ {
+ if ( x & 1)
+ {
+ result = m64Add(result, addend);
+ }
+ x >>= 1;
+ addend = m64ASL(addend);
+ }
+ return result;
+}
+
+
+/**************************************************************************
+ f i c l L o n g D i v
+** Portable versions of ficlLongMul and ficlLongDiv in C
+** Contributed by:
+** Michael A. Gauland gaulandm@mdhost.cse.tek.com
+**************************************************************************/
+UNSQR ficlLongDiv(DPUNS q, FICL_UNS y)
+{
+ UNSQR result;
+ DPUNS quotient;
+ DPUNS subtrahend;
+ DPUNS mask;
+
+ quotient.lo = 0;
+ quotient.hi = 0;
+
+ subtrahend.lo = y;
+ subtrahend.hi = 0;
+
+ mask.lo = 1;
+ mask.hi = 0;
+
+ while ((m64Compare(subtrahend, q) < 0) &&
+ (subtrahend.hi & CELL_HI_BIT) == 0)
+ {
+ mask = m64ASL(mask);
+ subtrahend = m64ASL(subtrahend);
+ }
+
+ while (mask.lo != 0 || mask.hi != 0)
+ {
+ if (m64Compare(subtrahend, q) <= 0)
+ {
+ q = m64Sub( q, subtrahend);
+ quotient = m64Or(quotient, mask);
+ }
+ mask = m64ASR(mask);
+ subtrahend = m64ASR(subtrahend);
+ }
+
+ result.quot = quotient.lo;
+ result.rem = q.lo;
+ return result;
+}
+
+#endif
+
OpenPOWER on IntegriCloud