math64.c

/*******************************************************************
** m a t h 6 4 . c
** Forth Inspired Command Language - 64 bit math support routines
** Author: John Sadler ([email protected])
** Created: 25 January 1998
** Rev 2.03: Support for 128 bit DP math. This file really ouught to
** be renamed!
** $Id: math64.c,v 1.9 2001/12/05 07:21:34 jsadler Exp $
*******************************************************************/
/*
** Copyright (c) 1997-2001 John Sadler ([email protected])
** 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"
#include "math64.h"


/**************************************************************************
                        m 6 4 A b s
** Returns the absolute value of an DPINT
**************************************************************************/
DPINT m64Abs(DPINT x)
{
    if (m64IsNegative(x))
        x = m64Negate(x);

    return x;
}


/**************************************************************************
                        m 6 4 F l o o r e d D i v I
** 
** 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
**************************************************************************/
INTQR m64FlooredDivI(DPINT num, FICL_INT den)
{
    INTQR qr;
    UNSQR uqr;
    int signRem = 1;
    int signQuot = 1;

    if (m64IsNegative(num))
    {
        num = m64Negate(num);
        signQuot = -signQuot;
    }

    if (den < 0)
    {
        den      = -den;
        signRem  = -signRem;
        signQuot = -signQuot;
    }

    uqr = ficlLongDiv(m64CastIU(num), (FICL_UNS)den);
    qr = m64CastQRUI(uqr);
    if (signQuot < 0)
    {
        qr.quot = -qr.quot;
        if (qr.rem != 0)
        {
            qr.quot--;
            qr.rem = den - qr.rem;
        }
    }

    if (signRem < 0)
        qr.rem = -qr.rem;

    return qr;
}


/**************************************************************************
                        m 6 4 I s N e g a t i v e
** Returns TRUE if the specified DPINT has its sign bit set.
**************************************************************************/
int m64IsNegative(DPINT x)
{
    return (x.hi < 0);
}


/**************************************************************************
                        m 6 4 M a c
** Mixed precision multiply and accumulate primitive for number building.
** Multiplies DPUNS u by FICL_UNS mul and adds FICL_UNS 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
**************************************************************************/
DPUNS m64Mac(DPUNS u, FICL_UNS mul, FICL_UNS add)
{
    DPUNS resultLo = ficlLongMul(u.lo, mul);
    DPUNS resultHi = ficlLongMul(u.hi, mul);
    resultLo.hi += resultHi.lo;
    resultHi.lo = resultLo.lo + add;

    if (resultHi.lo < resultLo.lo)
        resultLo.hi++;

    resultLo.lo = resultHi.lo;

    return resultLo;
}


/**************************************************************************
                        m 6 4 M u l I
** Multiplies a pair of FICL_INTs and returns an DPINT result.
**************************************************************************/
DPINT m64MulI(FICL_INT x, FICL_INT y)
{
    DPUNS prod;
    int sign = 1;

    if (x < 0)
    {
        sign = -sign;
        x = -x;
    }

    if (y < 0)
    {
        sign = -sign;
        y = -y;
    }

    prod = ficlLongMul(x, y);
    if (sign > 0)
        return m64CastUI(prod);
    else
        return m64Negate(m64CastUI(prod));
}


/**************************************************************************
                        m 6 4 N e g a t e
** Negates an DPINT by complementing and incrementing.
**************************************************************************/
DPINT m64Negate(DPINT x)
{
    x.hi = ~x.hi;
    x.lo = ~x.lo;
    x.lo ++;
    if (x.lo == 0)
        x.hi++;

    return x;
}


/**************************************************************************
                        m 6 4 P u s h
** Push an DPINT onto the specified stack in the order required
** by ANS Forth (most significant cell on top)
** These should probably be macros...
**************************************************************************/
void  i64Push(FICL_STACK *pStack, DPINT i64)
{
    stackPushINT(pStack, i64.lo);
    stackPushINT(pStack, i64.hi);
    return;
}

void  u64Push(FICL_STACK *pStack, DPUNS u64)
{
    stackPushINT(pStack, u64.lo);
    stackPushINT(pStack, u64.hi);
    return;
}


