Perfect Squares in C, No Math Library Required

Here’s a puzzle. Write this function:

bool is_square(uint64_t n);

The rules:

• No math.h. No sqrt, no pow, nothing from libm.
• No floats. Not even a quick cast to double “just to get close.”
• It has to be right for every uint64_t. All 18,446,744,073,709,551,616 of them.

Take a second and think about how you’d do it. I’ll wait.

Attempt 1: subtract odd numbers

My first instinct was a fact from grade school: the sum of the first k odd numbers is k².

1 = 1. 1 + 3 = 4. 1 + 3 + 5 = 9. 1 + 3 + 5 + 7 = 16.

There’s a nice picture behind it. To grow a k×k square into a (k+1)×(k+1) square, you add a strip along the right, a strip along the bottom, and one corner cell. That’s k + k + 1 cells, which is the next odd number.

So just run it backwards. Subtract 1, then 3, then 5, and keep going until there’s nothing left. Land on exactly zero and n was a square. Run out before you get there and it wasn’t.

bool is_square(uint64_t n)
{
    uint64_t odd = 1;
    while (n >= odd) {
        n -= odd;
        odd += 2;
    }
    return n == 0;
}

It’s correct and it’s cute. It’s also √n steps. For UINT64_MAX that’s 4.3 billion subtractions, which took 0.8 seconds on my machine. One call. For one number.

Rules satisfied, dignity not.

Attempt 2: binary search

The root of a uint64_t fits in 32 bits, so guess a root, square it, and see if you’re high or low. Halve the range and repeat.

bool is_square(uint64_t n)
{
    uint64_t lo = 0, hi = 0xFFFFFFFF;   /* not n: mid * mid would overflow */
    while (lo <= hi) {
        uint64_t mid = lo + (hi - lo) / 2;
        uint64_t sq = mid * mid;
        if (sq == n) return true;
        if (sq < n)  lo = mid + 1;
        else if (mid == 0) break;
        else hi = mid - 1;
    }
    return false;
}

Watch that hi. Start it at n and the first mid * mid overflows, so you get a wrong answer that looks perfectly reasonable.

Four billion steps become 32. About 85–91 ns per call, depending on the compiler. Most people would stop here.

But every round, binary search throws away what it learned. It picks a brand new mid, squares it from scratch, and then hits a branch that goes either way with 50/50 odds, so the CPU guesses wrong about half the time. We can do a lot better.

Attempt 3: subtract odds in bulk

Go back to the odds. Attempt 1 wasn’t wrong, it was just slow. So what if we subtracted them in blocks instead of one at a time?

Say we’ve already eaten the first R odds, so what’s left is $x = n - R^2$. The next d odds are $2R+1, 2R+3, \ldots, 2R+2d-1$. Each one is 2R plus one of the first d odds, so the whole block sums to:

So which d? Here’s the whole trick. That $d^2$ is the awkward part: squaring means multiplying, and the whole point was to stop squaring from scratch every round like binary search does. But if d is a power of two, $d = 2^j$, then $d^2 = 4^j$, and powers of four are about the friendliest numbers a CPU will ever see. In binary, every one of them is a single 1 followed by an even number of zeros:

 j    d = 2ʲ    d² = 4ʲ    in binary
 0         1          1            1
 1         2          4          100
 2         4         16        10000
 3         8         64      1000000
 4        16        256    100000000

Each step down in j just slides that 1 two places to the right. The square costs nothing. It’s a shift.

The other term gets the same treatment: $2Rd = R \cdot 2^{j+1}$ is just R shifted left. So the whole block, $R \cdot 2^{j+1} + 4^j$, is two shifts and an add. Geometrically, it’s the same picture as before, scaled up: two R×d strips plus a d×d corner, and the corner is itself a square (the first $2^j$ odds).

Now each step is a yes/no question. Walk j from high to low, and at each step ask whether the next $2^j$ odds fit in what’s left. If they do, subtract them and set bit j of the root. Picking bits greedily from the top lands on the largest R with $R^2 \le n$, and whatever’s left over is $n - R^2$. It’s zero exactly when n is a perfect square.

Try it. The dashed square is √n. Every block either fits inside it or doesn’t:

If that feels familiar, it should. It is a binary search on the root, one bit per round. The difference is it never starts over. It carries the remainder along, so each round is a compare and a subtract instead of a fresh multiply. And it only runs one round per pair of bits in n, so small numbers finish early.

The code

One more trick gets rid of the multiply entirely. Don’t store R. Store it pre-shifted, as $r = R \cdot 2^{j+1}$. Then the block sum is just r + bit, where bit $= 4^j$: that lone 1 from the table. It starts at 1ULL << 62, which is $4^{31}$, the biggest power of four a uint64_t can hold, and bit >>= 2 slides it to the next one. Moving to the next j halves r, and adds bit if that root bit was taken.

#include <stdint.h>
#include <stdbool.h>

