I'm working on modernizing Rosetta Code's infrastructure. Starting with communications. Please accept this time-limited open invite to RC's Slack.. --Michael Mol (talk) 20:59, 30 May 2020 (UTC)

Square form factorization

Square form factorization
You are encouraged to solve this task according to the task description, using any language you may know.

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 = b24ac. 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)24ac = a(a·k24c) = 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.

 A detailed analysis of SquFoF (2007)

C

Based on Wp

From the Wikipedia entry for Shanks's square forms factorization []

#include <math.h>
#include <stdio.h>

#define nelems(x) (sizeof(x) / sizeof((x)))

const unsigned long multiplier[] = {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};

unsigned long long gcd(unsigned long long a, unsigned long long b)
{
while (b != 0)
{
a %= b;
a ^= b;
b ^= a;
a ^= b;
}

return a;
}

unsigned long long SQUFOF( unsigned long long N )
{
unsigned long long D, Po, P, Pprev, Q, Qprev, q, b, r, s;
unsigned long L, B, i;
s = (unsigned long long)(sqrtl(N)+0.5);
if (s*s == N) return s;
for (int k = 0; k < nelems(multiplier) && N <= 0xffffffffffffffff/multiplier[k]; k++) {
D = multiplier[k]*N;
Po = Pprev = P = sqrtl(D);
Qprev = 1;
Q = D - Po*Po;
L = 2 * sqrtl( 2*s );
B = 3 * L;
for (i = 2 ; i < B ; i++) {
b = (unsigned long long)((Po + P)/Q);
P = b*Q - P;
q = Q;
Q = Qprev + b*(Pprev - P);
r = (unsigned long long)(sqrtl(Q)+0.5);
if (!(i & 1) && r*r == Q) break;
Qprev = q;
Pprev = P;
};
if (i >= B) continue;
b = (unsigned long long)((Po - P)/r);
Pprev = P = b*r + P;
Qprev = r;
Q = (D - Pprev*Pprev)/Qprev;
i = 0;
do {
b = (unsigned long long)((Po + P)/Q);
Pprev = P;
P = b*Q - P;
q = Q;
Q = Qprev + b*(Pprev - P);
Qprev = q;
i++;
} while (P != Pprev);
r = gcd(N, Qprev);
if (r != 1 && r != N) return r;
}
return 0;
}

int main(int argc, char *argv[]) {
int i;
const unsigned long long data[] = {
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};

for(int i = 0; i < nelems(data); i++) {
unsigned long long example, factor, quotient;
example = data[i];
factor = SQUFOF(example);
if(factor == 0) {
printf("%llu was not factored.\n", example);
}
else {
quotient = example / factor;
printf("Integer %llu has factors %llu and %llu\n",
example, factor, quotient);
}
}
}

Output:
Integer 2501 has factors 61 and 41
Integer 12851 has factors 71 and 181
Integer 13289 has factors 137 and 97
Integer 75301 has factors 293 and 257
Integer 120787 has factors 43 and 2809
Integer 967009 has factors 1609 and 601
Integer 997417 has factors 257 and 3881
Integer 7091569 has factors 2663 and 2663
Integer 13290059 has factors 3119 and 4261
Integer 42854447 has factors 9689 and 4423
Integer 223553581 has factors 11213 and 19937
Integer 2027651281 has factors 46061 and 44021
Integer 11111111111 has factors 21649 and 513239
Integer 100895598169 has factors 112303 and 898423
1002742628021 was not factored.
Integer 60012462237239 has factors 6862753 and 8744663
Integer 287129523414791 has factors 6059887 and 47381993
Integer 9007199254740931 has factors 10624181 and 847801751
Integer 11111111111111111 has factors 2071723 and 5363222357
Integer 314159265358979323 has factors 317213509 and 990371647
Integer 384307168202281507 has factors 415718707 and 924440401
Integer 419244183493398773 has factors 48009977 and 8732438749
Integer 658812288346769681 has factors 62222119 and 10588072199
Integer 922337203685477563 has factors 110075821 and 8379108103
1000000000000000127 was not factored.
Integer 1152921505680588799 has factors 139001459 and 8294312261
1537228672809128917 was not factored.
Integer 4611686018427387877 has factors 343242169 and 13435662733

Classical heuristic

See Discussion.

//SquFoF: minimalistic version without queue.
//Classical heuristic. Tested: tcc 0.9.27
#include <math.h>
#include <stdio.h>

//input maximum
#define MxN ((unsigned long long) 1 << 62)

//reduce indefinite form
#define rho(a, b, c) { \
t = c; c = a; a = t; t = b; \
q = (rN + b) / a; \
b = q * a - b; \
c += q * (t - b); }

//initialize
#define rhoin(a, b, c) { \
rho(a, b, c) h = b; \
c = (mN - h * h) / a; }

#define gcd(a, b) while (b) { \
t = a % b; a = b; b = t; }

//multipliers
const unsigned long m[] = {1, 3, 5, 7, 11, 0};

//square form factorization
unsigned long squfof( unsigned long long N ) {
unsigned long a, b, c, u, v, w, rN, q, t, r;
unsigned long long mN, h;
int i, ix, k = 0;

if ((N & 1)==0) return 2;

h = floor(sqrt(N)+ 0.5);
if (h * h == N) return h;

while (m[k]) {
if (k && N % m[k]==0) return m[k];
//check overflow m * N
if (N > MxN / m[k]) break;
mN = N * m[k++];

r = floor(sqrt(mN));
h = r; //float64 fix
if (h * h > mN) r -= 1;
rN = r;

//principal form
b = r; c = 1;
rhoin(a, b, c)

//iteration bound
ix = floor(sqrt(2*r)) * 4;

//search principal cycle
for (i = 2; i < ix; i += 2) {
rho(a, b, c)
//even step

r = floor(sqrt(c)+ 0.5);
if (r * r == c) {
//square form found

//inverse square root
v = -b; w = r;
rhoin(u, v, w)

//search ambiguous cycle
do { r = v;
rho(u, v, w)
} while (v != r);
//symmetry point

h = N; gcd(h, u)
if (h != 1) return h;
}
rho(a, b, c)
//odd step
}
}
return 1;
}

