/* * 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); }