Square form factorization
- 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 (bqf) is a polynomial f(x,y) = ax2 + bxy + cy2 with integer coefficients and discriminant D = b2 – 4ac. For each positive discriminant there are nearly always multiple forms (a, b, c).
The next form in a periodic sequence (cycle) of adjacent forms is found by applying a reduction operator rho. It is 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 works because there are cycles containing ambiguous forms, with the property that a divides b. They come in pairs of associated forms (c, b, a) and (a, b, c), easy to spot and obviously called symmetry points. If an ambiguous form is found (there is one for each divisor of D), write the discriminant as D = (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 a form on an ambiguous cycle is squared, that square form will always land in the principal cycle. Conversely, the square root of a form on 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 store 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 proper square forms. 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. (A short trial division loop acts as sentry). If N is prime or the cube of a prime, there are only trivial squares and the program will duly report failure.
- Reference.
[1] A detailed analysis of SquFoF (2007)
FreeBASIC
<lang freebasic> ' *********************************************** 'subject: Shanks's square form factorization: ' ambiguous forms of discriminant 4N ' give factors of N. 'tested : FreeBasic 1.07.1
const qx = 50
'queue size
'------------------------------------------------ const MxN = culngint(1) shl 62 'argument maximum
type arg
as ulong m, f as integer vb
end type
type bqf
declare sub rho (byval sw as integer) 'reduce indefinite form, set sw = 0 to initialize a declare function issqr (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 ulongint mN 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 integer t as ulong a(qx + 1), L, m
end type
'global variable
dim shared N as ulongint
dim shared flag as integer
'------------------------------------------------
sub bqf.rho (byval sw as integer)
dim as ulong q, t = b
swap a, c 'residue q = (rN + b) \ a b = q * a - b if sw then 'pseudo-square c += q * (t - b) else 'initialize c = (mN - culngint(b) * b) \ a end if
end sub
function bqf.issqr (byref r as ulong) as integer static as integer q64(63), q63(62), q55(54), sw = 0 if sw = 0 then
sw = -1 'tabulate quadratic residues dim i as ulong for i = 0 to 31 r = i * i q64(r and 63) =-1 q63(r mod 63) =-1 q55(r mod 55) =-1 next i
end if issqr = 0
if q64(c and 63) then
if q63(c mod 63) then if q55(c mod 55) then '>98% non-squares filtered r = culng(sqr(c)) if c = r * r then return -1 end if end if
end if end function
sub bqf.qform (byref g as string, byval t as integer) dim as longint d, u = a, v = b, w = c
if vb = 0 then exit sub
'{D/4 = mN} d = v * v + u * w if mN - d then print "fail:"; d: exit sub end if
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 > 1 then r = iif(r mod m, r, r \ m) end if
- endmacro
sub queue.enq (byref P as bqf) if t = qx then exit sub dim s as ulong
red(s, P.c) if s < L then t += 2 'enqueue P.b, P.c a(t) = P.b mod s a(t + 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 h 'principal and ambiguous cycles dim as bqf P, A dim Q as queue
h = culngint(sqr(N)) if N = h * h then
'N is square rp->f = culng(h) flag =-1: exit sub
end if
h = N rp->f = 1 m = rp->m if m > 1 then
'check overflow m * N if h > (MxN \ m) then exit sub h *= m
end if
P.mN = h A.mN = h r = int(sqr(h)) 'float64 fix if culngint(r) * r > h 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.t = -2: 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 P.rho(0) 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(1) if (i and 1) = 0 then
if P.issqr(r) then 'square form found
if P.c = 1 then exit for 'cycle too short
if Q.pro(P, r) then
P.qform("P", i) 'inverse square root A.b =-P.b: A.c = r A.rho(0): j = 1 A.qform("A", j)
do 'search ambiguous cycle t = A.b A.rho(1): 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 end if ' even indices
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 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
'main '------------------------------------------------ const tx = 4 dim as double tim = timer dim as ulongint i, f dim h(4) as any ptr dim a(4) as arg dim t as integer
width 64, 30 cls
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
do : read N loop until N < MxN if N < 2 then exit do
print "N = "; N
f = iif(N and 1, 1, 2) if f = 1 then for i = 3 to 37 step 2 if (N mod i) = 0 then f = i: exit for next i end if
if f = 1 then flag = 0
for t = 1 to tx h(t) = threadcreate(@squfof, @a(t)) next t
squfof(@a(0))
f = a(0).f for t = 1 to tx threadwait(h(t)) if f = 1 then f = a(t).f next t
'assert if N mod f then f = 1
elseif f = N then 'small prime N f = 1 end if
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 </lang>
- For reference only:
N = 2501 f = 41 N/f = 61 N = 12851 f = 181 N/f = 71 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 = 42854447 f = 4423 N/f = 9689 N = 223553581 f = 19937 N/f = 11213 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.0399536 s