bool is_square(uint64_t n)
{
    /* squares mod 16 are 0, 1, 4, 9: rejects 3/4 of inputs */
    if (!((0x0213u >> (n & 15)) & 1))
        return false;

    uint64_t x = n;             /* what's left: n - R^2        */
    uint64_t r = 0;             /* root so far: R * 2^(j+1)    */
    uint64_t bit = 1ULL << 62;  /* 4^j                         */

    while (bit > n)
        bit >>= 2;

    while (bit) {
        uint64_t block = r + bit;   /* sum of the next 2^j odds */
        if (x >= block) {
            x -= block;
            r = (r >> 1) + bit;
        } else {
            r >>= 1;
        }
        bit >>= 2;
    }
    return x == 0;              /* r now holds floor(sqrt(n)) */
}

Trace n = 100. The first block takes the first 8 odds (64), leaving 36. The next block of 4 odds (17 through 23, summing to 80) doesn’t fit. The one after that takes 17 + 19 = 36, leaving 0. Root is 10. Square.

No multiply, no divide, no floats, and at most 32 rounds. Zero doesn’t even need a special case: the loop never runs and x is already 0.

Here’s the part I didn’t expect. That if compiles to a pair of cmovs under both GCC and Clang, so the loop has no data-dependent branches at all. There’s nothing to mispredict. On squares, where every call runs the whole loop, it takes 13–15 ns against binary search’s 85–91, so 6× faster. On random inputs, where the filter below turns most of them away, it’s 15×. And attempt 1? That’s tens of millions of times slower.

The mod-16 line on top is a free bonus. Squares can only end in 0, 1, 4, or 9 in hex, so three out of four random inputs get turned away before the loop even starts.

Did I get it right?

“Every uint64_t” is a big promise, and I don’t trust fifteen lines of bit twiddling just because they look right, so I tested it against a reference. (The reference breaks the rules. You’ll see why at the end.)

Random 64-bit numbers are almost never squares, so random testing mostly exercises return false. The off-by-ones live right next to the squares, so that’s where I aimed:

for (int i = 0; i < 2000000; i++) {
    uint64_t k = rnd() >> 32;
    uint64_t cases[] = { k * k, k * k - 1, k * k + 1, rnd(), (uint64_t)i };
    for (int j = 0; j < 5; j++)
        mismatches += is_square(cases[j]) != sqrt_fixed(cases[j]);
}
/* plus the edges: 0 through 4, (2^32 - 1)^2, and UINT64_MAX */

10,000,007 inputs, 0 mismatches, under both GCC and Clang.

Okay, now the sane way

If you’re allowed floats, just use the hardware:

bool sqrt_fixed(uint64_t n)
{
    uint64_t r = (uint64_t)sqrt((double)n);
    if (r > 0xFFFFFFFFu) r = 0xFFFFFFFFu;
    while (r * r > n) r--;
    while (r < 0xFFFFFFFFu && (r + 1) * (r + 1) <= n) r++;
    return r * r == n;
}

sqrtsd gets you within one or two of the true root in a single instruction, and the two while loops nudge it onto the exact integer. Don’t skip them. A double only holds 53 bits, so the version you’ll find pasted everywhere, rounding and squaring back in floating point, starts calling non-squares squares once n passes 2⁵³. This is the reference I tested against, too.

It runs in about 2 ns. Our loop takes 13–15 ns once inputs get past the filter. It isn’t close.

Every number in this post comes from one benchmark file: best of 7 passes over 4 million inputs, pinned to one core of a Ryzen 9 7950X. Here’s the whole session:

gcc -O2 -march=native perfect-squares-bench.c -lm -o bench-gcc
taskset -c 2 ./bench-gcc > gcc.json
jq -c '{checked, mismatches, naiveMaxSeconds}' gcc.json
{"checked":10000007,"mismatches":0,"naiveMaxSeconds":0.798}
jq -r '.bench[] | select(.mix == "squares") | "\(.fn) \(.ns)"' gcc.json | column -t
sqrt_all_double     3.563
sqrt_fixed          2.042
bsearch_sq          85.394
is_square           13.333
is_square_nofilter  18.440
 
clang -O2 -march=native perfect-squares-bench.c -lm -o bench-clang
taskset -c 2 ./bench-clang > clang.json
jq -c '{checked, mismatches, naiveMaxSeconds}' clang.json
{"checked":10000007,"mismatches":0,"naiveMaxSeconds":0.795}
jq -r '.bench[] | select(.mix == "squares") | "\(.fn) \(.ns)"' clang.json | column -t
sqrt_all_double     1.750
sqrt_fixed          2.128
bsearch_sq          90.516
is_square           15.430
is_square_nofilter  21.800

And the hardware doesn’t care how big n is. Our loop does, one round per pair of bits:

So if you have an FPU, use it. The block loop is for when you don’t: microcontrollers, freestanding code, kernels, anywhere you want an answer that’s exact by construction and not by careful reasoning about float rounding.

And honestly? It’s for the fun of watching a fact from grade school turn into fifteen lines that beat binary search by 6×.