Issue #13219 has been updated by Akira Tanaka.


Jabari Zakiya wrote:

> 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?)

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


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

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


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

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

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

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

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

require "benchmark/ips"

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

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

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

 => true 
2.4.0 :179 > 

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

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



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

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