## Feature #13219

closed### bug in Math.sqrt(n).to_i, to compute integer squareroot, new word to accurately fix it

**Description**

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 >

#### Updated by jzakiya (Jabari Zakiya) almost 5 years ago

On my 64-bit Linux OS laptop with an I7 cpu, I tested that ** Math.sqrt(n).to_i** gives correct

answers from

`(0..9_999_899_999_899_999_322_536_673_279)`

. This was also the same on JRuby-9.1.7.0.Thus, I can create a hybrid method to maximize speed and accuracy like this.

```
def sqrt_i(n)
if n <= 9_999_899_999_899_999_322_536_673_279
Math.sqrt(n).to_i
else
intsqrt1(n)
end
end
def intsqrt1(n)
return n if n | 1 == 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
```

Test results below.

2.4.0 :415 > n = 9999899999899999322536673279; puts n, (a = intsqrt(n)), a*a, (b = intsqrt1(n)), b*b, (c = sqrt_i(n)), c*c, (d = Math.sqrt(n).to_i), d*d 9999899999899999322536673279 99999499998249 9999899999899801751003066001 99999499998249 9999899999899801751003066001 99999499998249 9999899999899801751003066001 99999499998249 9999899999899801751003066001 => nil 2.4.0 :416 > n = 9999899999899999322536673279 - 1; puts n, (a = intsqrt(n)), a*a, (b = intsqrt1(n)), b*b, (c = sqrt_i(n)), c*c, (d = Math.sqrt(n).to_i), d*d 9999899999899999322536673278 99999499998249 9999899999899801751003066001 99999499998249 9999899999899801751003066001 99999499998249 9999899999899801751003066001 99999499998249 9999899999899801751003066001 => nil 2.4.0 :417 > n = 9999899999899999322536673279 + 1; puts n, (a = intsqrt(n)), a*a, (b = intsqrt1(n)), b*b, (c = sqrt_i(n)), c*c, (d = Math.sqrt(n).to_i), d*d 9999899999899999322536673280 99999499998249 9999899999899801751003066001 99999499998249 9999899999899801751003066001 99999499998249 9999899999899801751003066001 99999499998250 9999899999900001750003062500 => nil

#### Updated by jzakiya (Jabari Zakiya) almost 5 years ago

Since this method is designed to operate on integers it would be nice to make it

a ** class Integer** method, then you can do this:

`4829391.sqrt_i`

(or whatever name it's given), which allows you to easily do method chaining.

#### Updated by jzakiya (Jabari Zakiya) almost 5 years ago

This is the C version of the algorithm.

I put the part before the start of the **while** loop of the ruby

version to make it applicable for arbitrary size integers.

http://www.codecodex.com/wiki/Calculate_an_integer_square_root

```
unsigned int sqrt32(unsigned long n)
{
unsigned int c = 0x8000;
unsigned int g = 0x8000;
for (;;) {
if (g*g > n) {
g ^= c;
}
c >>= 1;
if (c == 0) {
return g;
}
g |= c;
}
}
```

#### Updated by jzakiya (Jabari Zakiya) almost 5 years ago

I just realized ** class Integer** has a

**instance_method.**

`bit_length`

It makes the code a bit simpler, easier to understand, and faster.

```
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
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
def intsqrt2(n)
bits_shift = (n.bit_length)/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("intsqrt2(n)" ) { intsqrt2(n) }
x.report("Math.sqrt(n).to_i") { Math.sqrt(n).to_i }
x.compare!
end
```

It's a bit faster.

2.4.0 :299 > load 'intsqrttest.rb' integer squareroot tests for n = 10000000000000000000000000000000000000000 Warming up -------------------------------------- intsqrt(n) 5.495k i/100ms intsqrt1(n) 5.679k i/100ms intsqrt2(n) 5.759k i/100ms Math.sqrt(n).to_i 261.831k i/100ms Calculating ------------------------------------- intsqrt(n) 56.563k (± 5.9%) i/s - 285.740k in 5.069484s intsqrt1(n) 58.567k (± 6.1%) i/s - 295.308k in 5.062216s intsqrt2(n) 59.769k (± 6.1%) i/s - 299.468k in 5.029980s Math.sqrt(n).to_i 5.377M (± 6.7%) i/s - 26.969M in 5.039343s Comparison: Math.sqrt(n).to_i: 5377357.4 i/s intsqrt2(n): 59768.9 i/s - 89.97x slower intsqrt1(n): 58566.6 i/s - 91.82x slower intsqrt(n): 56563.4 i/s - 95.07x slower => true

#### Updated by jzakiya (Jabari Zakiya) almost 5 years ago

This works well in my code.

Easier to code and read with accurate results.

```
class Integer
def sqrt_i
self <= MAX_RANGE ? Math.sqrt(self).to_i : intsqrt(self)
end
private
MAX_RANGE = 9_999_899_999_899_999_322_536_673_279
def intsqrt(n)
bits_shift = (n.bit_length)/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
end
```

2.4.0 :020 > n = 10**111; puts n, (a = n.sqrt_i), a*a, (b = Math.sqrt(n).to_i), b*b 1000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000 31622776601683793319988935444327185337195551393252168268 999999999999999999999999999999999999999999999999999999963630737732541472910196882732895436793219183483386119824 31622776601683794595141492119969113562259413928748515328 1000000000000000080647728865639515055517417465817645783945119006952630720180836471433888096459369651964250947584

#### Updated by jzakiya (Jabari Zakiya) almost 5 years ago

A new version ** intsqrt3** is one line shorter, and more intuitive (to me)

but there is a neglible difference in speed versus

**in Ruby.**

`intsqrt2`

Maybe it's faster in C?

```
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
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
def intsqrt2(n)
bits_shift = (n.bit_length)/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 intsqrt3(n)
bits_shift = (n.bit_length)/2 + 2
root, bitn_mask = 0, (1 << bits_shift)
until (bitn_mask >>= 1) == 0
root |= bitn_mask
root ^= bitn_mask if (root * root) > n
end
root
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("intsqrt2(n)" ) { intsqrt2(n) }
x.report("intsqrt3(n)" ) { intsqrt3(n) }
# x.report("Math.sqrt(n).to_i") { Math.sqrt(n).to_i }
x.compare!
end
```

2.4.0 :391 > load 'intsqrttest.rb' integer squareroot tests for n = 10000000000000000000000000000000000000000 Warming up -------------------------------------- intsqrt(n) 5.621k i/100ms intsqrt1(n) 5.890k i/100ms intsqrt2(n) 5.879k i/100ms intsqrt3(n) 5.872k i/100ms Calculating ------------------------------------- intsqrt(n) 58.255k (± 6.9%) i/s - 292.292k in 5.041702s intsqrt1(n) 59.962k (± 6.8%) i/s - 300.390k in 5.033585s intsqrt2(n) 60.283k (± 6.4%) i/s - 305.708k in 5.092434s intsqrt3(n) 60.120k (± 6.7%) i/s - 305.344k in 5.102386s Comparison: intsqrt2(n): 60283.1 i/s intsqrt3(n): 60119.7 i/s - same-ish: difference falls within error intsqrt1(n): 59962.4 i/s - same-ish: difference falls within error intsqrt(n): 58254.7 i/s - same-ish: difference falls within error => true

#### Updated by shevegen (Robert A. Heiler) almost 5 years ago

I can not evaluate whether the above is correct or not (not saying that it is not,

I just simply do not know), but I believe that method names such as "`intsqrt2`

" or

"`intsqrt`

" may possibly not be an ideal choice, due to the name alone.

`Math.sqrt()`

is already a bit shortish, adding more abbreviations, in particular

something such as "int" as part of it, would be even more peculiar.

If it would be a bug, as behaviour, why should `Math.sqrt()`

itself not be able

to deal with it? Internally it could call any other method before it outputs

a result.

#### Updated by nobu (Nobuyoshi Nakada) almost 5 years ago

Robert A. Heiler wrote:

If it would be a bug, as behaviour, why should

`Math.sqrt()`

itself not be able

to deal with it? Internally it could call any other method before it outputs

a result.

Expected results differ.

`Math.sqrt()`

is expected to return a `Float`

with much precision as possible.

While the proposed method will return an `Integer`

which does not exceed the real square root.

#### Updated by jzakiya (Jabari Zakiya) almost 5 years ago

To be clear, I am not saying ** Math.sqrt** has a bug in it.

As stated, it produces a floating point result, which inherently has finite (not infinite) precision.

The "bug", or more accurately, the realization of the limitations of using ** Math.sqrt(n).to_i** to

produce correct results for integers past a certain size, led me to suggest creating a new method

which will produce correct results for arbitrary sized integers (i.e. the largest integer value for

`root`

such that `root*root <= n`

).This need emanated for me in numerical applications pertaining to generating very large prime numbers.

As I previously suggested, this method should most naturally be a ** class Integer** method, and after much

consideration I like the method name

**, versus**

`sqrt_i`

`intsqrt`

, `isqrt`

, etc. This makes its usage as -

`n.sqrt_i`

super simple, and idiomatic, and would give an error message if applied to a non-

`Integer`

number/object, to not confuse its use with

**, which inherently produces a**

`Math.sqrt(n)`

**float**value.

I think this name is more ruby idiomatic, because the name of the function is first, and the **x_i**

extension is consistent with other naming conventions to denote class Type, like ** to_i**,

**,**

`to_f`

**.**

`to_s`

The algorithm I provided (both C and Ruby) show the method is very simple to code, and would be very fast

if coded in C, using its native shift operators, etc. If I had any significant C chops I would try to do it myself. :-)

Also, providing a fast native implementation of this function makes Ruby more attractive for numerical applications,

especially for Number Theory and Cryptography, and also enhances the Ruby 3x3 goals of making it at least 2x faster

than Ruby 2.0.

#### Updated by jzakiya (Jabari Zakiya) almost 5 years ago

After further testing I found the same errors when using ** Math.sqrt(n).to_i** with

large number when using

`(n**(1.0/2)).to_i`

. This reinforces to me the need to provide **(by whatever name) so users won't fall prey to this undocumented anomaly.**

`sqrt_i`

2.4.0 :129 > n= 10**35; puts n, a = intsqrt4(n), a*a, b = Math.sqrt(n).to_i, b*b, c = (n**(1.0/2)).to_i, c*c 100000000000000000000000000000000000 316227766016837933 99999999999999999873578871987712489 316227766016837952 100000000000000011890233980627554304 316227766016837952 100000000000000011890233980627554304 => nil 2.4.0 :130 >

I also realized that the algorithm used to determine the integer squareroot can be easily

generalized to correctly generate integer nth roots for arbitrary size integers.

This capability will now allows Ruby to be used for many Number Theory and Cryptographic problems,

such as generating correct integer values for eliptic curves, real integer roots of polynomials,

public key, et al, encryption, for very large number sets.

Because of this, I am going to add this capability to my **roots** rubygem.

https://rubygems.org/gems/roots

It has two (2) methods ** root** and

**which generate all the n roots of an nth root for**

`roots`

all

**types (integers, floats, complex, rational).**

`Numeric`