void main(void) {
const unsigned long long data[] = {
2501,
12851,
13289,
75301,
120787,
967009,
997417,
7091569,

5214317,
20834839,
23515517,
33409583,
44524219,

13290059,
223553581,
2027651281,
11111111111,
100895598169,
1002742628021,
60012462237239,
287129523414791,
9007199254740931,
11111111111111111,
314159265358979323,
384307168202281507,
419244183493398773,
658812288346769681,
922337203685477563,
1000000000000000127,
1152921505680588799,
1537228672809128917,
4611686018427387877,
0};

unsigned long long N, f;
int i = 0;

while (1) {
N = data[i++];
//scanf("%llu", &N);
if (N < 2) break;

printf("N = %llu\n", N);

f = squfof(N);
if (N % f) f = 1;

if (f == 1) printf("fail\n\n");
else printf("f = %llu N/f = %llu\n\n", f, N/f);
}
}
Showing problem cases only:
...
N = 5214317
f = 73  N/f = 71429

N = 20834839
f = 3361  N/f = 6199

N = 23515517
f = 53  N/f = 443689

N = 33409583
f = 991  N/f = 33713

N = 44524219
f = 593  N/f = 75083
...
N = 1000000000000000127
f = 111756107  N/f = 8948056861
...
N = 1537228672809128917
f = 26675843  N/f = 57626245319
...

FreeBASIC

' ***********************************************
'subject: Shanks's square form factorization:
' ambiguous forms of discriminant 4N
' give factors of N.
'tested : FreeBasic 1.08.1

'------------------------------------------------
const MxN = culngint(1) shl 62
'input maximum

const qx = (1 shl 5) - 1
'queue size

type arg
'squfof arguments
as ulong m, f
as integer vb
end type

type bqf
declare sub rho ()
'reduce indefinite form
declare function issq (byref r as ulong) as integer
'return -1 if c is square, set r:= sqrt(c)
declare sub qform (byref g as string, byval t as integer)
'print binary quadratic form #t (a, 2b, c)

as ulong rN, a, b, c
as integer vb
end type

type queue
declare sub enq (byref P as bqf)
'enqueue P.c, P.b if appropriate
declare function pro (byref P as bqf, byval r as ulong) as integer
'return -1 if a proper square form is found

as ulong a(qx), L, m
as integer k, t
end type

'global variables
dim shared N as ulongint
'the number to split

dim shared flag as integer

dim shared as ubyte q1024(1023), q3465(3464)

'------------------------------------------------
sub bqf.rho ()
dim as ulong q, t
swap a, c
'residue
q = culng(rN + b) \ a
t = b: b = q * a - b
'pseudo-square
c += q * (t - b)
end sub

'initialize form
#macro rhoin(F)
F.rho : h = F.b
F.c = (mN - h * h) \ F.a
#endmacro

function bqf.issq (byref r as ulong) as integer
if q1024(c and 1023) andalso q3465(c mod 3465) then
'98.6% non-squares filtered
r = culng(sqr(c))
if r * r = c then return -1
end if
issq = 0
end function

sub bqf.qform (byref g as string, byval t as integer)
if vb = 0 then exit sub
dim as longint u = a, v = b, w = c
if t and 1 then
w = -w
else
u = -u
end if
v shl= 1
print g;str(t);" = (";u;",";v;",";w;")"
end sub

'------------------------------------------------
#macro red(r, a)
r = iif(a and 1, a, a shr 1)
if m > 2 then
r = iif(r mod m, r, r \ m)
end if
#endmacro

sub queue.enq (byref P as bqf)
dim s as ulong
red(s, P.c)
if s < L then
'circular queue
k = (k + 2) and qx
if k > t then t = k
'enqueue P.b, P.c
a(k) = P.b mod s
a(k + 1) = s
end if
end sub

function queue.pro (byref P as bqf, byval r as ulong) as integer
dim as integer i, sw
'skip improper square forms
for i = 0 to t step 2
sw = (P.b - a(i)) mod r = 0
sw and= a(i + 1) = r
if sw then return 0
next i
pro = -1
end function

'------------------------------------------------
sub squfof (byval ap as any ptr)
dim as arg ptr rp = cptr(arg ptr, ap)
dim as ulong L2, m, r, t, f = 1
dim as integer ix, i, j
dim as ulongint mN, h
'principal and ambiguous cycles
dim as bqf P, A
dim Q as queue

if (N and 1) = 0 then
rp->f = 2 ' even N
flag =-1: exit sub
end if

h = culngint(sqr(N))
if h * h = N then
'N is square
rp->f = culng(h)
flag =-1: exit sub
end if

rp->f = 1
'multiplier
m = rp->m
if m > 1 then
if (N mod m) = 0 then
rp->f = m ' m | N
flag =-1: exit sub
end if

'check overflow m * N
if N > (MxN \ m) then exit sub
end if
mN = N * m

r = int(sqr(mN))
'float64 fix
if culngint(r) * r > mN then r -= 1
P.rN = r
A.rN = r

P.vb = rp->vb
A.vb = rp->vb
'verbosity switch
if P.vb then print "r = "; r

Q.k = -2: Q.t = -1: Q.m = m
'Queue entry bounds
Q.L = int(sqr(r * 2))
L2 = Q.L * m shl 1

'principal form
P.b = r: P.c = 1
rhoin(P)
P.qform("P", 1)

ix = Q.L shl 2
for i = 2 to ix
'search principal cycle

if P.c < L2 then Q.enq(P)

P.rho
if (i and 1) = 0 andalso P.issq(r) then
'square form found

if Q.pro(P, r) then

P.qform("P", i)
'inverse square root
A.b =-P.b: A.c = r
rhoin(A): j = 1
A.qform("A", j)

do
'search ambiguous cycle
t = A.b
A.rho: j += 1

if A.b = t then
'symmetry point
A.qform("A", j)
red(f, A.a)
if f = 1 then exit do

flag = -1
'factor found
end if
loop until flag

end if ' proper square
end if ' square form

if flag then exit for
next i

rp->f = f
end sub

'------------------------------------------------
data 2501
data 12851
data 13289
data 75301
data 120787
data 967009
data 997417
data 7091569
data 13290059
data 23515517
data 42854447
data 223553581
data 2027651281
data 11111111111
data 100895598169
data 1002742628021
data 60012462237239
data 287129523414791
data 9007199254740931
data 11111111111111111
data 314159265358979323
data 384307168202281507
data 419244183493398773
data 658812288346769681
data 922337203685477563
data 1000000000000000127
data 1152921505680588799
data 1537228672809128917
data 4611686018427387877
data 0

'main
'------------------------------------------------
const tx = 4
dim as double tim = timer
dim h(4) as any ptr
dim a(4) as arg
dim as ulongint f
dim as integer s, t

