h_.hpp

/*
Copyright 2024 Owain Davies
SPDX-License-Identifier: Apache-2.0
*/
 
#ifndef BI_SRC_H__HPP_
#define BI_SRC_H__HPP_
 
#include <algorithm>
#include <array>
#include <cassert>
#include <cmath>
#include <limits>
#include <random>
#include <string>
#include <utility>
 
#include "bi.hpp"
#include "bi.inl"
#include "bi_exceptions.hpp"
#include "constants.hpp"
#include "uints.hpp"
 
 
namespace bi {
 
struct h_ {
  // increment/decrement
  static void increment_abs(bi_t& x);
  static void decrement_abs(bi_t& x);
 
  // additive
  static void add_abs(bi_t& result, const bi_t& a, const bi_t& b);
  static void add(bi_t& result, const bi_t& a, const bi_t& b);
  static void sub_abs_gt(bi_t& result, const bi_t& a, const bi_t& b);
  static void sub_abs(bi_t& result, const bi_t& a, const bi_t& b);
  static void sub(bi_t& result, const bi_t& a, const bi_t& b);
 
  // multiplicative
  static void imul1add1(bi_t& x, digit v, digit k);
  static void mul_algo_square(bi_t& result, const bi_t& a, const size_t m);
  static void mul_algo_knuth(bi_t& result, const bi_t& a, const bi_t& b,
                             const size_t m, const size_t n);
  static void mul_karatsuba(bi_t& result, const bi_t& a, const bi_t& b);
  static void mul_standard(bi_t& result, const bi_t& a, const bi_t& b);
  static void mul(bi_t& result, const bi_t& a, const bi_t& b);
  static digit div_algo_digit(bi_t& q, const bi_t& u, digit v) noexcept;
  static void div_algo_single(bi_t& q, bi_t& r, const bi_t& n,
                              const bi_t& d) noexcept;
  static void div_algo_binary(bi_t& q, bi_t& r, const bi_t& n, const bi_t& d);
  static void div_algo_knuth(bi_t& q, bi_t& r, const bi_t& n, const bi_t& d);
  static void divide(bi_t& q, bi_t& r, const bi_t& n, const bi_t& d);
 
  // bits
  static void left_shift(bi_t& result, const bi_t& a, bi_bitcount_t shift);
  static void right_shift(bi_t& result, const bi_t& a, bi_bitcount_t shift);
  enum class BitwiseOperation { AND, OR, XOR };
  template <BitwiseOperation Op>
  static void bitwise_operation_impl(bi_t& res, const bi_t& a, const bi_t& b);
  template <BitwiseOperation Op>
  static bi_t bitwise_operation(const bi_t& a, const bi_t& b);
  template <BitwiseOperation Op>
  static bi_t& ibitwise_operation(bi_t&, const bi_t& other);
 
  // comparisons
  static int cmp_abs(const bi_t&, const bi_t&) noexcept;
  static int cmp(const bi_t&, const bi_t&) noexcept;
  template <std::integral T>
  static int cmp(const bi_t& a, T b) noexcept;
 
  // to_string()
  static uint8_t idiv10(bi_t& x) noexcept;
  static size_t base_length(const bi_t& x, int base);
 
  // initializing
  template <std::integral T>
  static void init_one_digit(bi_t& x, T value);
  template <std::integral T>
  static void init_atleast_one_digit(bi_t& x, T value);
  // NOLINTNEXTLINE
  static void init_string(bi_t& x, const std::string& str, int base = 10);
 
  // misc.
  static dvector to_twos_complement(const dvector& vec);
  static void to_twos_complement_in_place(dvector& vec) noexcept;
  static void bisect(const bi_t&, bi_t&, bi_t&, size_t m);
 
  // double
  static void assign_from_double(bi_t&, double);
  static int cmp_abs(const bi_t&, double) noexcept;
  static int cmp(const bi_t&, double) noexcept;
 
  // exponentiation
  template <std::unsigned_integral T>
  static bi_t expo_left_to_right(const bi_t& base, T exp);
 
