Square form factorization: Difference between revisions

From Rosetta Code
Content added Content deleted
mNo edit summary
(→‎{{header|REXX}}: added the computer programming language REXX.)
Line 686: Line 686:
1537228672809128917 0 fail
1537228672809128917 0 fail
4611686018427387877 343242169 13435662733
4611686018427387877 343242169 13435662733
</pre>

=={{header|REXX}}==
<lang rexx>/*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); pprev= po
p= po; qprev= 1
QQ= d - po * po
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 & 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
p= b * QQ - p; q= QQ
QQ= qprev + b * (pprev - p)
qprev= q
end /*forever*/
r= gcd(n, qprev)
if r\==1 then if r\==n then return r
end /*#*/
return 0</lang>
{{out|output|text=&nbsp; when using the internal default input:}}
<pre>
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
</pre>
</pre>


Revision as of 20:19, 12 March 2021

Task
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 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 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)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 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

REXX

<lang rexx>/*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);           pprev= po
              p= po;                  qprev= 1
              QQ= d  -  po * po
              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  &  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
                                                 p= b * QQ  -  p;        q= QQ
                                                 QQ= qprev  +  b * (pprev - p)
                                                 qprev= q
                                                 end   /*forever*/
              r= gcd(n, qprev)
              if r\==1  then if  r\==n  then return r
              end   /*#*/
      return 0</lang>
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.

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
Retrieved from "https://rosettacode.org/wiki/Square_form_factorization?oldid=309024"