width 64, 30
cls

for t = 0 to 1540
s = t * t
q1024(s and 1023) =-1
q3465(s mod 3465) =-1
next t

a(0).vb = 0
'set one verbosity switch only

a(0).m = 1
'multipliers
a(1).m = 3
a(2).m = 5
a(3).m = 7
a(4).m = 11

do
print

loop until N < MxN
if N < 2 then exit do

print "N = "; N

flag = 0

for t = 1 to tx + 1 step 2
if t < tx then
end if

squfof(@a(t - 1))
f = a(t - 1).f

if t < tx then
if f = 1 then f = a(t).f
end if

if f > 1 then exit for
next t

'assert
if N mod f then f = 1

if f = 1 then
print "fail"
else
print "f = ";f;" N/f = ";N \ f
end if
loop

print "total time:"; csng(timer - tim); " s"
end
Examples:
N = 2501
f = 61  N/f = 41

N = 12851
f = 71  N/f = 181

N = 13289
f = 97  N/f = 137

N = 75301
f = 293  N/f = 257

N = 120787
f = 43  N/f = 2809

N = 967009
f = 1609  N/f = 601

N = 997417
f = 257  N/f = 3881

N = 7091569
f = 2663  N/f = 2663

N = 13290059
f = 3119  N/f = 4261

N = 23515517
f = 53  N/f = 443689

N = 42854447
f = 4423  N/f = 9689

N = 223553581
f = 11213  N/f = 19937

N = 2027651281
f = 44021  N/f = 46061

N = 11111111111
f = 21649  N/f = 513239

N = 100895598169
f = 112303  N/f = 898423

N = 1002742628021
fail

N = 60012462237239
f = 6862753  N/f = 8744663

N = 287129523414791
f = 6059887  N/f = 47381993

N = 9007199254740931
f = 10624181  N/f = 847801751

N = 11111111111111111
f = 2071723  N/f = 5363222357

N = 314159265358979323
f = 317213509  N/f = 990371647

N = 384307168202281507
f = 415718707  N/f = 924440401

N = 419244183493398773
f = 48009977  N/f = 8732438749

N = 658812288346769681
f = 62222119  N/f = 10588072199

N = 922337203685477563
f = 110075821  N/f = 8379108103

N = 1000000000000000127
f = 111756107  N/f = 8948056861

N = 1152921505680588799
f = 139001459  N/f = 8294312261

N = 1537228672809128917
f = 26675843  N/f = 57626245319

N = 4611686018427387877
f = 343242169  N/f = 13435662733

total time: 0.0170462 s

Go

Yet another solution based on the Wikipedia C code.

Rather than juggling with big.Int, I've just allowed the two 'awkward' cases to fail.

package main

import (
"fmt"
"math"
)

func isqrt(x uint64) uint64 {
x0 := x >> 1
x1 := (x0 + x/x0) >> 1
for x1 < x0 {
x0 = x1
x1 = (x0 + x/x0) >> 1
}
return x0
}

func gcd(x, y uint64) uint64 {
for y != 0 {
x, y = y, x%y
}
return x
}

var multiplier = []uint64{
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,
}

func squfof(N uint64) uint64 {
s := uint64(math.Sqrt(float64(N)) + 0.5)
if s*s == N {
return s
}
for k := 0; k < len(multiplier) && N <= math.MaxUint64/multiplier[k]; k++ {
D := multiplier[k] * N
P := isqrt(D)
Pprev := P
Po := Pprev
Qprev := uint64(1)
Q := D - Po*Po
L := uint32(isqrt(8 * s))
B := 3 * L
i := uint32(2)
var b, q, r uint64
for ; i < B; i++ {
b = uint64((Po + P) / Q)
P = b*Q - P
q = Q
Q = Qprev + b*(Pprev-P)
r = uint64(math.Sqrt(float64(Q)) + 0.5)
if (i&1) == 0 && r*r == Q {
break
}
Qprev = q
Pprev = P
}
if i >= B {
continue
}
b = uint64((Po - P) / r)
P = b*r + P
Pprev = P
Qprev = r
Q = (D - Pprev*Pprev) / Qprev
i = 0
for {
b = uint64((Po + P) / Q)
Pprev = P
P = b*Q - P
q = Q
Q = Qprev + b*(Pprev-P)
Qprev = q
i++
if P == Pprev {
break
}
}
r = gcd(N, Qprev)
if r != 1 && r != N {
return r
}
}
return 0
}

func main() {
examples := []uint64{
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,
}
fmt.Println("Integer Factor Quotient")
fmt.Println("------------------------------------------")
for _, N := range examples {
fact := squfof(N)
quot := "fail"
if fact > 0 {
quot = fmt.Sprintf("%d", N/fact)
}
fmt.Printf("%-20d %-10d %s\n", N, fact, quot)
}
}
Output:
Integer              Factor     Quotient
------------------------------------------
2501                 61         41
12851                71         181
13289                137        97
75301                293        257
120787               43         2809
967009               1609       601
997417               257        3881
7091569              2663       2663
13290059             3119       4261
42854447             9689       4423
223553581            11213      19937
2027651281           46061      44021
11111111111          21649      513239
100895598169         112303     898423
1002742628021        0          fail
60012462237239       6862753    8744663
287129523414791      6059887    47381993
9007199254740931     10624181   847801751
11111111111111111    2071723    5363222357
314159265358979323   317213509  990371647
384307168202281507   415718707  924440401
419244183493398773   48009977   8732438749
658812288346769681   62222119   10588072199
922337203685477563   110075821  8379108103
1000000000000000127  0          fail
1152921505680588799  139001459  8294312261
1537228672809128917  0          fail
4611686018427387877  343242169  13435662733

jq

Works with gojq, the Go implementation of jq

gojq has support for unbounded-precision integer arithmetic and accordingly the output shown below is from a run thereof; the C implementation of jq produces correct results up to and including [287129523414791,6059887,47381993].

Preliminaries

def gcd(a; b):
# subfunction expects [a,b] as input
# i.e. a ~ . and b ~ .
def rgcd: if . == 0 then .
else [., . % .] | rgcd
end;
[a,b] | rgcd;

# for infinite precision integer-arithmetic
def idivide(\$p; \$q): (\$p - (\$p % \$q)) / \$q ;
def idivide(\$q): (. - (. % \$q)) / \$q ;