/**************************************************************************
                        m 6 4 P o p
** Pops an DPINT off the stack in the order required by ANS Forth
** (most significant cell on top)
** These should probably be macros...
**************************************************************************/
DPINT i64Pop(FICL_STACK *pStack)
{
    DPINT ret;
    ret.hi = stackPopINT(pStack);
    ret.lo = stackPopINT(pStack);
    return ret;
}

DPUNS u64Pop(FICL_STACK *pStack)
{
    DPUNS ret;
    ret.hi = stackPopINT(pStack);
    ret.lo = stackPopINT(pStack);
    return ret;
}


/**************************************************************************
                        m 6 4 S y m m e t r i c D i v
** Divide an DPINT by a FICL_INT and return a FICL_INT quotient and a
** FICL_INT 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)
**************************************************************************/
INTQR m64SymmetricDivI(DPINT num, FICL_INT den)
{
    INTQR qr;
    UNSQR uqr;
    int signRem = 1;
    int signQuot = 1;

    if (m64IsNegative(num))
    {
        num = m64Negate(num);
        signRem  = -signRem;
        signQuot = -signQuot;
    }

    if (den < 0)
    {
        den      = -den;
        signQuot = -signQuot;
    }

    uqr = ficlLongDiv(m64CastIU(num), (FICL_UNS)den);
    qr = m64CastQRUI(uqr);
    if (signRem < 0)
        qr.rem = -qr.rem;

    if (signQuot < 0)
        qr.quot = -qr.quot;

    return qr;
}


/**************************************************************************
                        m 6 4 U M o d
** Divides a DPUNS by base (an UNS16) and returns an UNS16 remainder.
** Writes the quotient back to the original DPUNS as a side effect.
** This operation is typically used to convert an DPUNS to a text string
** in any base. See words.c:numberSignS, for example.
** Mechanics: performs 4 ficlLongDivs, each of which produces 16 bits
** of the quotient. C does not provide a way to divide an FICL_UNS by an
** UNS16 and get an FICL_UNS quotient (ldiv is closest, but it's signed,
** unfortunately), so I've used ficlLongDiv.
**************************************************************************/
#if (BITS_PER_CELL == 32)

#define UMOD_SHIFT 16
#define UMOD_MASK 0x0000ffff

#elif (BITS_PER_CELL == 64)

#define UMOD_SHIFT 32
#define UMOD_MASK 0x00000000ffffffff

#endif

UNS16 m64UMod(DPUNS *pUD, UNS16 base)
{
    DPUNS ud;
    UNSQR qr;
    DPUNS result;

    result.hi = result.lo = 0;

    ud.hi = 0;
    ud.lo = pUD->hi >> UMOD_SHIFT;
    qr = ficlLongDiv(ud, (FICL_UNS)base);
    result.hi = qr.quot << UMOD_SHIFT;

    ud.lo = (qr.rem << UMOD_SHIFT) | (pUD->hi & UMOD_MASK);
    qr = ficlLongDiv(ud, (FICL_UNS)base);
    result.hi |= qr.quot & UMOD_MASK;

    ud.lo = (qr.rem << UMOD_SHIFT) | (pUD->lo >> UMOD_SHIFT);
    qr = ficlLongDiv(ud, (FICL_UNS)base);
    result.lo = qr.quot << UMOD_SHIFT;

    ud.lo = (qr.rem << UMOD_SHIFT) | (pUD->lo & UMOD_MASK);
    qr = ficlLongDiv(ud, (FICL_UNS)base);
    result.lo |= qr.quot & UMOD_MASK;

    *pUD = result;

    return (UNS16)(qr.rem);
}


/**************************************************************************
** Contributed by
** Michael A. Gauland   [email protected]  
**************************************************************************/
#if PORTABLE_LONGMULDIV != 0
/**************************************************************************
                        m 6 4 A d d
** 
**************************************************************************/
DPUNS m64Add(DPUNS x, DPUNS y)
{
    DPUNS result;
    int carry;
    
    result.hi = x.hi + y.hi;
    result.lo = x.lo + y.lo;


    carry  = ((x.lo | y.lo) & CELL_HI_BIT) && !(result.lo & CELL_HI_BIT);
    carry |= ((x.lo & y.lo) & CELL_HI_BIT);

    if (carry)
    {
        result.hi++;
    }

    return result;
}


