src/bigint/div-barrett.cc

// Copyright 2021 the V8 project authors. All rights reserved.
// Use of this source code is governed by a BSD-style license that can be
// found in the LICENSE file.

// Barrett division, finding the inverse with Newton's method.
// Reference: "Fast Division of Large Integers" by Karl Hasselström,
// found at https://treskal.com/s/masters-thesis.pdf

// Many thanks to Karl Wiberg, k@w5.se, for both writing up an
// understandable theoretical description of the algorithm and privately
// providing a demo implementation, on which the implementation in this file is
// based.

#include <algorithm>

#include "src/bigint/bigint-internal.h"
#include "src/bigint/digit-arithmetic.h"
#include "src/bigint/div-helpers.h"
#include "src/bigint/vector-arithmetic.h"

namespace v8 {
namespace bigint {

namespace {

void DcheckIntegerPartRange(Digits X, digit_t min, digit_t max) {
#if DEBUG
  digit_t integer_part = X.msd();
  DCHECK(integer_part >= min);
  DCHECK(integer_part <= max);
#else
  USE(X);
  USE(min);
  USE(max);
#endif
}

}  // namespace

// Z := (the fractional part of) 1/V, via naive division.
// See comments at {Invert} and {InvertNewton} below for details.
void ProcessorImpl::InvertBasecase(RWDigits Z, Digits V, RWDigits scratch) {
  DCHECK(Z.len() > V.len());
  DCHECK(V.len() > 0);
  DCHECK(scratch.len() >= 2 * V.len());
  int n = V.len();
  RWDigits X(scratch, 0, 2 * n);
  digit_t borrow = 0;
  int i = 0;
  for (; i < n; i++) X[i] = 0;
  for (; i < 2 * n; i++) X[i] = digit_sub2(0, V[i - n], borrow, &borrow);
  DCHECK(borrow == 1);
  RWDigits R(nullptr, 0);  // We don't need the remainder.
  if (n < kBurnikelThreshold) {
    DivideSchoolbook(Z, R, X, V);
  } else {
    DivideBurnikelZiegler(Z, R, X, V);
  }
}

// This is Algorithm 4.2 from the paper.
// Computes the inverse of V, shifted by kDigitBits * 2 * V.len, accurate to
// V.len+1 digits. The V.len low digits of the result digits will be written
// to Z, plus there is an implicit top digit with value 1.
// Needs InvertNewtonScratchSpace(V.len) of scratch space.
// The result is either correct or off by one (about half the time it is
// correct, half the time it is one too much, and in the corner case where V is
// minimal and the implicit top digit would have to be 2 it is one too little).
// Barrett's division algorithm can handle that, so we don't care.
void ProcessorImpl::InvertNewton(RWDigits Z, Digits V, RWDigits scratch) {
  const int vn = V.len();
  DCHECK(Z.len() >= vn);
  DCHECK(scratch.len() >= InvertNewtonScratchSpace(vn));
  const int kSOffset = 0;
  const int kWOffset = 0;  // S and W can share their scratch space.
  const int kUOffset = vn + kInvertNewtonExtraSpace;

  // The base case won't work otherwise.
  DCHECK(V.len() >= 3);

  constexpr int kBasecasePrecision = kNewtonInversionThreshold - 1;
  // V must have more digits than the basecase.
  DCHECK(V.len() > kBasecasePrecision);
  DCHECK(IsBitNormalized(V));

  // Step (1): Setup.
  // Calculate precision required at each step.
  // {k} is the number of fraction bits for the current iteration.
  int k = vn * kDigitBits;
  int target_fraction_bits[8 * sizeof(vn)];  // "k_i" in the paper.
  int iteration = -1;  // "i" in the paper, except inverted to run downwards.
  while (k > kBasecasePrecision * kDigitBits) {
    iteration++;
    target_fraction_bits[iteration] = k;
    k = DIV_CEIL(k, 2);
  }
  // At this point, k <= kBasecasePrecision*kDigitBits is the number of
  // fraction bits to use in the base case. {iteration} is the highest index
  // in use for f[].

  // Step (2): Initial approximation.
  int initial_digits = DIV_CEIL(k + 1, kDigitBits);
  Digits top_part_of_v(V, vn - initial_digits, initial_digits);
  InvertBasecase(Z, top_part_of_v, scratch);
  Z[initial_digits] = Z[initial_digits] + 1;  // Implicit top digit.
  // From now on, we'll keep Z.len updated to the part that's already computed.
  Z.set_len(initial_digits + 1);

  // Step (3): Precision doubling loop.
  while (true) {
    DcheckIntegerPartRange(Z, 1, 2);

    // (3b): S = Z^2
    RWDigits S(scratch, kSOffset, 2 * Z.len());
    Multiply(S, Z, Z);
    if (should_terminate()) return;
    S.TrimOne();  // Top digit of S is unused.
    DcheckIntegerPartRange(S, 1, 4);

    // (3c): T = V, truncated so that at least 2k+3 fraction bits remain.
    int fraction_digits = DIV_CEIL(2 * k + 3, kDigitBits);
    int t_len = std::min(V.len(), fraction_digits);
    Digits T(V, V.len() - t_len, t_len);

    // (3d): U = T * S, truncated so that at least 2k+1 fraction bits remain
    // (U has one integer digit, which might be zero).
    fraction_digits = DIV_CEIL(2 * k + 1, kDigitBits);
    RWDigits U(scratch, kUOffset, S.len() + T.len());
    DCHECK(U.len() > fraction_digits);
    Multiply(U, S, T);
    if (should_terminate()) return;
    U = U + (U.len() - (1 + fraction_digits));
    DcheckIntegerPartRange(U, 0, 3);

    // (3e): W = 2 * Z, padded with "0" fraction bits so that it has the
    // same number of fraction bits as U.
    DCHECK(U.len() >= Z.len());
    RWDigits W(scratch, kWOffset, U.len());
    int padding_digits = U.len() - Z.len();
    for (int i = 0; i < padding_digits; i++) W[i] = 0;
    LeftShift(W + padding_digits, Z, 1);
    DcheckIntegerPartRange(W, 2, 4);

    // (3f): Z = W - U.
    // This check is '<=' instead of '<' because U's top digit is its
    // integer part, and we want vn fraction digits.
    if (U.len() <= vn) {
      // Normal subtraction.
      // This is not the last iteration.
      DCHECK(iteration > 0);
      Z.set_len(U.len());
      digit_t borrow = SubtractAndReturnBorrow(Z, W, U);
      DCHECK(borrow == 0);
      USE(borrow);
      DcheckIntegerPartRange(Z, 1, 2);
    } else {
      // Truncate some least significant digits so that we get vn
      // fraction digits, and compute the integer digit separately.
      // This is the last iteration.
      DCHECK(iteration == 0);
      Z.set_len(vn);
      Digits W_part(W, W.len() - vn - 1, vn);
      Digits U_part(U, U.len() - vn - 1, vn);
      digit_t borrow = SubtractAndReturnBorrow(Z, W_part, U_part);
      digit_t integer_part = W.msd() - U.msd() - borrow;
      DCHECK(integer_part == 1 || integer_part == 2);
      if (integer_part == 2) {
        // This is the rare case where the correct result would be 2.0, but
        // since we can't express that by returning only the fractional part
        // with an implicit 1-digit, we have to return [1.]9999... instead.
        for (int i = 0; i < Z.len(); i++) Z[i] = ~digit_t{0};
      }
      break;
    }
    // (3g, 3h): Update local variables and loop.
    k = target_fraction_bits[iteration];
    iteration--;
  }
}

// Computes the inverse of V, shifted by kDigitBits * 2 * V.len, accurate to
// V.len+1 digits. The V.len low digits of the result digits will be written
// to Z, plus there is an implicit top digit with value 1.
// (Corner case: if V is minimal, the implicit digit should be 2; in that case
// we return one less than the correct answer. DivideBarrett can handle that.)
// Needs InvertScratchSpace(V.len) digits of scratch space.
void ProcessorImpl::Invert(RWDigits Z, Digits V, RWDigits scratch) {
  DCHECK(Z.len() > V.len());
  DCHECK(V.len() >= 1);
  DCHECK(IsBitNormalized(V));
  DCHECK(scratch.len() >= InvertScratchSpace(V.len()));

  int vn = V.len();
  if (vn >= kNewtonInversionThreshold) {
    return InvertNewton(Z, V, scratch);
  }
  if (vn == 1) {
    digit_t d = V[0];
    digit_t dummy_remainder;
    Z[0] = digit_div(~d, ~digit_t{0}, d, &dummy_remainder);
    Z[1] = 0;
  } else {
    InvertBasecase(Z, V, scratch);
    if (Z[vn] == 1) {
      for (int i = 0; i < vn; i++) Z[i] = ~digit_t{0};
      Z[vn] = 0;
    }
  }
}

// This is algorithm 3.5 from the paper.
// Computes Q(uotient) and R(emainder) for A/B using I, which is a
// precomputed approximation of 1/B (e.g. with Invert() above).
// Needs DivideBarrettScratchSpace(A.len) scratch space.
void ProcessorImpl::DivideBarrett(RWDigits Q, RWDigits R, Digits A, Digits B,
                                  Digits I, RWDigits scratch) {
  DCHECK(Q.len() > A.len() - B.len());
  DCHECK(R.len() >= B.len());
  DCHECK(A.len() > B.len());  // Careful: This is *not* '>=' !
  DCHECK(A.len() <= 2 * B.len());
  DCHECK(B.len() > 0);
  DCHECK(IsBitNormalized(B));
  DCHECK(I.len() == A.len() - B.len());
  DCHECK(scratch.len() >= DivideBarrettScratchSpace(A.len()));

  int orig_q_len = Q.len();

  // (1): A1 = A with B.len fewer digits.
  Digits A1 = A + B.len();
  DCHECK(A1.len() == I.len());

  // (2): Q = A1*I with I.len fewer digits.
  // {I} has an implicit high digit with value 1, so we add {A1} to the high
  // part of the multiplication result.
  RWDigits K(scratch, 0, 2 * I.len());
  Multiply(K, A1, I);
  if (should_terminate()) return;
  Q.set_len(I.len() + 1);
  Add(Q, K + I.len(), A1);
  // K is no longer used, can re-use {scratch} for P.

  // (3): R = A - B*Q (approximate remainder).
  RWDigits P(scratch, 0, A.len() + 1);
  Multiply(P, B, Q);
  if (should_terminate()) return;
  digit_t borrow = SubtractAndReturnBorrow(R, A, Digits(P, 0, B.len()));
  // R may be allocated wider than B, zero out any extra digits if so.
  for (int i = B.len(); i < R.len(); i++) R[i] = 0;
  digit_t r_high = A[B.len()] - P[B.len()] - borrow;

  // Adjust R and Q so that they become the correct remainder and quotient.
  // The number of iterations is guaranteed to be at most some very small
  // constant, unless the caller gave us a bad approximate quotient.
  if (r_high >> (kDigitBits - 1) == 1) {
    // (5b): R < 0, so R += B
    digit_t q_sub = 0;
    do {
      r_high += AddAndReturnCarry(R, R, B);
      q_sub++;
      DCHECK(q_sub <= 5);
    } while (r_high != 0);
    Subtract(Q, q_sub);
  } else {
    digit_t q_add = 0;
    while (r_high != 0 || GreaterThanOrEqual(R, B)) {
      // (5c): R >= B, so R -= B
      r_high -= SubtractAndReturnBorrow(R, R, B);
      q_add++;
      DCHECK(q_add <= 5);
    }
    Add(Q, q_add);
  }
  // (5a): Return.
  int final_q_len = Q.len();
  Q.set_len(orig_q_len);
  for (int i = final_q_len; i < orig_q_len; i++) Q[i] = 0;
}

// Computes Q(uotient) and R(emainder) for A/B, using Barrett division.
void ProcessorImpl::DivideBarrett(RWDigits Q, RWDigits R, Digits A, Digits B) {
  DCHECK(Q.len() > A.len() - B.len());
  DCHECK(R.len() >= B.len());
  DCHECK(A.len() > B.len());  // Careful: This is *not* '>=' !
  DCHECK(B.len() > 0);

  // Normalize B, and shift A by the same amount.
  ShiftedDigits b_normalized(B);
  ShiftedDigits a_normalized(A, b_normalized.shift());
  // Keep the code below more concise.
  B = b_normalized;
  A = a_normalized;

  // The core DivideBarrett function above only supports A having at most
  // twice as many digits as B. We generalize this to arbitrary inputs
  // similar to Burnikel-Ziegler division by performing a t-by-1 division
  // of B-sized chunks. It's easy to special-case the situation where we
  // don't need to bother.
  int barrett_dividend_length = A.len() <= 2 * B.len() ? A.len() : 2 * B.len();
  int i_len = barrett_dividend_length - B.len();
  ScratchDigits I(i_len + 1);  // +1 is for temporary use by Invert().
  int scratch_len =
      std::max(InvertScratchSpace(i_len),
               DivideBarrettScratchSpace(barrett_dividend_length));
  ScratchDigits scratch(scratch_len);
  Invert(I, Digits(B, B.len() - i_len, i_len), scratch);
  if (should_terminate()) return;
  I.TrimOne();
  DCHECK(I.len() == i_len);
  if (A.len() > 2 * B.len()) {
    // This follows the variable names and and algorithmic steps of
    // DivideBurnikelZiegler().
    int n = B.len();  // Chunk length.
    // (5): {t} is the number of B-sized chunks of A.
    int t = DIV_CEIL(A.len(), n);
    DCHECK(t >= 3);
    // (6)/(7): Z is used for the current 2-chunk block to be divided by B,
    // initialized to the two topmost chunks of A.
    int z_len = n * 2;
    ScratchDigits Z(z_len);
    PutAt(Z, A + n * (t - 2), z_len);
    // (8): For i from t-2 downto 0 do
    int qi_len = n + 1;
    ScratchDigits Qi(qi_len);
    ScratchDigits Ri(n);
    // First iteration unrolled and specialized.
    {
      int i = t - 2;
      DivideBarrett(Qi, Ri, Z, B, I, scratch);
      if (should_terminate()) return;
      RWDigits target = Q + n * i;
      // In the first iteration, all qi_len = n + 1 digits may be used.
      int to_copy = std::min(qi_len, target.len());
      for (int j = 0; j < to_copy; j++) target[j] = Qi[j];
      for (int j = to_copy; j < target.len(); j++) target[j] = 0;
#if DEBUG
      for (int j = to_copy; j < Qi.len(); j++) {
        DCHECK(Qi[j] == 0);
      }
#endif
    }
    // Now loop over any remaining iterations.
    for (int i = t - 3; i >= 0; i--) {
      // (8b): If i > 0, set Z_(i-1) = [Ri, A_(i-1)].
      // (De-duped with unrolled first iteration, hence reading A_(i).)
      PutAt(Z + n, Ri, n);
      PutAt(Z, A + n * i, n);
      // (8a): Compute Qi, Ri such that Zi = B*Qi + Ri.
      DivideBarrett(Qi, Ri, Z, B, I, scratch);
      DCHECK(Qi[qi_len - 1] == 0);
      if (should_terminate()) return;
      // (9): Return Q = [Q_(t-2), ..., Q_0]...
      PutAt(Q + n * i, Qi, n);
    }
    Ri.Normalize();
    DCHECK(Ri.len() <= R.len());
    // (9): ...and R = R_0 * 2^(-leading_zeros).
    RightShift(R, Ri, b_normalized.shift());
  } else {
    DivideBarrett(Q, R, A, B, I, scratch);
    if (should_terminate()) return;
    RightShift(R, R, b_normalized.shift());
  }
}

}  // namespace bigint
}  // namespace v8