def isqrt:
def irt:
. as \$x
| 1 | until(. > \$x; . * 4) as \$q
| {\$q, \$x, r: 0}
| until( .q <= 1;
.q |= idivide(4)
| .t = .x - .r - .q
| .r |= idivide(2)
| if .t >= 0
then .x = .t
| .r += .q
else .
end)
| .r ;
if type == "number" and (isinfinite|not) and (isnan|not) and . >= 0
then irt
else "isqrt requires a non-negative integer for accuracy" | error
end ;

def multipliers:
[
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
];

# input should be a number
def squfof:
def toi : floor | tostring | tonumber;
. as \$N
| ((\$N|sqrt + 0.5)|toi) as \$s
| if (\$s*\$s == \$N) then \$s
else label \$out
| {}
| multipliers[] as \$multiplier
| (\$N * \$multiplier) as \$D
| .P = (\$D|isqrt)
| .Pprev = .P
| .Pprev as \$Po
| .Qprev = 1
| .Q = \$D - \$Po*\$Po
| ((\$s * 8)|isqrt) as \$L
| (3 * \$L) as \$B
| .i = 2
| .b = 0
| .q = 0
| .r = 0
| .stop = false
| until( (.i >= \$B) or .stop;
.b = idivide(\$Po + .P; .Q)
| .P = .b * .Q - .P
| .q = .Q
| .Q = .Qprev + .b * (.Pprev - .P)

| .r = (((.Q|isqrt) + 0.5)|toi)

| if ((.i % 2) == 0 and (.r*.r) == .Q) then .stop = true
else
.Qprev = .q
| .Pprev = .P
| .i += 1
end )
| if .i < \$B
then
.b = idivide(\$Po - .P; .r)
| .P = .b*.r + .P
| .Pprev = .P
| .Qprev = .r
| .Q = idivide(\$D - .Pprev*.Pprev; .Qprev)
| .i = 0
| .stop = false
| until (.stop;
.b = idivide(\$Po + .P; .Q)
| .Pprev = .P
| .P = .b * .Q - .P
| .q = .Q
| .Q = .Qprev + .b * (.Pprev - .P)
| .Qprev = .q
| .i += 1
| if (.P == .Pprev) then .stop = true else . end )
| .r = gcd(\$N; .Qprev)
| if .r != 1 and .r != \$N then .r, break \$out else empty end
else empty
end
end
// 0 ;

def examples: [
"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"
];

"[Integer, Factor, Quotient]"
"---------------------------",
(examples[] as \$example
| (\$example|tonumber) as \$N
| (\$N | squfof) as \$fact
| if \$fact == 0 then "fail"
else idivide(\$N; \$fact) as \$quot
| [\$N, \$fact, \$quot]
end
)
Output:
[Integer, Factor, Quotient]
---------------------------
[2501,61,41]
[12851,71,181]
[13289,137,97]
[75301,293,257]
[120787,43,2809]
[967009,1609,601]
[997417,257,3881]
[7091569,2663,2663]
[13290059,3119,4261]
[42854447,9689,4423]
[223553581,11213,19937]
[2027651281,46061,44021]
[11111111111,21649,513239]
[100895598169,112303,898423]
fail
[60012462237239,6862753,8744663]
[287129523414791,6059887,47381993]
[9007199254740931,10624181,847801751]
[11111111111111111,2071723,5363222357]
[314159265358979323,317213509,990371647]
[384307168202281507,415718707,924440401]
[419244183493398773,48009977,8732438749]
[658812288346769681,62222119,10588072199]
[922337203685477563,110075821,8379108103]
[1000000000000000127,111756107,8948056861]
[1152921505680588799,139001459,8294312261]
[1537228672809128917,26675843,57626245319]
[4611686018427387877,343242169,13435662733]

Julia

Modified from Wikipedia's article at []

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])
factr = square_form_factor(n)
println(factr == 0 ? "fail" : n ÷ factr)
end

Output:
Integer              Factor   Quotient
---------------------------------------------
2501                  61        41
12851                 71        181
13289                 137       97
75301                 293       257
120787                43        2809
967009                1609      601
997417                257       3881
7091569               2663      2663
13290059              3119      4261
42854447              9689      4423
223553581             11213     19937
2027651281            46061     44021
11111111111           21649     513239
100895598169          112303    898423
1002742628021         0         fail
60012462237239        6862753   8744663
287129523414791       6059887   47381993
9007199254740931      10624181  847801751
11111111111111111     2071723   5363222357
314159265358979323    317213509 990371647
384307168202281507    415718707 924440401
419244183493398773    48009977  8732438749
658812288346769681    62222119  10588072199
922337203685477563    110075821 8379108103
1000000000000000127   111756107 8948056861
1152921505680588799   139001459 8294312261
1537228672809128917   26675843  57626245319
4611686018427387877   343242169 13435662733
0.039027 seconds (698 allocations: 38.312 KiB)

Nim

Translation of: Phix
import math, strformat

const M = [uint64 1, 3, 5, 7, 11]

template isqrt(n: uint64): uint64 = uint64(sqrt(float(n)))
template isEven(n: uint64): bool = (n and 1) == 0

proc squfof(n: uint64): uint64 =

if n.isEven: return 2
var h = uint64(sqrt(float(n)) + 0.5)
if h * h == n: return h

for m in M:
if m > 1 and (n mod m == 0): return m
# Check overflow m * n.
if n > uint64.high div m: break
let mn = m * n
var r = isqrt(mn)
if r * r > mn: dec r
let rn = r

# Principal form.
var b = r
var a = 1u64
h = (rn + b) div a * a - b
var c = (mn - h * h) div a

for i in 2..<(4 * isqrt(2 * r)):
# Search principal cycle.
swap a, c
var q = (rn + b) div a
let t = b
b = q * a - b
c += q * (t - b)

if i.isEven:
r = uint64(sqrt(float(c)) + 0.5)
if r * r == c: # Square form found?

# Inverse square root.
q = (rn - b) div r
var v = q * r + b
var w = (mn - v * v) div r

# Search ambiguous cycle.
var u = r
while true:
swap w, u
r = v
q = (rn + v) div u
v = q * u - v
if v == r: break
w += q * (r - v)

# Symmetry point.
h = gcd(n, u)
if h != 1: return h

result = 1

const Data = [2501u64,
12851u64,
13289u64,
75301u64,
120787u64,
967009u64,
997417u64,
7091569u64,
13290059u64,
42854447u64,
223553581u64,
2027651281u64,
11111111111u64,
100895598169u64,
1002742628021u64,
60012462237239u64,
287129523414791u64,
9007199254740931u64,
11111111111111111u64,
314159265358979323u64,
384307168202281507u64,
419244183493398773u64,
658812288346769681u64,
922337203685477563u64,
1000000000000000127u64,
1152921505680588799u64,
1537228672809128917u64,
4611686018427387877u64]

