Identifying Squares, Revisited

May 7, 2019

Some algorithms require the ability to rapidly identify numbers such as 36 and 121 that are squares; for instance, in Fermat’s factoring algorithm, and variants derived from it, identifiying squares is the performance bottleneck of the procedure. A naive program computes the integer square root of the target number, squares it, and reports a square if it is equal to the original target:

def isqrt(n):
  if n < 1: return 0
  u, s = n, n+1
  while u < s:
    s = u
    t = s + n // s
    u = t // 2
  return s
def isSquare(n): # naive
  s = isqrt(n); return s*s == n

(We are using Python instead of Scheme for this exercise because bit-operators like & are easier to type than functions like bitwise-and.)

This is very slow because the calculation to compute the integer square root takes a long time. In a previous exercise, we looked at a method for reducing the number of integer square roots that need to be computed:

def isSquare(n): # fenderbender
  m = n & 127
  if ((m*0x8bc40d7d) & (m*0xa1e2f5d1) & 0x14020a): return False
  largeMod = n % (63*25*11*17*19*23*31)
  m = largeMod % 63
  if ((m*0x3d491df7) & (m*0xc824a9f9) & 0x10f14008): return False
  m = largeMod % 25
  if ((m*0x1929fc1b) & (m*0x4c9ea3b2) & 0x51001005): return False
  m = 0xd10d829a * (largeMod % 31) 
  if (m & (m+0x672a5354) & 0x21025115): return False
  m = largeMod % 23
  if ((m*0x7bd28629) & (m*0xe7180889) & 0xf8300): return False
  m = largeMod % 19
  if ((m*0x1b8bead3) & (m*0x4d75a124) & 0x4280082b): return False
  m = largeMod % 17 
  if ((m*0x6736f323) & (m*0x9b1d499) & 0xc0000300): return False
  m = largeMod % 11 
  if ((m*0xabf1a3a7) & (m*0x2612bf93) & 0x45854000): return False
  s = isqrt(n); return s*s == n

This is wicked fast because it computes residue classes based on small primes and quickly identifies numbers that cannot possibly be a square; fenderbender eliminates the expensive square root calculation in 99.92% of cases. But all those magic numbers make the function wicked, as well as wicked fast, and it is easy to make an editing error when copying the function to a new program (you will quickly and accurately guess how I know that).

Your task is to write a program that quickly identifies squares without all the magic numbers. When you are finished, you are welcome to read or run a suggested solution, or to post your own solution or discuss the exercise in the comments below.

Advertisements

Pages: 1 2

3 Responses to “Identifying Squares, Revisited”

  1. Paul said

    This method is using quadratic residues for a number of modular values (each resulting in a set of allowed numbers). The speed of the method depends on how fast these sets can be accessed. In accesstimings.py the access times are measured for the different possiblillities (a set, a list, a tuple, a string, an int, a dict). The fastest methods appear to be a set or a list. The winner depends on what Python interpreter is used. In is_square.py the methods of @programmingpraxis (speeded up a little), using a set and a list are compared.

    # Quadratic residues mod
    d64 = set((0, 1, 4, 9, 16, 17, 25, 33, 36, 41, 49, 57))
    d63 = set((0, 1, 4, 7, 9, 16, 18, 22, 25, 28, 36, 37, 43, 46, 49, 58))
    d25 = set((0, 1, 4, 6, 9, 11, 14, 16, 19, 21, 24))
    d31 = set((0, 1, 2, 4, 5, 7, 8, 9, 10, 14, 16, 18, 19, 20, 25, 28))
    d23 = set((0, 1, 2, 3, 4, 6, 8, 9, 12, 13, 16, 18))
    d19 = set((0, 1, 4, 5, 6, 7, 9, 11, 16, 17))
    d17 = set((0, 1, 2, 4, 8, 9, 13, 15, 16))
    d11 = set((0, 1, 3, 4, 5, 9))
    
    def is_square(n):
        """using sets"""
        return (n % 64 in d64 and n % 63 in d63 and n % 25 in d25 and
                n % 31 in d31 and n % 23 in d23 and n % 19 in d19 and
                n % 17 in d17 and n % 11 in d11 and isqrt(n)**2 == n)
    
    a64 = [1 if i in d64 else 0 for i in range(64)]
    a63 = [1 if i in d63 else 0 for i in range(63)]
    a25 = [1 if i in d25 else 0 for i in range(25)]
    a31 = [1 if i in d31 else 0 for i in range(31)]
    a23 = [1 if i in d23 else 0 for i in range(23)]
    a19 = [1 if i in d19 else 0 for i in range(19)]
    a17 = [1 if i in d17 else 0 for i in range(17)]
    a11 = [1 if i in d11 else 0 for i in range(11)]
    
    def is_square2(n):
        """using lists"""
        return (a64[n % 64] and a63[n % 63] and a25[n % 25] and a31[n % 31] and
                a23[n % 23] and a19[n % 19] and a17[n % 17] and a11[n % 11] and
                isqrt(n)**2 == n)
    
    """
    CPython3.6
    isSquare 7.860324823604874 50  (ProgrammingPraxis)
    isSquare2 11.401496343810134 50 (fenderbender)
    is_square 5.880065683130212 50
    is_square2 6.948259140474661 50
    
    PyPy3.6
    isSquare 2.482781270044901 50
    isSquare2 4.008553688261706 50
    is_square 2.052393072482361 50
    is_square2 1.6373788325849894 50
    """
    
  2. r. clayton said

    Some solutions in Racket.

  3. Milbrae said

    @r. clayton… sorry to tell you, but your contribution is off topic and certainly belongs here: https://programmingpraxis.com/2019/05/03/three-homework-problems-4/
    Nice work though ;)

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out /  Change )

Google photo

You are commenting using your Google account. Log Out /  Change )

Twitter picture

You are commenting using your Twitter account. Log Out /  Change )

Facebook photo

You are commenting using your Facebook account. Log Out /  Change )

Connecting to %s

%d bloggers like this: