Square form factorization

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

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 is a polynomial f(x,y) = ax² + bxy + cy² with integer coefficients and discriminant D = b² – 4ac. 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)² – 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 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. (A trial division loop acts as sentry). If N is prime or the cube of a prime, there are improper squares only 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

'------------------------------------------------

1. 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

1. 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

  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>

Examples:

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

Julia

Modified from Wikipedia's article at [[2]] <lang julia>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) 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, 0])
   print(rpad(n, 22))
   factr = square_form_factor(n)
   print(rpad(factr, 10))
   println(factr == 0 ? "fail" : n ÷ factr)

end

</lang>

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                     0         fail

Phix

Translation of: Julia – and wikipedia

An extra couple of the high 19-digit tests fail compared to other entries, but I can live with that. 0 moved up to mark the 32-bit limit. <lang Phix>--requires(64) -- (decided to limit 32-bit explicitly instead) constant 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},

        INT_MAX = power(2,iff(machine_bits()=32?53:64))

function square_form_factor(atom n)

   atom s = round(sqrt(n))
   if s * s == n then return s end if
   for mi=1 to length(multiplier) do
       atom k = multiplier[mi]
       if n > INT_MAX/k then exit end if
       atom d = k * n,
            p0 = floor(sqrt(d)),
            pprev = p0, p = p0,
            qprev = 1,
            Q = d - p0 * p0,
            l = floor(2 * sqrt(2 * s)),
            B = 3*l,
            i = 2, r, b, q
       while i < B do
           b = floor((p0 + p) / Q)
           p = b * Q - p
           q = Q
           Q = qprev + b * (pprev - p)
           r = round(sqrt(Q))
           if remainder(i,2)!=1 and r * r == Q then exit end if
           qprev = q
           pprev = p
           i += 1
       end while
       if i < B then
           b = floor((p0 - p) / r)
           p = b * r + p
           pprev = p
           qprev = r
           Q = floor((d - pprev * pprev) / qprev)
           i = 0
           while true do
               b = floor((p0 + p) / Q)
               pprev = p
               p = b * Q - p
               q = Q
               Q = qprev + b * (pprev - p)
               qprev = q
               i += 1
               if p == pprev then exit end if
           end while
           r = gcd(n, qprev)
           if r != 1 and r != n then return r end if
       end if
   end for
   return 0

end function

constant tests = {2501, 12851, 13289, 75301, 120787, 967009, 997417, 7091569, 13290059, 42854447, 223553581,

                 2027651281, 11111111111, 100895598169, 1002742628021, 60012462237239, 287129523414791,
                 0, -- (limit of 32-bit)
                 9007199254740931, 11111111111111111, 314159265358979323, 384307168202281507, 419244183493398773,
                 658812288346769681, 922337203685477563, 1000000000000000127, 1152921505680588799,
                 1537228672809128917, 4611686018427387877}

printf(1,"Integer Factor Quotient\n") printf(1,"======================================\n") for i=1 to length(tests) do

   atom n = tests[i],
        factr = square_form_factor(n)
   string f = iff(factr=0?"fail":sprintf("%d",n/factr))
   printf(1,"%-22d %-10d %s\n",{n,factr,f})
   if machine_bits()=32 and n=0 then exit end if

end for</lang>

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
0                      0          fail
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

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.

This is still 6 times faster than processing the whole lot using BigInt. <lang ecmascript>import "/long" for ULong import "/big" for BigInt import "/fmt" for Fmt

var 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

]

var squfof = Fn.new { |N|

   var s = N.isqrt
   if (s*s == N) return s
   for (k in 0...multiplier.count) {
       var T = ULong
       var n = N
       if (n > T.largest/multiplier[k]) {
           T = BigInt
           n = BigInt.new(n.toString)
       }
       var D = n * multiplier[k]
       var P = D.isqrt
       var Pprev = P
       var Po = Pprev
       var Qprev = T.one
       var Q = D - Po*Po
       var L = ((s * 2).isqrt * 2).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 = Q.isqrt
           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)

}</lang>

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