How well does this method of checking if an integer $N$ is a square perform?

Question

The method is based on the following observation: an integer $N=n^2$ is a square which can also be written as $N=n^2=(a+b)^2=a^2+2ab+b^2$ with:
$a=(n−1)/2$
$b=(n+1)/2$

We need to consider two cases, $N$ odd and $N$ even. The first case is simpler. Let's consider $N=n^2=9^2=81$. The value of $ab$ is given by $ab=(N−1)/4$.
Since $2ab=2(n−1)(n+1)/4=2(n^2−1)/4=2(N−1)/4=40$.
So $ab=20$ and $N=81=9^2=a^2+2∗20+b^2$.

At this point we don't know the values of a and b. However we know that their sum is $a^2+b^2=81−2∗20=41$. We also know that $a^2<ab<b^2$ since we assume $b>a$.

The important property is that $ab=2T_{k}$ with $T_{k}=k(k+1)/2$ a triangular number whose index $k$ needs to be determined.
This property is that of all numbers lying on the diagonal immediately below the main diagonal which contains all the integers square.

An important property of triangular numbers is $T_{k}^2= 1^3 + 2^3 + 3^3 + ... + k^3$, that is the square of a triangular number is the sum of the cubes from $1$ to $k$.
This property allows us to calculate the index of $T_{k}=ab/2=10$.

$T_{k}^2=10^2=100=1^3 + 2^3 + 3^3 + 4^3$. We can infer that $k=4$ and $T_{4}=4(4+1)/2=4*5/2=10$.

Knowing the numerical value of the index of $T_{k}$, we can immediately conclude that:

$a^2=4^2$ and $b^2=5^2$ and their sum is $s=4^2 + 5^2=41$.

So $N=81= 4^2 + 2*20 + 5^2= (4+5)^2=9^2$.

Basically this method allows to determine if an integer $N$ is a square without using the square root calculation and the complex algorithm of calculating the sum of two squares.

Since all numbers $ab$ on the diagonal below the main diagonal which contains the squares are always even we can use a simple test to discard odd integers that cannot be a square if the value of $ab$ is odd. For example $61$ cannot be a square since $(61-1)/4=15$. The same apply to $29$ and $85$.

For even numbers $N$, there is an additional step. For even numbers, we first calculate $M=N/4$. If $M$ is odd, then we decompose $M$ as was done before for $N$.
Example: $N=100$, then $M/4=25$. We try to decompose $M=25=a^2 + 2ab + b^2$ and we find $M=2^2 + 2*6 + 3^2=5^2$. Since $N=4M=4*5^2=2^2*5^2=10^2$.

In some cases with even numbers $N$, $N/4=M_1$ with even $M_1$. In this case, we calculate $M_2=M_1/4$. If $M_2$ is odd, we apply the method for odd potential squares.
An example is $N=144$, $N/4=M_1=36$ so we need another step. We calculate $M_2=36/4=9$. Since 9 is odd, we decompose 9 as done above. In the end, we just need to keep track of how many divisions by 4 we did so they can be incorporated in the final result.

How fast is this method?
How does it compare with other classical methods of checking for squareness of an integer?

Answer 1

Up to about 2**13, the rate at which Sum(i**3, (i, 1, k)) = (k*(k + 1)/2)**2 grows faster than 2**x. Past 13, it grows slower and takes longer and longer to determine the value of k for a given value of (ab/2)**2. This means that bisection for the bracketing square root of the number being tested will be faster; conversely, the number of Tk values needed to reach the desired value will be (too) large and costly to compute/store. e.g. try the method for n = primorial(30)**2 = 31610054640417607788145206291543662493274686990**2

The following is my implementation of the method in Python (using SymPy to compute the multiplicity of the factor of 2 (i.e. with trailing):

from bisect import bisect_left
from sympy import trailing

_tk2 = [1]
def tk2(i):
    # T(k)**2 = Sum(i**3, (i, 1, k)) = (k*(k+1)/2)**2
    k = len(_tk2)
    while i > _tk2[-1]:
        k += 1
        _tk2.append(k**3 + _tk2[-1])
    k = _tk2[bisect_left(_tk2, i)]
    if k == i:
        return k

def tk2_bin(i):
    if i == 1: return 1
    if i == 2: return
    if i == 3: return 2
    kold = k = 1
    j = 1
    while j<i:
        kold = k
        k*=2
        j = k*(k+1)//2
    if j == i:
        return k
    while k - kold > 1:
        dk = (k - kold)//2
        kk = kold + dk
        j = kk*(kk+1)//2
        if j > i:
            k = kk
        elif j < i:
            kold = kk
        else:
            return kk

def nr(n):
    x = n//2
    while 1:
        x2 = (x + n//x)//2
        if abs(x2 - x) < .5:
            return round(x)**2 == n
        x = x2

def isq(n):
    if not n:
        return True
    if n < 0:
        return False
    assert n == int(n)
    # remove powers of 4
    if not n % 2:
        t = trailing(n)
        if t % 2:
            return False
        n >>= t
    if n == 1:
        return True
    ab, r = divmod(n - 1, 4)
    if r:
        return False
    if ab % 2:
        return False
    k = tk2_bin(ab//2)  # tk2((ab//2)**2)
    return k is not None