In order to not clash with my proposed ** class Integer** method name of

**as a Ruby core method,**

`sqrt_i`

I will name two more

**methods**

`class Integer`

**and**

`iroot2`

**to be added to**

`irootn(n)`

**.**

`roots`

Since finding squareroots is much more frequently done than other roots I give it its own method name.

Below is the code for doing this (which I may refine).

These methods will return the correct nth integer real (not complex) root for positive integers, and

for odd roots of negative integers.

```
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
sign = self < 0 ? -1 : 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*sign
end
def iroot2; irootn(2) end
end
```

As with generating integer squareroots, doing `(num**(1.0/n)).to_i`

produces incorrect answers for

integers past some maximun value, as shown below.

2.4.0 :129 > n = 10**37; puts n, a = n.irootn(3), a**3, b = (n**(1.0/3)).to_i, b**3 10000000000000000000000000000000000000 2154434690031 9999999999987694380850132072511299791 2154434690031 9999999999987694380850132072511299791 => nil 2.4.0 :130 > n = 10**38; puts n, a = n.irootn(3), a**3, b = (n**(1.0/3)).to_i, b**3 100000000000000000000000000000000000000 4641588833612 99999999999949657815157877549034676928 4641588833612 99999999999949657815157877549034676928 => nil 2.4.0 :131 > n = 10**39; puts n, a = n.irootn(3), a**3, b = (n**(1.0/3)).to_i, b**3 1000000000000000000000000000000000000000 10000000000000 1000000000000000000000000000000000000000 9999999999999 999999999999700000000000029999999999999 => nil 2.4.0 :134 > n = 10**112; puts n, a = n.irootn(4), a**4, b = (n**(1.0/4)).to_i, b**4 10000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000 10000000000000000000000000000 10000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000 9999999999999999583119736832 9999999999999998332478947328000104273492291412559539763672807300805169073432752214389506884609933472640769458176 => nil 2.4.0 :135 >

These methods are also preferable because they work correctly for negative integers, while the exponential root

method throws errors.

2.4.0 :150 > n = (10**112); puts n, a = n.irootn(3), a**3 10000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000 21544346900318837217592935665193504952 9999999999999999999999999999999999999173635536969138804543654476122863556031297471619838179917400017798906449408 => nil 2.4.0 :151 > n = (10**112); puts n, b = (n**(1.0/3)).to_i, b**3 10000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000 21544346900318734919579678222552399872 9999999999999857552405189045342271937486394088257690663380520853991701500619736949767614488813401623055562702848 => nil 2.4.0 :152 > n = -(10**112); puts n, a = n.irootn(3), a**3 -10000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000 -21544346900318837217592935665193504952 -9999999999999999999999999999999999999173635536969138804543654476122863556031297471619838179917400017798906449408 => nil 2.4.0 :153 > n = -(10**112); puts n, b = (n**(1.0/3)).to_i, b**3 RangeError: can't convert 1.077217345015937e+37+1.865795172362055e+37i into Integer from (irb):153:in `to_i' from (irb):153 from /home/jzakiya/.rvm/rubies/ruby-2.4.0/bin/irb:11:in `<main>' 2.4.0 :154 >

Again, creating C coded core versions of these methods not only will be more efficient but faster.

#### Updated by stomar (Marcus Stollsteimer) almost 5 years ago

To clarify further the limitations of using Float here:

```
n = 10**35
Math.sqrt(n).to_i # => 316227766016837952
Math.sqrt(n).next_float.to_i # => 316227766016838016
```

BigDecimal / BigMath should be able to produce the correct results, but possibly slower(?) than a specialized integer method:

```
require "bigdecimal"
n = 10**35
r = BigDecimal(n).sqrt(20).to_i # => 316227766016837933
r ** 2 # => 99999999999999999873578871987712489
(r+1) ** 2 # => 100000000000000000506034404021388356
```

__Addendum:__

These methods are also preferable because they work correctly for negative integers, while the exponential root method throws errors.

To be precise: "correct" in __your__ sense (i.e. returning a real number); the "exponential root" itself __does__ work correctly -- returning a complex root is just fine --, the exception is raised by trying to convert a complex number into an integer. You could easily avoid this specific problem by calculating -(|n|^{1/3}).

#### Updated by jzakiya (Jabari Zakiya) almost 5 years ago

One if the really nice things about Ruby is its intent to make users/programmers happy,

and to adhere to the "principle of least surprises (pols)".

I don't think it was the intent of Matz to create math methods that give correct results

some of the time, but not tell users under what conditions they wouldn't give correct results.

Unless you had my use cases, it was probably never tested to cover giving correct results for them.

I have identified verified incorrect results that I have seen no documentation warning users of.

I have presented a verified method to always produce correct results, not only for the initial

case for integer squareroots for arbitrary size integers, but now also for any real nth integer

root for arbitrary sized integers.

As I single user, I have solved these inherent errors in using these native math functions for

my use cases. So I now know what I have to do so I will get correct answers in my applications.

I now will also include these methods in my **roots** gem so if others need/want them they are available.

I think it would be in the best interest of Ruby users for the language to provide them the methods

I have presented to create correct results for their use cases, which should be fairly easy to code

and include in the core language, and documented accordingly for users, versus messing with the

inherent floating point methods and figuring out how to coerce them into producing correct integer

results for every possible arbitrary sized integer, or forcing users to figure out how to do so.

#### Updated by jzakiya (Jabari Zakiya) almost 5 years ago

```
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
require "bigdecimal"
require "benchmark/ips"
Benchmark.ips do |x|
n = 10**35
puts "integer squareroot tests for n = #{n}"
x.report("iroot2" ) { n.iroot2 }
x.report("irootn(2)" ) { n.irootn(2) }
x.report("BigDecimal(n).sqrt(5 ).to_i") { BigDecimal(n).sqrt(5 ).to_i }
x.report("BigDecimal(n).sqrt(10).to_i") { BigDecimal(n).sqrt(10).to_i }
x.report("BigDecimal(n).sqrt(20).to_i") { BigDecimal(n).sqrt(20).to_i }
x.compare!
end
```

Yes, its much slower, even to the highlevel Ruby versions.

2.4.0 :201 > load 'irootstest.rb' integer squareroot tests for n = 100000000000000000000000000000000000 Warming up -------------------------------------- iroot2 5.681k i/100ms irootn(2) 5.714k i/100ms BigDecimal(n).sqrt(5 ).to_i 3.021k i/100ms BigDecimal(n).sqrt(10).to_i 2.953k i/100ms BigDecimal(n).sqrt(20).to_i 2.616k i/100ms Calculating ------------------------------------- iroot2 57.825k (± 3.3%) i/s - 289.731k in 5.016021s irootn(2) 57.462k (± 3.7%) i/s - 291.414k in 5.078940s BigDecimal(n).sqrt(5 ).to_i 30.543k (± 2.8%) i/s - 154.071k in 5.048265s BigDecimal(n).sqrt(10).to_i 30.709k (± 3.1%) i/s - 153.556k in 5.005239s BigDecimal(n).sqrt(20).to_i 26.725k (± 3.0%) i/s - 136.032k in 5.094723s Comparison: iroot2: 57825.2 i/s irootn(2): 57461.9 i/s - same-ish: difference falls within error BigDecimal(n).sqrt(10).to_i: 30708.9 i/s - 1.88x slower BigDecimal(n).sqrt(5 ).to_i: 30543.4 i/s - 1.89x slower BigDecimal(n).sqrt(20).to_i: 26725.0 i/s - 2.16x slower => true 2.4.0 :202

And you need to know beforehand the needed correct precision to display the correct results.

2.4.0 :214 > n = 10**35; puts n, n.iroot2, BigDecimal(n).sqrt(5).to_i, BigDecimal(n).sqrt(10).to_i,BigDecimal(n).sqrt(20).to_i 100000000000000000000000000000000000 316227766016837933 316227766016837933 316227766016837933 316227766016837933 => nil 2.4.0 :215 > n = 10**45; puts n, n.iroot2, BigDecimal(n).sqrt(5).to_i, BigDecimal(n).sqrt(10).to_i,BigDecimal(n).sqrt(20).to_i 1000000000000000000000000000000000000000000000 31622776601683793319988 31622776601683666666666 31622776601683666666666 31622776601683793319988 => nil 2.4.0 :216 > n = 10**55; puts n, n.iroot2, BigDecimal(n).sqrt(5).to_i, BigDecimal(n).sqrt(10).to_i,BigDecimal(n).sqrt(20).to_i 10000000000000000000000000000000000000000000000000000000 3162277660168379331998893544 3162277660168379331499021527 3162277660168379331499021527 3162277660168379331998893544 => nil 2.4.0 :217 > n = 10**65; puts n, n.iroot2, BigDecimal(n).sqrt(5).to_i, BigDecimal(n).sqrt(10).to_i,BigDecimal(n).sqrt(20).to_i 100000000000000000000000000000000000000000000000000000000000000000 316227766016837933199889354443271 316227766016837466536394723986322 316227766016837466536394723986322 316227766016837933199889000000000 => nil 2.4.0 :218 >

#### Updated by stomar (Marcus Stollsteimer) almost 5 years ago

@Jabari Note that my comments were not meant as arguments against your feature request, I only wanted to point out the reasons for the observed behavior. Personally, I always welcome improved math functionality; however, I'm not informed enough to know whether this should be in the stdlib or not.

There certainly is a problem with the docs for Math.sqrt, which is incomplete (return value types are not given), misleading (it promises to avoid rounding errors), and actually doesn't match the actual behavior:

$ ri Math::sqrt (from ruby site) ------------------------------------------------------------------------------ sqrt(a) ------------------------------------------------------------------------------ Computes the square root of a. It makes use of Complex and Rational to have no rounding errors if possible. Math.sqrt(4/9) # => 2/3 Math.sqrt(- 4/9) # => Complex(0, 2/3) Math.sqrt(4.0/9.0) # => 0.666666666666667

Based on this you would for instance expect the square root of an integer square number to be returned as an (exact) integer or rational; but actually it's a float:

```
Math.sqrt(9) # => 3.0 ("no rounding errors if possible" ?!?)
```

#### Updated by jzakiya (Jabari Zakiya) almost 5 years ago

No, no, I didn't take it as a negative comment per se,

I just hadn't tried to use **BigDecimal** before, so I just went ahead and

created tests to see for myself how they performed (accuracy and speed).

Ruby has a much nicer community compared to some others, which is why

I've tried to document the problem, and a solution for it, in great detail.

I want to see Ruby become more used in mathematical and numerical applications

like Python is now becoming. Thus, it needs to be as accurate as possible, and

then, of course, fast (enough).

For users (new and old) using Ruby for advanced mathematical/numerical applications,

and the new rave of "big data" things, I don't want it developing a reputation of

producing erroneous results, especially compared to Python, et al, which maybe don't

produce these errors (and I have no idea how inherently inaccurate Python is).

Here is a section of the README.md file for my new version of my **roots** gem I'll be releasing.

## Description Starting with Roots 2.0.0 (2017-2-20) the methods **iroot2** and **irootn** were added. They are instance_methods of **class Integer** that will accurately compute the real value squareroot and nth_root for arbitrary sized Integers. If you have an application where you need the exact correct integer real value for roots, especially for arbitrary sized (large) integers, use these methods instead or **root|s**. These methods work strictly in the Integer domain, and do not first create floating point approximations of these roots and then convert them to integers. They are useful for applications in pure math, Number Theory, Cryptography, etc, where you need to find the exact real integer roots of polynomials, encryption functions, etc. The methods **root** and **roots** work for all the **Numeric** classes (integers, floats, complex, rationals), and produce floating point results. They, thus, produce floating point approximations to the exact **Integer** values for real roots for arbitrary sized integers. Usage **iroot2** Use syntax: ival.iroot2 Return the largest Integer +root+ of Integer ival such that root**2 <= ival A negative ival will result in 'nil' being returned. 9.iroot2 => 3 -9.iroot2 => nil 120.iroot2 => 10 121.iroot2 => 11 121.iroot2.class => Integer (10**10).iroot2 => 100000 (10**10).root(2) => 100000.0 (10**10).roots(2)[0] = 100000.0+0.0i (10**11).iroot2 => 316227 (10**11).root(2) => 316227.76601684 (10**11).roots(2)[0] = 316227.76601684+0.0i (10**110).iroot2 => 10000000000000000000000000000000000000000000000000000000 (10**110).root(2) => 1.0e+55 (10**110).roots(2)[0] = 1.0e+55+0.0i (10**111).iroot2 => 31622776601683793319988935444327185337195551393252168268 (10**111).root(2) => 3.1622776601683795e+55 (10**111).roots(2)[0] = 3.1622776601683795e+55+0.0i

#### Updated by jzakiya (Jabari Zakiya) almost 5 years ago

Here's the README.md writeup on the general integer roots algorithm.

At the cpu level, all these bit operations -- >>, <<, |, ^ -- are one clock cycle instructions,

so you see this can be very fast if done in C. Maybe for arbitrary size numbers that can't

reside in a single cpu register they take longer, but they would still be faster than high level Ruby.

For integer numbers: For binary based (digital) computers (cpu), all integers are represented as binary numbers. To find the root "n" of an integer "num" we can use the following process to find the largest integer nth "root" such that root**n <= num. For an integer composed of 'b' bits an nth root 'n' will have at most (b/n + 1) bits. For squareroots, n = 2: For 9 = 0b1001, it has b = 4 bits, and its squareroot has at most (4/2 + 1) = 3 bits. The squareroot of 9 is 3 = 0b11, which is 2 bits. For 100 = 0b1100100, b = 7, and its squareroot has at most (7/2 + 1) = 4 bits. The squareroot of 100 is 10 = 0b1010, which is 4 bits. For cuberoots, n = 3: For 8 = 0b1000, it has b = 4 bits, and its cuberoot has at most (4/3 + 1) = 2 bits. The cuberoot of 8 is 2 = 0b10, which is 2 bits. For 125 = 0b1111101, b = 7, and its cuberoot has at most (7/3 + 1) = 3 bits. The cuberoot of 125 is 5 = 0b101, which is 3 bits. Algorithm: bits_shift = (num.bits.length)/n + 1 # determine max number of root bits bitn_mask = 1 << bits_shift # set value for max bit position of root root = 0 # initialize the value for root until bitn_mask == 0 # step through all the bit positions for root root |= bitn_mask # set current bit position to '1' in root root ^= bitn_mask if root**n > num # set it back to '0' if root too large bitn_mask >>= 1 # set bitn_mask to next smaller bit position end root # return exact integer value for root

#### Updated by akr (Akira Tanaka) almost 5 years ago

Is there many application for this method?

Concrete examples may help to persuade matz.

Another important problem is the method name.

If we use same name for same functionality in other languages/libraries,

it makes learning this method easier.

What's the name of this function in other languages/libraries?

I found isqrt in Common Lisp and Tcl, mpz_sqrt in GMP.

Others?

#### Updated by jzakiya (Jabari Zakiya) almost 5 years ago

I updated my **roots** rubygem to 2.0.0 to include **iroot2** and **irootn**.

https://rubygems.org/gems/roots

https://github.com/jzakiya/roots

#### Updated by Student (Nathan Zook) almost 5 years ago

I just noticed this thread. One bit per cycle methods are not fast by modern methods, they are quite slow. I will attempt to have something more concrete within 24 hours.

Nathan Zook

#### Updated by Student (Nathan Zook) almost 5 years ago

You might want to consider the following articles:

https://www.reddit.com/r/algorithms/comments/1zt63v/fast_algorithm_to_calculate_integer_square_root/

Which lead me to wikipedia:

https://en.wikipedia.org/wiki/Methods_of_computing_square_roots#Binary_numeral_system_.28base_2.29

(You might want to start a bit higher in the article to get the context.)

and also to Zimmerman's implementation:

https://gmplib.org/repo/gmp-5.1/file/c5010c039373/mpn/generic/sqrtrem.c

The wikipedia article implements a faster bit-by-bit computation. Zimmerman implements the standard solution, which converges quadratically.

Looking at your comments, however, you seem really to be interested in a high performance multiprecision library for ruby. Is there any particular reason that we should not take an existing C library and drop a wrapper around it?

#### Updated by Student (Nathan Zook) almost 5 years ago

I am confused by your solution in comment #3. By your own report, there is no issue for numbers much larger than a single word. But this code is clearly appropriate only for single-word inputs.

#### Updated by matz (Yukihiro Matsumoto) almost 5 years ago

It's OK to add the functionality. Is `Integer.sqrt(n)`

OK for you?

Matz.

#### Updated by nobu (Nobuyoshi Nakada) almost 5 years ago

Yukihiro Matsumoto wrote:

It's OK to add the functionality. Is

`Integer.sqrt(n)`

OK for you?

When `n`

is not an `Integer`

(but a `Numeric`

), is it truncated first then `sqrt`

?

#### Updated by tompng (tomoya ishida) almost 5 years ago

using Newton's method might be another good way to implement it.

```
def intsqrt_newton(n)
raise if n<0
return n if n<=1
r = 1<<((n.bit_length+1)/2)
# r*r >= n
while true # r*r-n >= 2*r
r2 = r-(r*r-n)/(2*r)
break if r2 == r
r=r2
end
# r*r-n < 2*r
# (r-1)**2 < n+1
# √n-1 < r-1 <= √n < √(n+1)
r-1
end
```

#### Updated by stomar (Marcus Stollsteimer) almost 5 years ago

__Edit: please ignore (misunderstanding)__

~~I think ~~`Integer.sqrt`

or `Numeric.sqrt`

is extremely problematic. I guess everybody would expect it to be an alias for `Math::sqrt`

and behave like this:

```
# (probably) expected behavior
2.sqrt(2) # => 1.4142135623730951
10.sqrt(2) # => 3.1622776601683795
(6.25).sqrt(2) # => 2.5
# proposed behavior for this feature
10.sqrt(2) # => 3
# _not_ part of this request
(6.25).sqrt(2) # 2 ??? 2.5 ???
```

~~Moreover, ~~`sqrt`

(__square__ root) implies `n = 2`

, so `sqrt(3)`

, `sqrt(4)`

etc. is contradictory.

#### Updated by stomar (Marcus Stollsteimer) almost 5 years ago

For square roots, `isqrt`

seems to be an option with some(?) spread, see https://en.wikipedia.org/wiki/Integer_square_root. It would look like this:

```
26.isqrt # => 5
30.isqrt # => 5
```

For higher roots, maybe `iroot(n)`

could be an option:

```
26.iroot(3) # => 2
30.iroot(3) # => 3
```

~~(However, I'm not sure whether the request goes beyond square roots.)~~

#### Updated by akr (Akira Tanaka) almost 5 years ago

Marcus Stollsteimer wrote:

I think

`Integer.sqrt`

or`Numeric.sqrt`

is extremely problematic. I guess everybody would expect it to be an alias for`Math::sqrt`

and behave like this:

No. Integer.sqrt is a class method.

```
Integer.sqrt(1000) #=> 31
```

It is clealy different from "Math.sqrt(n)" (or "sqrt(n)" after "include Math").

The benefit of Integer.sqrt is that it doesn't introduce a new word which expands

vocabulary of ruby.

If there is common name for this function, expanding the vocabulary is not a big problem, though.

Also, Integer.sqrt(n) is square root.

So, n-th root should be a different method:

Integer.root(m, n) for n-th root, for example.

#### Updated by matz (Yukihiro Matsumoto) almost 5 years ago

Marcus, probably you are confusing `Integer.sqrt()`

and `Integer#sqrt()`

