Rosetta Code

Square form factorization: Difference between revisions

Page Discussion

Square form factorization (view source)

Revision as of 16:34, 20 March 2021

2,989 bytes added ,  3 years ago
See Discussion
(Added Go)
(See Discussion)
Line 41:
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.
and the program will duly report failure.
 
;Reference.
Line 52 ⟶ 51:
 
=={{header|C}}==
===Based on Wikipedia===
From the Wikipedia entry for Shanks's square forms factorization [[//en.wikipedia.org/wiki/Shanks%27s_square_forms_factorization.]]
<lang c>#include <math.h>
Line 195:
</pre>
 
===Classical heuristic===
See Discussion page.
<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 = a; a = c; c = 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 (1) {
if (m[k] < 1) break;
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)
 
ix = floor(sqrt(2*r)) * 4;
//iteration bound
 
for (i = 2; i < ix; i++) {
//search principal cycle
 
rho(a, b, c)
if ((i & 1) == 0) {
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;
}
}
}
}
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>
 
{{out|Showing problem cases only}}
<pre>
...
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
...
</pre>
 
=={{header|FreeBASIC}}==
<lang freebasic>' ***********************************************
' ***********************************************
'subject: Shanks's square form factorization:
' ambiguous forms of discriminant 4N
Line 350 ⟶ 517:
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 N = h * h = N then
'N is square
rp->f = culng(h)
Line 363 ⟶ 535:
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 h > (MxN \ m) then exit sub
Line 402 ⟶ 579:
if P.issqr(r) then
'square form found
 
if P.c = 1 then exit for
'cycle too short
 
if Q.pro(P, r) then
Line 475 ⟶ 649:
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 as ulongint f
dim t as integer
 
Line 502 ⟶ 676:
print "N = "; N
 
fflag = iif(N and 1, 1, 2)0
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
 
iffor ft = 1 thento tx
flagh(t) = 0threadcreate(@squfof, @a(t))
next t
 
squfof(@a(0))
for t = 1 to tx
h(t) = threadcreate(@squfof, @a(t))
next t
 
f = squfof(@a(0)).f
for t = 1 to tx
threadwait(h(t))
if f = 1 then f = a(t).f
next t
 
'assert
f = a(0).f
if N mod forf tthen f = 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
Line 540 ⟶ 701:
 
print "total time:"; csng(timer - tim); " s"
end</lang>
</lang>
 
{{out|Examples}}
<pre>
N = 2501
f = 4161 N/f = 6141
 
N = 12851
f = 18171 N/f = 71181
 
N = 13289
Line 629 ⟶ 789:
f = 343242169 N/f = 13435662733
 
total time: 0.02295520227905 s
</pre>
 
Retrieved from "https://rosettacode.org/wiki/Special:MobileDiff/309034"