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 10¹⁰ and 10¹⁸, and will likely remain so.’

An integral binary quadratic form (bqf) is a polynomial f(x,y) = ax² + bxy + cy² with integer coefficients and discriminant D = b² – 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)² – 4ac = a(a·k² – 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 1 'higher if necessary 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

  print

  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 = 898423  N/f = 112303

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.0362443 s