newton_bignum_div_to_s.rb

def big_divmod(numerator, denominator)
  numerator_bits = numerator.bit_length
  denominator_bits = denominator.bit_length
  inv = 1
  accuracy_bits = 1
  target_bits = [numerator_bits - denominator_bits]
  while target_bits[0] > 1
    target_bits.unshift (target_bits[0] * 4 + 1)/ 7
  end
  accuracy_bits = target_bits.shift
  until target_bits.empty?
    # Caluculate X = 1 / denominator
    # Repeat X = X * (2 - denominator * X) until X is accurate enough
    # X is represented as X = inv / (1<<denominator_bits) / (1<<accuracy_bits)
    # t = Time.now
    next_accuracy_bits = target_bits.shift
    dshift = [denominator_bits - next_accuracy_bits, 0].max
    deno = denominator >> dshift
    mul = (2 << (denominator_bits + accuracy_bits - dshift)) - inv * deno
    inv = inv * mul >> (denominator_bits + 2 * accuracy_bits - next_accuracy_bits - dshift)
    accuracy_bits = next_accuracy_bits
    # p t: Time.now - t
  end
  # t = Time.now
  candidate = ((numerator >> denominator_bits) * inv) >> accuracy_bits
  # p candidate_t: Time.now - t
  # t = Time.now
  diff, mod = (numerator - candidate * denominator).divmod denominator
  # p final_diff: diff, t: Time.now - t
  [candidate + diff, mod]
end