Square form factorization: Difference between revisions
m (→{{header|REXX}}: changed a comment in the REXX section header.) |
m (timings) |
||
Line 535: | Line 535: | ||
println("Integer Factor Quotient\n", "-"^45) |
println("Integer Factor Quotient\n", "-"^45) |
||
for n in Int128.([ |
@time for n in Int128.([ |
||
2501, 12851, 13289, 75301, 120787, 967009, 997417, 7091569, 13290059, 42854447, 223553581, |
2501, 12851, 13289, 75301, 120787, 967009, 997417, 7091569, 13290059, 42854447, 223553581, |
||
2027651281, 11111111111, 100895598169, 1002742628021, 60012462237239, 287129523414791, |
2027651281, 11111111111, 100895598169, 1002742628021, 60012462237239, 287129523414791, |
||
9007199254740931, 11111111111111111, 314159265358979323, 384307168202281507, 419244183493398773, |
9007199254740931, 11111111111111111, 314159265358979323, 384307168202281507, 419244183493398773, |
||
658812288346769681, 922337203685477563, 1000000000000000127, 1152921505680588799, |
658812288346769681, 922337203685477563, 1000000000000000127, 1152921505680588799, |
||
1537228672809128917, 4611686018427387877 |
1537228672809128917, 4611686018427387877]) |
||
print(rpad(n, 22)) |
print(rpad(n, 22)) |
||
factr = square_form_factor(n) |
factr = square_form_factor(n) |
||
Line 578: | Line 578: | ||
1537228672809128917 26675843 57626245319 |
1537228672809128917 26675843 57626245319 |
||
4611686018427387877 343242169 13435662733 |
4611686018427387877 343242169 13435662733 |
||
0.039027 seconds (698 allocations: 38.312 KiB) |
|||
0 0 fail |
|||
</pre> |
</pre> |
||
Revision as of 23:25, 12 March 2021
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
'------------------------------------------------
- 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>
- 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) @time 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]) 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.039027 seconds (698 allocations: 38.312 KiB)
Phix
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
This REXX version is a modified version of the C code in the Wikipedia's article at Shanks' square forms factorization.
The only "failure" of this REXX version of Shanks' square forms factorization (when using the numbers supplied
within the
REXX program) is because the number being tested is a prime (1,002,742,628,021).
<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 then if 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 /*until*/ 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
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