echo "N f N/f"
echo "======================================"
for n in Data:
let f = squfof(n)
let res = if f == 1: "fail" else: &"{f:<10} {n div f}"
echo &"{n:<22} {res}"
Output:
N                      f          N/f
======================================
2501                   61         41
12851                  71         181
13289                  97         137
75301                  293        257
120787                 43         2809
967009                 1609       601
997417                 257        3881
7091569                2663       2663
13290059               3119       4261
42854447               4423       9689
223553581              11213      19937
2027651281             44021      46061
11111111111            21649      513239
100895598169           112303     898423
1002742628021          fail
60012462237239         6862753    8744663
287129523414791        6059887    47381993
9007199254740931       10624181   847801751
11111111111111111      2071723    5363222357
314159265358979323     317213509  990371647
384307168202281507     415718707  924440401
419244183493398773     48009977   8732438749
658812288346769681     62222119   10588072199
922337203685477563     110075821  8379108103
1000000000000000127    111756107  8948056861
1152921505680588799    139001459  8294312261
1537228672809128917    26675843   57626245319
4611686018427387877    343242169  13435662733

Perl

Translation of: Raku
Library: ntheory
use strict;
use warnings;
use feature 'say';
use ntheory <is_prime gcd forcomb vecprod>;

my @multiplier;
my @p = <3 5 7 11>;
forcomb { push @multiplier, vecprod @p[@_] } scalar @p;

sub sff {
my(\$N) = shift;
return 1 if is_prime \$N; # if n is prime
return sqrt \$N if sqrt(\$N) == int sqrt \$N; # if n is a perfect square

for my \$k (@multiplier) {
my \$P0 = int sqrt(\$k*\$N); # P=floor(sqrt(N)
my \$Q0 = 1; # Q=1
my \$Q = \$k*\$N - \$P0**2; # Q=N-P^2 & Q[i]
my \$P1 = \$P0; # P[i-1] = P
my \$Q1 = \$Q0; # Q[i-1] = Q
my \$P = 0; # P[i]
my \$Qn = 0; # \$P[\$i+1];
my \$b = 0; # b[i]

until (sqrt(\$Q) == int(sqrt(\$Q))) { # until Q[i] is a perfect square
\$b = int( int(sqrt(\$k*\$N) + \$P1 ) / \$Q); # floor(floor(sqrt(N+P[i-1])/Q[i])
\$P = \$b*\$Q - \$P1; # P[i]=b*Q[i]-P[i-1]
\$Qn = \$Q1 + \$b*(\$P1 - \$P); # Q[i+1]=Q[i-1]+b(P[i-1]-P[i])
(\$Q1, \$Q, \$P1) = (\$Q, \$Qn, \$P);
}

\$b = int( int( sqrt(\$k*\$N)+\$P ) / \$Q ); # b=floor((floor(sqrt(N)+P[i])/Q)
\$P1 = \$b*\$Q0 - \$P; # P[i-1]=b*Q-P[i]
\$Q = ( \$k*\$N - \$P1**2 )/\$Q0; # Q=(N-P^2)/Q & Q[i]
\$Q1 = \$Q0; # Q[i-1] = Q

while () {
\$b = int( int(sqrt(\$k*\$N)+\$P1 ) / \$Q ); # b=floor(floor(sqrt(N)+P[i-1])/Q[i])
\$P = \$b*\$Q - \$P1; # P[i]=b*Q[i]-P[i-1]
\$Qn = \$Q1 + \$b*(\$P1 - \$P); # Q[i+1]=Q[i-1]+b(P[i-1]-P[i])
last if \$P == \$P1; # until P[i+1]=P[i]
(\$Q1, \$Q, \$P1) = (\$Q, \$Qn, \$P);
}
for (gcd \$N, \$P) { return \$_ if \$_ != 1 and \$_ != \$N }
}
return 0
}

for my \$data (
11111, 2501, 12851, 13289, 75301, 120787, 967009, 997417, 4558849, 7091569, 13290059,
42854447, 223553581, 2027651281, 11111111111, 100895598169, 1002742628021, 60012462237239,
287129523414791, 11111111111111111, 384307168202281507, 1000000000000000127, 9007199254740931,
922337203685477563, 314159265358979323, 1152921505680588799, 658812288346769681,
419244183493398773, 1537228672809128917) {
my \$v = sff(\$data);
if (\$v == 0) { say 'The number ' . \$data . ' is not factored.' }
elsif (\$v == 1) { say 'The number ' . \$data . ' is a prime.' }
else { say "\$data = " . join ' * ', sort {\$a <=> \$b} \$v, int \$data/int(\$v) }
}
Output:
11111 = 41 * 271
2501 = 41 * 61
12851 = 71 * 181
13289 = 97 * 137
75301 = 257 * 293
120787 = 43 * 2809
967009 = 601 * 1609
997417 = 257 * 3881
4558849 = 383 * 11903
7091569 = 2663 * 2663
13290059 = 3119 * 4261
42854447 = 4423 * 9689
223553581 = 11213 * 19937
2027651281 = 44021 * 46061
11111111111 = 21649 * 513239
100895598169 = 112303 * 898423
The number 1002742628021 is a prime.
60012462237239 = 6862753 * 8744663
287129523414791 = 6059887 * 47381993
11111111111111111 = 2071723 * 5363222357
384307168202281507 = 415718707 * 924440401
1000000000000000127 = 111756107 * 8948056861
9007199254740931 = 10624181 * 847801751
922337203685477563 = 110075821 * 8379108103
314159265358979323 = 317213509 * 990371647
1152921505680588799 = 139001459 * 8294312261
658812288346769681 = 62222119 * 10588072199
419244183493398773 = 48009977 * 8732438749
1537228672809128917 = 26675843 * 57626245319

Phix

Translation of: C – (Classical heuristic - fixes the two incorrectly failing cases of the previous version)
--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, 120787, 967009, 997417, 7091569, 5214317, 20834839,
23515517, 33409583, 44524219, 13290059, 223553581, 42854447, 223553581, 2027651281,
11111111111, 100895598169, 1002742628021, -- (prime/expected to fail)
60012462237239, 287129523414791, 9007199254740931, 32, -- (limit of 32-bit)
11111111111111111, 314159265358979323, 384307168202281507, 419244183493398773,
658812288346769681, 922337203685477563, 1000000000000000127, 1152921505680588799,
1537228672809128917, 4611686018427387877}

