From: Xavier Noria <fxn@...>
Date: 2007-02-13T03:32:18+09:00
Subject: Re: What's the correct and fast way to determine if a (gig) number is a perfect square?

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

# 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