ref: 7d02e382d314d5bdde7978ccb7a64ea9201d03db
dir: /double.c/
/******************************************************************* ** 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.3 2010/11/01 14:10:27 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 <stdint.h> #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; int sign = 1; if (x < 0) { sign = -sign; x = -x; } if (y < 0) { sign = -sign; y = -y; } prod = ficl2UnsignedMultiply(x, y); if (sign > 0) return FICL_2UNSIGNED_TO_2INTEGER(prod); else return ficl2IntegerNegate(FICL_2UNSIGNED_TO_2INTEGER(prod)); } ficl2Integer ficl2IntegerDecrement(ficl2Integer x) { if (x.low == INT_MIN) x.high--; x.low--; return x; } ficl2Unsigned ficl2UnsignedAdd(ficl2Unsigned x, ficl2Unsigned y) { ficl2Unsigned result; int carry; result.high = x.high + y.high; result.low = x.low + y.low; carry = ((x.low | y.low) & FICL_CELL_HIGH_BIT) && !(result.low & FICL_CELL_HIGH_BIT); carry |= ((x.low & y.low) & FICL_CELL_HIGH_BIT); if (carry) { 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 */ /************************************************************************** ficl2IntegerAbsoluteValue ** Returns the absolute value of an ficl2Unsigned **************************************************************************/ ficl2Integer ficl2IntegerAbsoluteValue(ficl2Integer x) { if (ficl2IntegerIsNegative(x)) return ficl2IntegerNegate(x); return x; } /************************************************************************** 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; int signRem = 1; int signQuot = 1; if (ficl2IntegerIsNegative(num)) { num = ficl2IntegerNegate(num); signQuot = -signQuot; } if (den < 0) { den = -den; signRem = -signRem; signQuot = -signQuot; } uqr = ficl2UnsignedDivide(FICL_2INTEGER_TO_2UNSIGNED(num), (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; int signRem = 1; int signQuot = 1; if (ficl2IntegerIsNegative(num)) { num = ficl2IntegerNegate(num); signRem = -signRem; signQuot = -signQuot; } if (den < 0) { den = -den; signQuot = -signQuot; } uqr = ficl2UnsignedDivide(FICL_2INTEGER_TO_2UNSIGNED(num), (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; }