printf(1,"N                      f          N/f\n")
printf(1,"======================================\n")
for i=1 to length(tests) do
atom N = tests[i]
if N=32 then
if machine_bits()=32 then exit end if
else
atom f = squfof(N)
printf(1,"%-22d %s\n",{N,iff(f=1?"fail":sprintf("%-10d %d",{f,N/f}))})
end if
end for
Output:

(on 64-bit, whereas the last 10 entries, ie 11111111111111111 on, are simply omitted on 32-bit)

N                      f          N/f
======================================
2501                   61         41
12851                  71         181
13289                  97         137
75301                  293        257
120787                 43         2809
967009                 1609       601
997417                 257        3881
7091569                2663       2663
5214317                73         71429
20834839               3361       6199
23515517               53         443689
33409583               991        33713
44524219               593        75083
13290059               3119       4261
223553581              11213      19937
42854447               4423       9689
223553581              11213      19937
2027651281             44021      46061
11111111111            21649      513239
100895598169           112303     898423
1002742628021          fail
60012462237239         6862753    8744663
287129523414791        6059887    47381993
9007199254740931       10624181   847801751
11111111111111111      2071723    5363222357
314159265358979323     317213509  990371647
384307168202281507     415718707  924440401
419244183493398773     48009977   8732438749
658812288346769681     62222119   10588072199
922337203685477563     110075821  8379108103
1000000000000000127    111756107  8948056861
1152921505680588799    139001459  8294312261
1537228672809128917    26675843   57626245319
4611686018427387877    343242169  13435662733

Raku

A just for fun snail ..

References: Algorithm, C program example from the WP page and also the pseudocode from here.

# 20210325 Raku programming solution

my @multiplier = ( 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 );

sub circumfix:<⌊ ⌋>{ \$^n.floor }; sub prefix:<>{ \$^n.sqrt }; # just for fun

sub SQUFOF ( \𝑁 ) {

return 1 if 𝑁.is-prime; # if n is prime return 1
return √𝑁 if √𝑁 == Int(√𝑁); # if n is a perfect square return √𝑁

for @multiplier -> \𝑘 {
my \Pₒ = \$ = ⌊ √(𝑘*𝑁); # P=floor(√N)
my \Qₒ = \$ = 1 ; # Q=1
my \Q = \$ = 𝑘*𝑁 - Pₒ²; # Q=N-P^2 & Q[i]
my \Pₚᵣₑᵥ = \$ = Pₒ; # P[i-1] = P
my \Qₚᵣₑᵥ = \$ = Qₒ; # Q[i-1] = Q
my \P = \$ = 0; # P[i]
my \Qₙₑₓₜ = \$ = 0; # P[i+1]
my \b = \$ = 0; # b[i]
# i = 1
repeat until √Q == Int(√Q) { # until Q[i] is a perfect square
b = ⌊⌊ √(𝑘*𝑁) + Pₚᵣₑᵥ ⌋ / Q ⌋; # floor(floor(√N+P[i-1])/Q[i])
P = b*Q - Pₚᵣₑᵥ; # P[i]=b*Q[i]-P[i-1]
Qₙₑₓₜ = Qₚᵣₑᵥ + b*(Pₚᵣₑᵥ - P); # Q[i+1]=Q[i-1]+b(P[i-1]-P[i])
( Qₚᵣₑᵥ, Q, Pₚᵣₑᵥ ) = Q, Qₙₑₓₜ, P; # i++
}

b = ⌊ ⌊ √(𝑘*𝑁)+P ⌋ / Q ⌋; # b=floor((floor(√N)+P[i])/Q)
Pₚᵣₑᵥ = b*Qₒ - P; # P[i-1]=b*Q-P[i]
Q = ( 𝑘*𝑁 - Pₚᵣₑᵥ² )/Qₒ; # Q=(N-P^2)/Q & Q[i]
Qₚᵣₑᵥ = Qₒ; # Q[i-1] = Q
# i = 1
loop { # repeat
b = ⌊ ⌊ √(𝑘*𝑁)+Pₚᵣₑᵥ ⌋ / Q ⌋; # b=floor(floor(√N)+P[i-1])/Q[i])
P = b*Q - Pₚᵣₑᵥ; # P[i]=b*Q[i]-P[i-1]
Qₙₑₓₜ = Qₚᵣₑᵥ + b*(Pₚᵣₑᵥ - P); # Q[i+1]=Q[i-1]+b(P[i-1]-P[i])
last if (P == Pₚᵣₑᵥ); # until P[i+1]=P[i]
( Qₚᵣₑᵥ, Q, Pₚᵣₑᵥ ) = Q, Qₙₑₓₜ, P; # i++
}
given 𝑁 gcd P { return \$_ if \$_ != 1|𝑁 }
} # gcd(N,P[i]) (if != 1 or N) is a factor of N, otherwise try next k
return 0 # give up
}

race for (
11111, # wikipedia.org/wiki/Shanks%27s_square_forms_factorization#Example
4558849, # example from talk page
# all of the rest are taken from the FreeBASIC entry
2501,12851,13289,75301,120787,967009,997417,7091569,13290059,
42854447,223553581,2027651281,11111111111,100895598169,1002742628021,
# time hoarders
60012462237239, # = 6862753 * 8744663 15s
287129523414791, # = 6059887 * 47381993 80s
11111111111111111, # = 2071723 * 5363222357 2m
384307168202281507, # = 415718707 * 924440401 5m
1000000000000000127, # = 111756107 * 8948056861 12m
9007199254740931, # = 10624181 * 847801751 17m
922337203685477563, # = 110075821 * 8379108103 41m
314159265358979323, # = 317213509 * 990371647 61m
1152921505680588799, # = 139001459 * 8294312261 93m
658812288346769681, # = 62222119 * 10588072199 112m
419244183493398773, # = 48009977 * 8732438749 135m
1537228672809128917, # = 26675843 * 57626245319 254m
# don't know how to handle this one
# for 1e-323, 1e-324 { my \$*TOLERANCE = \$_ ;
# say 4611686018427387877.sqrt ≅ 4611686018427387877.sqrt.Int }
# skip the perfect square check and start k with 3 to get the following
# 4611686018427387877, # = 343242169 * 13435662733 217m
) -> \data {
given data.&SQUFOF {
when 0 { say "The number ", data, " is not factored." }
when 1 { say "The number ", data, " is a prime." }
default { say data, " = ", \$_, " * ", data div \$_.Int }
}
}