.

The former, you call `Integer.sqrt(4) # => 2`

and the latter, you call `4.sqrt # => 2`

.

We are talking about the former.

Matz.

#### Updated by Student (Nathan Zook) almost 5 years ago

Nobuyoshi Nakada wrote:

Yukihiro Matsumoto wrote:

It's OK to add the functionality. Is

`Integer.sqrt(n)`

OK for you?When

`n`

is not an`Integer`

(but a`Numeric`

), is it truncated first then`sqrt`

?

The beauty of this particular problem is that it does not matter. ;)

#### Updated by jzakiya (Jabari Zakiya) almost 5 years ago

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.

#### Updated by stomar (Marcus Stollsteimer) almost 5 years ago

Matz, akr:

sorry for the noise, it's of course completely clear. I only didn't read close enough and missed the switch from the proposed instance method and n-th root to a class method. Stupid.

Now it's also clear to me that you were only talking about square roots, not higher roots.

#### Updated by jzakiya (Jabari Zakiya) almost 5 years ago

No, I am sorry for any confusion.

My initial concern was for accurately computing exact integer squareroots

for my immediate use case needs.

After finding the **bbm** algorithm for squareroots, and modifying it to

compute any nth-root, I realized this should be the primitive to create,

because the squareroot (n=2) is just a particular instance, which is the

most frequently used root.

#### Updated by Student (Nathan Zook) almost 5 years ago

Jabari Zakiya wrote:

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

iroot2andirootnare now in myrootsgem and usebbm.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).

First, I love the idea of giving ruby the ability to play in the MP space. Anything to allow this beautiful tool to be used more broadly must be a good thing. I very much support adding such support. (Yes, this is me volunteering.)

Having said this, to do so involves a number of issues. For instance, I just attempted to time Integer#prime? on a 4096-bit number. Got a floating point overflow. Same for 2048 bits. So, I tried 1024. It's still running. In fact, I expect my computer to have a mechanical breakdown before the routine returns true.

The articles which I linked included both a faster bit-by-bit algorithm, and mentioned Newton-Raphson-derived methods, which clearly converge with better O(). The given benchmark is of a 133-bit number. I would like to see timing numbers for 128, 512, 1024, 2048, and 4096-bit integers before suggesting that we are going there. Serious work looks a much larger integers, btw.

I'll see what I can whip up for benchmarks...

#### Updated by jzakiya (Jabari Zakiya) almost 5 years ago

Hey that's great! I'd love to help work on this too, though

I'm not much of a C programmer. I can help with testing though.

