divmod_by_mult.rb

# returns (1<<(size + x.bit_length)) / x
def inv(x, size)
  x_size = x.bit_length
  bl = size.bit_length
  n = 2
  y = (1 << (n + 2)) / (x >> (x_size - 2))
  (0..bl).reverse_each do |i|
    n2 = [(size >> i) + 2, size].min
    # y = (1<<(x_size + n)) / x
    # y = 1/x * (1<<(x_size+n))
    # y2 = (1<<x_size+n2)*inv2
    x_shift = [x_size - n2, 0].max
    y = ((y << (x_size + n - x_shift)) + ((1 << (x_size + n - x_shift)) - (x >> x_shift) * y) * y) >> (x_size + 2 * n - n2 - x_shift)
    n = n2
  end
  y
end

def divmod_by_inv(x, y, inv, inv_bits)
  # inv = (1<<(inv_bits + y.bit_length)) / y
  div = ((x >> y.bit_length) * inv) >> inv_bits
  mod = x - div * y
  while mod < 0
    mod += y
    div -= 1
  end
  while mod >= y
    mod -= y
    div += 1
  end
  [div, mod]
end

def divmod(x, y)
  x_bits = x.bit_length
  y_bits = y.bit_length
  n = [(x_bits - y_bits).fdiv(y_bits).round, 1].max
  bits = (x_bits - y_bits + n - 1) / n
  yinv = inv(y, bits)
  ans = 0
  mod = x[n * bits, y_bits]
  (n - 1).downto(0) do |i|
    a = (mod << bits) | x[i * bits, bits]
    div, mod = divmod_by_inv(a, y, yinv, bits)
    ans <<= bits
    ans |= div
  end
  [ans, mod]
end

def split_by_bits(x, n)
  if x.bit_length >= 4 * n
    return split_by_bits(x, 2 * n).flat_map { a,b=split_by_bits(it, n); [a, b||0] }
  end
  a = []
  f = (1 << n) - 1
  while x > 0
    a << (x & f)
    x >>= n
  end
  a<<0 if a.empty?
  a
end

def join_by_bits(values, n)
  while values.size > 1
    values = values.each_slice(2).map do
      _2 ? (_2 << n) | _1 : _1
    end
    n *= 2
  end
  values.first
end

def bdmultshift(x, y, shift)
  x_digits = x.n_significant_digits
  y_digits = y.n_significant_digits
  z_digits = x_digits + y_digits - shift
  xs = x_digits - z_digits - 2
  ys = y_digits - z_digits - 2
  if xs > 0
    shift -= xs
    x = bdrshift(x, xs)
  end
  if ys
    shift -= ys
    y = bdrshift(y, ys)
  end
  bdrshift(x * y, shift)
end
def bdlshift(x,n) = (x * BigDecimal("1e#{n}")).fix
def bdrshift(x,n) = bdlshift(x, -n)

# returns 10**(size + x.exponent) / x where x is a BigDecimal integer

def bdinv(x, size)
  x_size = x.exponent
  bl = size.bit_length
  n = 2
  y = (bdlshift(1, n + 2) / bdrshift(x, x_size - 2)).fix
  (0..bl).reverse_each do |i|
    n2 = [(size >> i) + 2, size].min
    # y = (1<<(x_size + n)) / x
    # y = 1/x * (1<<(x_size+n))
    # y2 = (1<<x_size+n2)*inv2
    x_shift = [x_size - n2, 0].max

    y = bdlshift(y, n2 - n) + bdmultshift(y, (bdlshift(1, x_size + n - x_shift) - bdrshift(x, x_shift) * y), x_size + 2 * n - n2 - x_shift)
    n = n2
  end
  y
end

def bddivmod(x, y)
  inv_digits = x.exponent - y.exponent
  inv = bdinv(y, inv_digits)
  div = bdmultshift(bdrshift(x, y.exponent), inv, inv_digits)
  mod = x - div * y
  while mod < 0
    mod += y
    div -= 1
  end
  while mod >= y
    mod -= y
    div += 1
  end
  [div, mod]
end

require 'bigdecimal'
n=100000;x=BigDecimal('789'*n);y=(BigDecimal('123'*(n/2)+'1'))
binding.irb