/*
* m a t h 6 4 . c
* Forth Inspired Command Language - 64 bit math support routines
* Authors: Michael A. Gauland (gaulandm@mdhost.cse.tek.com)
* Larry Hastings (larry@hastings.org)
* 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: double.c,v 1.2 2010/09/12 15:18:07 asau 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.
*/
#include "ficl.h"
#if FICL_PLATFORM_HAS_2INTEGER
ficl2UnsignedQR
ficl2UnsignedDivide(ficl2Unsigned q, ficlUnsigned y)
{
ficl2UnsignedQR result;
result.quotient = q / y;
/*
* Once we have the quotient, it's cheaper to calculate the
* remainder this way than with % (mod). --lch
*/
result.remainder = (ficlInteger)(q - (result.quotient * y));
return (result);
}
#else /* FICL_PLATFORM_HAS_2INTEGER */
#define FICL_CELL_HIGH_BIT ((uintmax_t)1 << (FICL_BITS_PER_CELL-1))
#define UMOD_SHIFT (FICL_BITS_PER_CELL / 2)
#define UMOD_MASK ((1L << (FICL_BITS_PER_CELL / 2)) - 1)
/*
* ficl2IntegerIsNegative
* Returns TRUE if the specified ficl2Unsigned has its sign bit set.
*/
int
ficl2IntegerIsNegative(ficl2Integer x)
{
return (x.high < 0);
}
/*
* ficl2IntegerNegate
* Negates an ficl2Unsigned by complementing and incrementing.
*/
ficl2Integer
ficl2IntegerNegate(ficl2Integer x)
{
x.high = ~x.high;
x.low = ~x.low;
x.low ++;
if (x.low == 0)
x.high++;
return (x);
}
/*
* ficl2UnsignedMultiplyAccumulate
* Mixed precision multiply and accumulate primitive for number building.
* Multiplies ficl2Unsigned u by ficlUnsigned mul and adds ficlUnsigned 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
*/
ficl2Unsigned
ficl2UnsignedMultiplyAccumulate(ficl2Unsigned u, ficlUnsigned mul,
ficlUnsigned add)
{
ficl2Unsigned resultLo = ficl2UnsignedMultiply(u.low, mul);
ficl2Unsigned resultHi = ficl2UnsignedMultiply(u.high, mul);
resultLo.high += resultHi.low;
resultHi.low = resultLo.low + add;
if (resultHi.low < resultLo.low)
resultLo.high++;
resultLo.low = resultHi.low;
return (resultLo);
}
/*
* ficl2IntegerMultiply
* Multiplies a pair of ficlIntegers and returns an ficl2Integer result.
*/
ficl2Integer
ficl2IntegerMultiply(ficlInteger x, ficlInteger y)
{
ficl2Unsigned prod;
ficl2Integer result;
int sign = 1;
if (x < 0) {
sign = -sign;
x = -x;
}
if (y < 0) {
sign = -sign;
y = -y;
}
prod = ficl2UnsignedMultiply(x, y);
FICL_2INTEGER_SET(FICL_2UNSIGNED_GET_HIGH(prod),
FICL_2UNSIGNED_GET_LOW(prod), result);
if (sign > 0)
return (result);
else
return (ficl2IntegerNegate(result));
}
ficl2Integer
ficl2IntegerDecrement(ficl2Integer x)
{
if (x.low == INTMAX_MIN)
x.high--;
x.low--;
return (x);
}
ficl2Unsigned
ficl2UnsignedAdd(ficl2Unsigned x, ficl2Unsigned y)
{
ficl2Unsigned result;
result.high = x.high + y.high;
result.low = x.low + y.low;
if (result.low < y.low)
result.high++;
return (result);
}
/*
* ficl2UnsignedMultiply
* Contributed by:
* Michael A. Gauland gaulandm@mdhost.cse.tek.com
*/
ficl2Unsigned
ficl2UnsignedMultiply(ficlUnsigned x, ficlUnsigned y)
{
ficl2Unsigned result = { 0, 0 };
ficl2Unsigned addend;
addend.low = y;
addend.high = 0; /* No sign extension--arguments are unsigned */
while (x != 0) {
if (x & 1) {
result = ficl2UnsignedAdd(result, addend);
}
x >>= 1;
addend = ficl2UnsignedArithmeticShiftLeft(addend);
}
return (result);
}
/*
* ficl2UnsignedSubtract
*/
ficl2Unsigned
ficl2UnsignedSubtract(ficl2Unsigned x, ficl2Unsigned y)
{
ficl2Unsigned result;
result.high = x.high - y.high;
result.low = x.low - y.low;
if (x.low < y.low) {
result.high--;
}
return (result);
}
/*
* ficl2UnsignedArithmeticShiftLeft
* 64 bit left shift
*/
ficl2Unsigned
ficl2UnsignedArithmeticShiftLeft(ficl2Unsigned x)
{
ficl2Unsigned result;
result.high = x.high << 1;
if (x.low & FICL_CELL_HIGH_BIT) {
result.high++;
}
result.low = x.low << 1;
return (result);
}
/*
* ficl2UnsignedArithmeticShiftRight
* 64 bit right shift (unsigned - no sign extend)
*/
ficl2Unsigned
ficl2UnsignedArithmeticShiftRight(ficl2Unsigned x)
{
ficl2Unsigned result;
result.low = x.low >> 1;
if (x.high & 1) {
result.low |= FICL_CELL_HIGH_BIT;
}
result.high = x.high >> 1;
return (result);
}
/*
* ficl2UnsignedOr
* 64 bit bitwise OR
*/
ficl2Unsigned
ficl2UnsignedOr(ficl2Unsigned x, ficl2Unsigned y)
{
ficl2Unsigned result;
result.high = x.high | y.high;
result.low = x.low | y.low;
return (result);
}
/*
* ficl2UnsignedCompare
* Return -1 if x < y; 0 if x==y, and 1 if x > y.
*/
int
ficl2UnsignedCompare(ficl2Unsigned x, ficl2Unsigned y)
{
if (x.high > y.high)
return (1);
if (x.high < y.high)
return (-1);
/* High parts are equal */
if (x.low > y.low)
return (1);
else if (x.low < y.low)
return (-1);
return (0);
}
/*
* ficl2UnsignedDivide
* Portable versions of ficl2Multiply and ficl2Divide in C
* Contributed by:
* Michael A. Gauland gaulandm@mdhost.cse.tek.com
*/
ficl2UnsignedQR
ficl2UnsignedDivide(ficl2Unsigned q, ficlUnsigned y)
{
ficl2UnsignedQR result;
ficl2Unsigned quotient;
ficl2Unsigned subtrahend;
ficl2Unsigned mask;
quotient.low = 0;
quotient.high = 0;
subtrahend.low = y;
subtrahend.high = 0;
mask.low = 1;
mask.high = 0;
while ((ficl2UnsignedCompare(subtrahend, q) < 0) &&
(subtrahend.high & FICL_CELL_HIGH_BIT) == 0) {
mask = ficl2UnsignedArithmeticShiftLeft(mask);
subtrahend = ficl2UnsignedArithmeticShiftLeft(subtrahend);
}
while (mask.low != 0 || mask.high != 0) {
if (ficl2UnsignedCompare(subtrahend, q) <= 0) {
q = ficl2UnsignedSubtract(q, subtrahend);
quotient = ficl2UnsignedOr(quotient, mask);
}
mask = ficl2UnsignedArithmeticShiftRight(mask);
subtrahend = ficl2UnsignedArithmeticShiftRight(subtrahend);
}
result.quotient = quotient;
result.remainder = q.low;
return (result);
}
#endif /* !FICL_PLATFORM_HAS_2INTEGER */
/*
* ficl2IntegerDivideFloored
*
* 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
*/
ficl2IntegerQR
ficl2IntegerDivideFloored(ficl2Integer num, ficlInteger den)
{
ficl2IntegerQR qr;
ficl2UnsignedQR uqr;
ficl2Unsigned u;
int signRem = 1;
int signQuot = 1;
if (ficl2IntegerIsNegative(num)) {
num = ficl2IntegerNegate(num);
signQuot = -signQuot;
}
if (den < 0) {
den = -den;
signRem = -signRem;
signQuot = -signQuot;
}
FICL_2UNSIGNED_SET(FICL_2UNSIGNED_GET_HIGH(num),
FICL_2UNSIGNED_GET_LOW(num), u);
uqr = ficl2UnsignedDivide(u, (ficlUnsigned)den);
qr = FICL_2UNSIGNEDQR_TO_2INTEGERQR(uqr);
if (signQuot < 0) {
qr.quotient = ficl2IntegerNegate(qr.quotient);
if (qr.remainder != 0) {
qr.quotient = ficl2IntegerDecrement(qr.quotient);
qr.remainder = den - qr.remainder;
}
}
if (signRem < 0)
qr.remainder = -qr.remainder;
return (qr);
}
/*
* ficl2IntegerDivideSymmetric
* Divide an ficl2Unsigned by a ficlInteger and return a ficlInteger quotient
* and a ficlInteger 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)
*/
ficl2IntegerQR
ficl2IntegerDivideSymmetric(ficl2Integer num, ficlInteger den)
{
ficl2IntegerQR qr;
ficl2UnsignedQR uqr;
ficl2Unsigned u;
int signRem = 1;
int signQuot = 1;
if (ficl2IntegerIsNegative(num)) {
num = ficl2IntegerNegate(num);
signRem = -signRem;
signQuot = -signQuot;
}
if (den < 0) {
den = -den;
signQuot = -signQuot;
}
FICL_2UNSIGNED_SET(FICL_2UNSIGNED_GET_HIGH(num),
FICL_2UNSIGNED_GET_LOW(num), u);
uqr = ficl2UnsignedDivide(u, (ficlUnsigned)den);
qr = FICL_2UNSIGNEDQR_TO_2INTEGERQR(uqr);
if (signRem < 0)
qr.remainder = -qr.remainder;
if (signQuot < 0)
qr.quotient = ficl2IntegerNegate(qr.quotient);
return (qr);
}