Rosetta Code

Square form factorization: Difference between revisions

Page Discussion

Square form factorization (view source)

Revision as of 16:05, 9 February 2024

74,872 bytes added ,  3 months ago
m
→‎{{header|Wren}}: Minor tidy
(Added Go)
m (→‎{{header|Wren}}: Minor tidy)
 
(35 intermediate revisions by 14 users not shown)
Line 1:
{{task|Mathematics}}
[[Category:Prime Numbers]]
 
 
;Task.
 
[[wp:Daniel_Shanks|Daniel Shanks]]'s Square Form Factorization [[wp:Shanks%27s_square_forms_factorization|(SquFoF)]].
 
Line 41 ⟶ 40:
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 50 ⟶ 48:
 
__TOC__
 
=={{header|ALGOL 68}}==
{{works with|ALGOL 68G|Any - tested with release 2.8.3.win32}}
Assumes LONG INT is long enough - this is true in ALGOL 68G versioins 2 and 3.
{{Trans|Wren|...but only using a single size of integer}}
<syntaxhighlight lang="algol68">
BEGIN # Daniel Shanks's Square Form Factorization (SquFoF) - based on the Wren sample #
 
MODE INTEGER = LONG INT; # large enough INT type #
PROC(LONG REAL)LONG REAL size sqrt = long sqrt; # sqrt for INTEGER values #
 
[]INTEGER multipliers = ( 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
);
PROC gcd = ( INTEGER x, y )INTEGER: # iterative gcd #
BEGIN
INTEGER a := x, b := y;
WHILE b /= 0 DO
INTEGER next a = b;
b := a MOD b;
a := next a
OD;
ABS a
END # gcd # ;
 
PROC squfof = ( INTEGER n )INTEGER:
IF INTEGER s = ENTIER ( size sqrt( n ) + 0.5 );
s * s = n
THEN s
ELSE INTEGER result := 0;
FOR multiplier FROM LWB multipliers TO UPB multipliers WHILE result = 0 DO
INTEGER d = n * multipliers[ multiplier ];
INTEGER pp := ENTIER size sqrt( d );
INTEGER p prev := pp;
INTEGER po = p prev;
INTEGER q prev := 1;
INTEGER qq := d - ( po * po );
INTEGER l = ENTIER size sqrt( s * 8 );
INTEGER bb = 3 * l;
INTEGER i := 2;
INTEGER b := 0;
INTEGER q := 0;
INTEGER r := 0;
BOOL again := TRUE;
WHILE i < bb AND again DO
b := ( po + pp ) OVER qq;
pp := ( b * qq ) - pp;
q := qq;
qq := q prev + ( b * ( p prev - pp ) );
r := ENTIER ( size sqrt( qq ) + 0.5 );
IF i MOD 2 = 0 THEN again := r * r /= qq FI;
IF again THEN
q prev := q;
p prev := pp;
i +:= 1
FI
OD;
IF i < bb THEN
b := ( po - pp ) OVER r;
p prev := pp := ( b * r ) + pp;
q prev := r;
qq := ( d - ( p prev * p prev ) ) OVER q prev;
i := 0;
WHILE
b := ( po + pp ) OVER qq;
p prev := pp;
pp := ( b * qq ) - pp;
q := qq;
qq := q prev + ( b * ( p prev - pp ) );
q prev := q;
i +:= 1;
pp /= p prev
DO SKIP OD
FI;
r := gcd( n, q prev );
IF r /= 1 AND r /=n THEN result := r FI
OD;
result
FI # squfof # ;
 
[]INTEGER 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
);
 
print( ( "Integer Factor Quotient", newline ) );
print( ( "----------------------------------------", newline ) );
FOR example FROM LWB examples TO UPB examples DO
INTEGER n = examples[ example ];
INTEGER fact = squfof( n );
STRING quot = IF fact = 0 THEN "fail" ELSE whole( n OVER fact, 0 ) FI;
print( ( whole( n, -20 ), " ", whole( fact, -10 ), " ", quot, newline ) )
OD
END
</syntaxhighlight>
{{out}}
<pre>
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
</pre>
 
=={{header|C}}==
===Based on Wp===
From the Wikipedia entry for Shanks's square forms factorization [[//en.wikipedia.org/wiki/Shanks%27s_square_forms_factorization.]]
<langsyntaxhighlight lang="c">#include <math.h>
#include <stdio.h>
 
Line 163 ⟶ 302:
}
}
</langsyntaxhighlight>{{out}}
<pre>
Integer 2501 has factors 61 and 41
Line 195 ⟶ 334:
</pre>
 
===Classical heuristic===
See [[Talk:Square_Form_Factorization|Discussion]].
<syntaxhighlight 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);
}
}</syntaxhighlight>
 
{{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|C++}}==
<syntaxhighlight lang="c++">
#include <cmath>
#include <cstdint>
#include <iostream>
#include <numeric>
#include <random>
 
uint64_t test_value = 0;
uint64_t sqrt_test_value = 0;
 
class BQF { // Binary quadratic form
public:
BQF(const uint64_t& a, const uint64_t& b, const uint64_t& c) : a(a), b(b), c(c) {
q = ( sqrt_test_value + b ) / c;
bb = q * c - b;
}
 
BQF rho() {
return BQF(c, bb, a + q * ( b - bb ));
}
 
BQF rho_inverse() {
return BQF(c, bb, ( test_value - bb * bb ) / c);
}
 
uint64_t a, b, c;
private:
uint64_t q, bb;
};
 
uint64_t squfof(const uint64_t& number) {
const uint32_t sqrt = std::sqrt(number);
if ( sqrt * sqrt == number ) {
return sqrt;
}
 
test_value = number;
sqrt_test_value = std::sqrt(test_value);
 
// Principal form
BQF form(0, sqrt_test_value, 1);
form = form.rho_inverse();
 
// Search principal cycle
for ( uint32_t i = 0; i < 4 * std::sqrt(2 * sqrt_test_value); i += 2 ) {
// Even step
form = form.rho();
 
uint64_t sqrt_c = std::sqrt(form.c);
if ( sqrt_c * sqrt_c == form.c ) { // Square form found
// Inverse square root
BQF form_inverse(0, -form.b, sqrt_c);
form_inverse = form_inverse.rho_inverse();
 
// Search ambiguous cycle
uint64_t previous_b = 0;
do {
previous_b = form_inverse.b;
form_inverse = form_inverse.rho();
} while ( form_inverse.b != previous_b );
 
// Symmetry point
const uint64_t g = std::gcd(number, form_inverse.a);
if ( g != 1 ) {
return g;
}
}
 
// Odd step
form = form.rho();
}
 
if ( number % 2 == 0 ) {
return 2;
}
return 0; // Failed to factorise, possibly a prime number
}
 
int main() {
std::random_device random;
std::mt19937 generator(random());
const uint64_t lower_limit = 100'000'000'000'000'000;
std::uniform_int_distribution<uint64_t> distribution(lower_limit, 10 * lower_limit);
 
for ( uint32_t i = 0; i < 20; ++i ) {
uint64_t test = distribution(random);
uint64_t factor = squfof(test);
 
if ( factor == 0 ) {
std::cout << test << " - failed to factorise" << std::endl;
} else {
std::cout << test << " = " << factor << " * " << test / factor << std::endl;
}
std::cout << std::endl;
}
}
</syntaxhighlight>
{{ out }}
<pre>
822140815871714649 = 141 * 5830785928168189
 
473377979025428817 = 3 * 157792659675142939
 
482452941918160803 = 4410431 * 109389069213
 
165380937127655630 = 65438 * 2527292049385
 
191677853606692475 = 7589219 * 25256598025
 
480551815975206727 = 2843 * 169029833265989
 
178710207362206205 = 5 * 35742041472441241
 
484660189375949842 = 1094 * 443016626486243
 
758704390319635770 = 1605 * 472713015775474
 
820453356193182720 = 97280 * 8433936638499
 
706982627912630220 = 121273 * 5829678724140
 
614913973550671312 = 437204432 * 1406467841
 
601482456081568543 = 131 * 4591469130393653
 
610533314488947626 = 14 * 43609522463496259
 
336343281182924332 = 70108 * 4797502156429
 
308127213282933401 = 7 * 44018173326133343
 
582455924775519843 = 3 * 194151974925173281
 
694215100094443276 = 32070628 * 21646445467
 
398821795604697523 = 181 * 2203435334832583
 
477964959783291032 = 517029608 * 924444079
</pre>
 
