text: /*******************************************************************
  ** m a t h 6 4 . c
  ** Forth Inspired Command Language - 64 bit math support routines
  ** Author: 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!
  *******************************************************************/
  
  /* $FreeBSD: src/sys/boot/ficl/math64.c,v 1.2 1999/09/29 04:43:06 dcs Exp $ */
  
  #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	gaulandm@mdhost.cse.tek.com  
  **************************************************************************/
  #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	gaulandm@mdhost.cse.tek.com  
  **************************************************************************/
  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	gaulandm@mdhost.cse.tek.com  
  **************************************************************************/
  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
  
  