Below are test results comparing the **bbm** algorithm used in

**iroot2** and **irootn** versus the general Newton's method.

As you can see, for the range of values used, Newton becomes

much slower as the numbers gets bigger.

First **bbm**, for squareroot, is order O(log2(n)/2)), which I

think beats Newton. But also compared to Newton, all those

additions and divisions done in the inner loop become more

expensive (cpu intensive) as the numbers get bigger (taking care

of carries, etc), whereas the **bbm** merely does inplace bit operations.

You obviously can make Newton's mechanically better, but I suspect

it will never be faster on a real computer than **bbm** because of

its cpu work cost to implement. In fact my I7 cpu laptop's fan started

to kick in running these benchmarks, which I let cool down between runs.

Both methods did create the same correct results.

Test suite:

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 require "benchmark/ips" Benchmark.ips do |x| exp = 4000; n = 10**exp puts "integer squareroot tests for n = 10**#{exp}" x.report("iroot2" ) { n.iroot2 } x.report("irootn(2)" ) { n.irootn(2) } x.report("newtons" ) { intsqrt7(n) } x.compare! end

Benchmark results:

2.4.0 :042 > load 'irootstest.rb' integer squareroot tests for n = 10**500 Warming up -------------------------------------- iroot2 70.000 i/100ms irootn(2) 69.000 i/100ms newtons 43.000 i/100ms Calculating ------------------------------------- iroot2 636.901 (±14.9%) i/s - 3.150k in 5.089664s irootn(2) 676.828 (± 5.6%) i/s - 3.381k in 5.012279s newtons 427.424 (± 4.4%) i/s - 2.150k in 5.040686s Comparison: irootn(2): 676.8 i/s iroot2: 636.9 i/s - same-ish: difference falls within error newtons: 427.4 i/s - 1.58x slower => true 2.4.0 :043 > load 'irootstest.rb' integer squareroot tests for n = 10**1000 Warming up -------------------------------------- iroot2 22.000 i/100ms irootn(2) 22.000 i/100ms newtons 14.000 i/100ms Calculating ------------------------------------- iroot2 228.501 (± 3.9%) i/s - 1.144k in 5.014841s irootn(2) 227.442 (± 4.8%) i/s - 1.144k in 5.042017s newtons 139.466 (± 3.6%) i/s - 700.000 in 5.027213s Comparison: iroot2: 228.5 i/s irootn(2): 227.4 i/s - same-ish: difference falls within error newtons: 139.5 i/s - 1.64x slower => true 2.4.0 :044 > load 'irootstest.rb' integer squareroot tests for n = 10**2000 Warming up -------------------------------------- iroot2 6.000 i/100ms irootn(2) 6.000 i/100ms newtons 3.000 i/100ms Calculating ------------------------------------- iroot2 67.819 (± 2.9%) i/s - 342.000 in 5.049151s irootn(2) 68.311 (± 2.9%) i/s - 342.000 in 5.013356s newtons 38.863 (± 2.6%) i/s - 195.000 in 5.024956s Comparison: irootn(2): 68.3 i/s iroot2: 67.8 i/s - same-ish: difference falls within error newtons: 38.9 i/s - 1.76x slower => true 2.4.0 :045 > load 'irootstest.rb' integer squareroot tests for n = 10**4000 Warming up -------------------------------------- iroot2 1.000 i/100ms irootn(2) 1.000 i/100ms newtons 1.000 i/100ms Calculating ------------------------------------- iroot2 17.380 (± 5.8%) i/s - 87.000 in 5.011300s irootn(2) 17.417 (± 5.7%) i/s - 87.000 in 5.000415s newtons 9.512 (± 0.0%) i/s - 48.000 in 5.050706s Comparison: irootn(2): 17.4 i/s iroot2: 17.4 i/s - same-ish: difference falls within error newtons: 9.5 i/s - 1.83x slower => true 2.4.0 :046 >

#### Updated by Student (Nathan Zook) almost 5 years ago

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.

#### Updated by jzakiya (Jabari Zakiya) almost 5 years ago

Ah, these are good and interesting results Nathan. Thanks for posting.

It would be interesting to see how both faired as primitive implementations.

I will definitely study the code and play with it.

However **bbm** may still be the smallest/efficient/fastest way to exactly do

roots > 2, and can be used as the benchmark algorithm to test others against.

It certainly eliminates having to create specific optimal implementations for

individual roots > 2.

However, here's another reason this issue is important to the larger Ruby 3x3 goal.

In the library file **prime.rb** is the method **Integer#prime?**.

class Integer ... ... def prime? return self >= 2 if self <= 3 return true if self == 5 return false unless 30.gcd(self) == 1 (7..Math.sqrt(self).to_i).step(30) do |p| return false if self%(p) == 0 || self%(p+4) == 0 || self%(p+6) == 0 || self%(p+10) == 0 || self%(p+12) == 0 || self%(p+16) == 0 || self%(p+22) == 0 || self%(p+24) == 0 end true end ... ...

As in typical prime algorithms, you just check for the primes <= the number's squareroot.

Here **Math.sqrt(self).to_i** determines the integer squareroot, but as we know (and as I

found out, which initiated all this), as the numbers gets `> ~10**28`

, the approximation for

the squareroot becomes more and more larger than the actual value. This will cause **prime?**

to do increasingly more work than necessary as the numbers get bigger.

(Completely as an aside, for numbers this size you should probably look to use: n.to_bn.prime?)

In **prime.rb** I counted three places where **Math.sqrt** is used to find the integer squareroot.

Creating a primitive integer squareroot method would not only make these results exact, but smaller

as n gets bigger, increasing the accuracy and speed of any application which uses them.

It would be interesting to see what an audit of the entire Ruby code base would reveal where

this is used to provide the integer squareroot, versus the floating point one.

But also to Nathan's point, until I started looking for better alternatives to work with larger

numbers, I was totally unaware of the resources in the **openssl** library that can be applied

to arbitrary sized integers. Since these resource already exist within the Ruby ecosystem, maybe

a first step is to assess and test how well these existing resources meet these needs and

then standardize using them ubiquitously within the core library to optimize their performance.

From a users perspective, I first want accurate results, then speed. A fast (or slow) incorrect

result is of no value in doing mathematical or numerical analysis heavy applications.

#### Updated by tompng (tomoya ishida) almost 5 years ago

Newton seems to be faster than bbm, if the initial x is closer to √n.

when I use x=1<<((n.bit_length+1)/2) for the initial x in newton's method,

the iteration count (n=10**4000) became 11 (was 6656 when initial x is n).

I think the calculation cost is:

newton: O(log(k) * k**1.585)) Nathan Zook's isqrt: O(k**2) bbm: O(k * k**1.585)) k=log(n) bignum mult: O(k**1.585) # Bignum uses karatsuba algorithm

benchmark

```
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
def intsqrt7_with_good_initial_value(n) # Newton's method with √n < initial_x <= 2*√n
raise if n<0
return n if n<=1
x = 1<<((n.bit_length+1)/2)
y = (x + n/x)/2
while y < x
x = y
y = (x + n/x)/2
end
x
end
require "benchmark/ips"
raise unless 100000.times.all?{|i|intsqrt7_with_good_initial_value(i) == i.iroot2}
Benchmark.ips do |x|
exp = 4000; n = 10**exp
puts "integer squareroot tests for n = 10**#{exp}"
x.report("iroot2" ) { n.iroot2 }
x.report("newtons" ) { intsqrt7(n) }
x.report("newtons_fast" ) { intsqrt7_with_good_initial_value(n) }
x.compare!
end
```

result

integer squareroot tests for n = 10**4000 Warming up -------------------------------------- iroot2 1.000 i/100ms newtons 1.000 i/100ms newtons_fast 90.000 i/100ms Calculating ------------------------------------- iroot2 13.054 (±15.3%) i/s - 64.000 in 5.063850s newtons 2.455 (± 0.0%) i/s - 13.000 in 5.337181s newtons_fast 939.494 (± 7.3%) i/s - 4.680k in 5.022312s Comparison: newtons_fast: 939.5 i/s iroot2: 13.1 i/s - 71.97x slower newtons: 2.5 i/s - 382.74x slower

#### Updated by akr (Akira Tanaka) almost 5 years ago

Jabari Zakiya wrote:

In the library file

prime.rbis the methodInteger#prime?.class Integer ... ... def prime? return self >= 2 if self <= 3 return true if self == 5 return false unless 30.gcd(self) == 1 (7..Math.sqrt(self).to_i).step(30) do |p| return false if self%(p) == 0 || self%(p+4) == 0 || self%(p+6) == 0 || self%(p+10) == 0 || self%(p+12) == 0 || self%(p+16) == 0 || self%(p+22) == 0 || self%(p+24) == 0 end true end ... ...As in typical prime algorithms, you just check for the primes <= the number's squareroot.

HereMath.sqrt(self).to_idetermines the integer squareroot, but as we know (and as I

found out, which initiated all this), as the numbers gets`> ~10**28`

, the approximation for

the squareroot becomes more and more larger than the actual value. This will causeprime?

to do increasingly more work than necessary as the numbers get bigger.

(Completely as an aside, for numbers this size you should probably look to use: n.to_bn.prime?)

Math.sqrt(n) can be smaller than mathematical sqrt(n).

% ruby -ve '1.upto(50) {|i| j = 10**i; j2 = Math.sqrt(j*j).to_i; p [i, j2 - j] if j != j2 }' ruby 2.5.0dev (2017-02-21 trunk 57675) [x86_64-linux] [23, -8388608] [24, -16777216] [25, 905969664] [26, 4764729344] [27, 13287555072] [28, -416880263168] [29, -8566849142784] [30, 19884624838656] [31, -364103705034752] [32, 5366162204393472] [33, -54424769012957184] [34, -544247690129571840] [35, -3136633892082024448] [36, 42420637374017961984] [37, -461237341797878857728] [38, -2251190176543965970432] [39, -60290833628396821413888] [40, 303786028427003666890752] [41, 620008645040778319495168] [42, 44885712678075916785549312] [43, 139372116959414099130712064] [44, -10985679223259661757684121600] [45, -70242710975464448780069240832] [46, -68601809640529787052340805632] [47, 4384584304507619735463404765184] [48, 43845843045076197354634047651840] [49, -535097230524518206803127585210368] [50, 7629769841091887003294964970946560]

So, Integer#prime? can underestimate the limit to search factors.

Assume i is 1000000000000000031.

I think (i*i).prime? returns true.

1000000000000000031 is not searched because

Math.sqrt(i*i).to_i is 1000000000000000000.

