562 lines
15 KiB
C
562 lines
15 KiB
C
/*******************************************************************
|
|
** 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
|
|
|