=={{header|EasyLang}}==
{{trans|C}}
<syntaxhighlight lang=easylang>
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 ]
func gcd a b .
while b <> 0
a = a mod b
swap a b
.
return a
.
func squfof N .
s = floor (sqrt N + 0.5)
if s * s = N
return s
.
for multiplier in multiplier[]
if N > 9007199254740992 / multiplier
print "Number " & N & " is too big"
break 1
.
D = multiplier * N
P = floor sqrt D
Po = P
Pprev = P
Qprev = 1
Q = D - Po * Po
L = 2 * floor sqrt (2 * s)
B = 3 * L
for i = 2 to B - 1
b = (Po + P) div Q
P = b * Q - P
q = Q
Q = Qprev + b * (Pprev - P)
r = floor (sqrt Q + 0.5)
if i mod 2 = 0 and r * r = Q
break 1
.
Qprev = q
Pprev = P
.
if i < B
b = (Po - P) div r
P = b * r + P
Pprev = P
Qprev = r
Q = (D - Pprev * Pprev) / Qprev
i = 0
repeat
b = (Po + P) div Q
Pprev = P
P = b * Q - P
q = Q
Q = Qprev + b * (Pprev - P)
Qprev = q
i += 1
until P = Pprev
.
r = gcd N Qprev
if r <> 1 and r <> N
return r
.
.
.
return 0
.
data[] = [ 2501 12851 13289 75301 120787 967009 997417 7091569 13290059 42854447 223553581 2027651281 11111111111 100895598169 1002742628021 60012462237239 287129523414791 9007199254740931 ]
for example in data[]
factor = squfof example
if factor = 0
print example & " was not factored."
else
quotient = example / factor
print example & " has factors " & factor & " " & quotient
.
.
</syntaxhighlight>
 
=={{header|FreeBASIC}}==
<syntaxhighlight lang="freebasic">' ***********************************************
<lang freebasic>
' ***********************************************
'subject: Shanks's square form factorization:
' ambiguous forms of discriminant 4N
' give factors of N.
'tested : FreeBasic 1.0708.1
 
 
const qx = 50
'queue size
 
'------------------------------------------------
const MxN = culngint(1) shl 62
'input maximum
 
const qx = (1 shl 5) - 1
'queue size
 
type arg
Line 219 ⟶ 743:
 