(I couldn't complete (i*i).prime? because it takes too long time.)

% ruby -e 'i=1000000000000000031; p i, Math.sqrt(i**2).to_i' 1000000000000000031 1000000000000000000

I feel this is a good evidence that we need integer sqrt.

It would be interesting to see what an audit of the entire Ruby code base would reveal where

this is used to provide the integer squareroot, versus the floating point one.

Simple search result of 'sqrt(.*to_i' using gem-codesearch.

https://github.com/akr/gem-codesearch

Oompa-euler-1.0.3/lib/euler.rb: max = Math.sqrt(to_factor).to_i + 1 Oompa-euler-1.0.3/lib/euler.rb: 5.step(Math.sqrt(self).to_i, 2) do |x| Oompa-euler-1.0.3/lib/euler.rb: 3.step(Math.sqrt(max).to_i, 2) { |i| algebra-0.2.3/lib/algebra/polynomial-factor-int.rb: (Math.log(c * Math.sqrt(s))/Math.log(u)).to_i + 1 blackwinter-gsl-1.15.3.2/examples/random/randomwalk.rb:sigma = Math::sqrt(N).to_i boggle-0.0.2/bin/boggle:opts[:size] = Math.sqrt(opts[:board].size).to_i if !opts[:board].empty? boggler-0.0.1/lib/boggler/grid.rb: @size = Math.sqrt(@dice.length).to_i boqwij-0.1.0/lib/boqwij/integer_extensions.rb: 5.step(Math.sqrt(self).to_i, 2) do |x| boqwij-0.1.0/lib/boqwij/integer_extensions.rb: max = Math.sqrt(to_factor).to_i + 1 clusterer-0.1.9/lib/clusterer/clustering.rb: options[:no_of_clusters] ||= Math.sqrt(objects.size).to_i collective-0.2.1/lib/collective/checker.rb: @estimate > 0 ? Math.sqrt(@estimate).to_i : 1 cw-0.3.3/lib/cw/read.rb: # @magnitude = Math.sqrt(magnitude_squared).to_i euler-1.1.0/lib/euler.rb: max = Math.sqrt(to_factor).to_i + 1 euler-1.1.0/lib/euler.rb: 5.step(Math.sqrt(self).to_i, 2) do |x| euler-1.1.0/lib/euler.rb: 3.step(Math.sqrt(max).to_i, 2) { |i| euler-1.1.0/lib/rudoku.rb: @sqrt_n = Math.sqrt(@n).to_i gecoder-1.1.1/lib/gecoder/interface/constraints/int/arithmetic.rb: Gecode::Raw::sqrt(@model.active_space, int_op.to_int_var.bind, gecoder-with-gecode-1.1.1.1/lib/gecoder/interface/constraints/int/arithmetic.rb: Gecode::Raw::sqrt(@model.active_space, int_op.to_int_var.bind, gphys-1.5.1/lib/numru/gphys/interpolate.rb: ndiv[i] = Math.sqrt( (kx[i+1]-kx[i])**2 + (ky[i+1]-ky[i])**2).to_i green_shoes-1.1.374/samples/sample56.rb: .text(lambda {|d| d.class_name[0, Math.sqrt(d.node_value).to_i / 8]}) gsl-2.1.0.2/examples/random/randomwalk.rb:sigma = Math::sqrt(N).to_i gsl-nmatrix-1.15.3.2/examples/random/randomwalk.rb:sigma = Math::sqrt(N).to_i hotcocoa-0.6.3/lib/hotcocoa/graphics/elements/sandpainter.rb: #@grains = (sqrt((ox-x)*(ox-x)+(oy-y)*(oy-y))).to_i huge_enumerable-0.0.2/lib/huge_enumerable.rb: increment = next_prime(( 2 * domain_size / (1 + Math.sqrt(5)) ).to_i) ifmapper-2.0.4/lib/IFMapper/FXSpline.rb: num = Math.sqrt( x2 + y2 ).to_i / 16 jwthumbs-0.1.0/lib/shutter.rb: gridsize = Math.sqrt(@images.length).ceil.to_i lazylist-0.3.2/examples/examples.rb: not (2..Math.sqrt(x).to_i).find { |d| x % d == 0 } magic_cloud-0.0.4/lib/magic_cloud/layouter/place.rb: case (Math.sqrt(1 + 4 * sign * t) - sign).to_i & 3 mapbaidu-0.0.7/lib/map_baidu.rb: x = Math.sqrt(radius**2 - y**2).ceil.to_i mapcache-0.1.1/lib/mgmaps_export.rb: tpf_x = Math.sqrt(tpf).to_i mapcache-0.1.1/lib/mgmaps_export.rb: tpf_x = Math.sqrt(tpf * 2).to_i matrix_disp-0.0.3/lib/matrix_disp.rb: for i in (0...Math.sqrt(other[0].size).to_i) do matrix_disp-0.0.3/lib/matrix_disp.rb: for x in (0...Math.sqrt(other[0].size).to_i) do matrix_disp-0.0.3/lib/matrix_disp.rb: for y in (0...Math.sqrt(other[0].size).to_i) do matrix_disp-0.0.3/lib/matrix_disp.rb: adj[z] = other[0][x+(y*Math.sqrt(other[0].size).to_i)] matrix_disp-0.0.3/lib/matrix_disp.rb: resultado = resultado + tmp[i+(j*Math.sqrt(other[0].size).to_i)]*det(adj) matrix_disp-0.0.3/lib/matrix_disp.rb: resultado = resultado - tmp[i+(j*Math.sqrt(other[0].size).to_i)]*det(adj) minesweeper-console-1.0.1/test/minesweeper/console/prettyprinter/minefield_pretty_printer_test.rb: Math.sqrt(@string_representation.length).to_i mlanett-hive-0.4.0/lib/hive/checker.rb: @estimate > 0 ? Math.sqrt(@estimate).to_i : 1 mtah-ruby-treemap-0.0.3.2/lib/treemap/node.rb: (base_size * Math.sqrt(self.bounds.width * self.bounds.height) / 125).to_i naga_perfectsquare-1.0.1/lib/naga_perfectsquare.rb: Math.sqrt(self) == Math.sqrt(self).to_i narray-0.6.1.2/src/lib/narray/narray_ext.rb: m = ((n+Math::sqrt(n))*1.27).to_i narray-bigmem-0.0.0/lib/narray_ext.rb: m = ((n+Math::sqrt(n))*1.27).to_i nsudoku-0.2.0/lib/nsudoku/checker.rb: @width = Math.sqrt(@sudoku.length).to_i nsudoku-0.2.0/lib/nsudoku/checker.rb: @block_width = Math.sqrt(@width).to_i nsudoku-0.2.0/lib/nsudoku/solver.rb: @width = Math.sqrt(sudoku.length).to_i nsudoku-0.2.0/lib/nsudoku/solver.rb: @block_width = Math.sqrt(@width).to_i nsudoku-0.2.0/lib/nsudoku.rb: @width = Math.sqrt(@sudoku.length).to_i nsudoku-0.2.0/lib/nsudoku.rb: @block_width = Math.sqrt(@width).to_i numb-0.186.0/lib/numb/sum_of_squares.rb: s = Math.sqrt(i).to_i numeric_memoist-0.0.2/examples/isqrt_example.rb: $stderr.puts "calculate sqrt(#{self}).to_i (<= MaxIntFloat)" numeric_memoist-0.0.2/examples/isqrt_example.rb: return Math.sqrt(self).to_i numeric_memoist-0.0.2/examples/isqrt_example.rb: $stderr.puts "calculate sqrt(#{self}).to_i" numeric_memoist-0.0.2/examples/isqrt_example.rb: # STDERR> calculate sqrt(123456789).to_i (<= MaxIntFloat) numeric_memoist-0.0.2/examples/isqrt_example.rb: # STDERR> calculate sqrt(986960440108935861883449099987613301336842162813430234017512025).to_i numru-narray-1.0.3/lib/numru/narray_ext.rb: m = ((n+Math::sqrt(n))*1.27).to_i pixeldistance-0.2.1/lib/pixeldistance.rb: Math.sqrt((x1-x2)**2 + (y1-y2)**2).to_i >> (21 - zoom) prime-multiplication-1.0/lib/PrimeMultiplication.rb: (2..(Math.sqrt(number).to_i)).each do |i| prime_numbers-0.0.4/lib/prime_numbers/prime_int.rb: 2.upto(Math.sqrt(num).to_i) do |x| prime_table-0.0.2/lib/prime_table/prime.rb: sqrt = Math.sqrt(self).to_i primes-utils-2.7.0/lib/primes/utils.rb: sqrtN = Math.sqrt(n).to_i primes-utils-2.7.0/lib/primes/utils.rb: sqrtN = Math.sqrt(n).to_i primes-utils-2.7.0/lib/primes/utils.rb: factors << p; r -=1; n /= p; sqrtN = Math.sqrt(n).to_i primes-utils-2.7.0/lib/primes/utils.rb: sqrtN = Math.sqrt(num).to_i # sqrt of end_num (end of range) primes-utils-2.7.0/lib/primes/utils.rb: if start_num <= Math.sqrt(num).to_i # for one array of primes upto N print_square-0.0.2/lib/print_square/command_runner.rb: size = Math.sqrt(number).to_i quasi-0.0.6/lib/halton.rb: m=Math.sqrt(q).to_int rlsm-1.8.1/lib/rlsm/monoid.rb: @order = Math::sqrt(@table.size).to_i royw-drbman-0.0.6/examples/primes/lib/primes/sieve_of_eratosthenes.rb: sqr_primes = primes(Math.sqrt(maximum).to_i, drbman) ruby-compiler-0.1.1/vendor/ruby/lib/prime.rb: (7..Math.sqrt(self).to_i).step(30) do |p| ruby-spark-1.2.1/benchmark/comparison/ruby.rb: upper = Math.sqrt(x.to_f).to_i ruby-spark-1.2.1/lib/spark/rdd.rb: max_sample_size = Integer::MAX - (num_st_dev * Math.sqrt(Integer::MAX)).to_i ruby-treemap-0.0.4/lib/treemap/node.rb: (base_size * Math.sqrt(self.bounds.width * self.bounds.height) / 125).to_i ruby-treemap-fork-0.0.4/lib/treemap/node.rb: (base_size * Math.sqrt(self.bounds.width * self.bounds.height) / 125).to_i rubyvis-0.6.1/examples/5_pv_hierarchies/bubble_charts.rb: .text(lambda {|d| d.class_name[0, Math.sqrt(d.node_value).to_i / 8]}) rusby-0.1.2/lib/rusby/profiler.rb: SQRT_ITERATIONS = Math.sqrt(1000).to_i simstring_pure-1.1.1/lib/simstring_pure.rb: (alpha * Math.sqrt(query_size * y_size)).ceil.to_i sudoku_maker-0.0.2/lib/sudoku_maker/maker.rb: @decision_num = Math.sqrt(num).to_i sudoku_solver-1.0.0/lib/sudoku_solver/solver.rb: @dimension = Math.sqrt(@width).to_i tactical_tic_tac_toe-1.0.0/lib/board.rb: Math.sqrt(marked_spaces.count).to_i taskjuggler-3.6.0/lib/taskjuggler/reports/ChartPlotter.rb: rr = (r / Math.sqrt(2.0)).to_i tdiary-contrib-5.0.3/plugin/category_to_tagcloud.rb: level = ((Math.sqrt(count) - min) * factor).to_i tic_tac_toes-0.1.8/lib/tic_tac_toes/ui/serializer.rb: row_size = Math.sqrt(board_structure.count).to_i ul-wukong-4.1.1/lib/wukong/widget/reducers/bin.rb: self.num_bins = Math.sqrt(total_count).to_i whitestone-1.0.2/etc/extra_tests/output_examples_code.rb: (2..Math.sqrt(n).to_i).map { |i| n % i }.all? { |rem| rem != 0 } wukong-4.0.0/lib/wukong/widget/reducers/bin.rb: self.num_bins = Math.sqrt(total_count).to_i

#### Updated by jzakiya (Jabari Zakiya) almost 5 years ago

I did a simplifications in ** isqrt1**, replacing all the '+'s with '|' operations,

leaving only the math subtraction operation

`n -= t`

. I'll see if I can get rid of that.This makes it a bit faster, especially for numbers `< 10**400`

, where it starts being

faster than ** iroot2**.

This suggests, for optimal overall speed, that a hybrid method can be crafted, like I

proposed for using **Math.sqrt(n).to_i**, like:

```
def sqrt_i(n); n < UPPER_BOUND ? n.iroot2 : isqrt1(n) end
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 isqrt(n)
r = 0
x = 1 << ((n.bit_length >> 1 ) << 1)
while x != 0 do
t = r + x
if n >= t
n -= t
r = (r >> 1) + x
else
r >>= 1
end
x >>= 2
end
r
end
def isqrt1(n)
r = 0
x = 1 << ((n.bit_length >> 1 ) << 1)
until x == 0
t = r | x
r >>= 1
(n -= t; r |= x) if n >= t
x >>= 2
end
r
end
Benchmark.ips do |x|
exp = 400; n = 10**exp
puts "integer squareroot tests for n = 10**#{exp}"
x.report("iroot2" ) { n.iroot2 }
# x.report("irootn(2)" ) { n.irootn(2) }
x.report("isqrt(n)" ) { isqrt(n) }
x.report("isqrt1(n)" ) { isqrt1(n) }
x.compare!
end
```

2.4.0 :118 > load 'irootstest.rb' integer squareroot tests for n = 10**50 Warming up -------------------------------------- iroot2 4.157k i/100ms isqrt(n) 3.590k i/100ms isqrt1(n) 4.029k i/100ms Calculating ------------------------------------- iroot2 41.870k (± 3.7%) i/s - 212.007k in 5.070841s isqrt(n) 36.039k (± 3.6%) i/s - 183.090k in 5.087128s isqrt1(n) 40.580k (± 3.4%) i/s - 205.479k in 5.069789s Comparison: iroot2: 41870.0 i/s isqrt1(n): 40579.8 i/s - same-ish: difference falls within error isqrt(n): 36039.1 i/s - 1.16x slower => true 2.4.0 :119 > load 'irootstest.rb' integer squareroot tests for n = 10**100 Warming up -------------------------------------- iroot2 1.385k i/100ms isqrt(n) 879.000 i/100ms isqrt1(n) 883.000 i/100ms Calculating ------------------------------------- iroot2 14.004k (± 3.2%) i/s - 70.635k in 5.049521s isqrt(n) 8.558k (± 3.9%) i/s - 43.071k in 5.040986s isqrt1(n) 9.173k (± 3.8%) i/s - 45.916k in 5.012857s Comparison: iroot2: 14003.5 i/s isqrt1(n): 9173.2 i/s - 1.53x slower isqrt(n): 8557.9 i/s - 1.64x slower => true 2.4.0 :120 > load 'irootstest.rb' integer squareroot tests for n = 10**200 Warming up -------------------------------------- iroot2 425.000 i/100ms isqrt(n) 377.000 i/100ms isqrt1(n) 393.000 i/100ms Calculating ------------------------------------- iroot2 4.248k (± 4.1%) i/s - 21.250k in 5.010948s isqrt(n) 3.722k (± 4.5%) i/s - 18.850k in 5.075627s isqrt1(n) 3.876k (± 4.9%) i/s - 19.650k in 5.081749s Comparison: iroot2: 4248.2 i/s isqrt1(n): 3876.3 i/s - 1.10x slower isqrt(n): 3721.7 i/s - 1.14x slower => true 2.4.0 :121 > load 'irootstest.rb' integer squareroot tests for n = 10**300 Warming up -------------------------------------- iroot2 252.000 i/100ms isqrt(n) 219.000 i/100ms isqrt1(n) 238.000 i/100ms Calculating ------------------------------------- iroot2 2.477k (± 6.5%) i/s - 12.348k in 5.008024s isqrt(n) 2.263k (± 6.6%) i/s - 11.388k in 5.060793s isqrt1(n) 2.407k (± 4.0%) i/s - 12.138k in 5.050951s Comparison: iroot2: 2476.8 i/s isqrt1(n): 2407.2 i/s - same-ish: difference falls within error isqrt(n): 2262.7 i/s - same-ish: difference falls within error => true 2.4.0 :122 > load 'irootstest.rb' integer squareroot tests for n = 10**400 Warming up -------------------------------------- iroot2 101.000 i/100ms isqrt(n) 156.000 i/100ms isqrt1(n) 165.000 i/100ms Calculating ------------------------------------- iroot2 972.069 (± 6.0%) i/s - 4.848k in 5.006825s isqrt(n) 1.495k (± 7.3%) i/s - 7.488k in 5.038471s isqrt1(n) 1.611k (± 7.0%) i/s - 8.085k in 5.045823s Comparison: isqrt1(n): 1610.8 i/s isqrt(n): 1494.9 i/s - same-ish: difference falls within error iroot2: 972.1 i/s - 1.66x slower

#### Updated by jzakiya (Jabari Zakiya) almost 5 years ago

Thanks Akira for that list.

So, too small integer squareroots potentially can cause miss finding primes in some apps,

and being too large can create doing more work than necessary in some apps, and just being

wrong in other apps has unknown consequences. :(

I know how this affected my specific apps, and was diligent (lucky) enough to track it down.

How many other situations are out there being affected by this that people don't even know

are giving wrong answers, or know something is amiss and have no hint of what's causing it?

#### Updated by jzakiya (Jabari Zakiya) almost 5 years ago

So the clear winner, and new champeen is...** newtons_fast**, at least for squareroots.

Just goes to show, don't f**k with Newton!

```
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 isqrt1(n)
r = 0
x = 1 << ((n.bit_length >> 1 ) << 1)
until x == 0
t = r | x
r >>= 1
(n -= t; r |= x) if n >= t
x >>= 2
end
r
end
def newtons_fast(n) # Newton's method with √n < initial_x <= 2*√n
raise if n<0
return n if n<=1
x = 1<<((n.bit_length+1)/2)
y = (x + n/x)/2
while y < x
x = y
y = (x + n/x)/2
end
x
end
Benchmark.ips do |x|
exp = 400; n = 10**exp
puts "integer squareroot tests for n = 10**#{exp}"
x.report("iroot2" ) { n.iroot2 }
x.report("isqrt1(n)" ) { isqrt1(n) }
x.report("newtons_fast(n)") { newtons_fast(n)
x.compare!
end
```

2.4.0 :138 > load 'irootstest.rb' integer squareroot tests for n = 10**5 Warming up -------------------------------------- iroot2 85.000k i/100ms isqrt1(n) 87.646k i/100ms newtons_fast(n) 191.796k i/100ms Calculating ------------------------------------- iroot2 1.014M (± 2.9%) i/s - 5.100M in 5.031392s isqrt1(n) 1.049M (± 3.0%) i/s - 5.259M in 5.018645s newtons_fast(n) 2.846M (± 4.6%) i/s - 14.385M in 5.066760s Comparison: newtons_fast(n): 2845709.1 i/s isqrt1(n): 1048817.1 i/s - 2.71x slower iroot2: 1014498.5 i/s - 2.81x slower => true 2.4.0 :139 > load 'irootstest.rb' integer squareroot tests for n = 10**50 Warming up -------------------------------------- iroot2 4.077k i/100ms isqrt1(n) 3.992k i/100ms newtons_fast(n) 24.377k i/100ms Calculating ------------------------------------- iroot2 41.189k (± 6.1%) i/s - 207.927k in 5.071463s isqrt1(n) 40.293k (± 3.4%) i/s - 203.592k in 5.058807s newtons_fast(n) 260.960k (± 3.4%) i/s - 1.316M in 5.050191s Comparison: newtons_fast(n): 260959.7 i/s iroot2: 41188.8 i/s - 6.34x slower isqrt1(n): 40292.7 i/s - 6.48x slower => true 2.4.0 :140 > load 'irootstest.rb' integer squareroot tests for n = 10**500 Warming up -------------------------------------- iroot2 74.000 i/100ms isqrt1(n) 127.000 i/100ms newtons_fast(n) 3.775k i/100ms Calculating ------------------------------------- iroot2 741.359 (± 4.6%) i/s - 3.700k in 5.001895s isqrt1(n) 1.258k (± 7.0%) i/s - 6.350k in 5.077988s newtons_fast(n) 37.762k (± 4.5%) i/s - 188.750k in 5.008756s Comparison: newtons_fast(n): 37761.5 i/s isqrt1(n): 1257.6 i/s - 30.03x slower iroot2: 741.4 i/s - 50.94x slower => true 2.4.0 :141 > load 'irootstest.rb' integer squareroot tests for n = 10**5000 Warming up -------------------------------------- iroot2 1.000 i/100ms isqrt1(n) 5.000 i/100ms newtons_fast(n) 340.000 i/100ms Calculating ------------------------------------- iroot2 11.196 (± 8.9%) i/s - 56.000 in 5.013699s isqrt1(n) 53.056 (± 3.8%) i/s - 265.000 in 5.002572s newtons_fast(n) 3.402k (± 3.1%) i/s - 17.000k in 5.002400s Comparison: newtons_fast(n): 3401.9 i/s isqrt1(n): 53.1 i/s - 64.12x slower iroot2: 11.2 i/s - 303.85x slower

#### Updated by jzakiya (Jabari Zakiya) almost 5 years ago

Another reason this needs to change is because **Math.sqrt(n).to_i**

has an upper limit before crashing.

Here's something else to consider:

There is all that code out there using **Math.sqt(n).to_i** and I presume

most people expect it to return the correct result (which many probably

don't even realize what that should be).

Is there some possible way for the smart Ruby Gurus to fix the parser so

that when it sees the phrase **Math.sqrt(xyz).to_i** it uses whatever

chosen correct implementation decided upon, so all that old legacy code

out there, in the core, gems, and anywhere else, will just magically work?

I know it's a hack, but Ruby is designed to make programmers happy, and it

wasn't their fault this bug exists. I know it would make me happy to not

have to go back over everything I've done and replace it with a new method name.

Just something to consider.

2.4.0 :166 > load 'irootstest.rb' integer squareroot tests for n = 10**308 Warming up -------------------------------------- newtons_fast(n) 5.042k i/100ms Math.sqrt(n).to_i 167.539k i/100ms Calculating ------------------------------------- newtons_fast(n) 51.171k (± 4.3%) i/s - 257.142k in 5.034762s Math.sqrt(n).to_i 2.411M (± 6.1%) i/s - 12.063M in 5.022522s Comparison: Math.sqrt(n).to_i: 2411213.7 i/s newtons_fast(n): 51170.9 i/s - 47.12x slower => true 2.4.0 :167 > load 'irootstest.rb' integer squareroot tests for n = 10**309 Warming up -------------------------------------- newtons_fast(n) 4.190k i/100ms Math.sqrt(n).to_iFloatDomainError: Infinity from irootstest.rb:93:in `to_i' from irootstest.rb:93:in `block (2 levels) in <top (required)>' from /home/jzakiya/.rvm/gems/ruby-2.4.0/gems/benchmark-ips-2.7.2/lib/benchmark/ips/job/entry.rb:53:in `call_times' from /home/jzakiya/.rvm/gems/ruby-2.4.0/gems/benchmark-ips-2.7.2/lib/benchmark/ips/job.rb:220:in `block in run_warmup' from /home/jzakiya/.rvm/gems/ruby-2.4.0/gems/benchmark-ips-2.7.2/lib/benchmark/ips/job.rb:206:in `each' from /home/jzakiya/.rvm/gems/ruby-2.4.0/gems/benchmark-ips-2.7.2/lib/benchmark/ips/job.rb:206:in `run_warmup' from /home/jzakiya/.rvm/gems/ruby-2.4.0/gems/benchmark-ips-2.7.2/lib/benchmark/ips/job.rb:185:in `block in run' from /home/jzakiya/.rvm/gems/ruby-2.4.0/gems/benchmark-ips-2.7.2/lib/benchmark/ips/job.rb:184:in `times' from /home/jzakiya/.rvm/gems/ruby-2.4.0/gems/benchmark-ips-2.7.2/lib/benchmark/ips/job.rb:184:in `run' from /home/jzakiya/.rvm/gems/ruby-2.4.0/gems/benchmark-ips-2.7.2/lib/benchmark/ips.rb:59:in `ips' from irootstest.rb:85:in `<top (required)>' from (irb):167:in `load' from (irb):167 from /home/jzakiya/.rvm/rubies/ruby-2.4.0/bin/irb:11:in `<main>' 2.4.0 :168 >

#### Updated by Student (Nathan Zook) almost 5 years ago

First, please, PLEASE do not credit me with the optimized bit-by-bit method. I'm a formally trained mathematician, and I DO NOT want credit for someone else's work. It came from the above-referenced Reddit post.

Second, it appears that we must have a pretty efficient MP divide, or any direct version of Newton's method will perform poorly. Obviously, I now need to go ahead & see what I get with Zimmerman's method...

And probably look into Newton-Raphson for comparison. As I stated, this will take work. BTW, performance of Newton's method depends a lot on N. If I come up with a rival, we will likely need a lot of tests to validate speed claims.

Third, please don't get me started on Prime.rb when it comes to large numbers. :D The implementation there is simply not to be used even for full-word numbers. Try timing (1 << 56) + 81 to see what I mean. I've actually done experiments with Miller-Rabin, and...I'll probably be submitting some code soon.

I had not thought about the implications of Math.sqrt overestimating. It is a float, so yeah.

Actually, I need to leave prime.rb alone if I'm going to work here. There is some crazy work done for small (< 2^{64)} numbers...

#### Updated by Student (Nathan Zook) almost 5 years ago

# Core Algorithm by Paul Zimmerman, article entitled # Karatsuba Square Root # https://hal.inria.fr/inria-00072854/document # Presumably used in GMP. # Personal computations derived from implementation details (page 5) # n >= b**4 / 4 # b = 2**64**k # n * 4 >= b ** 4 # n.length + 2 >= b.length * 4 # (n.length + 2) >> 2 >= b.length # (n.length + 2) >> 2 >= 64 * k # ((n.length + 2) >> 2) / 64 = k def sqrtrem(n) nlength = n.bit_length if nlength >= (64 * 4 - 2) bbits = (n.bit_length + 2 >> 2) / 64 * 64 elsif nlength >= (32 * 4 - 2) bbits = (n.bit_length + 2 >> 2) / 32 * 32 elsif nlength >= (16 * 4 - 2) bbits = (n.bit_length + 2 >> 2) / 16 * 16 else # single word now -- my computer has 64-bit mantissas! s = Math.sqrt(n).to_i r = n - s * s return [s, r] end b = 1 << bbits bmask = b - 1 a3 = n >> bbits * 3 a2 = (n >> bbits * 2) & bmask a1 = (n >> bbits ) & bmask a0 = n & bmask i, j = sqrtrem((a3 << bbits) + a2) q, u = ((j << bbits) + a1).divmod(i << 1) s = (i << bbits) + q r = (u << bbits) + a0 - q ** 2 if r < 0 r += (s << 1) - 1 s -= 1 end [s, r] end def sqrt_z(n) s, r = sqrtrem(n) s end Benchmark.ips do |x| n = 10 ** 5000 x.report("iroot2") { n.iroot2 } x.report("sqrt_z") { sqrt_z(n) } end Warming up -------------------------------------- iroot2 1.000 i/100ms sqrt_z 987.000 i/100ms Calculating ------------------------------------- iroot2 8.665 (± 0.0%) i/s - 44.000 in 5.078551s sqrt_z 10.015k (± 1.6%) i/s - 50.337k in 5.027514s =>

I've been at this for way too long today.

A few notes:

1) Assuming no architectural weirdness, it looks like--for this PARTICULAR case--Zimmerman's algorithm is doing about 4x better than what you are calling newtons_fast. First, this is not a surprise, as no one uses Newton's method this way. What is interesting is that Zimmerman's algorithm is designed to be good on GMP, and I am using stock ruby.

2) That hack at the top can almost certainly be dramatically improved (It's my contribution). This is a demo of Zimmerman's work--which doesn't even really want to talk about numbers smaller than 2**254. Please don't bother looking at performance under 1000 bits until that is cleaned up.

3) In a note that Zimmerman included in his article, he indicated that Newton-Raphson would probably be faster. As I stated, I have been unable to find an implementation, and I am loathe to present anything without high confidence that it is correct. I will do some more looking tomorrow.

4) In some sense it is silly to do timing tests like this. We have no idea how much time is being spent in ruby and how much in C. It is entirely likely that we could have 10x speed improvements (or much, much more) across certain ranges for some of these methods by dropping into C. As I said, it is a question of how much work we want to do. At a minimum, we should do a GMP-linked build and see what Zimmerman's work can really do. Again, I don't know the legal considerations, and I would like guidance before attempting to reimplement an MP library from scratch.

#### Updated by Student (Nathan Zook) almost 5 years ago

Sorry to spam like this, but I really have to disagree regarding some of the implications drawn for Math.sqrt(n). This returns a float. If you don't know about the pointy edges of IEEE-754, but you do care, then shame on you. But allow me to calm you regarding Numeric#prime?. Let p be the largest prime < 2 ** 32. Run (p * p).prime?. When we fix that problem, we can take care of the rounding issue.

#### Updated by Student (Nathan Zook) almost 5 years ago

Following the discussion from the second answer at http://cs.stackexchange.com/questions/37596/arbitrary-precision-integer-square-root-algorithm, I implemented an integer Newton - Raphson on the inverse square root. This involved correcting an error and clarifying some things. The results were disappointing, as the asymptotic behavior is clearly beaten by the transcribed Zimmerman algorithm. The crossover was somewhat beyond 10 ** 1000. For reference, I made a faster Newton method by using Math.sqrt for the initial guess.

I also adjusted my transcription of the Zimmerman algorithm to reflect the fact that Matz implemented a simple , portable Math.sqrt. Before we use it, we would need to do serious testing, of course.

One more round of benchmarks. I have removed the warmup & comparison sections. I also implemented a newton_faster method, which uses Math.sqrt for the initial guess.

integer squareroot tests for n = 10**50 Calculating ------------------------------------- newtons_fast(n) 176.714k (± 1.3%) i/s - 891.495k in 5.045686s newton_faster(n) 389.709k (± 0.8%) i/s - 1.952M in 5.009508s sqrt_z(n) 166.894k (± 1.4%) i/s - 847.715k in 5.080312s inverse Newton(n) 250.008k (± 0.4%) i/s - 1.271M in 5.083042s integer squareroot tests for n = 10**500 Calculating ------------------------------------- newtons_fast(n) 31.530k (± 3.4%) i/s - 158.496k in 5.032977s newton_faster(n) 61.087k (± 3.1%) i/s - 306.020k in 5.014480s sqrt_z(n) 33.167k (± 0.5%) i/s - 167.128k in 5.039057s inverse Newton(n) 51.252k (± 3.6%) i/s - 257.972k in 5.041610s integer squareroot tests for n = 10**1000 Calculating ------------------------------------- newtons_fast(n) 14.085k (± 0.8%) i/s - 71.656k in 5.087601s newton_faster(n) 24.078k (± 3.3%) i/s - 121.074k in 5.034512s sqrt_z(n) 25.410k (± 2.3%) i/s - 127.500k in 5.020581s inverse Newton(n) 28.921k (± 2.7%) i/s - 145.350k in 5.030223s integer squareroot tests for n = 10**2000 Calculating ------------------------------------- newtons_fast(n) 3.994k (± 0.7%) i/s - 20.247k in 5.069925s newton_faster(n) 6.698k (± 0.7%) i/s - 34.017k in 5.079015s sqrt_z(n) 19.176k (± 2.3%) i/s - 96.951k in 5.058731s inverse Newton(n) 13.724k (± 1.7%) i/s - 68.952k in 5.025785s integer squareroot tests for n = 10**4000 Calculating ------------------------------------- newtons_fast(n) 1.035k (± 0.4%) i/s - 5.253k in 5.074204s newton_faster(n) 1.600k (± 0.9%) i/s - 8.000k in 5.001751s sqrt_z(n) 12.478k (± 0.7%) i/s - 62.475k in 5.006905s inverse Newton(n) 6.094k (± 0.7%) i/s - 30.957k in 5.079890s integer squareroot tests for n = 10**5000 Calculating ------------------------------------- newtons_fast(n) 628.308 (± 0.6%) i/s - 3.162k in 5.032783s newton_faster(n) 915.131 (± 0.4%) i/s - 4.641k in 5.071494s sqrt_z(n) 10.107k (± 0.6%) i/s - 50.796k in 5.026155s inverse Newton(n) 4.539k (± 3.0%) i/s - 22.850k in 5.038822s

It is possible that "inverse Newton" could be modified to advantage the particulars of our implementation of bignum, or GMP, should it be linked. I don't want spam any more performances after this until it is determined if we want to drop down to C for this.

One last thing to play with might be Halley's method. This converges cubically, but requires 4 multiplies. This should compute faster. But don't expect to see it very soon.

#### Updated by Student (Nathan Zook) almost 5 years ago

All of the algorithms are now at https://github.com/NathanZook/ruby_sqrt. Hope that's okay.

#### Updated by jzakiya (Jabari Zakiya) almost 5 years ago

Hey Nathan, Great Work!

I'm seeing an error in **inverse_newton_sqrt(n)** being one-off.

2.4.0 :224 > n = 10**10 ; puts sqrt_z(n), newtons_fast(n), newton_faster(n), inverse_newton_sqrt(n) 100000 100000 100000 100000 => nil 2.4.0 :225 > n = 10**11 ; puts sqrt_z(n), newtons_fast(n), newton_faster(n), inverse_newton_sqrt(n) 316227 316227 316227 316227 => nil 2.4.0 :226 > n = 10**12 ; puts sqrt_z(n), newtons_fast(n), newton_faster(n), inverse_newton_sqrt(n) 1000000 1000000 1000000 1000000 => nil 2.4.0 :227 > n = 10**13 ; puts sqrt_z(n), newtons_fast(n), newton_faster(n), inverse_newton_sqrt(n) 3162277 3162277 3162277 3162277 => nil 2.4.0 :228 > n = 10**14 ; puts sqrt_z(n), newtons_fast(n), newton_faster(n), inverse_newton_sqrt(n) 10000000 10000000 10000000 10000000 => nil 2.4.0 :229 > n = 10**15 ; puts sqrt_z(n), newtons_fast(n), newton_faster(n), inverse_newton_sqrt(n) 31622776 31622776 31622776 31622776 => nil 2.4.0 :230 > n = 10**16 ; puts sqrt_z(n), newtons_fast(n), newton_faster(n), inverse_newton_sqrt(n) 100000000 100000000 100000000 99999999 => nil 2.4.0 :231 > n = 10**17 ; puts sqrt_z(n), newtons_fast(n), newton_faster(n), inverse_newton_sqrt(n) 316227766 316227766 316227766 316227765 => nil 2.4.0 :232 > n = 10**20 ; puts sqrt_z(n), newtons_fast(n), newton_faster(n), inverse_newton_sqrt(n) 10000000000 10000000000 10000000000 9999999999 => nil 2.4.0 :233 > n = 10**30 ; puts sqrt_z(n), newtons_fast(n), newton_faster(n), inverse_newton_sqrt(n) 1000000000000000 1000000000000000 1000000000000000 1000000000000000 => nil 2.4.0 :234 > n = 10**40 ; puts sqrt_z(n), newtons_fast(n), newton_faster(n), inverse_newton_sqrt(n) 100000000000000000000 100000000000000000000 100000000000000000000 99999999999999999999 => nil 2.4.0 :235 > n = 10**50 ; puts sqrt_z(n), newtons_fast(n), newton_faster(n), inverse_newton_sqrt(n) 10000000000000000000000000 10000000000000000000000000 10000000000000000000000000 9999999999999999999999999 => nil 2.4.0 :236 >

#### Updated by jzakiya (Jabari Zakiya) almost 5 years ago

Oh I see where at least part of the error comes from.

In this

def inverse_newton_sqrt(n) raise if n < 0 return Math.sqrt(n).to_i if n < 1 << 53

See top of this thread at `#1`

I showed earlier I found that **Math.sqrt(n).to_i** starts given incorrect

results for n > 9_999_899_999_899_999_322_536_673_279

I replaced `if n < 1 << 53`

with `if n <= 9_999_899_999_899_999_322_536_673_279`

This eliminated the shown errors upto n = 10**40, when it started to be one-off again.

This is on a 64-bit Linux OS I7 cpu laptop.

#### Updated by stomar (Marcus Stollsteimer) almost 5 years ago

Jabari Zakiya wrote:

I showed earlier I found that

Math.sqrt(n).to_istarts given incorrect

results for n > 9_999_899_999_899_999_322_536_673_279I replaced

`if n < 1 << 53`

with`if n <= 9_999_899_999_899_999_322_536_673_279`

Is this guaranteed to be platform independent?

And a more general consideration:

I assume `Integer.sqrt()`

will be implemented in C, not in Ruby? How significant or useful are all those benchmarks then? Micro-optimizing Ruby code for speed and then reimplementing it in C seems to be the wrong approach.

#### Updated by Student (Nathan Zook) almost 5 years ago

Jabari Zakiya wrote:

Oh I see where at least part of the error comes from.

In this

def inverse_newton_sqrt(n) raise if n < 0 return Math.sqrt(n).to_i if n < 1 << 53See top of this thread at

`#1`

I showed earlier I found thatMath.sqrt(n).to_istarts given incorrect

results for n > 9_999_899_999_899_999_322_536_673_279

9999899999899999322536673279 < 1 << 53 => false

So this cannot be the source of this particular problem. I am not surprised that there are errors in inverse_newton_sqrt. It was quickly put together, designed to credibly determine if it is faster than the Zimmerman algorithm. Getting the last bit correct on perfect squares is actually hard. If it appeared that it would be a superior algorithm, I would have needed to have proven convergence. I had expected that the +2 in ins_find_initial_exponent would have sufficed, but it looks like I probably needed to add another bit because of the use of fixed point rather than floating point in the algorithm.

I will credit you with finding the errors.

The question now is, "where do we go from here?" Clearly, faster Newton and Zimmerman are the contenders. At this moment, faster Newton wins almost to 10 * 1000, which is about 3000 bits, on speed alone. It is also clearly correct. (It conceivably does not converge, but that is easy to fix.)

The Zimmerman algorithm only works if I got that hack at the top correct. Actually, if we switch to a different algorithm at that point, we can probably eliminate a full pass (or two??) through the algorithm, which will speed it up lot for smaller numbers. I guess I should reach out to him & see if he will release his work to the Ruby license, in which case I can rely upon his proofs & techniques... :D

Both will benefit if we drop to C, but Zimmerman will gain a lot more. (Besides all of the ruby that happens in the loop, consider also the tricks that can be played with extracting a3-a0, for instance.)

#### Updated by Student (Nathan Zook) almost 5 years ago

Marcus Stollsteimer wrote:

Jabari Zakiya wrote:

I showed earlier I found that

Math.sqrt(n).to_istarts given incorrect

results for n > 9_999_899_999_899_999_322_536_673_279I replaced

`if n < 1 << 53`

with`if n <= 9_999_899_999_899_999_322_536_673_279`

Is this guaranteed to be platform independent?

If this were the issue, it would be. Math.sqrt() runs through some IEEE-754-specificed functions in C. I know of no funny business in these particular functions.

And a more general consideration:

I assume

`Integer.sqrt()`

will be implemented in C, not in Ruby? How significant or useful are all those benchmarks then? Micro-optimizing Ruby code for speed and then reimplementing it in C seems to be the wrong approach.

I do expect that we want to go with C, but I don't have an answer. The purpose of these benchmarks has been to get some idea of which algorithms are credible. We have run large enough samples that it is clear which one have better O() performance based on the way that ruby handles large numbers. So long as we are smart, that does not depend on the ruby code.

And see my just-posted comment for more. ;)

#### Updated by Student (Nathan Zook) almost 5 years ago

It is now clear that Halley's method requires five multiplies, not the four previously reported by wikipedia. This makes it uninteresting as a candidate.

#### Updated by nobu (Nobuyoshi Nakada) almost 5 years ago

**Status**changed from*Open*to*Closed*

#### Updated by Student (Nathan Zook) almost 5 years ago

Jabari Zakiya wrote:

Oh I see where at least part of the error comes from.

In this

def inverse_newton_sqrt(n) raise if n < 0 return Math.sqrt(n).to_i if n < 1 << 53See top of this thread at

`#1`

I showed earlier I found thatMath.sqrt(n).to_istarts given incorrect

results for n > 9_999_899_999_899_999_322_536_673_279I replaced

`if n < 1 << 53`

with`if n <= 9_999_899_999_899_999_322_536_673_279`

This eliminated the shown errors upto n = 10**40, when it started to be one-off again.

This is on a 64-bit Linux OS I7 cpu laptop.

Actually, with some understanding of how floating point arithmetic works, the first even root past 1 << 53 is likely to be bad:

irb(main):1081:0> def msi(n) x = Math.sqrt(n).to_i big = x * x > n small = x * x + 2 * x < n puts [n, x, big, small] end irb(main):1082:0> msi(1<<53) 9007199254740992 94906265 false false => nil irb(main):1083:0> msi(94906266 ** 2 - 1) 9007199326062755 94906266 true false => nil

;)

#### Updated by jzakiya (Jabari Zakiya) almost 5 years ago

I want to raise for consideration, as a user, the **API behavior** for these methods.

I noticed in the other squareroot methods they raise an error for `n < 0`

whereas my version returns 'nil' when taking an even root of a negative integer.

Let me promote why this shouldn't be treated as a runtime error.

The reason I chose to return 'nil' is because in some applications where I'm

processing lots of data I don't want the program to bomb when encountering

negative integers, I want to continue and just choose how to process those conditions.

First, mathematically it is perfectly allowable to take the even root of a negative number,

it's just a complex number, so conceptually this shouldn't be an error, just unwanted for this method.

With **irootn|2** I know I will either get a real integer, which means n was a valid integer (and root)

or 'nil', which means it didn't have a real (as in complex part) root.

If I wanted to check for negativity in the first place I could always do something like:

```
def my_app_sqrt(n); n < 0 ? Complex(0, n.abs.iroot2) : n.iroot2 end
```

or using the **root** method from my **Roots** gem:

```
def my_app_sqrt(n); n < 0 ? n.root(2) : n.iroot2 end
```

As a user, I want to choose how to handle cases when my app may try to take an even root

of a negative integer/number and not automatically throw an error that crashes my program.

I think as long as this is clearly documented (either way) for the API a user will know

exactly how to deal with these conditions, but returning 'nil' gives the user total control.

#### Updated by stomar (Marcus Stollsteimer) almost 5 years ago

Jabari Zakiya wrote:

I noticed in the other squareroot methods they raise an error for

`n < 0`

whereas my version returns 'nil' when taking an even root of a negative integer.

I strongly agree with nobu's implementation. That's exactly the behavior I was about to suggest when I noticed it was already there.

The reason I chose to return 'nil' is because in some applications where I'm

processing lots of data I don't want the program to bomb when encountering

negative integers, I want to continue and just choose how to process those conditions.

It doesn't have to bomb.

- for typical use cases (number theory, etc.) you will have only positive values anyway
- you can check your data before using the method
- you can use a begin...rescue block and handle the exception appropriately

Returning `nil`

for a mathematical operation doesn't make much sense. And a silent failure you are not aware of will bite you much more than an exception.

First, mathematically it is perfectly allowable to take the even root of a negative number,

it's just a complex number, so conceptually this shouldn't be an error, just unwanted for this method.

"Not an error, but unwanted"???

The point is that by premise the method doesn't make __any__ sense in the complex domain. It's supposed to return the largest integer (natural number) less or equal to the square root. This wouldn't work for complex numbers at all: how would you truncate to an integer? how would you do the comparison? (Complex numbers are not ordered.)

Therefore, isqrt only does make sense for natural numbers. (Different from the "normal" sqrt.)

#### Updated by jzakiya (Jabari Zakiya) almost 5 years ago

I disagree, but as long as the API explicitly states what's happens for this case a user can adapt accordingly

#### Updated by nobu (Nobuyoshi Nakada) almost 5 years ago

**Description**updated (diff)