Task.
Daniel Shanks's Square Form Factorization (SquFoF).

Invented around 1975, ‘On a 32-bit computer, SquFoF is the clear champion factoring algorithm
for numbers between 1010 and 1018, and will likely remain so.’



An integral binary quadratic form is a polynomial
f(x,y) = ax2 + bxy + cy2
with integer coefficients and discriminant D = b2 – 4ac.
For each positive discriminant there are multiple forms (a, b, c).

The next form in a periodic sequence (cycle) of adjacent forms is found by applying a reduction operator
rho, essentially a variant of Euclid's algorithm for finding the continued fraction of a square root.
Using floor(√N), rho constructs a principal form
(1, b, c) with D = 4N.

SquFoF is based on the existence of cycles containing ambiguous forms, with the property that a divides b.
They come in pairs of associated forms (a, b, c) and (c, b, a) called symmetry points.
If an ambiguous form is found (there is one for each divisor of D), write the discriminant as
(ak)2 – 4ac = a(a·k2 – 4c) = 4N
and (if a is not equal to 1 or 2) N is split.

Shanks used square forms to jump to a random ambiguous cycle. Fact: if any form in an ambiguous cycle
is squared, that square form will always land in the principal cycle. Conversely, the square root of any
form in the principal cycle lies in an ambiguous cycle. (Possibly the principal cycle itself).

A square form is easy to find: the last coefficient c is a perfect square. This happens about once
every ∜N-th cycle step and for even indices only. Let rho compute the inverse square root form and track
the ambiguous cycle backward until the symmetry point is reached. (Taking the inverse reverses the cycle).
Then a or a/2 divides D and therefore N.

To avoid trivial factorizations, Shanks created a list (queue) to hold small coefficients appearing
early in the principal cycle, that may be roots of square forms found later on. If these forms are skipped,
no roots land in the principal cycle itself and cases a = 1 or a = 2 do not happen.

Sometimes the cycle length is too short to find a proper square form. This is fixed by running five instances
of SquFoF in parallel, with input N and 3, 5, 7, 11 times N; the discriminants then will have different periods.
If N is prime or the cube of a prime, there are improper squares only and the program will duly report failure.

Reference.
[1] A detailed analysis of SquFoF (2007)



function square_form_factor(n::T)::T where T <: Integer
    multiplier = T.([1, 3, 5, 7, 11, 3*5, 3*7, 3*11, 5*7, 5*11, 7*11, 3*5*7, 3*5*11, 3*7*11, 5*7*11, 3*5*7*11])
    s = T(round(sqrt(n)))
    s * s == n && return s
    for k in multiplier
        T != BigInt && n > typemax(T) ÷ k && break
        d = k * n
        p0 = pprev = p = isqrt(d)
        qprev = one(T)
        Q = d - p0 * p0
        l = T(floor(2 * sqrt(2 * s)))
        B, i = 3 * l, 2
        while i < B
            b = (p0 + p) ÷ Q
            p = b * Q - p
            q = Q
            Q = qprev + b * (pprev - p)
            r = T(round(sqrt(Q)))
            iseven(i) && r * r == Q && break
            qprev, pprev = q, p
            i += 1
        end
        i >= B && continue
        b = (p0 - p) ÷ r
        pprev = p = b * r + p
        qprev = r
        Q = (d - pprev * pprev) ÷ qprev
        i = 0
        while true
            b = (p0 + p) ÷ Q
            pprev = p
            p = b * Q - p
            q = Q
            Q = qprev + b * (pprev - p)
            qprev = q
            i += 1
            p == pprev && break
        end
        r = gcd(n, qprev)
        r != 1 && r != n && return r
    end
    return zero(T)
end
 
println("Integer              Factor   Quotient\n", "-"^45)
@time for n in Int128.([
    2501, 12851, 13289, 75301, 120787, 967009, 997417, 7091569, 13290059, 42854447, 223553581,
    2027651281, 11111111111, 100895598169, 1002742628021, 60012462237239, 287129523414791,
    9007199254740931, 11111111111111111, 314159265358979323, 384307168202281507, 419244183493398773,
    658812288346769681, 922337203685477563, 1000000000000000127, 1152921505680588799,
    1537228672809128917, 4611686018427387877])
    print(rpad(n, 22))
    factr = square_form_factor(n)
    print(rpad(factr, 10))
    println(factr == 0 ? "fail" : n ÷ factr)
end
 