Integer sqrt in stdlib

Crystal Contrib
8

I used: n.to_s(2).size to find the bit length because it was
the most reliable that worked for all number sizes for n.

Compile the following code as:

crystal run bitslengthtest.cr -Ddisable_overflow --release

def bit_length(n : BigInt); Math.log2(n).ceil.to_i end

puts
puts "for: def bit_length(n : BigInt); Math.log2(n).ceil.to_i"
puts
exp = 308; n= 10.to_big_i ** exp; puts "for n = (10**#{exp}), n.to_s(2).size = #{n.to_s(2).size}, bit_length(n) = #{bit_length(n)}, Math.log2(n) = #{Math.log2(n)}"
puts
exp = 309; n= 10.to_big_i ** exp; puts "for n = (10**#{exp}), n.to_s(2).size = #{n.to_s(2).size}, bit_length(n) = #{bit_length(n)}, Math.log2(n) = #{Math.log2(n)}"
puts
exp = 310; n= 10.to_big_i ** exp; puts "for n = (10**#{exp}), n.to_s(2).size = #{n.to_s(2).size}, bit_length(n) = #{bit_length(n)}, Math.log2(n) = #{Math.log2(n)}"
puts
exp = 3010; n= 10.to_big_i ** exp; puts "for n = (10**#{exp}), n.to_s(2).size = #{n.to_s(2).size}, bit_length(n) = #{bit_length(n)}, Math.log2(n) = #{Math.log2(n)}"

It produces this output:

for: def bit_length(n : BigInt); Math.log2(n).ceil.to_i end

for n = (10**308), n.to_s(2).size = 1024, bit_length(n) = 1024, Math.log2(n) = 1023.1538532253076

for n = (10**309), n.to_s(2).size = 1027, bit_length(n) = -2147483648, Math.log2(n) = Infinity

for n = (10**310), n.to_s(2).size = 1030, bit_length(n) = -2147483648, Math.log2(n) = Infinity

for n = (10**3010), n.to_s(2).size = 10000, bit_length(n) = -2147483648, Math.log2(n) = Infinity

Again, any operation using floats will beome inaccurate past a certain (large) size for n.

Ruby has method bit_length.

(10**308).bit_length  => 1024 
(10**309).bit_length  => 1027 
(10**310).bit_length  => 1030 
(10**3010).bit_length => 10000

See its implemenation in https://github.com/ruby/ruby/blob/master/numeric.c starting at line 4894.

Since any representable integer value has a fixed binary representation its bit_length can be directly assessed from it.

Actually, the same issue arises for taking the nth root of any integer.
The sqrt of the Newton method is a specific case for any Integer nth roots.
A little bit more code can provide for users a huge amount of functionality not in other languages.

require "big"

module IntRoots

  def isqrt(n = self)   # Newton's method for integer sqrt
    raise "ivalid negative integer" if n < 0
    return n if n < 2
    b = n.to_s(2).size
    one = typeof(self).new(1)  # value 1 of type self
    x = one << ((b-1) >> 1) | n >> ((b >> 1) + 1)
    while (t = n // x) < x; x = (x + t) >> 1 end
    x
  end

  def irootn(n)        # Newton's method for integer nth root
    return nil if self < 0 && n.even?
    raise "root n not an Integer >= 2" unless n.is_a?(Int) && n > 1
    return self if self == 0 || (self == -1 && n.odd?)
    num = self.abs
    one = typeof(self).new(1)  # value 1 of type self
    e, x = n-1, one << (num.to_s(2).size-1) // n + 1
    while (t = (e * x + num // x ** e) // n) < x
      x = t 
    end
    x *= self < 0 ? -1 : 1
  end
end

struct Int; include IntRoots end

(10.to_big_i ** 517).iroot(13) => 5878016072274912830666108305036555474358

I don’t know if GMP library provides these methods.
But here, all we need is BigInt and no other dependencies.

I’ll create an issues topic in github to discuss all the details around this.