type bqf
declare sub rho (byval sw as integer)
'reduce indefinite form, set sw = 0 to initialize
declare function issqrissq (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
Line 237 ⟶ 760:
'return -1 if a proper square form is found
 
as integerulong ta(qx), L, m
as ulonginteger a(qx + 1), Lk, mt
end type
 
Line 247 ⟶ 770:
dim shared flag as integer
'signal to end all threads
 
dim shared as ubyte q1024(1023), q3465(3464)
'quadratic residue tables
 
 
'------------------------------------------------
sub bqf.rho (byval sw as integer)
dim as ulong q, t = b
swap a, c
'residue
q = culng(rN + b) \ a
t = b: b = q * a - b
'pseudo-square
if sw then
c += q * (t 'pseudo-square b)
c += q * (t - b)
else
'initialize
c = (mN - culngint(b) * b) \ a
end if
end sub
 
'initialize form
function bqf.issqr (byref r as ulong) as integer
#macro rhoin(F)
static as integer q64(63), q63(62), q55(54), sw = 0
F.rho : h = F.b
if sw = 0 then
swF.c = (mN -1 h * h) \ F.a
#endmacro
'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
 
function bqf.issq (byref r as ulong) as integer
if q64(c and 63) then
if q1024(c and if1023) andalso q63q3465(c mod 633465) then
'98.6% non-squares filtered
if q55(c mod 55) then
r = culng(sqr(c))
'>98% non-squares filtered
if r * r = culng(sqr(c)) then return -1
if c = r * r then return -1
end if
end if
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 d, u = a, v = b, w = c
dim as longint ifu vb= a, v = 0b, thenw exit= subc
 
'{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
Line 319 ⟶ 822:
 
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'circular += 2queue
k = (k + 2) and qx
if k > t then t = k
'enqueue P.b, P.c
a(tk) = P.b mod s
a(tk + 1) = s
end if
end sub
Line 346 ⟶ 850:
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 N = h * h = N then
'N is square
rp->f = culng(h)
Line 358 ⟶ 867:
end if
 
h = N
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 hN > (MxN \ m) then exit sub
h *= m
end if
mN = N * m
 
r = int(sqr(mN))
P.mN = h
A.mN = h
r = int(sqr(h))
'float64 fix
if culngint(r) * r > hmN then r -= 1
P.rN = r
A.rN = r
Line 381 ⟶ 892:
if P.vb then print "r = "; r
 
Q.tk = -2: Q.t = -1: Q.m = m
'Queue entry bounds
Q.L = int(sqr(r * 2))
Line 388 ⟶ 899:
'principal form
P.b = r: P.c = 1
P.rhorhoin(0P)
P.qform("P", 1)
 
Line 397 ⟶ 908:
if P.c < L2 then Q.enq(P)
 
P.rho(1)
if (i and 1) = 0 andalso P.issq(r) then
'square form found
 
if PQ.issqrpro(P, r) then
'square form found
 
if P.c = 1 then exitqform("P", fori)
'cycleinverse toosquare shortroot
A.b =-P.b: A.c = r
rhoin(A): j = 1
A.qform("A", j)
 
if Q.pro(P, r) thendo
'search ambiguous cycle
t = A.b
A.rho: j += 1
 
Pif A.qform("P",b = t i)then
'inverse square root 'symmetry point
A.b =-P.b: A.c =qform("A", rj)
A.rho(0): j = 1red(f, A.a)
A.qform("A", j) if f = 1 then exit do
 
do flag = -1
'searchfactor ambiguous cyclefound
end t = A.bif
loop until A.rho(1): j += 1flag
 
end if A.b =' tproper thensquare
end if ' square form
'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
Line 450 ⟶ 955:
data 7091569
data 13290059
data 23515517
data 42854447
data 223553581
Line 475 ⟶ 981:
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 integerulongint 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
Line 502 ⟶ 1,015:
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 + 1 step 2
flagif =t 0< tx then
 
for t = 1 to tx
h(t) = threadcreate(@squfof, @a(t))
nextend tif
 
squfof(@a(0t - 1))
f = a(t - 1).f
 
fif =t a(0).f< tx then
for t = 1 to tx
threadwait(h(t))
if f = 1 then f = a(t).f
nextend tif
 
'assertif f > 1 then exit for
next t
if N mod f then f = 1
 
'assert
elseif f = N then
if N mod 'smallf primethen Nf = 1
f = 1
end if
 
if f = 1 then
Line 540 ⟶ 1,044:
 
print "total time:"; csng(timer - tim); " s"
end</syntaxhighlight>
end
</lang>
 
{{out|Examples}}
<pre>
N = 2501
f = 4161 N/f = 6141
 
N = 12851
f = 18171 N/f = 71181
 
N = 13289
Line 571 ⟶ 1,074:
N = 13290059
f = 3119 N/f = 4261
 
N = 23515517
f = 53 N/f = 443689
 
N = 42854447
Line 629 ⟶ 1,135:
f = 343242169 N/f = 13435662733
 
total time: 0.02295520170462 s
</pre>
 
Line 636 ⟶ 1,142:
 
Rather than juggling with big.Int, I've just allowed the two 'awkward' cases to fail.
<langsyntaxhighlight lang="go">package main
 
import (
Line 762 ⟶ 1,268:
fmt.Printf("%-20d %-10d %s\n", N, fact, quot)
}
}</langsyntaxhighlight>
 
{{out}}
Line 796 ⟶ 1,302:
1537228672809128917 0 fail
4611686018427387877 343242169 13435662733
</pre>
 
=={{header|J}}==
{{eff note|J|q:}}
J does not have an unsigned fixed width integer type, which is one of the reasons that (in J) this algorithm is less optimal than advertised:<syntaxhighlight lang="j">sqff=: {{
s=. <.%:y
if. y=*:s do. s return. end.
for_D. (x:y)*/:~*/@>,{1,each}.p:i.5 do.
if. -.'integer'-:datatype D=. x:inv D do. break. end.
P=. <.%:D
Q=. 1, D-P*P
lim=. <:6*<.%:2*s
for_i. }.i.lim do.
b=. <.(+/0 _1{P)%{:Q
P=. P,|(b*{:Q)-{:P
Q=. Q,|(_2{Q)+b*-/_2{.P
if. 2|i do. if. (=<.&.%:){:Q do. break. end. end.
end.
if. i>:lim do. continue. end.
Q=. <.%:{:Q
b=. <.(-/0 _1{P)%Q
P=. ,(b*Q)+{:P
Q=. Q, <.|(D-*:P)%Q
whilst. ~:/_2{.P do.
b=. <.(+/0 _1{P)%{:Q
P=. P,|(b*{:Q)-{:P
Q=. Q,|(_2{Q)+b*-/_2{.P
end.
f=. y+.x:_2{Q
if. -. f e. 1,y do. f return. end.
end.
1
}}</syntaxhighlight>
 
Task examples:<syntaxhighlight lang="j"> task ''
2501: 61 * 41
12851: 71 * 181
13289: 137 * 97
75301: 293 * 257
120787: 43 * 2809
967009: 601 * 1609
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 was not factored
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 was not factored
1152921505680588799: 139001459 * 8294312261
1537228672809128917 was not factored
4611686018427387877 was not factored</syntaxhighlight> where <syntaxhighlight lang="j">task=: {{
for_num. nums do.
factor=. x:sqff num
if. 1=factor do. echo num,&":' was not factored'
else. echo num,&":': ',factor,&":' * ',":x:num%factor
end.
end.
}}
 
nums=: ".{{)n
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
4611686018427387877x
}}-.LF</syntaxhighlight>
 
=={{header|Java}}==
<syntaxhighlight lang ="java">
import java.math.BigInteger;
import java.util.List;
import java.util.concurrent.ThreadLocalRandom;
 
public final class SquareFormFactorization {
 
public static void main(String[] args) {
ThreadLocalRandom random = ThreadLocalRandom.current();
final long lowerLimit = 10_000_000_000_000_000L;
final List<Long> tests = random.longs(20, lowerLimit, 10 * lowerLimit).boxed().toList();
 
for ( long test : tests ) {
long factor = squfof(test);
 
if ( factor == 0 ) {
System.out.println(test + " - failed to factorise");
} else if ( factor == 1 ) {
System.out.println(test + " is a prime number");
} else {
System.out.println(test + " = " + factor + " * " + test / factor);
}
System.out.println();
}
}
private static long squfof(long number) {
if ( BigInteger.valueOf(number).isProbablePrime(15) ) {
return 1; // Prime number
}
final int sqrt = (int) Math.sqrt(number);
if ( sqrt * sqrt == number ) {
return sqrt;
}
testValue = number;
sqrtTestValue = (long) Math.sqrt(testValue);
 
// Principal form
BQF form = new BQF(0, sqrtTestValue, 1);
form = form.rhoInverse();
 
// Search principal cycle
for ( int i = 0; i < 4 * (long) Math.sqrt(2 * sqrtTestValue); i += 2 ) {
// Even step
form = form.rho();
 
long sqrtC = (long) Math.sqrt(form.c);
if ( sqrtC * sqrtC == form.c ) { // Square form found
// Inverse square root
BQF formInverse = new BQF(0, -form.b, sqrtC);
formInverse = formInverse.rhoInverse();
 
// Search ambiguous cycle
long previousB = 0;
do {
previousB = formInverse.b;
formInverse = formInverse.rho();
} while ( formInverse.b != previousB );
// Symmetry point
final long gcd = gcd(number, formInverse.a);
if ( gcd != 1 ) {
return gcd;
}
}
// Odd step
form = form.rho();
}
if ( number % 2 == 0 ) {
return 2;
}
return 0; // Failed to factorise
}
 
private static long gcd(long a, long b) {
while ( b != 0 ) {
long temp = a; a = b; b = temp % b;
}
return a;
}
private static class BQF { // Binary quadratic form
public BQF(long aA, long aB, long aC) {
a = aA; b = aB; c = aC;
q = ( sqrtTestValue + b ) / c;
bb = q * c - b;
}
public BQF rho() {
return new BQF(c, bb, a + q * ( b - bb ));
}
public BQF rhoInverse() {
return new BQF(c, bb, ( testValue - bb * bb ) / c);
}
private long a, b, c;
private long q, bb;
}
private static long testValue, sqrtTestValue;
 
}
</syntaxhighlight>
{{ out }}
<pre>
20096060843736547 = 433 * 46411225967059
 
24628423963378844 = 7 * 3518346280482692
 
68276045265502398 = 37 * 1845298520689254
 
61072103663732497 = 8477 * 7204447760261
 
63462639942509072 = 16 * 3966414996406817
 
60313009405143787 = 89288189 * 675486983
 
76093594148871700 = 377900 * 201359074223
 
31796652636180617 is a prime number
 
87047981623879461 = 243 * 358222146600327
 
71567116631895554 = 73 * 980371460710898
 
50852012325831410 = 2 * 25426006162915705
 
65816967116185802 = 131280559 * 501345878
 
89627452852493643 = 31 * 2891208156532053
 
41735751565855318 = 10004047 * 4171886794
 
97291513005945602 = 2 * 48645756502972801
 
88974788272758998 = 59 * 1508047258860322
 
53903340306287681 = 21727 * 2480938017503
 
10811459482792395 = 546427 * 19785734385
 
95115727966103864 = 26105228 * 3643550938
 
11340988571009785 = 5 * 2268197714201957
</pre>
 
=={{header|jq}}==
'''Adapted from [[#Wren|Wren]]'''
'''Works with gojq, the Go implementation of jq'''
 
gojq has support for unbounded-precision
integer arithmetic and accordingly the output shown below is from a run thereof;
the C implementation of jq produces correct results up to and including
[287129523414791,6059887,47381993].
 
'''Preliminaries'''
<syntaxhighlight lang="jq">def gcd(a; b):
# subfunction expects [a,b] as input
# i.e. a ~ .[0] and b ~ .[1]
def rgcd: if .[1] == 0 then .[0]
else [.[1], .[0] % .[1]] | rgcd
end;
[a,b] | rgcd;
 
# for infinite precision integer-arithmetic
def idivide($p; $q): ($p - ($p % $q)) / $q ;
def idivide($q): (. - (. % $q)) / $q ;
 
def isqrt:
def irt:
. as $x
| 1 | until(. > $x; . * 4) as $q
| {$q, $x, r: 0}
| until( .q <= 1;
.q |= idivide(4)
| .t = .x - .r - .q
| .r |= idivide(2)
| if .t >= 0
then .x = .t
| .r += .q
else .
end)
| .r ;
if type == "number" and (isinfinite|not) and (isnan|not) and . >= 0
then irt
else "isqrt requires a non-negative integer for accuracy" | error
end ;</syntaxhighlight>
'''The Tasks'''
<syntaxhighlight lang="jq">def multipliers:
[
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
];
 
# input should be a number
def squfof:
def toi : floor | tostring | tonumber;
. as $N
| (($N|sqrt + 0.5)|toi) as $s
| if ($s*$s == $N) then $s
else label $out
| {}
| multipliers[] as $multiplier
| ($N * $multiplier) as $D
| .P = ($D|isqrt)
| .Pprev = .P
| .Pprev as $Po
| .Qprev = 1
| .Q = $D - $Po*$Po
| (($s * 8)|isqrt) as $L
| (3 * $L) as $B
| .i = 2
| .b = 0
| .q = 0
| .r = 0
| .stop = false
| until( (.i >= $B) or .stop;
.b = idivide($Po + .P; .Q)
| .P = .b * .Q - .P
| .q = .Q
| .Q = .Qprev + .b * (.Pprev - .P)
 
| .r = (((.Q|isqrt) + 0.5)|toi)
 
| if ((.i % 2) == 0 and (.r*.r) == .Q) then .stop = true
else
.Qprev = .q
| .Pprev = .P
| .i += 1
end )
| if .i < $B
then
.b = idivide($Po - .P; .r)
| .P = .b*.r + .P
| .Pprev = .P
| .Qprev = .r
| .Q = idivide($D - .Pprev*.Pprev; .Qprev)
| .i = 0
| .stop = false
| until (.stop;
.b = idivide($Po + .P; .Q)
| .Pprev = .P
| .P = .b * .Q - .P
| .q = .Q
| .Q = .Qprev + .b * (.Pprev - .P)
| .Qprev = .q
| .i += 1
| if (.P == .Pprev) then .stop = true else . end )
| .r = gcd($N; .Qprev)
| if .r != 1 and .r != $N then .r, break $out else empty end
else empty
end
end
// 0 ;
 
def 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"
];
"[Integer, Factor, Quotient]"
"---------------------------",
(examples[] as $example
| ($example|tonumber) as $N
| ($N | squfof) as $fact
| if $fact == 0 then "fail"
else idivide($N; $fact) as $quot
| [$N, $fact, $quot]
end
)</syntaxhighlight>
{{out}}
<pre>
[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]
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]
</pre>
 
=={{header|Julia}}==
Modified from Wikipedia's article at [[https://en.wikipedia.org/wiki/Shanks%27s_square_forms_factorization]]
<langsyntaxhighlight 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)))
Line 856 ⟶ 1,800:
println(factr == 0 ? "fail" : n ÷ factr)
end
</langsyntaxhighlight>{{out}}
<pre>
Integer Factor Quotient
Line 891 ⟶ 1,835:
</pre>
 
=={{header|PhixNim}}==
{{trans|Julia|and wikipediaPhix}}
<syntaxhighlight lang="nim">import math, strformat
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))
 
const M = [uint64 1, 3, 5, 7, 11]
function square_form_factor(atom n)
 
atom s = round(sqrt(n))
template isqrt(n: uint64): uint64 = uint64(sqrt(float(n)))
if s * s == n then return s end if
template isEven(n: uint64): bool = (n and 1) == 0
for mi=1 to length(multiplier) do
 
atom k = multiplier[mi]
proc squfof(n: uint64): uint64 =
if n > INT_MAX/k then exit end if
 
atom d = k * n,
if n.isEven: return 2
p0 = floor(sqrt(d)),
var h = uint64(sqrt(float(n)) + 0.5)
pprev = p0, p = p0,
if h * h == n: return qprev = 1,h
 
Q = d - p0 * p0,
for m in M:
l = floor(2 * sqrt(2 * s)),
if m > 1 and (n mod m == B0): =return 3*l,m
# Check overflow m * i = 2, r, b, qn.
if n > uint64.high whilediv im: < B dobreak
let bmn = floor((p0 + p)m /* Q)n
var pr = b * Q - pisqrt(mn)
if r * r > mn: dec q = Qr
let rn = r
Q = qprev + b * (pprev - p)
 
r = round(sqrt(Q))
# Principal form.
if remainder(i,2)!=1 and r * r == Q then exit end if
var qprevb = qr
var ppreva = p1u64
h = (rn + b) div a * ia +=- 1b
var c = (mn end- whileh * h) div a
 
if i < B then
for i in 2..<(4 * b = floorisqrt((p02 -* pr) / r):
# Search principal p = b * r + pcycle.
swap a, pprev = pc
var q = (rn + b) qprev =div ra
let t = b
Q = floor((d - pprev * pprev) / qprev)
b = q * a - i = 0b
c += q * (t - while true dob)
 
b = floor((p0 + p) / Q)
if pprev = pi.isEven:
pr = b * Quint64(sqrt(float(c)) -+ p0.5)
if r * r == c: # Square q =form Qfound?
 
Q = qprev + b * (pprev - p)
# Inverse square qprev = qroot.
q = (rn - b) div i += 1r
var v if p == pprevq then* exitr end+ ifb
var w end= while(mn - v * v) div r
 
r = gcd(n, qprev)
# Search ambiguous cycle.
if r != 1 and r != n then return r end if
end if var u = r
end for while true:
swap w, u
return 0
r = v
end function
q = (rn + v) div u
v = q * u - v
constant tests = {2501, 12851, 13289, 75301, 120787, 967009, 997417, 7091569, 13290059, 42854447, 223553581,
if v == r: break
2027651281, 11111111111, 100895598169, 1002742628021, 60012462237239, 287129523414791,
w += q * 0,(r -- (limit of 32-bitv)
 
9007199254740931, 11111111111111111, 314159265358979323, 384307168202281507, 419244183493398773,
# Symmetry point.
658812288346769681, 922337203685477563, 1000000000000000127, 1152921505680588799,
h = 1537228672809128917gcd(n, 4611686018427387877}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}"</syntaxhighlight>
 
{{out}}
printf(1,"Integer Factor Quotient\n")
<pre>N f N/f
printf(1,"======================================\n")
======================================
for i=1 to length(tests) do
2501 61 41
atom n = tests[i],
12851 factr = square_form_factor(n) 71 181
13289 97 137
string f = iff(factr=0?"fail":sprintf("%d",n/factr))
75301 293 257
printf(1,"%-22d %-10d %s\n",{n,factr,f})
120787 43 2809
if machine_bits()=32 and n=0 then exit end if
967009 1609 601
end for</lang>
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</pre>
 
=={{header|Pascal}}==
{{works with|Free Pascal}}
A console program written in Free Pascal. The code is based on:
Jason E. Gower and Samuel S. Wagstaff, jr., "Square form factorization",
Mathematics of Computation Volume 77, Number 261, January 2008, Pages 551–588
S 0025-5718(07)02010-8
Article electronically published on May 14, 2007
 
I'm not sure whether this is the same as reference [1] in the task description; the words "a detailed analysis of SquFoF" appear in the abstract.
 
The Pascal program includes some small changes to the Gower-Wagstaff algorithm, as noted in the comments. The output shows the successful multiplier (if any) and whether the factor was found in the main or preliminary part of the algorithm.
 
Further to the advice (on the Discussion page) not to use the Wikipedia version of the algorithm, I tested the Gower-Wagstaff version for all odd composite numbers less than 10^9, and it found a factor in every case. The Wikipedia version failed in 790 cases.
<syntaxhighlight lang="pascal">
program SquFoF_console;
 
{$mode objfpc}{$H+}
 
uses SquFoF_utils;
 
type TResultKind =
(rkPrelim, // a factor was found by the preliminary routine
rkMain, // a factor was found by the main algorithm
rkFail); // no factor was found
type TAlgoResult = record
Kind : TResultKind;
Mult : word;
Factor : UInt64;
end;
 
// Preliminary to G-W algorithm. Returns D and S of the algorithm.
// Also returns a non-trivial factor if found, else returns factor = 1.
procedure GWPrelim( N : UInt64; // number to be factorized
m : word; // multiplier
out D, S, factor : UInt64);
var
sqflag : boolean;
begin
D := m*N;
sqflag := SquFoF_utils.IsSquare( D, S);
if m = 1 then
if sqflag then factor := S
else factor := 1
else
factor := GCD( N,m);
end;
 
// Tries to factorize N by applying Gower-Wagstaff alsorithm to m*N.
function GW_with_multiplier( N : UInt64;
m : word) : TAlgoResult;
var
D, S, P, P_prev, Q, L, B: Uint64;
r : UInt64;
i, j, k : integer;
f, g : UInt64;
sqrtD : double;
endCycle : boolean;
// Queue is not much used, so make it a simple array.
type TQueueItem = record
Left, Right : UInt64;
end;
const QUEUE_CAPACITY = 50; // as suggested by Gower & Wagstaff
var
queue : array [0..QUEUE_CAPACITY - 1] of TQueueItem;
queueCount : integer;
begin
result.Mult := m;
 
// Filter out special cases (differs from Gower & Wagstaff). Note:
// (1) multiplier m is assumed to be squarefree;
// (2) if we proceed to the main algorithm, mN must not be square
// (otherwise Q = 0 and division by Q causes an error).
GWPrelim( N, m, {out} D, S, f);
if f > 1 then begin
result.Kind := rkPrelim;
result.Factor := f;
exit;
end;
// Not special, proceed to main algorithm
result.Kind := rkMain;
result.Mult := m;
result.Factor := 1;
queueCount := 0; // Clear queue
P := S;
Q := 1;
i := -1; // keep i same as in G & W; algo fails if i > B
sqrtD := SquFoF_utils.FSqrt( D);
// L := Trunc( 2.0*Sqrt( 2.0*sqrtD));
L := Trunc( Sqrt( 2.0*sqrtD)); // as in Section 5.2 of G&W paper
B := 2*L;
 
// Start forward cycle
endCycle := false;
while not endCycle do begin
// We update Q here, P at the end of the loop
Q := (D - P*P) div Q;
if (not Odd(i)) and SquFoF_utils.IsSquare( Q, r) then begin
// Q is square for even i.
// Possibly (?probably?) ends the forward cycle,
// but we need to inspect the queue first.
endCycle := true;
j := queueCount; // working backwards down the queue
if r = 1 then begin // the method may fail
while (j > 0) and (result.Kind <> rkFail) do begin
dec( j);
if queue[j].Left = 1 then result.Kind := rkFail;
end;
if result.Kind = rkFail then exit;
end
else begin // if r > 1
while (j > 0) and (endCycle) do begin
dec( j);
if (queue[j].Left = r)
and ((P - queue[j].Right) mod r = 0) then begin
// Deleting up to the *last* item in the list that
// satisfies the condition, Should it be the *first* ?
// Delete queue items 0..j inclusive
inc(j); k := 0;
while j < queueCount do begin
queue[k] := queue[j];
inc(j); inc(k);
end;
queueCount := k;
endCycle := false;
end; // if
end;
end;
end; // if i even and Q square
if not endCycle then begin
g := Q div SquFoF_utils.GCD( Q, 2*m);
if g <= L then begin
if queueCount < QUEUE_CAPACITY then begin
with queue[queueCount] do begin
Left := g; Right := P mod g;
end;
inc( queueCount);
end
else begin // queue overflow, fail
result.Kind := rkFail;
exit;
end;
end;
inc(i);
if i > B then begin
result.Kind := rkFail;
exit;
end;
end;
P := S - ((S + P) mod Q);
end; // while not endCycle
Assert( (D - P*P) mod r = 0); // optional check
P := S - ((S + P) mod r);
Q := r;
// Start backward cycle
endCycle := false;
while not endCycle do begin
P_prev := P;
Q := (D - P*P) div Q;
P := S - ((S + P) mod q);
endCycle := (P = P_prev);
end; // while not endCycle
// Finished
result.Factor := Q div SquFoF_utils.GCD( Q, 2*m);
end;
 
const NR_RC_VALUES = 28;
RC_VALUES : array [0..NR_RC_VALUES - 1] of 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);
 
type TMultAndMaxN = record
Mult : word; // small multiplier
MaxN : UInt64; // maximum N for that multiplier (N*multiplier < 2^64)
end;
const NR_MULTIPLIERS = 16;
const MULTIPLIERS : array [0..NR_MULTIPLIERS - 1] of TMultAndMaxN =
((Mult: 1; MaxN: 18446744073709551615),
(Mult: 3; MaxN: 6148914691236517205),
(Mult: 5; MaxN: 3689348814741910323),
(Mult: 7; MaxN: 2635249153387078802),
(Mult: 11; MaxN: 1676976733973595601),
(Mult: 15; MaxN: 1229782938247303441),
(Mult: 21; MaxN: 878416384462359600),
(Mult: 33; MaxN: 558992244657865200),
(Mult: 35; MaxN: 527049830677415760),
(Mult: 55; MaxN: 335395346794719120),
(Mult: 77; MaxN: 239568104853370800),
(Mult: 105; MaxN: 175683276892471920),
(Mult: 165; MaxN: 111798448931573040),
(Mult: 231; MaxN: 79856034951123600),
(Mult: 385; MaxN: 47913620970674160),
(Mult: 1155; MaxN: 15971206990224720));
 
function GowerWagstaff( N : UInt64) : TAlgoResult;
var
j : integer;
begin
j := 0;
result.Kind := rkFail;
while (result.Kind = rkFail)
and (j < NR_MULTIPLIERS)
and (N <= MULTIPLIERS[j].MaxN) do
begin
result := GW_with_multiplier( N, MULTIPLIERS[j].Mult);
if result.Kind = rkFail then inc(j);
end;
end;
 
// Main program
var
j : integer;
ar : TAlgoResult;
kindStr : string;
N, cofactor : UInt64;
begin
WriteLn( ' Number Mult M/P Factorization');
for j := 0 to NR_RC_VALUES - 1 do begin
N := RC_VALUES[j];
ar := GowerWagstaff( N);
if ar.Kind = rkFail then
WriteLn( N:20, ' No factor found')
else begin
case ar.Kind of
rkPrelim: kindStr := 'Prelim';
rkMain : kindStr := 'Main ';
end;
cofactor := N div ar.Factor;
Assert( cofactor * ar.Factor = N); // check that all has gone well
WriteLn( N:20, ar.Mult:5, ' ',
kindStr:6, ' ', ar.Factor, ' * ', cofactor);
end;
end;
end.
 
unit SquFoF_utils;
 
{$mode objfpc}{$H+}
 
interface
 
// Returns floating-point square root of 64-bit unsigned integer.
function FSqrt( x : UInt64) : double;
 
// Returns whether a 64-bit unsigned integer is a perfect square.
// In either case, returns floor(sqrt(x)) in the out parameter.
function IsSquare( x : UInt64; out iroot : UInt64) : boolean;
 
// Returns g.c.d. of 64-bit and 16-bit unsigned integer.
function GCD( u : UInt64; x : word) : word;
 
implementation
 
function FSqrt( x : UInt64) : double;
// Both Free Pascal and Delphi 7 seem unreliable when casting
// a 64-bit integer to floating point. We use a workaround.
type TSplitUint64 = packed record case boolean of
true: (All : UInt64);
false: (Lo, Hi : longword); // longword is 32-bit unsigned
end;
var
temp : TSplitUInt64;
begin
temp.All := x;
result := Sqrt( 1.0*temp.Lo + 4294967296.0*temp.Hi);
end;
 
// Based on Rosetta Code ISqrt, solution for Modula-2.
// Trunc of the f.p. square root won't do, bacause of rounding errors..
function IsSquare( x : UInt64; out iroot : UInt64) : boolean;
var
Xdiv4, q, r, s, z : UInt64;
begin
Xdiv4 := X shr 2;
q := 1;
while q <= Xdiv4 do q := q shl 2;
z := x;
r := 0;
repeat
s := q + r;
r := r shr 1;
if z >= s then begin
z := z - s;
r := r + q;
end;
q := q shr 2;
until q = 0;
iroot := r;
result := (z = 0);
end;
 
function GCD( u : UInt64; x : word) : word;
var
y, t : word;
begin
y := u mod x;
while y <> 0 do begin
t := x mod y;
x := y;
y := t;
end;
result := x;
end;
end.
</syntaxhighlight>
{{out}}
<pre>
Integer FactorNumber Mult QuotientM/P Factorization
2501 3 Main 61 * 41
12851 1 Main 71 * 181
13289 1 Main 97 * 137
75301 3 Main 293 * 257
120787 1 Main 43 * 2809
967009 7 Main 1609 * 601
997417 1 Main 257 * 3881
7091569 1 Prelim 2663 * 2663
13290059 1 Main 3119 * 4261
42854447 3 Main 4423 * 9689
223553581 1 Main 11213 * 19937
2027651281 1 Main 44021 * 46061
11111111111 3 Main 21649 * 513239
100895598169 11 Main 898423 * 112303
1002742628021 No factor found
60012462237239 1 Main 6862753 * 8744663
287129523414791 5 Main 6059887 * 47381993
9007199254740931 1 Main 10624181 * 847801751
11111111111111111 1 Main 2071723 * 5363222357
314159265358979323 1 Main 317213509 * 990371647
384307168202281507 1 Main 415718707 * 924440401
419244183493398773 1 Main 48009977 * 8732438749
658812288346769681 1 Main 62222119 * 10588072199
922337203685477563 1 Main 110075821 * 8379108103
1000000000000000127 1 Main 111756107 * 8948056861
1152921505680588799 3 Main 139001459 * 8294312261
1537228672809128917 3 Main 26675843 * 57626245319
4611686018427387877 1 Main 343242169 * 13435662733
</pre>
 
=={{header|Perl}}==
{{trans|Raku}}
{{libheader|ntheory}}
<syntaxhighlight 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) }
}</syntaxhighlight>
{{out}}
<pre>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</pre>
 
=={{header|Phix}}==
{{trans|C|<small>(Classical heuristic - fixes the two incorrectly failing cases of the previous version)</small>}}
<!--<syntaxhighlight lang="phix">(phixonline)-->
<span style="color: #000080;font-style:italic;">--requires(64) -- (decided to limit 32-bit explicitly instead)</span>
<span style="color: #008080;">constant</span> <span style="color: #000000;">MxN</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">power<span style="color: #0000FF;">(<span style="color: #000000;">2<span style="color: #0000FF;">,<span style="color: #008080;">iff<span style="color: #0000FF;">(<span style="color: #7060A8;">machine_bits<span style="color: #0000FF;">(<span style="color: #0000FF;">)<span style="color: #0000FF;">=<span style="color: #000000;">32<span style="color: #0000FF;">?<span style="color: #000000;">53<span style="color: #0000FF;">:<span style="color: #000000;">63<span style="color: #0000FF;">)<span style="color: #0000FF;">)<span style="color: #0000FF;">,</span>
<span style="color: #000000;">m</span> <span style="color: #0000FF;">=</span> <span style="color: #0000FF;">{<span style="color: #000000;">1<span style="color: #0000FF;">,</span> <span style="color: #000000;">3<span style="color: #0000FF;">,</span> <span style="color: #000000;">5<span style="color: #0000FF;">,</span> <span style="color: #000000;">7<span style="color: #0000FF;">,</span> <span style="color: #000000;">11<span style="color: #0000FF;">}</span>
<span style="color: #008080;">function</span> <span style="color: #000000;">squfof<span style="color: #0000FF;">(<span style="color: #004080;">atom</span> <span style="color: #000000;">N<span style="color: #0000FF;">)</span>
<span style="color: #000080;font-style:italic;">-- square form factorization</span>
<span style="color: #004080;">integer</span> <span style="color: #000000;">h<span style="color: #0000FF;">,</span> <span style="color: #000000;">a<span style="color: #0000FF;">=<span style="color: #000000;">0<span style="color: #0000FF;">,</span> <span style="color: #000000;">b<span style="color: #0000FF;">,</span> <span style="color: #000000;">c<span style="color: #0000FF;">,</span> <span style="color: #000000;">u<span style="color: #0000FF;">=<span style="color: #000000;">0<span style="color: #0000FF;">,</span> <span style="color: #000000;">v<span style="color: #0000FF;">,</span> <span style="color: #000000;">w<span style="color: #0000FF;">,</span> <span style="color: #000000;">rN<span style="color: #0000FF;">,</span> <span style="color: #000000;">q<span style="color: #0000FF;">,</span> <span style="color: #000000;">r<span style="color: #0000FF;">,</span> <span style="color: #000000;">t</span>
<span style="color: #008080;">if</span> <span style="color: #7060A8;">remainder<span style="color: #0000FF;">(<span style="color: #000000;">N<span style="color: #0000FF;">,<span style="color: #000000;">2<span style="color: #0000FF;">)<span style="color: #0000FF;">==<span style="color: #000000;">0</span> <span style="color: #008080;">then</span> <span style="color: #008080;">return</span> <span style="color: #000000;">2</span> <span style="color: #008080;">end</span> <span style="color: #008080;">if</span>
<span style="color: #000000;">h</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">floor<span style="color: #0000FF;">(<span style="color: #7060A8;">sqrt<span style="color: #0000FF;">(<span style="color: #000000;">N<span style="color: #0000FF;">)</span> <span style="color: #0000FF;">+</span> <span style="color: #000000;">0.5<span style="color: #0000FF;">)</span>
<span style="color: #008080;">if</span> <span style="color: #000000;">h<span style="color: #0000FF;">*<span style="color: #000000;">h<span style="color: #0000FF;">==<span style="color: #000000;">N</span> <span style="color: #008080;">then</span> <span style="color: #008080;">return</span> <span style="color: #000000;">h</span> <span style="color: #008080;">end</span> <span style="color: #008080;">if</span>
<span style="color: #008080;">for</span> <span style="color: #000000;">k<span style="color: #0000FF;">=<span style="color: #000000;">1</span> <span style="color: #008080;">to</span> <span style="color: #7060A8;">length<span style="color: #0000FF;">(<span style="color: #000000;">m<span style="color: #0000FF;">)</span> <span style="color: #008080;">do</span>
<span style="color: #004080;">integer</span> <span style="color: #000000;">mk</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">m<span style="color: #0000FF;">[<span style="color: #000000;">k<span style="color: #0000FF;">]</span>
<span style="color: #008080;">if</span> <span style="color: #000000;">mk<span style="color: #0000FF;">><span style="color: #000000;">1</span> <span style="color: #008080;">and</span> <span style="color: #7060A8;">remainder<span style="color: #0000FF;">(<span style="color: #000000;">N<span style="color: #0000FF;">,<span style="color: #000000;">mk<span style="color: #0000FF;">)<span style="color: #0000FF;">==<span style="color: #000000;">0</span> <span style="color: #008080;">then</span> <span style="color: #008080;">return</span> <span style="color: #000000;">mk</span> <span style="color: #008080;">end</span> <span style="color: #008080;">if</span>
<span style="color: #000080;font-style:italic;">//check overflow m * N</span>
<span style="color: #008080;">if</span> <span style="color: #000000;">N<span style="color: #0000FF;">><span style="color: #000000;">MxN<span style="color: #0000FF;">/<span style="color: #000000;">mk</span> <span style="color: #008080;">then</span> <span style="color: #008080;">exit</span> <span style="color: #008080;">end</span> <span style="color: #008080;">if</span>
<span style="color: #004080;">atom</span> <span style="color: #000000;">mN</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">N<span style="color: #0000FF;">*<span style="color: #000000;">mk</span>
<span style="color: #000000;">r</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">floor<span style="color: #0000FF;">(<span style="color: #7060A8;">sqrt<span style="color: #0000FF;">(<span style="color: #000000;">mN<span style="color: #0000FF;">)<span style="color: #0000FF;">)</span>
<span style="color: #008080;">if</span> <span style="color: #000000;">r<span style="color: #0000FF;">*<span style="color: #000000;">r<span style="color: #0000FF;">><span style="color: #000000;">mN</span> <span style="color: #008080;">then</span> <span style="color: #000000;">r</span> <span style="color: #0000FF;">-=</span> <span style="color: #000000;">1</span> <span style="color: #008080;">end</span> <span style="color: #008080;">if</span>
<span style="color: #000000;">rN</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">r</span>
<span style="color: #000080;font-style:italic;">//principal form</span>
<span style="color: #0000FF;">{<span style="color: #000000;">b<span style="color: #0000FF;">,<span style="color: #000000;">a<span style="color: #0000FF;">}</span> <span style="color: #0000FF;">=</span> <span style="color: #0000FF;">{<span style="color: #000000;">r<span style="color: #0000FF;">,<span style="color: #000000;">1<span style="color: #0000FF;">}</span>
<span style="color: #000000;">h</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">floor<span style="color: #0000FF;">(<span style="color: #0000FF;">(<span style="color: #000000;">rN<span style="color: #0000FF;">+<span style="color: #000000;">b<span style="color: #0000FF;">)<span style="color: #0000FF;">/<span style="color: #000000;">a<span style="color: #0000FF;">)<span style="color: #0000FF;">*<span style="color: #000000;">a<span style="color: #0000FF;">-<span style="color: #000000;">b</span>
<span style="color: #000000;">c</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">floor<span style="color: #0000FF;">(<span style="color: #0000FF;">(<span style="color: #000000;">mN<span style="color: #0000FF;">-<span style="color: #000000;">h<span style="color: #0000FF;">*<span style="color: #000000;">h<span style="color: #0000FF;">)<span style="color: #0000FF;">/<span style="color: #000000;">a<span style="color: #0000FF;">)</span>
<span style="color: #008080;">for</span> <span style="color: #000000;">i<span style="color: #0000FF;">=<span style="color: #000000;">2</span> <span style="color: #008080;">to</span> <span style="color: #7060A8;">floor<span style="color: #0000FF;">(<span style="color: #7060A8;">sqrt<span style="color: #0000FF;">(<span style="color: #000000;">2<span style="color: #0000FF;">*<span style="color: #000000;">r<span style="color: #0000FF;">)<span style="color: #0000FF;">)</span> <span style="color: #0000FF;">*</span> <span style="color: #000000;">4<span style="color: #0000FF;">-<span style="color: #000000;">1</span> <span style="color: #008080;">do</span>
<span style="color: #000080;font-style:italic;">//search principal cycle</span>
<span style="color: #0000FF;">{<span style="color: #000000;">a<span style="color: #0000FF;">,<span style="color: #000000;">c<span style="color: #0000FF;">,<span style="color: #000000;">t<span style="color: #0000FF;">}</span> <span style="color: #0000FF;">=</span> <span style="color: #0000FF;">{<span style="color: #000000;">c<span style="color: #0000FF;">,<span style="color: #000000;">a<span style="color: #0000FF;">,<span style="color: #000000;">b<span style="color: #0000FF;">}</span>
<span style="color: #000000;">q</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">floor<span style="color: #0000FF;">(<span style="color: #0000FF;">(<span style="color: #000000;">rN<span style="color: #0000FF;">+<span style="color: #000000;">b<span style="color: #0000FF;">)<span style="color: #0000FF;">/<span style="color: #000000;">a<span style="color: #0000FF;">)</span>
<span style="color: #000000;">b</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">q<span style="color: #0000FF;">*<span style="color: #000000;">a<span style="color: #0000FF;">-<span style="color: #000000;">b</span>
<span style="color: #000000;">c</span> <span style="color: #0000FF;">+=</span> <span style="color: #000000;">q<span style="color: #0000FF;">*<span style="color: #0000FF;">(<span style="color: #000000;">t<span style="color: #0000FF;">-<span style="color: #000000;">b<span style="color: #0000FF;">)</span>
<span style="color: #008080;">if</span> <span style="color: #7060A8;">remainder<span style="color: #0000FF;">(<span style="color: #000000;">i<span style="color: #0000FF;">,<span style="color: #000000;">2<span style="color: #0000FF;">)<span style="color: #0000FF;">==<span style="color: #000000;">0</span> <span style="color: #008080;">then</span>
<span style="color: #000000;">r</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">floor<span style="color: #0000FF;">(<span style="color: #7060A8;">sqrt<span style="color: #0000FF;">(<span style="color: #000000;">c<span style="color: #0000FF;">)<span style="color: #0000FF;">+<span style="color: #000000;">0.5<span style="color: #0000FF;">)</span>
<span style="color: #008080;">if</span> <span style="color: #000000;">r<span style="color: #0000FF;">*<span style="color: #000000;">r<span style="color: #0000FF;">==<span style="color: #000000;">c</span> <span style="color: #008080;">then</span>
<span style="color: #000080;font-style:italic;">//square form found
//inverse square root</span>
<span style="color: #000000;">q</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">floor<span style="color: #0000FF;">(<span style="color: #0000FF;">(<span style="color: #000000;">rN<span style="color: #0000FF;">-<span style="color: #000000;">b<span style="color: #0000FF;">)<span style="color: #0000FF;">/<span style="color: #000000;">r<span style="color: #0000FF;">)</span>
<span style="color: #000000;">v</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">q<span style="color: #0000FF;">*<span style="color: #000000;">r<span style="color: #0000FF;">+<span style="color: #000000;">b</span>
<span style="color: #000000;">w</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">floor<span style="color: #0000FF;">(<span style="color: #0000FF;">(<span style="color: #000000;">mN<span style="color: #0000FF;">-<span style="color: #000000;">v<span style="color: #0000FF;">*<span style="color: #000000;">v<span style="color: #0000FF;">)<span style="color: #0000FF;">/<span style="color: #000000;">r<span style="color: #0000FF;">)</span>
<span style="color: #000080;font-style:italic;">//search ambiguous cycle</span>
<span style="color: #000000;">u</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">r</span>
<span style="color: #008080;">while</span> <span style="color: #004600;">true</span> <span style="color: #008080;">do</span>
<span style="color: #0000FF;">{<span style="color: #000000;">u<span style="color: #0000FF;">,<span style="color: #000000;">w<span style="color: #0000FF;">,<span style="color: #000000;">r<span style="color: #0000FF;">}</span> <span style="color: #0000FF;">=</span> <span style="color: #0000FF;">{<span style="color: #000000;">w<span style="color: #0000FF;">,<span style="color: #000000;">u<span style="color: #0000FF;">,<span style="color: #000000;">v<span style="color: #0000FF;">}</span>
<span style="color: #000000;">q</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">floor<span style="color: #0000FF;">(<span style="color: #0000FF;">(<span style="color: #000000;">rN<span style="color: #0000FF;">+<span style="color: #000000;">v<span style="color: #0000FF;">)<span style="color: #0000FF;">/<span style="color: #000000;">u<span style="color: #0000FF;">)</span>
<span style="color: #000000;">v</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">q<span style="color: #0000FF;">*<span style="color: #000000;">u<span style="color: #0000FF;">-<span style="color: #000000;">v</span>
<span style="color: #008080;">if</span> <span style="color: #000000;">v<span style="color: #0000FF;">==<span style="color: #000000;">r</span> <span style="color: #008080;">then</span> <span style="color: #008080;">exit</span> <span style="color: #008080;">end</span> <span style="color: #008080;">if</span>
<span style="color: #000000;">w</span> <span style="color: #0000FF;">+=</span> <span style="color: #000000;">q<span style="color: #0000FF;">*<span style="color: #0000FF;">(<span style="color: #000000;">r<span style="color: #0000FF;">-<span style="color: #000000;">v<span style="color: #0000FF;">)</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">while</span>
<span style="color: #000080;font-style:italic;">//symmetry point</span>
<span style="color: #000000;">h</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">gcd<span style="color: #0000FF;">(<span style="color: #000000;">N<span style="color: #0000FF;">,<span style="color: #000000;">u<span style="color: #0000FF;">)</span>
<span style="color: #008080;">if</span> <span style="color: #000000;">h<span style="color: #0000FF;">!=<span style="color: #000000;">1</span> <span style="color: #008080;">then</span> <span style="color: #008080;">return</span> <span style="color: #000000;">h</span> <span style="color: #008080;">end</span> <span style="color: #008080;">if</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">if</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">if</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span>
<span style="color: #008080;">return</span> <span style="color: #000000;">1</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">function</span>
<span style="color: #008080;">constant</span> <span style="color: #000000;">tests</span> <span style="color: #0000FF;">=</span> <span style="color: #0000FF;">{<span style="color: #000000;">2501<span style="color: #0000FF;">,</span> <span style="color: #000000;">12851<span style="color: #0000FF;">,</span> <span style="color: #000000;">13289<span style="color: #0000FF;">,</span> <span style="color: #000000;">75301<span style="color: #0000FF;">,</span> <span style="color: #000000;">120787<span style="color: #0000FF;">,</span> <span style="color: #000000;">967009<span style="color: #0000FF;">,</span> <span style="color: #000000;">997417<span style="color: #0000FF;">,</span> <span style="color: #000000;">7091569<span style="color: #0000FF;">,</span> <span style="color: #000000;">5214317<span style="color: #0000FF;">,</span> <span style="color: #000000;">20834839<span style="color: #0000FF;">,</span>
<span style="color: #000000;">23515517<span style="color: #0000FF;">,</span> <span style="color: #000000;">33409583<span style="color: #0000FF;">,</span> <span style="color: #000000;">44524219<span style="color: #0000FF;">,</span> <span style="color: #000000;">13290059<span style="color: #0000FF;">,</span> <span style="color: #000000;">223553581<span style="color: #0000FF;">,</span> <span style="color: #000000;">42854447<span style="color: #0000FF;">,</span> <span style="color: #000000;">223553581<span style="color: #0000FF;">,</span> <span style="color: #000000;">2027651281<span style="color: #0000FF;">,</span>
<span style="color: #000000;">11111111111<span style="color: #0000FF;">,</span> <span style="color: #000000;">100895598169<span style="color: #0000FF;">,</span> <span style="color: #000000;">1002742628021<span style="color: #0000FF;">,</span> <span style="color: #000080;font-style:italic;">-- (prime/expected to fail)</span>
<span style="color: #000000;">60012462237239<span style="color: #0000FF;">,</span> <span style="color: #000000;">287129523414791<span style="color: #0000FF;">,</span> <span style="color: #000000;">9007199254740931<span style="color: #0000FF;">,</span> <span style="color: #000000;">32<span style="color: #0000FF;">,</span> <span style="color: #000080;font-style:italic;">-- (limit of 32-bit)</span>
<span style="color: #000000;">11111111111111111<span style="color: #0000FF;">,</span> <span style="color: #000000;">314159265358979323<span style="color: #0000FF;">,</span> <span style="color: #000000;">384307168202281507<span style="color: #0000FF;">,</span> <span style="color: #000000;">419244183493398773<span style="color: #0000FF;">,</span>
<span style="color: #000000;">658812288346769681<span style="color: #0000FF;">,</span> <span style="color: #000000;">922337203685477563<span style="color: #0000FF;">,</span> <span style="color: #000000;">1000000000000000127<span style="color: #0000FF;">,</span> <span style="color: #000000;">1152921505680588799<span style="color: #0000FF;">,</span>
<span style="color: #000000;">1537228672809128917<span style="color: #0000FF;">,</span> <span style="color: #000000;">4611686018427387877<span style="color: #0000FF;">}</span>
<span style="color: #7060A8;">printf<span style="color: #0000FF;">(<span style="color: #000000;">1<span style="color: #0000FF;">,<span style="color: #008000;">"N f N/f\n"<span style="color: #0000FF;">)</span>
<span style="color: #7060A8;">printf<span style="color: #0000FF;">(<span style="color: #000000;">1<span style="color: #0000FF;">,<span style="color: #008000;">"======================================\n"<span style="color: #0000FF;">)</span>
<span style="color: #008080;">for</span> <span style="color: #000000;">i<span style="color: #0000FF;">=<span style="color: #000000;">1</span> <span style="color: #008080;">to</span> <span style="color: #7060A8;">length<span style="color: #0000FF;">(<span style="color: #000000;">tests<span style="color: #0000FF;">)</span> <span style="color: #008080;">do</span>
<span style="color: #004080;">atom</span> <span style="color: #000000;">N</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">tests<span style="color: #0000FF;">[<span style="color: #000000;">i<span style="color: #0000FF;">]</span>
<span style="color: #008080;">if</span> <span style="color: #000000;">N<span style="color: #0000FF;">=<span style="color: #000000;">32</span> <span style="color: #008080;">then</span>
<span style="color: #008080;">if</span> <span style="color: #7060A8;">machine_bits<span style="color: #0000FF;">(<span style="color: #0000FF;">)<span style="color: #0000FF;">=<span style="color: #000000;">32</span> <span style="color: #008080;">then</span> <span style="color: #008080;">exit</span> <span style="color: #008080;">end</span> <span style="color: #008080;">if</span>
<span style="color: #008080;">else</span>
<span style="color: #004080;">atom</span> <span style="color: #000000;">f</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">squfof<span style="color: #0000FF;">(<span style="color: #000000;">N<span style="color: #0000FF;">)</span>
<span style="color: #7060A8;">printf<span style="color: #0000FF;">(<span style="color: #000000;">1<span style="color: #0000FF;">,<span style="color: #008000;">"%-22d %s\n"<span style="color: #0000FF;">,<span style="color: #0000FF;">{<span style="color: #000000;">N<span style="color: #0000FF;">,<span style="color: #008080;">iff<span style="color: #0000FF;">(<span style="color: #000000;">f<span style="color: #0000FF;">=<span style="color: #000000;">1<span style="color: #0000FF;">?<span style="color: #008000;">"fail"<span style="color: #0000FF;">:<span style="color: #7060A8;">sprintf<span style="color: #0000FF;">(<span style="color: #008000;">"%-10d %d"<span style="color: #0000FF;">,<span style="color: #0000FF;">{<span style="color: #000000;">f<span style="color: #0000FF;">,<span style="color: #000000;">N<span style="color: #0000FF;">/<span style="color: #000000;">f<span style="color: #0000FF;">}<span style="color: #0000FF;">)<span style="color: #0000FF;">)<span style="color: #0000FF;">}<span style="color: #0000FF;">)</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">if</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">for
<!--</syntaxhighlight>-->
{{out}}
<small>(on 64-bit, whereas the last 10 entries, ie 11111111111111111 on, are simply omitted on 32-bit)</small>
<pre>
N f N/f
======================================
2501 61 41
12851 71 181
13289 13797 97 137
75301 293 257
120787 43 2809
Line 975 ⟶ 2,503:
997417 257 3881
7091569 2663 2663
5214317 73 71429
20834839 3361 6199
23515517 53 443689
33409583 991 33713
44524219 593 75083
13290059 3119 4261
42854447 9689 4423
223553581 11213 19937
202765128142854447 46061 4423 44021 9689
223553581 11213 19937
2027651281 44021 46061
11111111111 21649 513239
100895598169 112303 898423
1002742628021 0 fail
60012462237239 6862753 8744663
287129523414791 6059887 47381993
0 0 fail
9007199254740931 10624181 847801751
11111111111111111 2071723 5363222357
Line 992 ⟶ 2,525:
658812288346769681 62222119 10588072199
922337203685477563 110075821 8379108103
1000000000000000127 0111756107 fail8948056861
1152921505680588799 139001459 8294312261
1537228672809128917 026675843 fail57626245319
4611686018427387877 343242169 13435662733
</pre>
 
=={{header|Raku}}==
===A just for fun snail ..===
References: [https://en.wikipedia.org/wiki/Shanks%27s_square_forms_factorization#Algorithm Algorithm], [https://en.wikipedia.org/wiki/Shanks%27s_square_forms_factorization#Example_implementations C program example] from the WP page and also the pseudocode from [http://colin.barker.pagesperso-orange.fr/lpa/big_squf.htm here].
<syntaxhighlight lang="raku" line># 20210325 Raku programming solution
 
my @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 );
 
sub circumfix:<⌊ ⌋>{ $^n.floor }; sub prefix:<√>{ $^n.sqrt }; # just for fun
 
sub SQUFOF ( \𝑁 ) {
 
return 1 if 𝑁.is-prime; # if n is prime return 1
return √𝑁 if √𝑁 == Int(√𝑁); # if n is a perfect square return √𝑁
 
for @multiplier -> \𝑘 {
my \Pₒ = $ = ⌊ √(𝑘*𝑁) ⌋; # P[0]=floor(√N)
my \Qₒ = $ = 1 ; # Q[0]=1
my \Q = $ = 𝑘*𝑁 - Pₒ²; # Q[1]=N-P[0]^2 & Q[i]
my \Pₚᵣₑᵥ = $ = Pₒ; # P[i-1] = P[0]
my \Qₚᵣₑᵥ = $ = Qₒ; # Q[i-1] = Q[0]
my \P = $ = 0; # P[i]
my \Qₙₑₓₜ = $ = 0; # P[i+1]
my \b = $ = 0; # b[i]
# i = 1
repeat until √Q == Int(√Q) { # until Q[i] is a perfect square
b = ⌊⌊ √(𝑘*𝑁) + Pₚᵣₑᵥ ⌋ / Q ⌋; # floor(floor(√N+P[i-1])/Q[i])
P = b*Q - Pₚᵣₑᵥ; # P[i]=b*Q[i]-P[i-1]
Qₙₑₓₜ = Qₚᵣₑᵥ + b*(Pₚᵣₑᵥ - P); # Q[i+1]=Q[i-1]+b(P[i-1]-P[i])
( Qₚᵣₑᵥ, Q, Pₚᵣₑᵥ ) = Q, Qₙₑₓₜ, P; # i++
}
 
b = ⌊ ⌊ √(𝑘*𝑁)+P ⌋ / Q ⌋; # b=floor((floor(√N)+P[i])/Q[0])
Pₚᵣₑᵥ = b*Qₒ - P; # P[i-1]=b*Q[0]-P[i]
Q = ( 𝑘*𝑁 - Pₚᵣₑᵥ² )/Qₒ; # Q[1]=(N-P[0]^2)/Q[0] & Q[i]
Qₚᵣₑᵥ = Qₒ; # Q[i-1] = Q[0]
# i = 1
loop { # repeat
b = ⌊ ⌊ √(𝑘*𝑁)+Pₚᵣₑᵥ ⌋ / Q ⌋; # b=floor(floor(√N)+P[i-1])/Q[i])
P = b*Q - Pₚᵣₑᵥ; # P[i]=b*Q[i]-P[i-1]
Qₙₑₓₜ = Qₚᵣₑᵥ + b*(Pₚᵣₑᵥ - P); # Q[i+1]=Q[i-1]+b(P[i-1]-P[i])
last if (P == Pₚᵣₑᵥ); # until P[i+1]=P[i]
( Qₚᵣₑᵥ, Q, Pₚᵣₑᵥ ) = Q, Qₙₑₓₜ, P; # i++
}
given 𝑁 gcd P { return $_ if $_ != 1|𝑁 }
} # gcd(N,P[i]) (if != 1 or N) is a factor of N, otherwise try next k
return 0 # give up
}
 
race for (
11111, # wikipedia.org/wiki/Shanks%27s_square_forms_factorization#Example
4558849, # example from talk page
# all of the rest are taken from the FreeBASIC entry
2501,12851,13289,75301,120787,967009,997417,7091569,13290059,
42854447,223553581,2027651281,11111111111,100895598169,1002742628021,
# time hoarders
60012462237239, # = 6862753 * 8744663 15s
287129523414791, # = 6059887 * 47381993 80s
11111111111111111, # = 2071723 * 5363222357 2m
384307168202281507, # = 415718707 * 924440401 5m
1000000000000000127, # = 111756107 * 8948056861 12m
9007199254740931, # = 10624181 * 847801751 17m
922337203685477563, # = 110075821 * 8379108103 41m
314159265358979323, # = 317213509 * 990371647 61m
1152921505680588799, # = 139001459 * 8294312261 93m
658812288346769681, # = 62222119 * 10588072199 112m
419244183493398773, # = 48009977 * 8732438749 135m
1537228672809128917, # = 26675843 * 57626245319 254m
# don't know how to handle this one
# for 1e-323, 1e-324 { my $*TOLERANCE = $_ ;
# say 4611686018427387877.sqrt ≅ 4611686018427387877.sqrt.Int }
# skip the perfect square check and start k with 3 to get the following
# 4611686018427387877, # = 343242169 * 13435662733 217m
) -> \data {
given data.&SQUFOF {
when 0 { say "The number ", data, " is not factored." }
when 1 { say "The number ", data, " is a prime." }
default { say data, " = ", $_, " * ", data div $_.Int }
}
}
</syntaxhighlight>
{{out}}
<pre>
11111 = 41 * 271
4558849 = 383 * 11903
2501 = 61 * 41
12851 = 71 * 181
13289 = 97 * 137
75301 = 293 * 257
120787 = 43 * 2809
967009 = 601 * 1609
997417 = 257 * 3881
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
</pre>
 
===Use NativeCall===
Use idea similar to the second approach from [https://rosettacode.org/wiki/Machine_code#Raku here], by compiling the C (Classical heuristic) entry to a shared library and then make use of the squfof routine.
 
squfof.raku
<syntaxhighlight lang="raku" line># 20210326 Raku programming solution
 
use NativeCall;
 
constant LIBSQUFOF = '/home/user/LibSQUFOF.so';
 
sub squfof(uint64 $n) returns uint64 is native(LIBSQUFOF) { * };
 
race for (
11111, # wikipedia.org/wiki/Shanks%27s_square_forms_factorization#Example
4558849, # example from talk page
# all of the rest are taken from the FreeBASIC entry
2501,12851,13289,75301,120787,967009,997417,7091569,13290059,
42854447,223553581,2027651281,11111111111,100895598169,1002742628021,
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
4611686018427387877, # = 343242169 * 13435662733
) -> \data {
given squfof(data) { say data, " = ", $_, " * ", data div $_ }
}</syntaxhighlight>
{{out}}<pre>gcc -Wall -fPIC -shared -o LibSQUFOF.so squfof.c
file LibSQUFOF.so
LibSQUFOF.so: ELF 64-bit LSB shared object, x86-64, version 1 (SYSV), dynamically linked, BuildID[sha1]=f7f9864e5508ba1de80cbc6e2f4d789f4ec448e9, not stripped
raku -c squfof.raku && time ./squfof.raku
Syntax OK
11111 = 41 * 271
4558849 = 383 * 11903
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 = 1 * 1002742628021
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
4611686018427387877 = 343242169 * 13435662733
 
real 0m0.784s
user 0m0.716s
sys 0m0.056s
</pre>
 
Line 1,005 ⟶ 2,726:
within the
<br>REXX program) &nbsp; is because the number being tested is a prime &nbsp; (1,002,742,628,021).
<langsyntaxhighlight 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
Line 1,038 ⟶ 2,759:
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); p= po
popprev= iSqrt(d)po; QQ= d - pprev= po*po
pqprev= po1; qprevBB= 1iSqrt(s+s)*6
QQ do i=2 d while -i<BB; po * b= (po+p)%QQ
p= b*QQ - p; q= QQ
BB= iSqrt(s+s) * 6
do iQQ=2 qprev while+ i<BBb*(pprev-p); br= iSqrt(po + pQQ) % 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
ppprev= bp; * r + p qprev= r
pprevQQ= p; (d - pprev*pprev)%qprev= r
QQ= (d - do until p==pprev*pprev); % qprev pprev= p
do untilb= (po+p==pprev)%QQ; q= pprevQQ; p= b*QQ - p
bQQ= (poqprev + b*(pprev-p); % QQ qprev= q
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</langsyntaxhighlight>
{{out|output|text=&nbsp; when using the internal default input:}}
<pre>
Line 1,095 ⟶ 2,810:
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>
 
=={{header|Sidef}}==
{{trans|Perl}}
<syntaxhighlight lang="ruby">const multipliers = divisors(3*5*7*11).grep { _ > 1 }
 
func sff(N) {
 
N.is_prime && return 1 # n is prime
N.is_square && return N.isqrt # n is square
 
multipliers.each {|k|
 
var P0 = isqrt(k*N) # P[0]=floor(sqrt(N)
var Q0 = 1 # Q[0]=1
var Q = (k*N - P0*P0) # Q[1]=N-P[0]^2 & Q[i]
var P1 = P0 # P[i-1] = P[0]
var Q1 = Q0 # Q[i-1] = Q[0]
var P = 0 # P[i]
var Qn = 0 # P[i+1]
var b = 0 # b[i]
 
while (!Q.is_square) { # until Q[i] is a perfect square
b = idiv(isqrt(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 = idiv(isqrt(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*P1)/Q0 # Q[1]=(N-P[0]^2)/Q[0] & Q[i]
Q1 = Q0 # Q[i-1] = Q[0]
 
loop {
b = idiv(isqrt(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])
break if (P == P1) # until P[i+1]=P[i]
(Q1, Q, P1) = (Q, Qn, P)
}
with (gcd(N,P)) {|g|
return g if g.is_ntf(N)
}
}
 
return 0
}
 
[ 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
].each {|n|
var v = sff(n)
if (v == 0) { say "The number #{n} is not factored." }
elsif (v == 1) { say "The number #{n} is a prime." }
else { say "#{n} = #{[n/v, v].sort.join(' * ')}" }
}</syntaxhighlight>
 
{{out}}
<pre>
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
^C
</pre>
 
Line 1,108 ⟶ 2,906:
 
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.
<langsyntaxhighlight ecmascriptlang="wren">import "./long" for ULong
import "./big" for BigInt
import "./fmt" for Fmt
 
var multipliers = [
Line 1,210 ⟶ 3,008:
var quot = (fact.isZero) ? "fail" : (N / fact).toString
Fmt.print("$-20s $-10s $s", N, fact, quot)
}</langsyntaxhighlight>
 
{{out}}
PureFox
9,476

edits

Retrieved from "https://rosettacode.org/wiki/Special:MobileDiff/309033...354627"