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