xref: /illumos-gate/usr/src/common/ficl/double.c (revision afc2ba1d)

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