Square form factorization: Difference between revisions
mNo edit summary |
mNo edit summary |
||
Line 1: | Line 1: | ||
{{draft task|mathematics}} |
{{draft task|mathematics}} |
||
;Task. |
;Task. |
||
Daniel Shanks's Square Form Factorization [[wp:Shanks%27s_square_forms_factorization|(SquFoF)]]. |
[[wp:Daniel_Shanks|Daniel Shanks]]'s Square Form Factorization [[wp:Shanks%27s_square_forms_factorization|(SquFoF)]]. |
||
Line 9: | Line 10: | ||
''for numbers between 10<sup>10</sup> and 10<sup>18</sup>, and will likely remain so.’'' |
''for numbers between 10<sup>10</sup> and 10<sup>18</sup>, and will likely remain so.’'' |
||
An integral |
An integral [[wp:Binary_quadratic_form|binary quadratic form]] is a polynomial |
||
{{math|''f''(''x,y'') = ''ax<sup>2</sup>'' + ''bxy'' + ''cy<sup>2</sup>''}} |
{{math|''f''(''x,y'') = ''ax<sup>2</sup>'' + ''bxy'' + ''cy<sup>2</sup>''}} |
||
with integer coefficients and discriminant {{math|''D'' = ''b<sup>2</sup>'' – ''4ac''}}. |
with integer coefficients and discriminant {{math|''D'' = ''b<sup>2</sup>'' – ''4ac''}}. |
||
Line 20: | Line 21: | ||
SquFoF works because there are cycles containing ''ambiguous forms'', with the property that ''a'' divides ''b''. |
SquFoF works because there are cycles containing ''ambiguous forms'', with the property that ''a'' divides ''b''. |
||
They come in pairs of associated forms {{math|('' |
They come in pairs of associated forms {{math|(''a, b, c'') and (''c, b, a'')}}, 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 |
symmetry points. If an ambiguous form is found (there is one for each divisor of D), write the discriminant as |
||
{{math| |
{{math|(''ak'')''<sup>2</sup>'' – ''4ac'' = ''a''(''a·k<sup>2</sup>'' – ''4c'') = ''4N''}} |
||
and (if a is not equal to 1 or 2) N is split. |
and (if a is not equal to 1 or 2) N is split. |
||
Line 34: | Line 35: | ||
Then ''a'' or ''a/2'' divides D and therefore N. |
Then ''a'' or ''a/2'' divides D and therefore N. |
||
To avoid trivial factorizations, Shanks created a list (queue) to |
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, |
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. |
no roots land in the principal cycle itself and cases a = 1 or a = 2 do not happen. |
||
Line 45: | Line 46: | ||
;Reference. |
;Reference. |
||
⚫ | |||
[https://www.ams.org/journals/mcom/2008-77-261/S0025-5718-07-02010-8/S0025-5718-07-02010-8.pdf] |
|||
⚫ | |||
Line 66: | Line 66: | ||
'------------------------------------------------ |
'------------------------------------------------ |
||
const MxN = culngint(1) shl 62 |
const MxN = culngint(1) shl 62 |
||
' |
'input maximum |
||
type arg |
type arg |
||
'SquFoF arguments |
|||
as ulong m, f |
as ulong m, f |
||
as integer vb |
as integer vb |
||
Line 75: | Line 76: | ||
type bqf |
type bqf |
||
declare sub rho (byval sw as integer) |
declare sub rho (byval sw as integer) |
||
'reduce indefinite form, set sw = 0 to initialize |
'reduce indefinite form, set sw = 0 to initialize |
||
declare function issqr (byref r as ulong) as integer |
declare function issqr (byref r as ulong) as integer |
||
'return -1 if c is square, set r:= sqrt(c) |
'return -1 if c is square, set r:= sqrt(c) |
||
Line 96: | Line 97: | ||
end type |
end type |
||
⚫ | |||
⚫ | |||
dim shared N as ulongint |
dim shared N as ulongint |
||
'the number to split |
|||
dim shared flag as integer |
dim shared flag as integer |
||
'signal to end all threads |
|||
Line 167: | Line 168: | ||
#macro red(r, a) |
#macro red(r, a) |
||
r = iif(a and 1, a, a shr 1) |
r = iif(a and 1, a, a shr 1) |
||
if m > |
if m > 2 then |
||
r = iif(r mod m, r, r \ m) |
r = iif(r mod m, r, r \ m) |
||
end if |
end if |
||
Line 214: | Line 215: | ||
h = N |
h = N |
||
rp->f = 1 |
rp->f = 1 |
||
'multiplier |
|||
m = rp->m |
m = rp->m |
||
if m > 1 then |
if m > 1 then |
||
Line 322: | Line 324: | ||
data 1537228672809128917 |
data 1537228672809128917 |
||
data 4611686018427387877 |
data 4611686018427387877 |
||
data 0 |
|||
'main |
'main |
||
Line 395: | Line 398: | ||
</lang> |
</lang> |
||
{{out| |
{{out|(No need to replicate all)}} |
||
<pre> |
<pre> |
||
N = 2501 |
N = 2501 |
||
Line 428: | Line 431: | ||
N = 223553581 |
N = 223553581 |
||
f = |
f = 11213 N/f = 19937 |
||
N = 2027651281 |
N = 2027651281 |
||
Line 481: | Line 484: | ||
f = 343242169 N/f = 13435662733 |
f = 343242169 N/f = 13435662733 |
||
total time: 0. |
total time: 0.0229552 s |
||
</pre> |
</pre> |
Revision as of 16:07, 11 March 2021
- 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 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 (a, b, c) and (c, b, a), 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 (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 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 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 'input maximum
type arg
'SquFoF arguments 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 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 variables dim shared N as ulongint 'the number to split dim shared flag as integer 'signal to end all threads
'------------------------------------------------
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 > 2 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 'multiplier 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 data 0
'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>
- (No need to replicate all):
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 = 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.0229552 s