lpi (Lode's Programming Interface) Code

 
/*
Copyright (c) 2005-2007 Lode Vandevenne
All rights reserved.

Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met:

  * Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer.
  * 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.
  * Neither the name of Lode Vandevenne nor the names of his contributors may be used to endorse or promote products derived from this software without specific prior written permission.

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS 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 COPYRIGHT OWNER 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 "lpi_largenums.h"

namespace lpi
{

fixed128& fixed128::operator+=(const fixed128& rhs)
{
  uint32 carry = 0U;
  for(int i = 0; i < 4; i++)
  {
    data[i] += rhs.data[i] + carry;
    if(!carry && data[i] < rhs.data[i]) carry = 1U; //if it has overflown
    else if(carry && data[i] <= rhs.data[i]) carry = 1U;
    else carry = 0U;
  }
  return *this;
}

fixed128& fixed128::operator-=(const fixed128& rhs)
{
  uint32 borrow = 0U;
  for(int i = 0; i < 4; i++)
  {
    if(borrow)
    {
      borrow = rhs.data[i] >= data[i] ? 1U : 0U;
      data[i] -= rhs.data[i] + 1U;
    }
    else
    {
      borrow = rhs.data[i] > data[i] ? 1U : 0U;
      data[i] -= rhs.data[i];
    }
  }
  return *this;
}

fixed128& fixed128::operator+=(uint32 rhs)
{
  data[1] += rhs; //data[0], the fractional part, is left untouched
  if(data[1] < rhs) //if it has overflown
  for(int i = 2; i < 4; i++)
  {
    data[i]++;
    if(data[i] != 0U) break; //only if data[i] is 0 (4294967295 + 1), it has overflown and we should continue
  }
  return *this;
}

fixed128& fixed128::operator-=(uint32 rhs)
{
  if(rhs <= data[1])
  {
    data[1] -= rhs; //data[0], the fractional part, is left untouched
  }
  else //overflow case
  {
    data[1] -= rhs; //data[0], the fractional part, is left untouched
    for(int i = 2; i < 4; i++)
    {
      data[i]--;
      if(data[i] != 4294967295U) break; //only if data[i] is 4294967295 (0 - 1), it has overflown and we should continue
    }
  }
  return *this;
}

void fixed128::addFractional(uint32 rhs)
{
  data[0] += rhs; //data[0], the fractional part, is left untouched
  if(data[0] < rhs) //if it has overflown
  for(int i = 1; i < 4; i++)
  {
    data[i]++;
    if(data[i] != 0U) break; //only if data[i] is 0 (4294967295 + 1), it has overflown and we should continue
  }
}

void fixed128::subtractFractional(uint32 rhs)
{
  if(rhs <= data[0])
  {
    data[0] -= rhs; //data[0], the fractional part, is left untouched
  }
  else //overflow case
  {
    data[0] -= rhs; //data[0], the fractional part, is left untouched
    for(int i = 1; i < 4; i++)
    {
      data[i]--;
      if(data[i] != 4294967295U) break; //only if data[i] is 4294967295 (0 - 1), it has overflown and we should continue
    }
  }
}

fixed128::operator double() const //this only returns a correct result if the double can hold this value
{
  double result = 0.0;
  
  if(isNegative())
  {
    double multiplier = 1.0 / 4294967296.0;
    for(int i = 0; i < 4; i++)
    {
      result += multiplier * (~data[i]);
      multiplier *= 4294967296.0;
    }
    return -result - 1.0 / 4294967296.0;
  }
  else
  {
    double multiplier = 1.0 / 4294967296.0;
    for(int i = 0; i < 4; i++)
    {
      result += multiplier * data[i];
      multiplier *= 4294967296.0;
    }
    return result;
  }
}

void fixed128::operator=(double d)
{
  bool negative = d < 0.0;
  if(negative) d = -d;
  double multiplier = 4294967296.0;
  for(int i = 0; i < 4; i++)
  {
    data[i] = uint32(d * multiplier);
    multiplier /= 4294967296.0;
  }
  if(negative) negate();
}

fixed128::operator uint32() const //this only returns a correct result if the double can hold this value
{
  return data[1];
}

void fixed128::operator=(uint32 d)
{
  data[0] = 0;
  data[1] = d;
  data[2] = 0;
  data[3] = 0;
}

fixed128& fixed128::operator++()
{
  data[1]++;
  if(data[1] == 0)
  {
    data[2]++;
    if(data[2] == 0)
    {
      data[3]++;
    }
  }
  return *this;
}

fixed128 fixed128::operator++(int)
{
  fixed128 result = *this;
  (*this)++;
  return result;
}

fixed128& fixed128::operator--()
{
  data[1]--;
  if(data[1] == 4294967295U)
  {
    data[2]--;
    if(data[2] == 4294967295U)
    {
      data[3]--;
    }
  }
  return *this;
}

fixed128 fixed128::operator--(int)
{
  fixed128 result = *this;
  (*this)--;
  return result;
}

void fixed128::addLSB()
{
  data[0]++;
  if(data[0] == 0)
  {
    data[1]++;
    if(data[1] == 0)
    {
      data[2]++;
      if(data[2] == 0)
      {
        data[3]++;
      }
    }
  }
}

void fixed128::subtractLSB()
{
  data[0]--;
  if(data[0] == 4294967295U)
  {
    data[1]--;
    if(data[1] == 4294967295U)
    {
      data[2]--;
      if(data[2] == 4294967295U)
      {
        data[3]--;
      }
    }
  }
}

void fixed128::bit_invert()
{
  for(int i = 0; i < 4; i++) data[i] = ~data[i];
}

void fixed128::negate()
{
  bit_invert();
  addLSB();
}

bool fixed128::isNegative() const
{
  return data[3] > 2147483647U;
}

fixed128& fixed128::operator>>=(uint32 shift)
{
  //will make this more efficient if needed
  if(shift >= 128)
  {
    if(isNegative()) data[0] = data[1] = data[2] = data[3] = 4294967295U;
    else data[0] = data[1] = data[2] = data[3] = 0U;
  }
  else
  {
    while(shift >= 32)
    {
      shift -= 32;
      data[0] = data[1];
      data[1] = data[2];
      data[2] = data[3];
      data[3] = isNegative() ? 4294967295U : 0U;
    }
    while(shift >= 8)
    {
      void* v = &data[0];
      unsigned char* c = (unsigned char*)v;
      shift -= 8;
      for(int i = 0; i < 15; i++) c[i] = c[i + 1];
      c[15] = isNegative() ? 255U : 0U;
    }
    while(shift != 0)
    {
      bool negative = isNegative();
      shift--;
      data[0] >>= 1; if(data[1] & 1) data[0] |= 2147483648U; else data[0] &= 2147483647U;
      data[1] >>= 1; if(data[2] & 1) data[1] |= 2147483648U; else data[1] &= 2147483647U;
      data[2] >>= 1; if(data[3] & 1) data[2] |= 2147483648U; else data[2] &= 2147483647U;
      data[3] >>= 1; if(negative) data[3] |= 2147483648U; else data[3] &= 2147483647U;
    }
  }
  
  return *this;
}

fixed128& fixed128::operator<<=(uint32 shift)
{
  //will make this more efficient if needed
  if(shift >= 128)
  {
    data[0] = data[1] = data[2] = data[3] = 0U;
  }
  else
  {
    while(shift >= 32)
    {
      shift -= 32;
      data[3] = data[2];
      data[2] = data[1];
      data[1] = data[0];
      data[0] = 0;
    }
    while(shift >= 8)
    {
      void* v = &data[0];
      unsigned char* c = (unsigned char*)v;
      shift -= 8;
      for(int i = 0; i < 15; i++) c[i + 1] = c[i];
      c[0] = 0;
    }
    while(shift != 0)
    {
      shift--;
      data[3] <<= 1; if(data[2] & 2147483648U) data[3] |= 1; else data[3] &= 4294967294U;
      data[2] <<= 1; if(data[1] & 2147483648U) data[2] |= 1; else data[2] &= 4294967294U;
      data[1] <<= 1; if(data[0] & 2147483648U) data[1] |= 1; else data[1] &= 4294967294U;
      data[0] <<= 1;
    }
  }
  
  return *this;
}

void fixed128::makeZero()
{
  data[0] = data[1] = data[2] = data[3] = 0U;
}

fixed128& fixed128::operator*=(uint32 rhs)
{
  fixed128 a = *this;
  
  makeZero();
  
  for(int i = 0; i < 32; i++)
  {
    if((rhs >> i) & 1) (*this) += a;
    a <<= 1;
  }

  return *this;
}

fixed128& fixed128::operator*=(const fixed128& rhs)
{
  //160 bit precision for multiplication: one 32-bit group extra, because at the end we need to do a right shift
  uint32 a[5] = { data[0], data[1], data[2], data[3], isNegative() ? 4294967295U : 0U };
  uint32 extra = 0U; //similar: 5th 32-bit group for *this, representing "data[4]"
  
  fixed128 b = rhs;
  
  makeZero(); //result will be stored in *this, start at zero
  
  for(int j = 0; j < 128; j++)
  {
    if(b.getLSB()) //the code from the += operator implemented again, to fill the 5th uint32 (*this += a)
    {
      uint32 carry = 0U;
      for(int i = 0; i < 4; i++)
      {
        data[i] += a[i] + carry;
        if(!carry && data[i] < a[i]) carry = 1U; //if it has overflown
        else if(carry && data[i] <= a[i]) carry = 1U;
        else carry = 0U;
      }
      extra += a[4] + carry;
    }
    //the code from left shift with 1 bit implemented again, to fill the 5th uint32 (a <<= 1)
    a[4] <<= 1; if(a[3] & 2147483648U) a[4] |= 1; else a[4] &= 4294967294U;
    a[3] <<= 1; if(a[2] & 2147483648U) a[3] |= 1; else a[3] &= 4294967294U;
    a[2] <<= 1; if(a[1] & 2147483648U) a[2] |= 1; else a[2] &= 4294967294U;
    a[1] <<= 1; if(a[0] & 2147483648U) a[1] |= 1; else a[1] &= 4294967294U;
    a[0] <<= 1;
      
    b >>= 1;
  }
  
  data[0] = data[1];
  data[1] = data[2];
  data[2] = data[3];
  data[3] = extra;

  return *this;
}

fixed128& fixed128::operator*=(double rhs)
{
  return operator*=((fixed128)rhs);
}

//this one loses the leftmost 32 bits, but shows a simpler version of the multiplication implementation for documentation
void fixed128::crudeMultiply(const fixed128& rhs)
{
  fixed128 a = *this;
  fixed128 b = rhs;
  makeZero();
  for(int i = 0; i < 128; i++)
  {
    if(b.getLSB()) (*this) += a;
    a <<= 1;
    b >>= 1;
  }
  
  (*this) >>= 32; //the part behind the comma is smaller
}

bool fixed128::operator==(const fixed128& rhs) const
{
  return data[0] == rhs.data[0]
      && data[1] == rhs.data[1]
      && data[2] == rhs.data[2]
      && data[3] == rhs.data[3];
}

bool fixed128::operator==(fixed128::uint32 rhs) const
{
  return data[0] == 0U
      && data[1] == rhs
      && data[2] == 0U
      && data[3] == 0U;
}

bool fixed128::operator>(const fixed128& rhs) const
{
  if(isNegative() && !rhs.isNegative()) return false;
  if(!isNegative() && rhs.isNegative()) return true;
  
  if(data[3] == rhs.data[3])
  {
    if(data[2] == rhs.data[2])
    {
      if(data[1] == rhs.data[1])
      {
        return data[0] > rhs.data[0];
      }
      else return data[1] > rhs.data[1];
    }
    else return data[2] > rhs.data[2];
  }
  else return data[3] > rhs.data[3];
}

bool fixed128::operator<(const fixed128& rhs) const
{
  if(isNegative() && !rhs.isNegative()) return true;
  if(!isNegative() && rhs.isNegative()) return false;
  
  if(data[3] == rhs.data[3])
  {
    if(data[2] == rhs.data[2])
    {
      if(data[1] == rhs.data[1])
      {
        return data[0] < rhs.data[0];
      }
      else return data[1] < rhs.data[1];
    }
    else return data[2] < rhs.data[2];
  }
  else return data[3] < rhs.data[3];
}

bool fixed128::operator>=(const fixed128& rhs) const
{
  return !((*this) > rhs);
}

bool fixed128::operator<=(const fixed128& rhs) const
{
  return !((*this) < rhs);
}

fixed128& fixed128::operator&=(const fixed128& rhs)
{
  for(int i = 0; i < 4; i++)
  {
    data[i] &= rhs.data[i];
  }
  
  return *this;
}

fixed128& fixed128::operator|=(const fixed128& rhs)
{
  for(int i = 0; i < 4; i++)
  {
    data[i] |= rhs.data[i];
  }
  
  return *this;
}

fixed128& fixed128::operator^=(const fixed128& rhs)
{
  for(int i = 0; i < 4; i++)
  {
    data[i] ^= rhs.data[i];
  }
  
  return *this;
}

////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////

fixed128 operator-(const fixed128& a)
{
  fixed128 result = a;
  result.negate();
  return result;
}

fixed128 operator~(const fixed128& a)
{
  fixed128 result = a;
  result.bit_invert();
  return result;
}


fixed128 operator+(const fixed128& a, const fixed128& b)
{
  fixed128 result = a;
  result += b;
  return result;
}

fixed128 operator-(const fixed128& a, const fixed128& b)
{
  fixed128 result = a;
  result -= b;
  return result;
}

fixed128 operator+(const fixed128& a, fixed128::uint32& b)
{
  fixed128 result = a;
  result += b;
  return result;
}

fixed128 operator+(fixed128::uint32& a, const fixed128& b)
{
  fixed128 result = b;
  result += a;
  return result;
}

fixed128 operator-(const fixed128& a, fixed128::uint32& b)
{
  fixed128 result = a;
  result -= b;
  return result;
}

fixed128 operator>>(const fixed128& a, fixed128::uint32 b)
{
  fixed128 result = a;
  result >>= b;
  return result;
}

fixed128 operator<<(const fixed128& a, fixed128::uint32 b)
{
  fixed128 result = a;
  result <<= b;
  return result;
}

fixed128 operator&(const fixed128& a, const fixed128& b)
{
  fixed128 result = a;
  result &= b;
  return result;
}

fixed128 operator|(const fixed128& a, const fixed128& b)
{
  fixed128 result = a;
  result |= b;
  return result;
}

fixed128 operator^(const fixed128& a, const fixed128& b)
{
  fixed128 result = a;
  result ^= b;
  return result;
}

fixed128 operator*(const fixed128& a, const fixed128& b)
{
  fixed128 result = a;
  result *= b;
  return result;
}

fixed128 operator*(const fixed128& a, fixed128::uint32 b)
{
  fixed128 result = a;
  result *= b;
  return result;
}

fixed128 operator*(fixed128::uint32 a, const fixed128& b)
{
  fixed128 result = b;
  result *= a;
  return result;
}


fixed128 operator*(const fixed128& a, double b)
{
  fixed128 result = a;
  result *= b;
  return result;
}

fixed128 operator*(double a, const fixed128& b)
{
  fixed128 result = b;
  result *= a;
  return result;
}

bool operator==(fixed128::uint32 a, const fixed128& b)
{
  return b == a;
}


} //namespace lpi

//test
#if 0
#include <iostream>

class Testfixed128
{
  public:
  
  void testaddsub()
  {
    using namespace lpi;
    std::cout.precision(20);
    
    std::cout << "testing addition and subtraction" << std::endl;
    
    fixed128 n = -5.0;
    fixed128 p = +5.0;
    
    std::cout<<(double)(p + p)<<std::endl;
    std::cout<<(double)(p + n)<<std::endl;
    std::cout<<(double)(n + p)<<std::endl;
    std::cout<<(double)(n + n)<<std::endl;
    std::cout<<(double)(p - p)<<std::endl;
    std::cout<<(double)(p - n)<<std::endl;
    std::cout<<(double)(n - p)<<std::endl;
    std::cout<<(double)(n - n)<<std::endl;
    
    std::cout << std::endl;
    
    fixed128 f = 5.5;
    fixed128 g = -5.5;
    fixed128::uint32 i = 5;
    fixed128::uint32 j = 6;

    std::cout<<(double)(f + i)<<std::endl;
    std::cout<<(double)(f - i)<<std::endl;
    std::cout<<(double)(g + i)<<std::endl;
    std::cout<<(double)(g - i)<<std::endl;
    std::cout<<(double)(f + j)<<std::endl;
    std::cout<<(double)(f - j)<<std::endl;
    std::cout<<(double)(g + j)<<std::endl;
    std::cout<<(double)(g - j)<<std::endl;
    
    std::cout<<std::endl;
  }
  
  void testsmallmul()
  {
    using namespace lpi;
    std::cout.precision(20);
    
    std::cout << "testing small multiplication" << std::endl;
    
    fixed128 a = 0.3;
    fixed128 b = 0.01;
    
    std::cout<<(double)a<<" "<<(double)b<<std::endl;
    
    a *= b;
    
    std::cout<<(double)a<<std::endl;
    
    std::cout<<std::endl;
  }
  
  void testlargemul()
  {
    using namespace lpi;
    std::cout.precision(30);
    
    std::cout << "testing large multiplication" << std::endl;
    
    fixed128 a = 8000000000.0;
    fixed128 b = 9000000000.0;
    fixed128 c = -8000000000.0;
    fixed128 d = -9000000000.0;
    
    std::cout<<(double)(a*b)<<std::endl;
    std::cout<<(double)(a*d)<<std::endl;
    std::cout<<(double)(c*b)<<std::endl;
    std::cout<<(double)(c*d)<<std::endl;
    
    std::cout<<std::endl;
  }
  
  void testintmul()
  {
    using namespace lpi;
    std::cout.precision(30);
    
    std::cout << "testing multiplication with integer" << std::endl;
    
    fixed128 a = 8000000000.0;
    fixed128::uint32 b = 5;
    
    std::cout<<(double)a<<" "<<b<<std::endl;
    
    a *= b;
    
    std::cout<<(double)a<<std::endl;
    
    std::cout<<std::endl;
  }
  
  void testdoublemul()
  {
    using namespace lpi;
    std::cout.precision(30);
    
    std::cout << "testing multiplication with doubles" << std::endl;
    
    fixed128 a = 8000000000.0;
    
    std::cout<<(double)(a * 0.5)<<std::endl;
    std::cout<<(double)(a * -1.5)<<std::endl;
    std::cout<<(double)(a * 1000.0)<<std::endl;
    std::cout<<(double)(a * -1.0)<<std::endl;
    std::cout<<(double)(a * 0.0)<<std::endl;
    
    std::cout << std::endl;
    fixed128 b = -8000000000.0;
    
    std::cout<<(double)(b * 0.5)<<std::endl;
    std::cout<<(double)(b * -1.5)<<std::endl;
    std::cout<<(double)(b * 1000.0)<<std::endl;
    std::cout<<(double)(b * -1.0)<<std::endl;
    std::cout<<(double)(b * 0.0)<<std::endl;
  }
  
  void testltgt()
  {
    using namespace lpi;
    std::cout.precision(30);
    
    std::cout << "testing lesser than and greater than" << std::endl;
    
    fixed128 a = +5.0;
    fixed128 b = +6.0;
    fixed128 c = -5.0;
    fixed128 d = -6.0;
    
    std::cout<<(a < b)<<std::endl;
    std::cout<<(a < c)<<std::endl;
    std::cout<<(c < d)<<std::endl;
    std::cout<<(b < a)<<std::endl;
    std::cout<<(c < a)<<std::endl;
    std::cout<<(d < c)<<std::endl;
    std::cout << std::endl;
    std::cout<<(a > b)<<std::endl;
    std::cout<<(a > c)<<std::endl;
    std::cout<<(c > d)<<std::endl;
    std::cout<<(b > a)<<std::endl;
    std::cout<<(c > a)<<std::endl;
    std::cout<<(d > c)<<std::endl;
  }
  
  Testfixed128()
  {
    testaddsub();
    testsmallmul();
    testlargemul();
    testintmul();
    testdoublemul();
    testltgt();
  }
} testfixed128;

int main(){}
#endif