From: Sam Kong <sam.s.kong@...>
Date: 2007-02-13T05:02:33+09:00
Subject: Re: What's the correct and fast way to determine if a (gig) number is a perfect square?

Hi Xavier,

I really appreciate your kind help.
I put my comment inline.

On Feb 12, 10:32 am, Xavier Noria <f...@hashref.com> wrote:
> On Feb 11, 2007, at 10:48 PM, Xavier Noria wrote:
>
> > There are faster algorithms based on the fact that most non-squares
> > aren't quadratic residues modulo some integers. I remember some is
> > explained in Henri Cohen's "A Course in Computational Algebraic
> > Number Theory", but do not have the book at hand.
>
> I've been able to consult the book today. There is a specialised
> algorithm to compute integer roots using integer arithmetic, and the
> integer square test itself. Thus, they guarantee a correct result. I
> attach them below translated to Ruby.
>
> Some operations would be written using bitwise stuff, but anyway the
> code performs very poorly (compared to Math::sqrt) according to the
> benchmarks. If performance is important a solution in C would be
> worth exploring.
>
> -- fxn
>
> # From Henri Cohen's "A Course in Computational Algebraic Number
> # Theory".
>
> # Algorithm 1.7.1 (Integer Square Root) Given a positive integer n,
> # this algorithm computes the integer part of the square root of n,
> # i.e. the number m such that m^2 <= n < (m + 1)^2.
> def isqrt(n)
>    x = n
>    loop do
>      y = ((x + n/x)/2)
>      if y < x
>        x = y
>      else
>        return x
>      end
>    end
> end

Actually you posted this message, I implemented an algorithm myself.
It's just the way we calculate sqrt with pen and paper.
It might be slow but correct even for big numbers.

def max_digit n, m
  result = 0
  (0..9).each do |i|
    break if ((n * 10) + i) * i > m
    result = i
  end
  result
end

def perfect_square? n
  arr = []
  a = n
  while true
    a, b = a / 100, a % 100
    arr.unshift b
    break if a == 0
  end
  remain = 0
  carried = 0
  arr.each do |i|
    num = remain * 100 + i
    digit = max_digit(carried, num)
    remain = num - (carried * 10 + digit) * digit
    carried = (carried * 10 + digit) + digit
  end
  remain == 0
end

I felt ashamed to see the code you posted.OTL

>
> # Cache the squares modulus 11, 63, 64, and 65. This is used to check
> # for non-squares, since a square is a square mod k for all k. The
> # choice of these numbers is based on the probability that a non-square
> # is a square mod any of them, which is 6/715, less than a 1%.
> $squares = {}
> [11, 63, 64, 65].each do |m|
>    $squares[m] = [false] * m
>    (0...(m/2)).each {|i| $squares[m][i**2 % m] = true}
> end
>
> # Algorithm 1.7.3 (Square Test). Given a positive integer n,
> # this algorithm determines whether n is a square or not,
> # and if it is, outputs the square root of n. We assume the
> # precomputations made above.
> def issquare(n)
>    return false unless $squares[64][n % 64]
>
>    r = n % 45045 # 45045 = 63*65*11
>    return false unless $squares[63][r % 63]
>    return false unless $squares[65][r % 65]
>    return false unless $squares[11][r % 11]
>
>    q = isqrt(n)
>    return q**2 == n ? q : false
> end
>
> require 'benchmark'
>
> $r = 1000
>
> # square of 32248581868698698768678697647823648238462348
> $s  =
> 103997103254216245831009973053627303273419242413513150637645310539211244
> 0875299413673104
>
> # non-square, the previous number minus 1
> $ns =
> 103997103254216245831009973053627303273419242413513150637645310539211244
> 0875299413673103
>
> # Just for the sake of curiosity, since the code based on Math::sqrt
> is not correct.
> Benchmark.benchmark do |x|
>    x.report("builtin is square (true)") do
>      1.upto($r) do
>        sqrt = Math::sqrt($s)
>        $s == sqrt.floor**2
>      end
>    end
>    x.report("modular is square (true)") do
>      1.upto($r) do
>        issquare($s)
>      end
>    end
>    x.report("builtin is square (false)") do
>      1.upto($r) do
>        sqrt = Math::sqrt($ns)
>        $ns == sqrt.floor**2
>      end
>    end
>    x.report("modular is square (false)") do
>      1.upto($r) do
>        issquare($ns)
>      end
>    end
> end

Thank you for your support.

Sam