Determine if a number is a perfect square

Verifying whether a number is a perfect square becomes a trivial task with the availability of a square root function; squaring the floored square root should yield the original number:

\[\sqrt{n}\in\mathbb{N}\Leftrightarrow\lfloor\sqrt{n}\rfloor^2 = n\]

However, on many systems, a calculation of a square root is costly, so the easiest way is an \(O(\sqrt{n})\) iteration like

function isPerfectSquare(n) {

  for (let i = 0; i * i <= n; i++) { 

    if (i * i === n) { 
      return true; 
    } 
  }
  return false;
}

Babylonian Method

The intuitive algorithm is already quite fast in practice, but we can utilize the Babylonian Method for square root finding to check if a number is a perfect square. The idea is that we can express our number \(n\) as \(m\cdot\frac{n}{m}\) and on the other hand also as \(\sqrt{n}\cdot\sqrt{n}\), from which follows that if \(m<\sqrt{n}\), then \(\frac{n}{m}>\sqrt{n}\) and if \(m>\sqrt{n}\), then \(\frac{n}{m}<\sqrt{n}\). This way we have one upper estimate that is too big and one lower estimate that is too small for our real \(\sqrt{n}\), so we simply take the average of those and form an iterative process:

\[m_{k+1} = \frac{m_k + n / m_k}{2}\]

It is obvious that by halving on each iteration, the error towards the real square root gets also halved as well, forming an \(O(\log(n))\) algorithm:

function isPerfectSquare(n) {

  let m = n;
  while (m * m > n) {
    m = (m + n / m) >> 1;
  }
  return m * m === n;
}

Please note that the division is actually an integer division, but we use JavaScript's implicit cast via the right shift to avoid an additional call. We can go even further, because many numbers follow patterns that can be exploited to make the best-case performance \(O(1)\) by excluding whole ranges of numbers. The \(O(\log{n})\) step is then only needed as a worst-case fallback.

Check the last digit

It is well known that in base 10, a perfect square can have the least significant digit to be only 0, 1, 4, 5, 6 or 9, which forbids 2, 3, 7 and 8 under modulo 10.

The proof is pretty straightforward. Let's say the least significant digit of our number \(n\in\mathbb{Z}\) is \(a\) and the remaining digits are called \(b\), we can then say \(n=10b+a\).

Now when squared, \(n^2=(10b+a)^2=100b^2+20ab+a^2\). It is clear that the first part \(100b^2+2ab\cdot 10\) does not affect the least significant digit, but only \(a\) does.

Listing all possible squares results in

\[\{a^2 : a\in[0, 9]\} = \{0,1,4,9,16,25,36,49,64,81\}\]

which ultimately shows that the only possible least significant digits of a perfect square are 0,1,4,5,6 and 9.

Check digital root

The digital root of a number \(n\) is the sum of all digits of that number. If the sum has more than one digit, the procedure is repeated until only one digit remains. Instead of implementing a recursive solution, we can use a trick. When we work in base \(b := 10\), it is clear that for any \(k\) the following congruence holds: \(b^k\equiv 1\pmod{b-1}\). As the number \(n\) consists of its digits \(n_i\) times \(10^i\), the same idea lets us reduce the sum modulo 9:

\[ \begin{array}{rl} \text{dr}(n) &= \sum\limits_{i=0}10^i\cdot n_i \pmod{9}\\ &= \sum\limits_{i=0}1\cdot n_i \pmod{9}\\ &= \sum\limits_{i=0}n_i \pmod{9}\\ &= \begin{cases} n\pmod{9} & n\not\equiv 0\pmod{9}\\ 9 & n\equiv 0\pmod{9} \end{cases}\\ &= 1 + (n-1\pmod{9}) \end{array} \]

Using the digital root function, we can now check the behavior when getting a squared number:

\[\begin{array}{rl} \text{dr}(n^2) &= 1 + (n^2-1\pmod{9})\\ &= 1 + (n^2 + 8\pmod{9})\\ &= \underbrace{1 + \big(\underbrace{(n^2\pmod{9}}_{\in\{0, 1, 4, 7\}}) +8\pmod{9}\big)}_{\in\{9, 1, 4, 7\}}\\ \end{array}\]

That is interesting: a perfect square has only four possible digital roots, namely 1, 4, 7 and 9. Equivalently, it has only four possible residues modulo 9: 0, 1, 4 and 7.

Improve findings with other moduli

An interesting question that arises from these two findings is whether we can find a better modulus that covers a larger range of numbers, in other words a better modulus \(m\) that increases the coverage and reduces the number of checks.

\[n^2 \equiv k \pmod{m}\]

Since the possible digits repeat quickly, a simple brute force enumeration approach can give us possible candidates:

const LIMIT = 500;

const res = [];
for (let m = 2; m < LIMIT; m++) {
  const residues = new Set();
  for (let i = 0; i < m; i++)
    residues.add((i * i) % m);

  if (residues.size < 10)
    res.push({module: m, coverage: (m - residues.size) / m, candidates: [...residues]});
}

Sorting the result by coverage descending leads to the following table:

The most interesting result is modulus \(16\), which has only four possible residue classes and covers \(\frac{3}{4}\) of all numbers, giving a best-case performance of \(O(1)\).

A much better implementation to check for perfect squares is then

function isPerfectSquare(n) {
  let m = n & 15;

  if (m !== 0 && m !== 1 && m !== 4 && m !== 9) {
    return false;
  }

  for (m = n; m * m > n; ) {
    m = (m + n / m) >> 1;
  }
  return m * m === n;
}

