Issue #13219 has been updated by Nathan Zook.


~~~
def isqrt(n)
    r = 0
    o = 1 << ((n.bit_length >> 1 ) << 1)
    while o != 0 do
        t = r + o
        if n >= t
            n -= t
            r = (r >> 1) + o
        else
            r >>= 1
        end
        o >>= 2
    end
    r
end

def bench(n)
  Benchmark.ips do |x|
    puts "integer squareroot tests for n = #{n}"
    x.report("iroot2"     ) { n.iroot2    }
    x.report("isqrt"       ) { isqrt(n) }
    x.compare!
  end
  nil
end

bench(41)
integer squareroot tests for n = 10 ** 41
Warming up --------------------------------------
              iroot2     3.264k i/100ms
               isqrt     3.088k i/100ms
Calculating -------------------------------------
              iroot2     33.749k ( 0.6%) i/s -    169.728k in   5.029342s
               isqrt     31.469k ( 1.4%) i/s -    157.488k in   5.005489s

Comparison:
              iroot2:    33748.7 i/s
               isqrt:    31469.0 i/s - 1.07x  slower

=> nil
irb(main):245:0> bench(500)
integer squareroot tests for n = 10 ** 500
Warming up --------------------------------------
              iroot2   100.000  i/100ms
               isqrt    97.000  i/100ms
Calculating -------------------------------------
              iroot2      1.002k ( 0.7%) i/s -      5.100k in   5.092307s
               isqrt    976.542  ( 0.8%) i/s -      4.947k in   5.066141s

Comparison:
              iroot2:     1001.6 i/s
               isqrt:      976.5 i/s - 1.03x  slower

=> nil
irb(main):246:0> bench(1000)
integer squareroot tests for n = 10 ** 1000
Warming up --------------------------------------
              iroot2    31.000  i/100ms
               isqrt    38.000  i/100ms
Calculating -------------------------------------
              iroot2    322.487  ( 1.6%) i/s -      1.643k in   5.095980s
               isqrt    392.904  ( 1.3%) i/s -      1.976k in   5.029970s

Comparison:
               isqrt:      392.9 i/s
              iroot2:      322.5 i/s - 1.22x  slower

=> nil
irb(main):247:0> bench(2000)
integer squareroot tests for n = 10 ** 2000
Warming up --------------------------------------
              iroot2     7.000  i/100ms
               isqrt    15.000  i/100ms
Calculating -------------------------------------
              iroot2     74.351  ( 1.3%) i/s -    378.000  in   5.084339s
               isqrt    154.343  ( 1.3%) i/s -    780.000  in   5.054739s

Comparison:
               isqrt:      154.3 i/s
              iroot2:       74.4 i/s - 2.08x  slower

=> nil
irb(main):248:0> bench(4000)
integer squareroot tests for n = 10 ** 4000
Warming up --------------------------------------
              iroot2     1.000  i/100ms
               isqrt     5.000  i/100ms
Calculating -------------------------------------
              iroot2     15.017  ( 0.0%) i/s -     76.000  in   5.061323s
               isqrt     56.217  ( 1.8%) i/s -    285.000  in   5.070483s

Comparison:
               isqrt:       56.2 i/s
              iroot2:       15.0 i/s - 3.74x  slower

=> nil

~~~

It is interesting to note that isqrt is actually slower in the low end--this is because the ruby interpreter is doing more per loop. This in spite of the fact that this method does not require multiplies!

But back to my earlier question: if we want folks to use Ruby for multi-precision work, why not encourage them to link to GMP? (And use GMP's sqrt?) We've supported this since 2.1. I know that the square root there will be much, much faster. If we want to implement a non-GMP solution (because of the license issue), then we have to decide how much work to put into it.

I've read Zimmerman's paper (he appears to be the author of the GMP code), https://hal.inria.fr/inria-00072854/document, and it does not look too hard to implement, assuming that we have the primitives already in place. Implementing the primitives would not be overly difficult, but again, we need to start asking just how far to go. Any serious MP work must drop down to asm, because of things like the carry bit and full-word multiplication.

I'm sorry to have confused you about Newton's method. The fast way to compute square roots using Newton's method is to compute 1/sqrt(n), then multiply by n. (This can be made to always give the correct answer.) This is called "Newton-Raphson" in the literature. It has been amazingly hard to find an mp integer implementation anywhere, and I'm good enough to know better than to trust something that I just put together.


----------------------------------------
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-63122

* 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>