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. 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)
C
Based on Wp
From the Wikipedia entry for Shanks's square forms factorization [[2]] <lang c>#include <math.h>
- include <stdio.h>
- define nelems(x) (sizeof(x) / sizeof((x)[0]))
const unsigned long 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};
unsigned long long gcd(unsigned long long a, unsigned long long b) {
while (b != 0)
{
a %= b;
a ^= b;
b ^= a;
a ^= b;
}
return a;
}
unsigned long long SQUFOF( unsigned long long N ) {
unsigned long long D, Po, P, Pprev, Q, Qprev, q, b, r, s;
unsigned long L, B, i;
s = (unsigned long long)(sqrtl(N)+0.5);
if (s*s == N) return s;
for (int k = 0; k < nelems(multiplier) && N <= 0xffffffffffffffff/multiplier[k]; k++) {
D = multiplier[k]*N;
Po = Pprev = P = sqrtl(D);
Qprev = 1;
Q = D - Po*Po;
L = 2 * sqrtl( 2*s );
B = 3 * L;
for (i = 2 ; i < B ; i++) {
b = (unsigned long long)((Po + P)/Q);
P = b*Q - P;
q = Q;
Q = Qprev + b*(Pprev - P);
r = (unsigned long long)(sqrtl(Q)+0.5);
if (!(i & 1) && r*r == Q) break;
Qprev = q;
Pprev = P;
};
if (i >= B) continue;
b = (unsigned long long)((Po - P)/r);
Pprev = P = b*r + P;
Qprev = r;
Q = (D - Pprev*Pprev)/Qprev;
i = 0;
do {
b = (unsigned long long)((Po + P)/Q);
Pprev = P;
P = b*Q - P;
q = Q;
Q = Qprev + b*(Pprev - P);
Qprev = q;
i++;
} while (P != Pprev);
r = gcd(N, Qprev);
if (r != 1 && r != N) return r;
}
return 0;
}
int main(int argc, char *argv[]) {
int i;
const unsigned long long data[] = {
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};
for(int i = 0; i < nelems(data); i++) {
unsigned long long example, factor, quotient;
example = data[i];
factor = SQUFOF(example);
if(factor == 0) {
printf("%llu was not factored.\n", example);
}
else {
quotient = example / factor;
printf("Integer %llu has factors %llu and %llu\n",
example, factor, quotient);
}
}
}
</lang>
- Output:
Integer 2501 has factors 61 and 41 Integer 12851 has factors 71 and 181 Integer 13289 has factors 137 and 97 Integer 75301 has factors 293 and 257 Integer 120787 has factors 43 and 2809 Integer 967009 has factors 1609 and 601 Integer 997417 has factors 257 and 3881 Integer 7091569 has factors 2663 and 2663 Integer 13290059 has factors 3119 and 4261 Integer 42854447 has factors 9689 and 4423 Integer 223553581 has factors 11213 and 19937 Integer 2027651281 has factors 46061 and 44021 Integer 11111111111 has factors 21649 and 513239 Integer 100895598169 has factors 112303 and 898423 1002742628021 was not factored. Integer 60012462237239 has factors 6862753 and 8744663 Integer 287129523414791 has factors 6059887 and 47381993 Integer 9007199254740931 has factors 10624181 and 847801751 Integer 11111111111111111 has factors 2071723 and 5363222357 Integer 314159265358979323 has factors 317213509 and 990371647 Integer 384307168202281507 has factors 415718707 and 924440401 Integer 419244183493398773 has factors 48009977 and 8732438749 Integer 658812288346769681 has factors 62222119 and 10588072199 Integer 922337203685477563 has factors 110075821 and 8379108103 1000000000000000127 was not factored. Integer 1152921505680588799 has factors 139001459 and 8294312261 1537228672809128917 was not factored. Integer 4611686018427387877 has factors 343242169 and 13435662733
Classical heuristic
See Discussion. <lang c>//SquFoF: minimalistic version without queue. //Classical heuristic. Tested: tcc 0.9.27
- include <math.h>
- include <stdio.h>
//input maximum
- define MxN ((unsigned long long) 1 << 62)
//reduce indefinite form
- define rho(a, b, c) { \
t = c; c = a; a = t; t = b; \ q = (rN + b) / a; \ b = q * a - b; \ c += q * (t - b); }
//initialize
- define rhoin(a, b, c) { \
rho(a, b, c) h = b; \ c = (mN - h * h) / a; }
- define gcd(a, b) while (b) { \
t = a % b; a = b; b = t; }
//multipliers const unsigned long m[] = {1, 3, 5, 7, 11, 0};
//square form factorization unsigned long squfof( unsigned long long N ) { unsigned long a, b, c, u, v, w, rN, q, t, r; unsigned long long mN, h; int i, ix, k = 0;
if ((N & 1)==0) return 2;
h = floor(sqrt(N)+ 0.5); if (h * h == N) return h;
while (m[k]) {
if (k && N % m[k]==0) return m[k];
//check overflow m * N
if (N > MxN / m[k]) break;
mN = N * m[k++];
r = floor(sqrt(mN));
h = r; //float64 fix
if (h * h > mN) r -= 1;
rN = r;
//principal form
b = r; c = 1;
rhoin(a, b, c)
//iteration bound
ix = floor(sqrt(2*r)) * 4;
//search principal cycle
for (i = 2; i < ix; i += 2) {
rho(a, b, c)
//even step
r = floor(sqrt(c)+ 0.5);
if (r * r == c) {
//square form found
//inverse square root
v = -b; w = r;
rhoin(u, v, w)
//search ambiguous cycle
do { r = v;
rho(u, v, w)
} while (v != r);
//symmetry point
h = N; gcd(h, u)
if (h != 1) return h;
}
rho(a, b, c)
//odd step
}
}
return 1;
}
void main(void) { const unsigned long long data[] = {
2501, 12851, 13289, 75301, 120787, 967009, 997417, 7091569,
5214317, 20834839, 23515517, 33409583, 44524219,
13290059, 223553581, 2027651281, 11111111111, 100895598169, 1002742628021, 60012462237239, 287129523414791, 9007199254740931, 11111111111111111, 314159265358979323, 384307168202281507, 419244183493398773, 658812288346769681, 922337203685477563, 1000000000000000127, 1152921505680588799, 1537228672809128917, 4611686018427387877, 0};
unsigned long long N, f; int i = 0;
while (1) {
N = data[i++];
//scanf("%llu", &N);
if (N < 2) break;
printf("N = %llu\n", N);
f = squfof(N);
if (N % f) f = 1;
if (f == 1) printf("fail\n\n");
else printf("f = %llu N/f = %llu\n\n", f, N/f);
}
}</lang>
- Showing problem cases only:
... N = 5214317 f = 73 N/f = 71429 N = 20834839 f = 3361 N/f = 6199 N = 23515517 f = 53 N/f = 443689 N = 33409583 f = 991 N/f = 33713 N = 44524219 f = 593 N/f = 75083 ... N = 1000000000000000127 f = 111756107 N/f = 8948056861 ... N = 1537228672809128917 f = 26675843 N/f = 57626245319 ...
FreeBASIC
<lang freebasic>' *********************************************** 'subject: Shanks's square form factorization: ' ambiguous forms of discriminant 4N ' give factors of N. 'tested : FreeBasic 1.07.1
'------------------------------------------------
const MxN = culngint(1) shl 62
'input maximum
const qx = (1 shl 5) - 1 'queue size
type arg
'squfof arguments as ulong m, f as integer vb
end type
type bqf
declare sub rho () 'reduce indefinite form declare function issq (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 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 ulong a(qx), L, m as integer k, t
end type
'global variables dim shared N as ulongint 'the number to split
dim shared flag as integer 'signal to end all threads
dim shared as ubyte q1024(1023), q3465(3464) 'quadratic residue tables
'------------------------------------------------
sub bqf.rho ()
dim as ulong q, t
swap a, c 'residue q = (rN + b) \ a t = b: b = q * a - b 'pseudo-square c += q * (t - b)
end sub
'initialize form
- macro rhoin(F)
F.rho : h = F.b F.c = (mN - h * h) \ F.a
- endmacro
function bqf.issq (byref r as ulong) as integer if q1024(c and 1023) andalso q3465(c mod 3465) then
'98.6% non-squares filtered r = culng(sqr(c)) if r * r = c then return -1
end if issq = 0 end function
sub bqf.qform (byref g as string, byval t as integer) if vb = 0 then exit sub dim as longint u = a, v = b, w = c
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) dim s as ulong
red(s, P.c)
if s < L then
'circular queue
k = (k + 2) and qx
if k > t then t = k
'enqueue P.b, P.c
a(k) = P.b mod s
a(k + 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 mN, h 'principal and ambiguous cycles dim as bqf P, A dim Q as queue
if (N and 1) = 0 then
rp->f = 2 ' even N flag =-1: exit sub
end if
h = culngint(sqr(N)) if h * h = N then
'N is square rp->f = culng(h) flag =-1: exit sub
end if
rp->f = 1 'multiplier m = rp->m if m > 1 then
if (N mod m) = 0 then
rp->f = m ' m | N
flag =-1: exit sub
end if
'check overflow m * N if N > (MxN \ m) then exit sub
end if mN = N * m
r = int(sqr(mN)) 'float64 fix if culngint(r) * r > mN 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.k = -2: Q.t = -1: 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 rhoin(P) 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
if (i and 1) = 0 andalso P.issq(r) then
'square form found
if Q.pro(P, r) then
P.qform("P", i)
'inverse square root
A.b =-P.b: A.c = r
rhoin(A): j = 1
A.qform("A", j)
do
'search ambiguous cycle
t = A.b
A.rho: 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
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 23515517 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 h(4) as any ptr dim a(4) as arg dim as ulongint f dim as integer s, t
width 64, 30 cls
'tabulate quadratic residues for t = 0 to 1540
s = t * t q1024(s and 1023) =-1 q3465(s mod 3465) =-1
next t
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
flag = 0
for t = 1 to tx + 1 step 2
if t < tx then
h(t) = threadcreate(@squfof, @a(t))
end if
squfof(@a(t - 1))
f = a(t - 1).f
if t < tx then
threadwait(h(t))
if f = 1 then f = a(t).f
end if
if f > 1 then exit for next t
'assert if N mod f then f = 1
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 = 61 N/f = 41 N = 12851 f = 71 N/f = 181 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 = 23515517 f = 53 N/f = 443689 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.0170462 s
Go
Yet another solution based on the Wikipedia C code.
Rather than juggling with big.Int, I've just allowed the two 'awkward' cases to fail. <lang go>package main
import (
"fmt" "math"
)
func isqrt(x uint64) uint64 {
x0 := x >> 1
x1 := (x0 + x/x0) >> 1
for x1 < x0 {
x0 = x1
x1 = (x0 + x/x0) >> 1
}
return x0
}
func gcd(x, y uint64) uint64 {
for y != 0 {
x, y = y, x%y
}
return x
}
var multiplier = []uint64{
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,
}
func squfof(N uint64) uint64 {
s := uint64(math.Sqrt(float64(N)) + 0.5)
if s*s == N {
return s
}
for k := 0; k < len(multiplier) && N <= math.MaxUint64/multiplier[k]; k++ {
D := multiplier[k] * N
P := isqrt(D)
Pprev := P
Po := Pprev
Qprev := uint64(1)
Q := D - Po*Po
L := uint32(isqrt(8 * s))
B := 3 * L
i := uint32(2)
var b, q, r uint64
for ; i < B; i++ {
b = uint64((Po + P) / Q)
P = b*Q - P
q = Q
Q = Qprev + b*(Pprev-P)
r = uint64(math.Sqrt(float64(Q)) + 0.5)
if (i&1) == 0 && r*r == Q {
break
}
Qprev = q
Pprev = P
}
if i >= B {
continue
}
b = uint64((Po - P) / r)
P = b*r + P
Pprev = P
Qprev = r
Q = (D - Pprev*Pprev) / Qprev
i = 0
for {
b = uint64((Po + P) / Q)
Pprev = P
P = b*Q - P
q = Q
Q = Qprev + b*(Pprev-P)
Qprev = q
i++
if P == Pprev {
break
}
}
r = gcd(N, Qprev)
if r != 1 && r != N {
return r
}
}
return 0
}
func main() {
examples := []uint64{
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,
}
fmt.Println("Integer Factor Quotient")
fmt.Println("------------------------------------------")
for _, N := range examples {
fact := squfof(N)
quot := "fail"
if fact > 0 {
quot = fmt.Sprintf("%d", N/fact)
}
fmt.Printf("%-20d %-10d %s\n", 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 0 fail 1152921505680588799 139001459 8294312261 1537228672809128917 0 fail 4611686018427387877 343242169 13435662733
Julia
Modified from Wikipedia's article at [[3]] <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)
Nim
<lang Nim>import math, strformat
const M = [uint64 1, 3, 5, 7, 11]
template isqrt(n: uint64): uint64 = uint64(sqrt(float(n))) template isEven(n: uint64): bool = (n and 1) == 0
proc squfof(n: uint64): uint64 =
if n.isEven: return 2 var h = uint64(sqrt(float(n)) + 0.5) if h * h == n: return h
for m in M: if m > 1 and (n mod m == 0): return m # Check overflow m * n. if n > uint64.high div m: break let mn = m * n var r = isqrt(mn) if r * r > mn: dec r let rn = r
# Principal form. var b = r var a = 1u64 h = (rn + b) div a * a - b var c = (mn - h * h) div a
for i in 2..<(4 * isqrt(2 * r)):
# Search principal cycle.
swap a, c
var q = (rn + b) div a
let t = b
b = q * a - b
c += q * (t - b)
if i.isEven:
r = uint64(sqrt(float(c)) + 0.5)
if r * r == c: # Square form found?
# Inverse square root.
q = (rn - b) div r
var v = q * r + b
var w = (mn - v * v) div r
# Search ambiguous cycle.
var u = r
while true:
swap w, u
r = v
q = (rn + v) div u
v = q * u - v
if v == r: break
w += q * (r - v)
# Symmetry point.
h = gcd(n, u)
if h != 1: return h
result = 1
const Data = [2501u64,
12851u64,
13289u64,
75301u64,
120787u64,
967009u64,
997417u64,
7091569u64,
13290059u64,
42854447u64,
223553581u64,
2027651281u64,
11111111111u64,
100895598169u64,
1002742628021u64,
60012462237239u64,
287129523414791u64,
9007199254740931u64,
11111111111111111u64,
314159265358979323u64,
384307168202281507u64,
419244183493398773u64,
658812288346769681u64,
922337203685477563u64,
1000000000000000127u64,
1152921505680588799u64,
1537228672809128917u64,
4611686018427387877u64]
echo "N f N/f" echo "======================================" for n in Data:
let f = squfof(n)
let res = if f == 1: "fail" else: &"{f:<10} {n div f}"
echo &"{n:<22} {res}"</lang>
- Output:
N f N/f ====================================== 2501 61 41 12851 71 181 13289 97 137 75301 293 257 120787 43 2809 967009 1609 601 997417 257 3881 7091569 2663 2663 13290059 3119 4261 42854447 4423 9689 223553581 11213 19937 2027651281 44021 46061 11111111111 21649 513239 100895598169 112303 898423 1002742628021 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
Perl
<lang perl>use strict; use warnings; use feature 'say'; use ntheory <is_prime gcd forcomb vecprod>;
my @multiplier; my @p = <3 5 7 11>; forcomb { push @multiplier, vecprod @p[@_] } scalar @p;
sub sff {
my($N) = shift; return 1 if is_prime $N; # if n is prime return sqrt $N if sqrt($N) == int sqrt $N; # if n is a perfect square
for my $k (@multiplier) {
my $P0 = int sqrt($k*$N); # P[0]=floor(sqrt(N)
my $Q0 = 1; # Q[0]=1
my $Q = $k*$N - $P0**2; # Q[1]=N-P[0]^2 & Q[i]
my $P1 = $P0; # P[i-1] = P[0]
my $Q1 = $Q0; # Q[i-1] = Q[0]
my $P = 0; # P[i]
my $Qn = 0; # $P[$i+1];
my $b = 0; # b[i]
until (sqrt($Q) == int(sqrt($Q))) { # until Q[i] is a perfect square
$b = int( int(sqrt($k*$N) + $P1 ) / $Q); # floor(floor(sqrt(N+P[i-1])/Q[i])
$P = $b*$Q - $P1; # P[i]=b*Q[i]-P[i-1]
$Qn = $Q1 + $b*($P1 - $P); # Q[i+1]=Q[i-1]+b(P[i-1]-P[i])
($Q1, $Q, $P1) = ($Q, $Qn, $P);
}
$b = int( int( sqrt($k*$N)+$P ) / $Q ); # b=floor((floor(sqrt(N)+P[i])/Q[0])
$P1 = $b*$Q0 - $P; # P[i-1]=b*Q[0]-P[i]
$Q = ( $k*$N - $P1**2 )/$Q0; # Q[1]=(N-P[0]^2)/Q[0] & Q[i]
$Q1 = $Q0; # Q[i-1] = Q[0]
while () {
$b = int( int(sqrt($k*$N)+$P1 ) / $Q ); # b=floor(floor(sqrt(N)+P[i-1])/Q[i])
$P = $b*$Q - $P1; # P[i]=b*Q[i]-P[i-1]
$Qn = $Q1 + $b*($P1 - $P); # Q[i+1]=Q[i-1]+b(P[i-1]-P[i])
last if $P == $P1; # until P[i+1]=P[i]
($Q1, $Q, $P1) = ($Q, $Qn, $P);
}
for (gcd $N, $P) { return $_ if $_ != 1 and $_ != $N }
}
return 0
}
for my $data (
11111, 2501, 12851, 13289, 75301, 120787, 967009, 997417, 4558849, 7091569, 13290059,
42854447, 223553581, 2027651281, 11111111111, 100895598169, 1002742628021, 60012462237239,
287129523414791, 11111111111111111, 384307168202281507, 1000000000000000127, 9007199254740931,
922337203685477563, 314159265358979323, 1152921505680588799, 658812288346769681,
419244183493398773, 1537228672809128917) {
my $v = sff($data);
if ($v == 0) { say 'The number ' . $data . ' is not factored.' }
elsif ($v == 1) { say 'The number ' . $data . ' is a prime.' }
else { say "$data = " . join ' * ', sort {$a <=> $b} $v, int $data/int($v) }
}</lang>
- Output:
11111 = 41 * 271 2501 = 41 * 61 12851 = 71 * 181 13289 = 97 * 137 75301 = 257 * 293 120787 = 43 * 2809 967009 = 601 * 1609 997417 = 257 * 3881 4558849 = 383 * 11903 7091569 = 2663 * 2663 13290059 = 3119 * 4261 42854447 = 4423 * 9689 223553581 = 11213 * 19937 2027651281 = 44021 * 46061 11111111111 = 21649 * 513239 100895598169 = 112303 * 898423 The number 1002742628021 is a prime. 60012462237239 = 6862753 * 8744663 287129523414791 = 6059887 * 47381993 11111111111111111 = 2071723 * 5363222357 384307168202281507 = 415718707 * 924440401 1000000000000000127 = 111756107 * 8948056861 9007199254740931 = 10624181 * 847801751 922337203685477563 = 110075821 * 8379108103 314159265358979323 = 317213509 * 990371647 1152921505680588799 = 139001459 * 8294312261 658812288346769681 = 62222119 * 10588072199 419244183493398773 = 48009977 * 8732438749 1537228672809128917 = 26675843 * 57626245319
Phix
--requires(64) -- (decided to limit 32-bit explicitly instead) constant MxN = power(2,iff(machine_bits()=32?53:63)), m = {1, 3, 5, 7, 11} function squfof(atom N) -- square form factorization integer h, a=0, b, c, u=0, v, w, rN, q, r, t if remainder(N,2)==0 then return 2 end if h = floor(sqrt(N) + 0.5) if h*h==N then return h end if for k=1 to length(m) do integer mk = m[k] if mk>1 and remainder(N,mk)==0 then return mk end if //check overflow m * N if N>MxN/mk then exit end if atom mN = N*mk r = floor(sqrt(mN)) if r*r>mN then r -= 1 end if rN = r //principal form {b,a} = {r,1} h = floor((rN+b)/a)*a-b c = floor((mN-h*h)/a) for i=2 to floor(sqrt(2*r)) * 4-1 do //search principal cycle {a,c,t} = {c,a,b} q = floor((rN+b)/a) b = q*a-b c += q*(t-b) if remainder(i,2)==0 then r = floor(sqrt(c)+0.5) if r*r==c then //square form found //inverse square root q = floor((rN-b)/r) v = q*r+b w = floor((mN-v*v)/r) //search ambiguous cycle u = r while true do {u,w,r} = {w,u,v} q = floor((rN+v)/u) v = q*u-v if v==r then exit end if w += q*(r-v) end while //symmetry point h = gcd(N,u