ruby-core

Issue #13219 has been updated by Jabari Zakiya.


Let me attempt to address some of the issues/questions that have recently been presented.
Excuse me if any of this redundant, at the time I am writing.

1) Fastest algorithm

After extensive searching and testing, the fastest and most flexible algorithm
to exactly compute not only the integer squareroot, but also any root, is the
one I'll just name the "binary bit method" (**bbm**). See benchmarks below of code.

The methods **iroot2** and **irootn** are now in my **roots** gem and use **bbm**.

I encourage people to take this code and include any other proposed methods to
create a reference test suite for assessing performance (assuming the other methods
correctly compute the exact root values).

```
class Integer
  def irootn(n)
    return nil if self < 0 && n.even?
    raise "root n is < 2 or not an Integer" unless n.is_a?(Integer) && n > 1
    num  = self.abs
    bits_shift = (num.bit_length)/n + 2
    root, bitn_mask = 0, (1 << bits_shift)
    until (bitn_mask >>= 1) == 0
      root |= bitn_mask
      root ^= bitn_mask if root**n > num
    end
    root *= (self < 0 ? -1 : 1)
  end

  def iroot2; irootn(2) end
end


def intsqrt7(n)  # Newton's method
  x = n
  y = (n + 1)/2
  while y < x
    x = y
    y = (x + n/x)/2
  end
  x
end

# from http://www.codecodex.com/wiki/Calculate_an_integer_square_root
def intsqrt8(square)  
    # Check our input  
    square = square.to_i # force integer  
    return 0 if square == 0 # quick exit  
    raise RangeError if square < 0  
    # Actual computation  
    n = iter(1, square)  
    n1 = iter(n, square)  
    n1, n = iter(n1, square), n1 until n1 >= n - 1  
    n1 = n1 - 1 until n1*n1 <= square  
    return n1  
end  
def iter(n, square) (n+(square/n))/2 end

require "bigdecimal"
require "benchmark/ips"

Benchmark.ips do |x|
  n = 10**40
  puts "integer squareroot tests for n = #{n}"
  x.report("iroot2"     ) { n.iroot2    }
  x.report("irootn(2)"  ) { n.irootn(2) }
  x.report("newtons"    ) { intsqrt7(n) }
  x.report("codecodex"  ) { intsqrt8(n) }
  x.compare!
end
``` 
All will correctly compute the exact integer squareroot, but bbm is clearly faster.

```
2.4.0 :024 > load 'irootstest.rb'
integer squareroot tests for n = 10000000000000000000000000000000000000000
Warming up --------------------------------------
              iroot2     4.368k i/100ms
           irootn(2)     4.375k i/100ms
             newtons     3.491k i/100ms
           codecodex     2.880k i/100ms
Calculating -------------------------------------
              iroot2     43.969k (賊 3.0%) i/s -    222.768k in   5.071034s
           irootn(2)     43.990k (賊 3.1%) i/s -    223.125k in   5.077264s
             newtons     35.369k (賊 2.9%) i/s -    178.041k in   5.038327s
           codecodex     29.093k (賊 3.0%) i/s -    146.880k in   5.053165s

Comparison:
           irootn(2):    43990.3 i/s
              iroot2:    43969.0 i/s - same-ish: difference falls within error
             newtons:    35369.0 i/s - 1.24x  slower
           codecodex:    29093.5 i/s - 1.51x  slower

 => true 
```

2) It is well understood in numerical analysis, Number Theory, cryptography, etc,
an integer nth_root n is the largest integer such that root*n <= num. It is not a
floating point value, nor just the truncation of a floating point value. A 
calculator "root" is, by definition, a floating point value, and is an approximation
of some possible exact value, to the precision of the system to produce it.

Other languages have at least the equivalent of an **isqrt** function that work
strictly on integers, not floats, for the reasons given above.

3) Usefulness in advanced mathematical and numerical analysis fields.

Most people/programmers probably first (and mostly) think of Ruby as a Rails (web)
language. Whereas Python is being promoted and used widely now in a number of
mathematical and numerical analysis based projects.

https://jupyter.org/
https://ipython.org/
https://github.com/jupyterhub/jupyterhub
http://www.sagemath.org/
https://en.wikipedia.org/wiki/SageMath

I think Ruby is way easier to use (in general) for these applications, but it doesn't
have anywhere near the breadth of libraries and focus from its user base to be applied
to them. So Python is being used, basically by default, because it's easier to program
in than C/C++, and being made fast enough to satisfy more users. I think as a focus for 
Ruby 3x3 it can start to create more necessary primitive methods to give it more creed 
to be applicable to these fields.

4) roots as primitive (machine level) methods.

Implementing the **bbm** algorithm as a primitive method can be at least 1-3 orders
faster than the high level implementations shown in **irooot2** and **irootn**.

The following is why, with computing squareroot, but conceptually applies to all roots too.

For a cpu with n even bits (32, 64) an integer squareroot will only be composed of at most
n/2 bits. So a 64 bit cpu can retain in one register the squareroot of a 128 bit number,
i.e. a value of 2**128 - 1 = 340,282,366,920,938,463,463,374,607,431,768,211,455. For an
arbitrary sized integer it will take take (log2(num).ceil)/64 +1 cpu mem locations to 
represent the number.

So at the primitive level, it would first determine how many register length chunks the
maximum root value would need. If it needed only one register (a num < 2**128), then all
the **bbm** process is super simple and fast. If the root value takes more than one register,
the primitive algorithm would start with the MSchunk, do the bbm, and switch in the lower
chunks until finished. For roots > 2, this becomes even faster because any nth_root n requires at
most (num.bit_length/n + 1) bit to represent it. So the cube (3) root can be fit in 64 bits for
any integer < 2**(64*3) = 2**192 = 62,77,101,735,386,680,763,835,789,423,207,666,416,102,355,444,464,034,512,896

No other algorithm has the flexibility and efficiency of the **bbm** because it requires no
additions, subtractions, or divisions, only a multiplication via an exponentiation.  And that
can be done efficiently for large numbers with an appropriate algorithm (better than brute force).
All the other math (besides the > and == comparisons) involve in place bit operations -- ``>>, <<, ^, |``.
which at the cpu are one cycle (or better) operations.

I hope this has sufficiently helped to address these issues/questions.


----------------------------------------
Feature #13219: bug in Math.sqrt(n).to_i, to compute integer squareroot,  new word to accurately fix it
https://bugs.ruby-lang.org/issues/13219#change-63115

* Author: Jabari Zakiya
* Status: Open
* Priority: Normal
* Assignee: 
* Target version: 
----------------------------------------
In doing a math application using **Math.sqrt(n).to_i** to compute integer squareroots 
of integers I started noticing errors for numbers > 10**28.


I coded an algorithm that accurately computes the integer squareroot for arbirary sized numbers
but its significantly slower than the math library written in C.

Thus, I request a new method **Math.intsqrt(n)** be created, that is coded in C and part of the
core library, that will compute the integer squareroots of integers accurately and fast.

Here is working highlevel code to accurately compute the integer squareroot.

```
def intsqrt(n)
  bits_shift = (n.to_s(2).size)/2 + 1
  bitn_mask = root = 1 << bits_shift
  while true
    root ^= bitn_mask if (root * root) > n
    bitn_mask >>= 1
    return root if bitn_mask == 0
    root |= bitn_mask
  end
end

def intsqrt1(n)
  return n if n | 1 == 1   # if n is 0 or 1
  bits_shift = (Math.log2(n).ceil)/2 + 1
  bitn_mask = root = 1 << bits_shift
  while true
    root ^= bitn_mask if (root * root) > n
    bitn_mask >>= 1
    return root if bitn_mask == 0
    root |= bitn_mask
  end
end

require "benchmark/ips"

Benchmark.ips do |x|
  n = 10**40
  puts "integer squareroot tests for n = #{n}"
  x.report("intsqrt(n)"       ) { intsqrt(n)  }
  x.report("intsqrt1(n)"      ) { intsqrt1(n) }
  x.report("Math.sqrt(n).to_i") { Math.sqrt(n).to_i }
  x.compare!
end
```
Here's why it needs to be done in C.

```
2.4.0 :178 > load 'intsqrttest.rb'
integer squareroot tests for n = 10000000000000000000000000000000000000000
Warming up --------------------------------------
          intsqrt(n)     5.318k i/100ms
         intsqrt1(n)     5.445k i/100ms
   Math.sqrt(n).to_i   268.281k i/100ms
Calculating -------------------------------------
          intsqrt(n)     54.219k (賊 5.5%) i/s -    271.218k in   5.017552s
         intsqrt1(n)     55.872k (賊 5.4%) i/s -    283.140k in   5.082953s
   Math.sqrt(n).to_i      5.278M (賊 6.1%) i/s -     26.560M in   5.050707s

Comparison:
   Math.sqrt(n).to_i:  5278477.8 i/s
         intsqrt1(n):    55871.7 i/s - 94.47x  slower
          intsqrt(n):    54219.4 i/s - 97.35x  slower

 => true 
2.4.0 :179 > 

```
Here are examples of math errors using **Math.sqrt(n).to_i** run on Ruby-2.4.0.

```
2.4.0 :101 > n = 10**27; puts n, (a = intsqrt(n)), a*a, (b = intsqrt1(n)), b*b, (c = Math.sqrt(n).to_i), c*c   
1000000000000000000000000000
31622776601683
999999999999949826038432489
31622776601683
999999999999949826038432489
31622776601683
999999999999949826038432489
 => nil 
2.4.0 :102 > n = 10**28; puts n, (a = intsqrt(n)), a*a, (b = intsqrt1(n)), b*b, (c = Math.sqrt(n).to_i), c*c 
10000000000000000000000000000
100000000000000
10000000000000000000000000000
100000000000000
10000000000000000000000000000
100000000000000
10000000000000000000000000000
 => nil 
2.4.0 :103 > n = 10**29; puts n, (a = intsqrt(n)), a*a, (b = intsqrt1(n)), b*b, (c = Math.sqrt(n).to_i), c*c 
100000000000000000000000000000
316227766016837
99999999999999409792567484569
316227766016837
99999999999999409792567484569
316227766016837
99999999999999409792567484569
 => nil 
2.4.0 :104 > n = 10**30; puts n, (a = intsqrt(n)), a*a, (b = intsqrt1(n)), b*b, (c = Math.sqrt(n).to_i), c*c  
1000000000000000000000000000000
1000000000000000
1000000000000000000000000000000
1000000000000000
1000000000000000000000000000000
1000000000000000
1000000000000000000000000000000
 => nil 
2.4.0 :105 > n = 10**31; puts n, (a = intsqrt(n)), a*a, (b = intsqrt1(n)), b*b, (c = Math.sqrt(n).to_i), c*c 
10000000000000000000000000000000
3162277660168379
9999999999999997900254631487641
3162277660168379
9999999999999997900254631487641
3162277660168379
9999999999999997900254631487641
 => nil 
2.4.0 :106 > n = 10**32; puts n, (a = intsqrt(n)), a*a, (b = intsqrt1(n)), b*b, (c = Math.sqrt(n).to_i), c*c 
100000000000000000000000000000000
10000000000000000
100000000000000000000000000000000
10000000000000000
100000000000000000000000000000000
10000000000000000
100000000000000000000000000000000
 => nil 
2.4.0 :107 > n = 10**33; puts n, (a = intsqrt(n)), a*a, (b = intsqrt1(n)), b*b, (c = Math.sqrt(n).to_i), c*c
1000000000000000000000000000000000
31622776601683793
999999999999999979762122758866849
31622776601683793
999999999999999979762122758866849
31622776601683792
999999999999999916516569555499264
 => nil 
2.4.0 :108 > n = 10**34; puts n, (a = intsqrt(n)), a*a, (b = intsqrt1(n)), b*b, (c = Math.sqrt(n).to_i), c*c
10000000000000000000000000000000000
100000000000000000
10000000000000000000000000000000000
100000000000000000
10000000000000000000000000000000000
100000000000000000
10000000000000000000000000000000000
 => nil 
2.4.0 :109 > n = 10**35; puts n, (a = intsqrt(n)), a*a, (b = intsqrt1(n)), b*b, (c = Math.sqrt(n).to_i), c*c
100000000000000000000000000000000000
316227766016837933
99999999999999999873578871987712489
316227766016837933
99999999999999999873578871987712489
316227766016837952
100000000000000011890233980627554304
 => nil 
2.4.0 :110 > n = 10**36; puts n, (a = intsqrt(n)), a*a, (b = intsqrt1(n)), b*b, (c = Math.sqrt(n).to_i), c*c
1000000000000000000000000000000000000
1000000000000000000
1000000000000000000000000000000000000
1000000000000000000
1000000000000000000000000000000000000
1000000000000000000
1000000000000000000000000000000000000
 => nil 
2.4.0 :111 > n = 10**37; puts n, (a = intsqrt(n)), a*a, (b = intsqrt1(n)), b*b, (c = Math.sqrt(n).to_i), c*c
10000000000000000000000000000000000000
3162277660168379331
9999999999999999993682442519108007561
3162277660168379331
9999999999999999993682442519108007561
3162277660168379392
10000000000000000379480317059650289664
 => nil 
2.4.0 :112 > 
```



-- 
https://bugs.ruby-lang.org/

Unsubscribe: <mailto:ruby-core-request@ruby-lang.org?subject=unsubscribe>
<http://lists.ruby-lang.org/cgi-bin/mailman/options/ruby-core>