Interestingly, the two previous ideas of checking the last digit and the digital root, which motivated the whole approach, do not cover very much and require many more individual checks.

Combining Multiple Moduli

The good thing about this method is that we can add additional checks, for example for \(m=7\) and \(m=9\), to further increase the \(O(1)\) test coverage. Each additional modular check removes numbers that cannot be perfect squares before the Babylonian method is needed.

For a modulus \(m\), let

\[ QR(m) = \{x^2 \bmod m : x \in \mathbb{Z}\} \]

be the set of quadratic residue classes modulo \(m\). Then a number \(n\) can only be a perfect square if

\[ n \bmod m \in QR(m). \]

If we use two moduli \(a\) and \(b\), then both conditions must hold:

\[ n \bmod a \in QR(a) \quad\text{and}\quad n \bmod b \in QR(b). \]

The combined test still runs in \(O(1)\). The interesting question is therefore which pair of moduli rejects the largest number of non-square candidates.

If \(a\) and \(b\) are coprime, the Chinese Remainder Theorem tells us that the two modular conditions are independent over the combined period \(a \cdot b\). In that case, the fraction of numbers that pass both tests is simply

\[ \frac{|QR(a)|}{a} \cdot \frac{|QR(b)|}{b}. \]

Therefore the coverage is

\[ 1 - \frac{|QR(a)|}{a} \cdot \frac{|QR(b)|}{b}. \]

For example, for \(m=48\), we found

\[ QR(48)=\{0,1,4,9,16,25,33,36\} \]

so only \( \frac{8}{48} \) of all numbers pass this filter. Combining it with \(m=7\), where

\[ QR(7)=\{0,1,2,4\} \]

gives

\[ \frac{8}{48}\cdot\frac{4}{7} = \frac{2}{21}. \]

So the coverage becomes

\[ 1-\frac{2}{21}=\frac{19}{21}\approx 90.48%. \]

However, \(m=48\) is not necessarily the best first choice when two moduli are allowed.

We can again find the best pair by brute force. For non-coprime moduli, the conditions are no longer independent, so we count over the combined period \(\operatorname{lcm}(a,b)\):

const gcd = (a, b) => { while (b) [a, b] = [b, a % b]; return a; };
const lcm = (a, b) => (a / gcd(a, b)) * b;

const LIMIT = 500;
const MAX_RESIDUE_CLASSES = 32;
const residueCache = new Map();

for (let m = 2; m < LIMIT; m++) {
  const residues = new Set();
  for (let i = 0; i < m; i++) {
    residues.add((i * i) % m);
  }
  residueCache.set(m, residues);
}

const moduli = [];
for (let m = 2; m < LIMIT; m++) {
  // Keep only moduli whose residue sets stay compact enough for bit-flag variants.
  if (residueCache.get(m).size <= MAX_RESIDUE_CLASSES) {
    moduli.push(m);
  }
}

const res = [];
for (let i = 0; i < moduli.length; i++) {
  for (let j = i + 1; j < moduli.length; j++) {
    const a = moduli[i], b = moduli[j];
    const A = residueCache.get(a), B = residueCache.get(b);
    const period = lcm(a, b);
    let passed = 0;

    for (let n = 0; n < period; n++) {
      if (A.has(n % a) && B.has(n % b)) {
        passed++;
      }
    }

    res.push({
      modulusA: a,
      modulusB: b,
      sizeA: A.size,
      sizeB: B.size,
      bitflagFriendly: a <= 32 && b <= 32,
      coverage: (period - passed) / period
    });
  }
}

Sorting by coverage yields:

The strongest pairs by pure coverage use larger moduli. The highlighted row \((a=21, b=32)\) is the best pair that uses only few residues and where \(a\) and \(b\) are \(\leq 32\). The resulting implementation is then:

function isPerfectSquare(n) {
  let m = n & 31;

  if (
    m !== 0 &&
    m !== 1 &&
    m !== 4 &&
    m !== 9 &&
    m !== 16 &&
    m !== 17 &&
    m !== 25
  ) {
    return false;
  }

  m = n % 21;

  if (
    m !== 0 &&
    m !== 1 &&
    m !== 4 &&
    m !== 7 &&
    m !== 9 &&
    m !== 15 &&
    m !== 16 &&
    m !== 18
  ) {
    return false;
  }

  for (m = n; m * m > n; ) {
    m = (m + n / m) >> 1;
  }

  return m * m === n;
}

Binary Flag Testing

If we restrict the residues to be \(\leq 32\), the residue checks can be written more compactly by storing the allowed residue classes as bit flags. If bit \(k\) is set, then \(k\) is an allowed residue class.

For \(m=32\) and \(m=21\), this gives

const QR32 = 0x2030213;
const QR21 = 0x58293;

and the implementation becomes this compact function that rejects \(\frac{11}{12}\approx 91.67\%\) of all numbers before the Babylonian method is required:

function isPerfectSquare(n) {

  const QR32 = 0x2030213;
  const QR21 = 0x58293;

  if (((QR32 >>> (n & 31)) & 1) === 0) { // 78.1% coverage
    return false;
  }

  if (((QR21 >>> (n % 21)) & 1) === 0) { // 61.9% coverage
    return false;
  }

  let m = n;
  for (; m * m > n; ) {
    m = (m + n / m) >> 1;
  }
  return m * m === n;
}