cpp-bignumber/src/bignum_helper.cpp

//
//  bignum_helper.cpp
//  bigdecimal
//
//  Created by Sam Jaffe on 7/4/17.
//

#include "bignum_helper.h"

#include <cmath>

using namespace math::detail;

namespace {
  constexpr int32_t const MAX_SEG { 999999999};
  constexpr int32_t const OVER_SEG{1000000000};
  constexpr int32_t const SEG_DIGITS{9};
}

namespace math { namespace detail {
#define IMPL_COMPARE(expr) \
  if (rhs expr < lhs expr) return -1; \
  else if (lhs expr < rhs expr) return 1
  int compare(data_type const & rhs, data_type const & lhs) {
    IMPL_COMPARE(.size());
    for (size_t i = rhs.size(); i > 0; --i) {
      IMPL_COMPARE([i-1]);
    }
    return 0;
  }
#undef IMPL_COMPARE
  
  void add(data_type & rhs, data_type const & lhs, size_t offset) {
    rhs.resize(std::max(rhs.size(), lhs.size()+offset)+1);
    auto const lbnd = lhs.size(), ubnd = rhs.size() - 1;
    // Add
    for (size_t i = 0; i < lbnd; ++i) {
      rhs[i+offset] += lhs[i];
    }
    // Carry
    for (size_t i = 0; i < ubnd; ++i) {
      if (rhs[i] > MAX_SEG) {
        rhs[i] -= OVER_SEG;
        rhs[i+1] += 1;
      }
    }
    if (rhs[ubnd] == 0) { rhs.pop_back(); }
  }
  
  void subtract_nounderflow(data_type & rhs, data_type const & lhs, size_t offset) {
    size_t const rbnd = rhs.size(), lbnd = lhs.size();
    // Subtract
    for (size_t i = 0; i < lbnd; ++i) {
      rhs[i+offset] -= lhs[i];
    }
    // Borrow
    for (size_t i = 0; i < rbnd; ++i) {
      if (rhs[i] < 0) {
        rhs[i] += OVER_SEG;
        rhs[i+1] -= 1;
      }
    }
    if (rhs[rbnd-1] == 0 && rbnd > 1) { rhs.pop_back(); }
  }
  
  void multiply(data_type & rhs, data_type const & lhs) {
    size_t const rbnd = rhs.size(), lbnd = lhs.size();
    size_t const ubnd = rbnd + lbnd;
    rhs.resize(ubnd);
    // Multiply
    for (size_t i = rbnd; i > 0; --i) {
      int32_t const value = rhs[i-1];
      for (size_t j = lbnd; j > 0; --j) {
        // Max input              999,999,999
        // Max output 999,999,998,000,000,001
        int64_t product = static_cast<int64_t>(value) * static_cast<int64_t>(lhs[j-1]);
        int64_t overflow = product / OVER_SEG;
        rhs[i+j-2] += static_cast<int32_t>(product - (overflow * OVER_SEG));
        rhs[i+j-1] += static_cast<int32_t>(overflow);
      }
      rhs[i-1] -= value;
    }
    // Carry
    for (size_t i = 0; i < ubnd-1; ++i) {
      if (rhs[i] > MAX_SEG) {
        int32_t overflow = rhs[i] / OVER_SEG;
        rhs[i] -= (overflow * OVER_SEG);
        rhs[i+1] += overflow;
      }
    }
    while (rhs.back() == 0) { rhs.pop_back(); }
  }
  
  data_type shift10(data_type const & data, int32_t pow, size_t shift) {
    size_t const bnd = data.size();
    data_type rval(bnd + shift + 1);
    for (size_t i = 0; i < bnd; ++i) {
      int64_t product = static_cast<int64_t>(data[i]) * static_cast<int64_t>(pow);
      int64_t overflow = product / OVER_SEG;
      rval[i+shift] += static_cast<int32_t>(product - (overflow * OVER_SEG));
      rval[i+shift+1] += static_cast<int32_t>(overflow);
    }
    if (rval.back() == 0 && rval.size() > 1) { rval.pop_back(); }
    return rval;
  }
  
  size_t digits(int32_t val) {
    return val == 0 ? 1 : static_cast<size_t>(floor(log10(val))) + 1;
  }
  
  size_t digits(data_type const & data) {
    return SEG_DIGITS * (data.size()-1) + digits(data.back());
  }
  
  int32_t powers[] = { 1, 10, 100, 1000, 10000, 100000, 1000000, 10000000, 100000000 };
  
  data_type divide(data_type & remainder, data_type const & divisor) {
    data_type accum{0};
    auto const dig = digits(divisor);
    do {
      auto const diff = std::max(1UL, digits(remainder) - dig) - 1;
      auto const shift = diff / SEG_DIGITS;
      auto const ipow = diff - (shift * SEG_DIGITS);
      data_type step{shift10({1}, powers[ipow], 0)};
      data_type value{shift10(divisor, powers[ipow], 0)};
      do {
        subtract_nounderflow(remainder, value, shift);
        add(accum, step, shift);
      } while (compare(remainder, value) >= 0);
    } while (compare(remainder, divisor) >= 0);
    return accum;
  }
} }