Output:
11111 = 41 * 271
4558849 = 383 * 11903
2501 = 61 * 41
12851 = 71 * 181
13289 = 97 * 137
75301 = 293 * 257
120787 = 43 * 2809
967009 = 601 * 1609
997417 = 257 * 3881
7091569 = 2663 * 2663
13290059 = 3119 * 4261
42854447 = 4423 * 9689
223553581 = 11213 * 19937
2027651281 = 44021 * 46061
11111111111 = 21649 * 513239
100895598169 = 112303 * 898423
The number 1002742628021 is a prime.
60012462237239 = 6862753 * 8744663
287129523414791 = 6059887 * 47381993
11111111111111111 = 2071723 * 5363222357
384307168202281507 = 415718707 * 924440401
1000000000000000127 = 111756107 * 8948056861
9007199254740931 = 10624181 * 847801751
922337203685477563 = 110075821 * 8379108103
314159265358979323 = 317213509 * 990371647
1152921505680588799 = 139001459 * 8294312261
658812288346769681 = 62222119 * 10588072199
419244183493398773 = 48009977 * 8732438749
1537228672809128917 = 26675843 * 57626245319

Use NativeCall

Use idea similar to the second approach from here, by compiling the C (Classical heuristic) entry to a shared library and then make use of the squfof routine.

squfof.raku

# 20210326 Raku programming solution

use NativeCall;

constant LIBSQUFOF = '/home/user/LibSQUFOF.so';

sub squfof(uint64 \$n) returns uint64 is native(LIBSQUFOF) { * };

race for (
11111, # wikipedia.org/wiki/Shanks%27s_square_forms_factorization#Example
4558849, # example from talk page
# all of the rest are taken from the FreeBASIC entry
2501,12851,13289,75301,120787,967009,997417,7091569,13290059,
42854447,223553581,2027651281,11111111111,100895598169,1002742628021,
60012462237239, # = 6862753 * 8744663
287129523414791, # = 6059887 * 47381993
11111111111111111, # = 2071723 * 5363222357
384307168202281507, # = 415718707 * 924440401
1000000000000000127, # = 111756107 * 8948056861
9007199254740931, # = 10624181 * 847801751
922337203685477563, # = 110075821 * 8379108103
314159265358979323, # = 317213509 * 990371647
1152921505680588799, # = 139001459 * 8294312261
658812288346769681, # = 62222119 * 10588072199
419244183493398773, # = 48009977 * 8732438749
1537228672809128917, # = 26675843 * 57626245319
4611686018427387877, # = 343242169 * 13435662733
) -> \data {
given squfof(data) { say data, " = ", \$_, " * ", data div \$_ }
}
Output:
gcc -Wall -fPIC -shared -o LibSQUFOF.so squfof.c

file LibSQUFOF.so LibSQUFOF.so: ELF 64-bit LSB shared object, x86-64, version 1 (SYSV), dynamically linked, BuildID[sha1]=f7f9864e5508ba1de80cbc6e2f4d789f4ec448e9, not stripped raku -c squfof.raku && time ./squfof.raku Syntax OK 11111 = 41 * 271 4558849 = 383 * 11903 2501 = 61 * 41 12851 = 71 * 181 13289 = 97 * 137 75301 = 293 * 257 120787 = 43 * 2809 967009 = 1609 * 601 997417 = 257 * 3881 7091569 = 2663 * 2663 13290059 = 3119 * 4261 42854447 = 4423 * 9689 223553581 = 11213 * 19937 2027651281 = 44021 * 46061 11111111111 = 21649 * 513239 100895598169 = 112303 * 898423 1002742628021 = 1 * 1002742628021 60012462237239 = 6862753 * 8744663 287129523414791 = 6059887 * 47381993 11111111111111111 = 2071723 * 5363222357 384307168202281507 = 415718707 * 924440401 1000000000000000127 = 111756107 * 8948056861 9007199254740931 = 10624181 * 847801751 922337203685477563 = 110075821 * 8379108103 314159265358979323 = 317213509 * 990371647 1152921505680588799 = 139001459 * 8294312261 658812288346769681 = 62222119 * 10588072199 419244183493398773 = 48009977 * 8732438749 1537228672809128917 = 26675843 * 57626245319 4611686018427387877 = 343242169 * 13435662733

real 0m0.784s user 0m0.716s sys 0m0.056s

REXX

This REXX version is a modified version of the   C   code in the Wikipedia's article at   Shanks' square forms factorization.

The only "failure" of this REXX version of Shanks' square forms factorization   (when using the numbers supplied within the
REXX program)   is because the number being tested is a prime   (1,002,742,628,021).