  // NOLINTNEXTLINE(cppcoreguidelines-avoid-non-const-global-variables)
  static thread_local std::mt19937 rng_;
  static bi_t random_(bi_bitcount_t z);
};
 
thread_local std::mt19937 h_::rng_{std::random_device{}()};  // NOLINT
 
void h_::increment_abs(bi_t& x) {
  if (x.size() == 0 || x[x.size() - 1] == std::numeric_limits<digit>::max()) {
    x.reserve_(x.size() + 1);
  }
 
  if (x.size() == 0) {
    x.resize_(1);
    x[0] = 1;
    return;
  }
 
  size_t i = 0;
  bool carry = true;
 
  while (carry && i < x.size()) {
    if (x[i] == std::numeric_limits<digit>::max()) {
      x[i] = 0;
    } else {
      ++x[i];
      carry = false;
    }
    ++i;
  }
 
  if (carry) {
    x.vec_.push_back(1);
  }
}
 
void h_::decrement_abs(bi_t& x) {
  if (x.size() == 0) {
    x.resize_(1);
    x[0] = 1;
    x.negative_ = true;
    return;
  }
 
  size_t i = 0;
  bool borrow = true;
 
  while (borrow && i < x.size()) {
    if (x[i] == 0) {
      // Borrow from the next more significant digit
      x[i] = std::numeric_limits<digit>::max();
    } else {
      --x[i];
      borrow = false;
    }
    ++i;
  }
 
  x.trim();
}
 
 
void h_::add_abs(bi_t& r, const bi_t& x, const bi_t& y) {
  // Algorithm assumes x.size() >= y.size()
  const bi_t& large = x.size() >= y.size() ? x : y;
  const bi_t& small = x.size() >= y.size() ? y : x;
 
  r.reserve_(large.size() + 1);
 
  digit temp{};
  bool carry = 0;
  size_t i = 0;
  for (; i < small.size(); ++i) {
    uints::uaddc(temp, large[i], small[i], carry);
    r[i] = temp;
  }
 
  while (i < large.size()) {
    temp = large[i] + carry;
    carry = temp < static_cast<uint8_t>(carry);  // MSVC warning C4804
    r[i++] = temp;
  }
 
  r[i] = carry;
  r.resize_(large.size() + 1);
  r.trim();
  r.negative_ = false;
}
 
void h_::add(bi_t& r, const bi_t& x, const bi_t& y) {
  if (x.negative()) {
    if (y.negative()) {
      // x < 0, y < 0 ==> x = -|x|, y = -|y|. ==> x + y = -(|x| + |y|)
      add_abs(r, x, y);
      r.negative_ = true;
    } else {
      // x < 0, y >= 0 ==> x = -|x|, y = |y| ==> x + y = |y| - |x|
      sub_abs(r, y, x);
    }
  } else {
    if (y.negative()) {
      // x >= 0, y < 0 ==> x = |x|, y = -|y| ==> x + y = |x| - |y|
      sub_abs(r, x, y);
    } else {
      // x >= 0, y >= 0 ==> x = |x|, y = |y| ==> x + y = |x| + |y|
      add_abs(r, x, y);
    }
  }
}
 
void h_::sub_abs_gt(bi_t& r, const bi_t& x, const bi_t& y) {
  assert(x.size() >= y.size());
 
  r.reserve_(x.size());
 
  digit temp{};
  bool borrow = 0;
  size_t i = 0;
  for (; i < y.size(); ++i) {
    uints::usubb(temp, x[i], y[i], borrow);
    r[i] = temp;
  }
  while (i < x.size()) {
    temp = x[i] - borrow;
    borrow = temp > x[i];
    r[i++] = temp;
  }
 
  r.resize_(x.size());
  r.trim();
  r.negative_ = false;
}
 
void h_::sub_abs(bi_t& r, const bi_t& x, const bi_t& y) {
  if (x.size() == y.size()) {
    const auto [it_x, it_y] = std::mismatch(x.rbegin(), x.rend(), y.rbegin());
 
    if (it_x == x.rend()) {
      // x, y have the same digits (or both are zero)
      r.resize_(0);
      r.negative_ = false;
      return;
    }
 
    if (*it_y > *it_x) {
      sub_abs_gt(r, y, x);
      r.negative_ = true;
    } else {
      sub_abs_gt(r, x, y);
    }
  } else if (x.size() > y.size()) {
    sub_abs_gt(r, x, y);
  } else {
    sub_abs_gt(r, y, x);
    r.negative_ = true;
  }
}
 
void h_::sub(bi_t& r, const bi_t& x, const bi_t& y) {
  if (x.negative()) {
    if (y.negative()) {
      // x < 0, y < 0 ==> x = -|x|, y = -|y| ==> x - y = |y| - |x|
      sub_abs(r, y, x);
    } else {
      // x < 0, y >= 0 ==> x = -|x|, y = |y| ==> x - y = -(|x| + |y|)
      add_abs(r, x, y);
      r.negative_ = true;
    }
  } else {
    if (y.negative()) {
      // x >= 0, y < 0 ==> x = |x|, y = -|y| ==> x - y = |x| + |y|
      add_abs(r, x, y);
    } else {
      // x >= 0, y >= 0 ==> x = |x|, y = |y| ==> x - y = |x| - |y|
      sub_abs(r, x, y);
    }
  }
}
 
 
inline void h_::imul1add1(bi_t& x, digit v, digit k = 0) {
  for (auto& d : x.vec_) {
    const ddigit t = static_cast<ddigit>(d) * v + k;
    k = t >> bi_dwidth;         // floor(t / 2^{bi_dwidth})
    d = static_cast<digit>(t);  // t mod 2^{bi_dwidth}
  }
 
  if (k) {
    x.vec_.push_back(k);
  }
}
 
void h_::mul_algo_square(bi_t& w, const bi_t& a, const size_t t) {
  std::fill(w.begin(), w.end(), 0);
 
  for (size_t i = 0; i < t; ++i) {
    ddigit uv = w[2 * i] + static_cast<ddigit>(a[i]) * a[i];
    w[2 * i] = static_cast<digit>(uv);  // set w[2 * i] <- v
    ddigit c = uv >> bi_dwidth;         // set c <- u
 
    for (size_t j = i + 1; j < t; ++j) {
      // (hi, lo) represents a quadruple digit
      ddigit hi{0};
      ddigit lo = static_cast<ddigit>(a[j]) * a[i];
 
      // Multiply product by two
      hi = lo >> (std::numeric_limits<ddigit>::digits - 1);
      lo <<= 1;
 
      // Now add w[i + j] and c to (hi, lo)
      hi += uints::uadd_overflow(lo, lo, static_cast<ddigit>(w[i + j]));
      hi += uints::uadd_overflow(lo, lo, c);
 
      w[i + j] = static_cast<digit>(lo);          // set w[i + j] <- v
      c = (hi << bi_dwidth) | (lo >> bi_dwidth);  // set c <- u
    }
 
    for (size_t k = i + t; c > 0; ++k) {
      c += w[k];
      w[k] = static_cast<digit>(c);
      c >>= bi_dwidth;
    }
  }
}
 
void h_::mul_algo_knuth(bi_t& w, const bi_t& a, const bi_t& b, const size_t m,
                        const size_t n) {
  // NOLINTNEXTLINE(cppcoreguidelines-pro-bounds-pointer-arithmetic)
  std::fill(w.begin(), w.begin() + m, 0);
 
  for (size_t j = 0; j < n; ++j) {
    digit k = 0;
    for (size_t i = 0; i < m; ++i) {
      const ddigit t =
          static_cast<ddigit>(a[i]) * b[j] + static_cast<ddigit>(w[i + j]) + k;
      k = t >> bi_dwidth;                // floor(t / 2^{bi_dwidth})
      w[i + j] = static_cast<digit>(t);  // t mod 2^{bi_dwidth}
    }
    w[j + m] = k;
  }
}
 
void h_::bisect(const bi::bi_t& x, bi::bi_t& lower, bi::bi_t& upper, size_t n) {
  const size_t bisect_point = std::min(n, x.size());
 
  lower.vec_ = dvector(x.vec_.begin(), x.vec_.begin() + bisect_point);
  upper.vec_ = dvector(x.vec_.begin() + bisect_point, x.vec_.end());
 
  lower.trim();
  upper.trim();
}
 
void h_::mul_karatsuba(bi_t& w, const bi_t& u, const bi_t& v) {
  const size_t n = std::min(u.size(), v.size()) >> 1;
 
  bi_t u0, u1, v0, v1;
  bisect(u, u0, u1, n);
  bisect(v, v0, v1, n);
 
  bi_t a, b, c;
  mul(a, u1, v1);  // a = u1 * v1
  mul(b, u0, v0);  // b = u0 * v0
 
  u1 += u0;
  v1 += v0;
  mul(c, u1, v1);
  c -= a + b;  // c = (u1 + u0)(v1 + v0) - (u0v0 + u1v1)
 
  a <<= (2 * n * bi_dbits);
  c <<= (n * bi_dbits);
  w = a + c + b;  // w = (bi_base ** 2n) * a + (bi_base ** n) * c + b
}
 
void h_::mul_standard(bi_t& result, const bi_t& a, const bi_t& b) {
  if (a.size() == 0 || b.size() == 0) {
    result.resize_(0);
    result.negative_ = false;
    return;
  }
 
  const size_t m = a.size();
  const size_t n = b.size();
 
  const bool overlap = &result == &a || &result == &b;
 
  bi_t temp;
  bi_t& target = overlap ? temp : result;
 
  target.resize_(m + n);
 
  if (&a != &b) {
    h_::mul_algo_knuth(target, a, b, m, n);
  } else {
    h_::mul_algo_square(target, a, m);
  }
 
  target.trim();
 
  if (overlap) {
    result.swap(target);
  }
}
 
void h_::mul(bi_t& w, const bi_t& u, const bi_t& v) {
  if (u.size() < karatsuba_threshold || v.size() < karatsuba_threshold) {
    h_::mul_standard(w, u, v);
  } else {
    h_::mul_karatsuba(w, u, v);
  }
 
  w.negative_ = u.negative() != v.negative();
}
 
digit h_::div_algo_digit(bi_t& q, const bi_t& u, digit v) noexcept {
  ddigit rem{0};
 
  for (size_t j = u.size() - 1; j < std::numeric_limits<size_t>::max(); --j) {
    const ddigit temp = (rem << bi_dwidth) | u[j];  // rem * bi_base + u[j]
    q[j] = temp / v;
    rem = temp % v;
  }
 
  q.trim();
  return static_cast<digit>(rem);
}
 
void h_::div_algo_single(bi_t& q, bi_t& r, const bi_t& u,
                         const bi_t& v) noexcept {
  ddigit rem{0};
 
  for (size_t j = u.size() - 1; j < std::numeric_limits<size_t>::max(); --j) {
    const ddigit temp = (rem << bi_dwidth) | u[j];  // rem * bi_base + u[j]
    q[j] = temp / v[0];
    rem = temp % v[0];
  }
 
  r[0] = static_cast<digit>(rem);
 
  q.trim();
  r.trim();
}
 
void h_::div_algo_binary(bi_t& q, bi_t& r, const bi_t& u, const bi_t& v) {
  // (2)
  q = 0;
  r = 0;
 
  // (3)
  for (unsigned long i = u.bit_length() - 1; i < ULONG_MAX; --i) {
    // (3)(I)
    r <<= 1;
 
    // (3)(II)
    if (u.test_bit(i)) {
      r.set_bit(0);
    }
 
    // (3)(III)
    if (h_::cmp_abs(r, v) >= 0) {
      // (3)(III)(a)
      h_::sub_abs(r, r, v);
      // (3)(III)(b)
      q.set_bit(i);
    }
  }
}
 
void h_::div_algo_knuth(bi_t& q, bi_t& r, const bi_t& u, const bi_t& v) {
#if defined(BI_DIGIT_64_BIT)
  using sddigit = __int128;
#else
  using sddigit = int64_t;
#endif
 
  const size_t m = u.size();
  const size_t n = v.size();
 
  /* (1) Normalize */
  // Find e such that b > 2^{e} * v_{n-1} >= b / 2
  constexpr digit b_half = bi_base >> 1;
  const digit v_msd = v[n - 1];
  int e = 0;
  while ((v_msd << e) < b_half) {
    ++e;
  }
  const unsigned compl_e = bi_dwidth - e;
 
  bi_t u_norm, v_norm;
  u_norm.resize_(m + 1);
  v_norm.resize_(n);
 
  // u_norm set to u * 2^{e}. Cast to ddigit as e may be 0 which would be UB
  u_norm[0] = u[0] << e;
  for (size_t i = 1; i < m; ++i) {
    u_norm[i] = (u[i] << e) | (static_cast<ddigit>(u[i - 1]) >> compl_e);
  }
  u_norm[m] = static_cast<ddigit>(u[m - 1]) >> compl_e;
 
  // v_norm set to v * 2^{e}
  v_norm[0] = v[0] << e;
  for (size_t i = 1; i < n; ++i) {
    v_norm[i] = (v[i] << e) | (static_cast<ddigit>(v[i - 1]) >> compl_e);
  }
 
  const digit vp = v_norm[n - 1], vpp = v_norm[n - 2];
  /* (2) Initialize j. Also (7) Loop on j */
  for (size_t j = m - n; j < std::numeric_limits<size_t>::max(); --j) {
    /* (3) Calculate q_hat */
    const ddigit tmp = u_norm[j + n] * bi_base + u_norm[j + n - 1];
    ddigit q_hat = tmp / vp;
    ddigit r_hat = tmp % vp;
 
    while (q_hat == bi_base ||
           q_hat * vpp > bi_base * r_hat + u_norm[j + n - 2]) {
      --q_hat;
      r_hat += vp;
      if (r_hat >= bi_base)
        break;
    }
 
    /* (4) Multiply and subtract */
    sddigit b{};
    for (size_t i = 0; i < n; ++i) {
      const ddigit product = q_hat * v_norm[i];
      const sddigit stmp =
          static_cast<sddigit>(u_norm[j + i]) - static_cast<digit>(product) - b;
      u_norm[j + i] = static_cast<digit>(stmp);
      b = static_cast<sddigit>(product >> bi_dbits) - (stmp >> bi_dbits);
    }
    const bool neg{u_norm[j + n] < b};
    u_norm[j + n] = static_cast<digit>(u_norm[j + n] - b);
 
    /* (5) Test remainder */
    q[j] = q_hat;
 
    if (neg) {
      /* (6) Add back */
      // Decrease q_{j} by 1
      --q[j];
 
      // Add (0v_{n-1}...v_{0})_{b} to (u_{j+n}...u_{j})_{b}
      digit c = 0;
      for (size_t i = 0; i < n; ++i) {
        const ddigit tmp = (static_cast<ddigit>(u_norm[j + i]) + v_norm[i]) + c;
        u_norm[j + i] = tmp;
        c = tmp >> bi_dwidth;
      }
      // A carry will occur to the left of u_{j+n} and it should be ignored
      u_norm[j + n] += c;
      assert(u_norm[j + n] == 0);
    }
  }
 
  /* (8) Unnormalize */
  // (u_{n-1}...u_{0})_{b} / 2^e
  const size_t last_idx = n - 1;
  r[last_idx] = u_norm[last_idx] >> e;
  for (size_t i = last_idx; i-- > 0;) {
    r[i] = (static_cast<ddigit>(u_norm[i + 1]) << compl_e) | (u_norm[i] >> e);
  }
 
  q.trim();
  r.trim();
}
 
void h_::divide(bi_t& Q, bi_t& R, const bi_t& N, const bi_t& D) {
  if (D.size() == 0) {
    throw division_by_zero("Division by zero attempt.");
  }
 
  const size_t size_N = N.size();
  const size_t size_D = D.size();
 
  // |N| < |D| case
  if ((size_N == size_D && N[size_N - 1] < D[size_D - 1]) || size_N < size_D) {
    Q = 0;  // |N| < |D| ==> |N| / |D| < 1 ==> Q computes to 0
    R = N;  // Computation (a/b)*b + a%b should equal a ==> a%b is a
    return;
  }
 
  // TRUE: size_N >= size_D > 0
 
  // Ensure enough space in Q and R
  const size_t size_Q = size_N - size_D + 1;
  Q.resize_(size_Q);
  R.resize_(size_D);
 
  // Unsigned integer division algorithms
  if (size_D == 1) {
    // Knuth (Vol. 2, 4.3.1, p. 272) recommends using the algorithm used in
    // div_algo_single() when size_D is 1.
    div_algo_single(Q, R, N, D);
  } else {
    div_algo_knuth(Q, R, N, D);
  }
 
  Q.negative_ = D.negative() ^ N.negative();
  R.negative_ = R.size() > 0 && N.negative();
}
 
 
void h_::left_shift(bi_t& result, const bi_t& x, bi_bitcount_t n_bits) {
  if (x.size() == 0) {
    result.resize_unsafe_(0);
    result.negative_ = false;
    return;
  }
 
  if (n_bits == 0) {
    result = x;  // Copy assignment avoids copying if &result == &x
    return;
  }
 
  const size_t size_x = x.size();
  const bi_bitcount_t digit_shift = n_bits / bi_dwidth;
  const unsigned bit_shift = n_bits % bi_dwidth;
 
  const size_t size_result = size_x + digit_shift + (bit_shift ? 1 : 0);
  if (size_result < size_x) {
    throw overflow_error("Result size exceeds SIZE_MAX.");
  }
 
  // Create a temporary vector if `result` and `x` overlap
  dvector temp_vec;
  dvector& target_vec = (&result == &x) ? temp_vec : result.vec_;
 
  target_vec.resize(size_result);
 
  size_t i = 0;
  // Zero-fill vacated digits
  for (; i < digit_shift; ++i) {
    target_vec[i] = 0;
  }
 
  const auto compl_bit_shift = bi_dwidth - bit_shift;
 
  digit carry = 0;
  for (size_t j = 0; j < size_x; ++j, ++i) {
    const digit current = x[j];
    target_vec[i] = (current << bit_shift) | carry;
    carry = bit_shift ? (current >> compl_bit_shift) : 0;
  }
 
  if (bit_shift) {
    target_vec[i] = carry;
  }
 
  if (&result == &x) {
    result.vec_ = std::move(temp_vec);
  }
 
  result.trim();
  result.negative_ = x.negative();
}
 
void h_::right_shift(bi_t& result, const bi_t& x, bi_bitcount_t n_bits) {
  // TODO: maybe set `result = x` for special cases x.size() == 0 || n_bits == 0
  const size_t size_x = x.size();
 
  const bi_bitcount_t digit_shift = n_bits / bi_dwidth;
  const unsigned bit_shift = n_bits % bi_dwidth;
 
  if (size_x <= digit_shift) {
    if (x.negative()) {
      // Set result to -1
      result.resize_(1);
      result[0] = 1;
      result.negative_ = true;
    } else {
      // Set result to 0
      result.resize_unsafe_(0);
      result.negative_ = false;
    }
    return;
  }
 
  result.resize_(size_x - digit_shift);
 
  if (bit_shift == 0) {
    for (size_t i = 0; i < result.size(); ++i) {
      result[i] = x[i + digit_shift];
    }
  } else {
    const unsigned compl_bit_shift = bi_dwidth - bit_shift;
 
    for (size_t i = 0, k = digit_shift + 1; i < result.size(); ++i, ++k) {
      result[i] = x[k - 1] >> bit_shift;
      if (k < size_x) {
        result[i] |= x[k] << compl_bit_shift;
      }
    }
  }
 
  result.trim();
  result.negative_ = x.negative();
 
  // At this point, the result is trunc(x / 2^{n_bits}) for all x
 
  // Adjust for floor division for negative numbers
  if (x.negative()) {
    bool subtract_one = false;
    if (bit_shift > 0) {
      const digit mask = (static_cast<digit>(1) << bit_shift) - 1;
      if ((x[digit_shift] & mask) != 0) {
        subtract_one = true;
      }
    }
    if (!subtract_one && digit_shift > 0) {
      for (size_t i = 0; i < digit_shift; ++i) {
        if (x[i] != 0) {
          subtract_one = true;
          break;
        }
      }
    }
    if (subtract_one) {
      --result;
    }
  }
}
 
 
template <h_::BitwiseOperation Op>
void h_::bitwise_operation_impl(bi_t& result, const bi_t& x, const bi_t& y) {
  size_t max_size = std::max(x.size(), y.size());
 
  const dvector& lhs_digits = x.negative_ ? to_twos_complement(x.vec_) : x.vec_;
  const dvector& rhs_digits = y.negative_ ? to_twos_complement(y.vec_) : y.vec_;
 
  bool result_negative;  // NOLINT(cppcoreguidelines-init-variables)
  if constexpr (Op == BitwiseOperation::AND) {
    result_negative = x.negative_ && y.negative_;
  } else if constexpr (Op == BitwiseOperation::OR) {
    result_negative = x.negative_ || y.negative_;
  } else if constexpr (Op == BitwiseOperation::XOR) {
    result_negative = x.negative_ != y.negative_;
    ++max_size;
  }
 
  result.resize_(max_size);
 
  for (size_t i = 0; i < result.size(); ++i) {
    const digit lhs_digit =
        i < lhs_digits.size() ? lhs_digits[i] : (x.negative_ ? bi_dmax : 0);
    const digit rhs_digit =
        i < rhs_digits.size() ? rhs_digits[i] : (y.negative_ ? bi_dmax : 0);
 
    if constexpr (Op == BitwiseOperation::AND) {
      result[i] = lhs_digit & rhs_digit;
    } else if constexpr (Op == BitwiseOperation::OR) {
      result[i] = lhs_digit | rhs_digit;
    } else if constexpr (Op == BitwiseOperation::XOR) {
      result[i] = lhs_digit ^ rhs_digit;
    }
  }
 
  result.negative_ = result_negative;
 
  if (result.negative_) {
    to_twos_complement_in_place(result.vec_);
  }
 
  result.trim();
}
 
template <h_::BitwiseOperation Op>
bi_t h_::bitwise_operation(const bi_t& a, const bi_t& b) {
  bi_t result;
  bitwise_operation_impl<Op>(result, a, b);
  return result;
}
 
template <h_::BitwiseOperation Op>
bi_t& h_::ibitwise_operation(bi_t& x, const bi_t& other) {
  bitwise_operation_impl<Op>(x, x, other);
  return x;
}
 
// Assumes x, y >= 0
int h_::cmp_abs(const bi_t& x, const bi_t& y) noexcept {
  if (x.size() > y.size()) {
    return 1;
  } else if (x.size() < y.size()) {
    return -1;
  } else {
    for (size_t i = x.size(); i-- > 0;) {
      if (x.vec_[i] > y.vec_[i]) {
        return 1;
      } else if (x.vec_[i] < y.vec_[i]) {
        return -1;
      }
    }
    return 0;
  }
}
 
int h_::cmp(const bi_t& x, const bi_t& y) noexcept {
  if (!x.negative() && y.negative()) {
    // x >= 0, y < 0
    return 1;
  } else if (x.negative() && !y.negative()) {
    // x < 0, y >= 0
    return -1;
  } else if (x.negative() && y.negative()) {
    // x, y < 0
    if (x.size() > y.size()) {
      return -1;
    } else if (x.size() < y.size()) {
      return 1;
    } else {
      for (size_t i = x.size(); i-- > 0;) {
        if (x.vec_[i] < y.vec_[i]) {
          return 1;
        } else if (x.vec_[i] > y.vec_[i]) {
          return -1;
        }
      }
      return 0;
    }
  } else {
    // x, y >= 0
    return h_::cmp_abs(x, y);
  }
}
 
template <std::integral T>
int h_::cmp(const bi_t& a, T b) noexcept {
  if (b == 0) {
    if (a.size() == 0) {
      return 0;  // a
    }
    return a.negative() ? -1 : 1;  // b
  }
 
  using UnsignedT = typename std::make_unsigned<T>::type;
  const bool b_negative = b < 0;
  const UnsignedT unsigned_b =
      b_negative ? -static_cast<UnsignedT>(b) : static_cast<UnsignedT>(b);
 
  if (a.negative() && !b_negative) {
    return -1;  // c
  }
  if (!a.negative() && b_negative) {
    return 1;  // d
  }
 
  size_t n_b_digits = 0;
  if constexpr (std::numeric_limits<UnsignedT>::max() <= bi_dmax) {
    n_b_digits = (unsigned_b != 0) ? 1 : 0;
  } else {
    UnsignedT temp_b = unsigned_b;
    while (temp_b != 0) {
      temp_b >>= bi_dwidth;
      n_b_digits++;
    }
  }
 
  const size_t a_size = a.size();
  if (a_size < n_b_digits) {
    return a.negative() ? 1 : -1;  // e
  }
 
  if (a_size > n_b_digits) {
    return a.negative() ? -1 : 1;  // f
  }
 
  for (size_t i = n_b_digits; i-- > 0;) {
    const digit a_digit = a[i];
    const digit b_digit = static_cast<digit>(unsigned_b >> (bi_dwidth * i));
 
    if (a_digit < b_digit) {
      return a.negative() ? 1 : -1;  // g
    }
    if (a_digit > b_digit) {
      return a.negative() ? -1 : 1;  // h
    }
  }
  return 0;  // i
}
 
 
uint8_t h_::idiv10(bi_t& x) noexcept {
  constexpr ddigit base = static_cast<ddigit>(bi_dmax) + 1;
  uint8_t carry = 0;
 
  // NOLINTBEGIN(cppcoreguidelines-avoid-magic-numbers)
  for (size_t i = x.size(); i-- > 0;) {
    const ddigit current = (base * carry) + x.vec_[i];
    x.vec_[i] = static_cast<digit>(current / 10);
    carry = static_cast<uint8_t>(current % 10);
  }
  // NOLINTEND(cppcoreguidelines-avoid-magic-numbers)
 
  x.trim();
  return carry;
}
 
// log_base_2[base] gives log_{base}(2) (with some rounding up)
constexpr std::array<double, 37> log_base_2 = {
    0,     0,     1.0,   0.631, 0.5,   0.431, 0.387, 0.357, 0.334, 0.316,
    0.302, 0.290, 0.279, 0.271, 0.263, 0.256, 0.25,  0.245, 0.240, 0.236,
    0.232, 0.228, 0.225, 0.222, 0.219, 0.216, 0.213, 0.211, 0.209, 0.206,
    0.204, 0.202, 0.200, 0.199, 0.197, 0.195, 0.194};
 
size_t h_::base_length(const bi_t& x, int base) {
  // TODO: a different exception is probably more appropriate in this func
  constexpr const char* s = "Digit estimation exceeds practical limit.";
 
  const bi_bitcount_t bitlen = x.bit_length();
  if (bitlen > dbl_max_int) {
    throw overflow_error(s);
  }
 
  const bi_bitcount_t r = static_cast<bi_bitcount_t>(
      // NOLINTNEXTLINE(cppcoreguidelines-pro-bounds-constant-array-index)
      std::floor(static_cast<double>(bitlen) * log_base_2[base]) + 1);
  if (r > std::numeric_limits<size_t>::max() - 2) {
    throw overflow_error(s);
  }
 
  return r;
}
 
 
template <std::integral T>
void h_::init_one_digit(bi_t& x, T value) {
  x.resize_(1);
 
  if (value >= 0) {
    x[0] = static_cast<digit>(value);
    x.negative_ = false;
  } else {
    x[0] = static_cast<digit>(0) - static_cast<digit>(value);
    x.negative_ = true;
  }
}
 
template <std::integral T>
void h_::init_atleast_one_digit(bi_t& x, T value) {
  using UnsignedT = typename std::make_unsigned<T>::type;
 
  size_t size = 1;
  UnsignedT temp, uvalue;
 
  if (value >= 0) {
    uvalue = value;
    x.negative_ = false;
  } else {
    uvalue = static_cast<UnsignedT>(0) - static_cast<UnsignedT>(value);
    x.negative_ = true;
  }
 
  temp = uvalue;
 
  while (temp >>= bi_dwidth) {
    ++size;
  }
 
  x.resize_(size);
 
  for (size_t i = 0; uvalue != 0; ++i) {
    x[i] = static_cast<digit>(uvalue);
    uvalue >>= bi_dwidth;
  }
}
 
// TODO: move this somewhere else.
 
 
constexpr auto pow2_8 = 1 << 8;
 
constexpr std::array<uint8_t, pow2_8> create_char_to_b36_map() {
  // NOLINTBEGIN(cppcoreguidelines-avoid-magic-numbers)
  constexpr auto n_digits = 10;
  constexpr auto n_letters = 26;
  constexpr std::array<char, n_letters> lowercase{
      'a', 'b', 'c', 'd', 'e', 'f', 'g', 'h', 'i', 'j', 'k', 'l', 'm',
      'n', 'o', 'p', 'q', 'r', 's', 't', 'u', 'v', 'w', 'x', 'y', 'z'};
  constexpr std::array<char, n_letters> uppercase{
      'A', 'B', 'C', 'D', 'E', 'F', 'G', 'H', 'I', 'J', 'K', 'L', 'M',
      'N', 'O', 'P', 'Q', 'R', 'S', 'T', 'U', 'V', 'W', 'X', 'Y', 'Z'};
 
  std::array<uint8_t, pow2_8> map{};
  map.fill(0xFF);
 
  for (uint8_t i = 0; i < n_digits; ++i) {
    map.at(static_cast<uint8_t>('0') + i) = i;
  }
 
  for (uint8_t i = 0; i < n_letters; ++i) {
    map.at(static_cast<uint8_t>(lowercase.at(i))) = 10 + i;
    map.at(static_cast<uint8_t>(uppercase.at(i))) = 10 + i;
  }
 
  return map;
  // NOLINTEND(cppcoreguidelines-avoid-magic-numbers)
}
 
constexpr auto char_to_int_map = create_char_to_b36_map();
 
constexpr uint8_t char_to_b36(char ch) {
  // NOLINTNEXTLINE(cppcoreguidelines-pro-bounds-constant-array-index)
  return char_to_int_map[static_cast<uint8_t>(ch)];
}
 
// Calculate base ** n
constexpr digit pow(int base, size_t n) {
  digit result = 1;
  for (size_t i = 0; i < n; ++i) {
    result *= base;
  }
  return result;
}
 
constexpr int calculate_max_batch_size(digit base) {
  int n = 0;
 
  ddigit b = base;
  while (b < bi_dmax) {
    b *= base;
    ++n;
  }
 
  return n;
}
 
struct BaseMBS {
  // max batch size
  unsigned mbs;
  // base ** mbs
  digit base_pow_mbs;
};
 
// For example, if digit <==> uint32_t (uint64_t), 10^{9} (10^{19}) is the
// highest power of 10 that fits in it.
// NOLINTBEGIN(cppcoreguidelines-avoid-magic-numbers)
constexpr std::array<BaseMBS, 37> create_base_mbs_array() {
  std::array<BaseMBS, 37> base_mbs{};
  for (int base = 2; base <= 36; ++base) {
    unsigned max_batch_size = calculate_max_batch_size(base);
    base_mbs.at(base) = {max_batch_size, pow(base, max_batch_size)};
  }
  return base_mbs;
}
// NOLINTEND(cppcoreguidelines-avoid-magic-numbers)
 
constexpr auto base_mbs = create_base_mbs_array();
 
void h_::init_string(bi_t& x, const std::string& s, int base) {
  // NOLINTNEXTLINE(cppcoreguidelines-avoid-magic-numbers)
  if (base <= 1 || base > 36) {
    throw std::invalid_argument("base argument must be in [2, 36]");
  }
 
  // NOLINTNEXTLINE(cppcoreguidelines-pro-bounds-constant-array-index)
  const auto [max_batch_size, base_pow_max_batch_size] = base_mbs[base];
 
  // std::stoi and company allow leading whitespace and a plus/minus sign. We
  // follow suit.
 
  // Allow leading whitespace
  auto it = std::find_if_not(s.begin(), s.end(),
                             [](char ch) { return std::isspace(ch); });
 
  // Allow plus/minus sign to precede first base-`base` digit
  x.negative_ = false;
  if (it != s.end()) {
    if (*it == '-') {
      x.negative_ = true;
      ++it;
    } else if (*it == '+') {
      ++it;
    }
  }
 
  const auto start_digit = it;
  // TODO: isalnum() is too broad!
  it = std::find_if_not(start_digit, s.end(),
                        [](char ch) { return std::isalnum(ch); });
 
  if (start_digit == it) {
    throw std::invalid_argument("Invalid string format.");
  }
 
  size_t n_base = std::distance(start_digit, it);  // it - start_digit
  const size_t n_digits = uints::div_ceil(n_base, max_batch_size);
 
  x.reserve_(n_digits);
  x.resize_unsafe_(0);
 
  auto dec_it = start_digit;
  const size_t rem_batch_size = n_base % max_batch_size;
 
  // Initialize batch value
  digit batch = 0;
  // Convert batch substring to integer value
  for (size_t j = 0; j < rem_batch_size; ++j, ++dec_it) {
    // NOLINTNEXTLINE(cppcoreguidelines-pro-bounds-constant-array-index)
    batch = char_to_int_map[*dec_it] + batch * base;
  }
  // x is initially zero, so multiplication is zero in h_::imul1add1()
  if (batch) {
    x.vec_.push_back(batch);
  }
 
  while (dec_it < it) {
    // Initialize batch value
    batch = 0;
    // Convert batch substring to integer value
    for (size_t j = 0; j < max_batch_size; ++j, ++dec_it) {
      // NOLINTNEXTLINE(cppcoreguidelines-pro-bounds-constant-array-index)
      batch = char_to_int_map[*dec_it] + batch * base;
    }
    h_::imul1add1(x, base_pow_max_batch_size, batch);
  }
 
  x.trim();
}
 
 
inline dvector h_::to_twos_complement(const dvector& vec) {
  const size_t vec_size = vec.size();
  dvector result;
  result.resize(vec_size);
 
  bool carry = true;
  for (size_t i = 0; i < vec_size; ++i) {
    const digit sum = ~vec[i] + carry;
    carry = sum < static_cast<uint8_t>(carry);
    result[i] = sum;
  }
 
  return result;
}
 
inline void h_::to_twos_complement_in_place(dvector& vec) noexcept {
  const size_t vec_size = vec.size();
 
  bool carry = true;
  for (size_t i = 0; i < vec_size; ++i) {
    const digit sum = ~vec[i] + carry;
    carry = sum < static_cast<uint8_t>(carry);
    vec[i] = sum;
  }
}
 
void h_::assign_from_double(bi_t& x, double d) {
  if (std::isnan(d) || std::isinf(d)) {
    throw from_float(
        "Conversion error: NaN or infinity cannot be converted to an integer.");
  }
 
  if (d > -1 && d < 1) {
    x = 0;
    return;
  }
 
  const bool neg = d < 0;
  if (neg) {
    d = -d;
  }
 
  size_t n_digits{1};
  while (bi_base_dbl <= d) {
    // Equiv. to `d /= bi_base_dbl;`, but multiplication is generally cheaper
    d *= bi_base_dbl_reciprocal;
    ++n_digits;
  }
 
  x.resize_(n_digits);
 
  for (size_t i = n_digits; i-- > 0;) {
    x[i] = static_cast<digit>(d);
    d = (d - static_cast<double>(x[i])) * bi_base_dbl;
  }
 
  x.negative_ = neg;
  x.trim();
}
 
int h_::cmp_abs(const bi_t& z, double dbl) noexcept {
  if (z.size() == 0) {
    return dbl > 0 ? -1 : 0;
  }
 
  if (z.size() >= bi_cmp_dbl_size_upper) {
    return 1;
  }
 
  // Same as dbl = std::ldexp(dbl, -static_cast<int>(bi_dbits) *
  //                                static_cast<int>(z.size() - 1));
  for (size_t j = 1; j < z.size(); ++j) {
    // Equiv. to `dbl /= bi_base_dbl;`, but multiplication is generally cheaper
    dbl *= bi_base_dbl_reciprocal;
  }
 
  if (bi_base_dbl <= dbl) {
    return -1;
  }
 
  for (auto it = z.rbegin(); it != z.rend(); ++it) {
    const digit cur = static_cast<digit>(dbl);
 
    if (*it < cur) {
      return -1;
    } else if (*it > cur) {
      return 1;
    } else {
      dbl = (dbl - static_cast<double>(cur)) * bi_base_dbl;
    }
  }
 
  return dbl > 0 ? -1 : 0;
}
 
int h_::cmp(const bi_t& z, double dbl) noexcept {
  if (!z.negative()) {
    return dbl < 0 ? 1 : h_::cmp_abs(z, dbl);
  } else {
    return dbl < 0 ? -h_::cmp_abs(z, -dbl) : -1;
  }
}
 
 
template <std::unsigned_integral T>
bi_t h_::expo_left_to_right(const bi_t& base, T exp) {
  assert(exp > 0);
 
  // (1)
  bi_t ret{1};
 
  // (2)
  for (int j = uints::bit_length(exp); j-- > 0;) {
    ret *= ret;
    // Test if j-th bit of exp is set
    if (exp & (static_cast<T>(1) << j)) {
      ret *= base;
    }
  }
 
  // (3)
  return ret;
}
 
 
bi_t h_::random_(bi_bitcount_t z) {
  bi_t result{};
 
  if (z == 0) {
    return result;
  }
 
  size_t num_digits = z / bi_dbits;
  size_t remainder = z % bi_dbits;
 
  if (remainder != 0) {
    num_digits += 1;
  }
 
  result.resize_(num_digits);
 
  std::uniform_int_distribution<digit> dist(0, bi_dmax);
 
  for (size_t i = 0; i < num_digits - 1; ++i) {
    result[i] = dist(rng_);
  }
 
  if (remainder != 0) {
    const digit mask = (static_cast<digit>(1) << remainder) - 1;
    result[num_digits - 1] = dist(rng_) & mask;
  } else {
    result[num_digits - 1] = dist(rng_);
  }
 
  result.trim();
  return result;
}
 
}  // namespace bi
 
#endif  // BI_SRC_H__HPP_