/**************************************************************************
                        m 6 4 S u b
** 
**************************************************************************/
DPUNS m64Sub(DPUNS x, DPUNS y)
{
    DPUNS result;
    
    result.hi = x.hi - y.hi;
    result.lo = x.lo - y.lo;

    if (x.lo < y.lo) 
    {
        result.hi--;
    }

    return result;
}


/**************************************************************************
                        m 6 4 A S L
** 64 bit left shift
**************************************************************************/
DPUNS m64ASL( DPUNS x )
{
    DPUNS result;
    
    result.hi = x.hi << 1;
    if (x.lo & CELL_HI_BIT) 
    {
        result.hi++;
    }

    result.lo = x.lo << 1;

    return result;
}


/**************************************************************************
                        m 6 4 A S R
** 64 bit right shift (unsigned - no sign extend)
**************************************************************************/
DPUNS m64ASR( DPUNS x )
{
    DPUNS result;
    
    result.lo = x.lo >> 1;
    if (x.hi & 1) 
    {
        result.lo |= CELL_HI_BIT;
    }

    result.hi = x.hi >> 1;
    return result;
}


/**************************************************************************
                        m 6 4 O r
** 64 bit bitwise OR
**************************************************************************/
DPUNS m64Or( DPUNS x, DPUNS y )
{
    DPUNS result;
    
    result.hi = x.hi | y.hi;
    result.lo = x.lo | y.lo;
    
    return result;
}


/**************************************************************************
                        m 6 4 C o m p a r e
** Return -1 if x < y; 0 if x==y, and 1 if x > y.
**************************************************************************/
int m64Compare(DPUNS x, DPUNS y)
{
    int result;
    
    if (x.hi > y.hi) 
    {
        result = +1;
    } 
    else if (x.hi < y.hi) 
    {
        result = -1;
    } 
    else 
    {
        /* High parts are equal */
        if (x.lo > y.lo) 
        {
            result = +1;
        } 
        else if (x.lo < y.lo) 
        {
            result = -1;
        } 
        else 
        {
            result = 0;
        }
    }
    
    return result;
}


/**************************************************************************
                        f i c l L o n g M u l
** Portable versions of ficlLongMul and ficlLongDiv in C
** Contributed by:
** Michael A. Gauland   [email protected]  
**************************************************************************/
DPUNS ficlLongMul(FICL_UNS x, FICL_UNS y)
{
    DPUNS result = { 0, 0 };
    DPUNS addend;
    
    addend.lo = y;
    addend.hi = 0; /* No sign extension--arguments are unsigned */
    
    while (x != 0) 
    {
        if ( x & 1) 
        {
            result = m64Add(result, addend);
        }
        x >>= 1;
        addend = m64ASL(addend);
    }
    return result;
}


/**************************************************************************
                        f i c l L o n g D i v
** Portable versions of ficlLongMul and ficlLongDiv in C
** Contributed by:
** Michael A. Gauland   [email protected]  
**************************************************************************/
UNSQR ficlLongDiv(DPUNS q, FICL_UNS y)
{
    UNSQR result;
    DPUNS quotient;
    DPUNS subtrahend;
    DPUNS mask;

    quotient.lo = 0;
    quotient.hi = 0;
    
    subtrahend.lo = y;
    subtrahend.hi = 0;
    
    mask.lo = 1;
    mask.hi = 0;
    
    while ((m64Compare(subtrahend, q) < 0) &&
           (subtrahend.hi & CELL_HI_BIT) == 0)
    {
        mask = m64ASL(mask);
        subtrahend = m64ASL(subtrahend);
    }
    
    while (mask.lo != 0 || mask.hi != 0) 
    {
        if (m64Compare(subtrahend, q) <= 0) 
        {
            q = m64Sub( q, subtrahend);
            quotient = m64Or(quotient, mask);
        }
        mask = m64ASR(mask);
        subtrahend = m64ASR(subtrahend);
    }
    
    result.quot = quotient.lo;
    result.rem = q.lo;
    return result;
}

#endif