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