/*REXX pgm factors an integer using Daniel Shanks' (1917-1996) square form factorization*/
numeric digits 100 /*ensure enough decimal digits.*/
call dMults 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
call dTests 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
w= length( commas(!.\$) ) /*the max width of test numbers*/
do tests=1 for !.0; n= !.tests; nc= commas(n)
f= ssff(n); fc= commas(f); wf= length(fc); if f\==0 then nf= commas(n%f)
if f\==0 then do; nfc= commas(n%f); wnfc= length(nfc); end
if f ==0 then _= " (Shank's square form factor failed.)"
else _= ' factors are: ' right( fc, max(w%2 , wf ) ) " and " ,
right(nfc, max(w%2+4, wnfc) )
say right(nc, w+5) _
end /*tests*/
exit 0 /*stick a fork in it, we're all done. */
/*──────────────────────────────────────────────────────────────────────────────────────*/
commas: parse arg ?; do jc=length(?)-3 to 1 by -3; ?=insert(',', ?, jc); end; return ?
dMults: @.\$= 0; do j=1 for arg(); @.j= arg(j); @.\$=max(@.\$, @.j); end; @.0=j-1; return
dTests: !.\$= 0; do j=1 for arg(); !.j= arg(j); !.\$=max(!.\$, !.j); end; !.0=j-1; return
gcd: procedure; parse arg x,y; do until _==0; _= x // y; x= y; y= _; end; return x
/*──────────────────────────────────────────────────────────────────────────────────────*/
iSqrt: procedure; parse arg x; r=0; q=1; do while q<=x; q=q*4; end
do while q>1; q=q%4; _=x-r-q; r=r%2; if _>=0 then do;x=_;r=r+q; end; end
return r
/*──────────────────────────────────────────────────────────────────────────────────────*/
ssff: procedure expose @.; parse arg n; n= abs(n); er= '***error***'
s= iSqrt(n); if s**2==n then return s; big= 2**digits()
do #=1 for @.0; k= @.# /*get a # from the list of low factors*/
if n>big/k then do; say er 'number is too large: ' commas(k); exit 8; end
d= n*k; po= iSqrt(d); p= po
pprev= po; QQ= d - po*po
qprev= 1; BB= iSqrt(s+s)*6
do i=2 while i<BB; b= (po+p)%QQ
p= b*QQ - p; q= QQ
QQ= qprev + b*(pprev-p); r= iSqrt(QQ)
if i//2==0 then if r*r==QQ then leave
qprev= q; pprev= p
end /*i*/
if i>=BB then iterate
b= (po-p)%r; p= b*r + p
pprev= p; qprev= r
QQ= (d - pprev*pprev)%qprev
do until p==pprev; pprev= p
b= (po+p)%QQ; q= QQ; p= b*QQ - p
QQ= qprev + b*(pprev-p); qprev= q
end /*until*/
r= gcd(n, qprev)
if r\==1 then if r\==n then return r
end /*#*/
return 0
output   when using the internal default input:
2,501  factors are:            61   and                 41
12,851  factors are:            71   and                181
13,289  factors are:           137   and                 97
75,301  factors are:           293   and                257
120,787  factors are:            43   and              2,809
967,009  factors are:         1,609   and                601
997,417  factors are:           257   and              3,881
7,091,569  factors are:         2,663   and              2,663
13,290,059  factors are:         3,119   and              4,261
42,854,447  factors are:         9,689   and              4,423
223,553,581  factors are:        11,213   and             19,937
2,027,651,281  factors are:        46,061   and             44,021
11,111,111,111  factors are:        21,649   and            513,239
100,895,598,169  factors are:       112,303   and            898,423
1,002,742,628,021    (Shank's square form factor failed.)
60,012,462,237,239  factors are:     6,862,753   and          8,744,663
287,129,523,414,791  factors are:     6,059,887   and         47,381,993
9,007,199,254,740,931  factors are:    10,624,181   and        847,801,751
11,111,111,111,111,111  factors are:     2,071,723   and      5,363,222,357
314,159,265,358,979,323  factors are:   317,213,509   and        990,371,647
384,307,168,202,281,507  factors are:   415,718,707   and        924,440,401
419,244,183,493,398,773  factors are:    48,009,977   and      8,732,438,749
658,812,288,346,769,681  factors are:    62,222,119   and     10,588,072,199
922,337,203,685,477,563  factors are:   110,075,821   and      8,379,108,103
1,000,000,000,000,000,127  factors are:   111,756,107   and      8,948,056,861
1,152,921,505,680,588,799  factors are:   139,001,459   and      8,294,312,261
1,537,228,672,809,128,917  factors are:    26,675,843   and     57,626,245,319
4,611,686,018,427,387,877  factors are:   343,242,169   and     13,435,662,733

Wren

Library: Wren-long
Library: Wren-big
Library: Wren-fmt

This is based on the C code in the Wikipedia article and the FreeBASIC entry examples.

'0' is not actually an example here but is used by FreeBASIC to mark the end of the 'data' statement list so I've ignored that.

As Wren doesn't natively support unsigned 64-bit arithmetic, I've used the ULong class in the first above named module for this task.

Even so, there are two examples which fail (1000000000000000127 and 1537228672809128917) because the code is unable to process enough 'multipliers' before an overflow situation is encountered. To deal with this, the program automatically switches to BigInt to process such cases.

import "/long" for ULong
import "/big" for BigInt
import "/fmt" for Fmt

var multipliers = [
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
]

var squfof = Fn.new { |N|
var s = ULong.new((N.toNum.sqrt + 0.5).floor)
if (s*s == N) return s
for (multiplier in multipliers) {
var T = ULong
var n = N
if (n > ULong.largest/multiplier) {
T = BigInt
n = BigInt.new(n.toString)
}
var D = n * multiplier
var P = D.isqrt
var Pprev = P
var Po = Pprev
var Qprev = T.one
var Q = D - Po*Po
var L = (s * 8).isqrt.toSmall
var B = 3 * L
var i = 2
var b = T.zero
var q = T.zero
var r = T.zero
while (i < B) {
b = (Po + P) / Q
P = b * Q - P
q = Q
Q = Qprev + b * (Pprev - P)
r = T.new((Q.toNum.sqrt + 0.5).floor)
if ((i & 1) == 0 && r*r == Q) break
Qprev = q
Pprev = P
i = i + 1
}
if (i < B) {
b = (Po - P) / r
Pprev = P = b*r + P
Qprev = r
Q = (D - Pprev*Pprev) / Qprev
i = 0
while (true) {
b = (Po + P) / Q
Pprev = P
P = b * Q - P
q = Q
Q = Qprev + b * (Pprev - P)
Qprev = q
i = i + 1
if (P == Pprev) break
}
r = T.gcd(n, Qprev)
if (r != T.one && r != n) return (r is ULong) ? r : ULong.new(r.toString)
}
}
return ULong.zero
}

var examples = [
"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"
]

System.print("Integer Factor Quotient")
System.print("------------------------------------------")
for (example in examples) {
var N = ULong.new(example)
var fact = squfof.call(N)
var quot = (fact.isZero) ? "fail" : (N / fact).toString
Fmt.print("\$-20s \$-10s \$s", N, fact, quot)
}
Output:
Integer              Factor     Quotient
------------------------------------------
2501                 61         41
12851                71         181
13289                137        97
75301                293        257
120787               43         2809
967009               1609       601
997417               257        3881
7091569              2663       2663
13290059             3119       4261
42854447             9689       4423
223553581            11213      19937
2027651281           46061      44021
11111111111          21649      513239
100895598169         112303     898423
1002742628021        0          fail
60012462237239       6862753    8744663
287129523414791      6059887    47381993
9007199254740931     10624181   847801751
11111111111111111    2071723    5363222357
314159265358979323   317213509  990371647
384307168202281507   415718707  924440401
419244183493398773   48009977   8732438749
658812288346769681   62222119   10588072199
922337203685477563   110075821  8379108103
1000000000000000127  111756107  8948056861
1152921505680588799  139001459  8294312261
1537228672809128917  26675843   57626245319
4611686018427387877  343242169  13435662733