(phixonline)-->
--requires(64) -- (decided to limit 32-bit explicitly instead)
constant MxN = power(2,iff(machine_bits()=32?53:63)),
m = {1, 3, 5, 7, 11}
function squfof(atom N)
-- square form factorization
integer h, a=0, b, c, u=0, v, w, rN, q, r, t
if remainder(N,2)==0 then return 2 end if
h = floor(sqrt(N) + 0.5)
if h*h==N then return h end if
for k=1 to length(m) do
integer mk = m[k]
if mk>1 and remainder(N,mk)==0 then return mk end if
//check overflow m * N
if N>MxN/mk then exit end if
atom mN = N*mk
r = floor(sqrt(mN))
if r*r>mN then r -= 1 end if
rN = r
//principal form
{b,a} = {r,1}
h = floor((rN+b)/a)*a-b
c = floor((mN-h*h)/a)
for i=2 to floor(sqrt(2*r)) * 4-1 do
//search principal cycle
{a,c,t} = {c,a,b}
q = floor((rN+b)/a)
b = q*a-b
c += q*(t-b)
if remainder(i,2)==0 then
r = floor(sqrt(c)+0.5)
if r*r==c then
//square form found
//inverse square root
q = floor((rN-b)/r)
v = q*r+b
w = floor((mN-v*v)/r)
//search ambiguous cycle
u = r
while true do
{u,w,r} = {w,u,v}
q = floor((rN+v)/u)
v = q*u-v
if v==r then exit end if
w += q*(r-v)
end while
//symmetry point
h = gcd(N,u)
if h!=1 then return h end if
end if
end if
end for
end for
return 1
end function
constant tests = {2501, 12851, 13289, 75301