Square form factorization: Difference between revisions
A classic O(N^1/4) factorization algorithm |
m →{{header|Wren}}: Minor tidy |
||
(60 intermediate revisions by 16 users not shown) | |||
Line 1: | Line 1: | ||
{{ |
{{task|Mathematics}} |
||
[[Category:Prime Numbers]] |
|||
;Task. |
;Task. |
||
[[wp:Daniel_Shanks|Daniel Shanks]]'s Square Form Factorization [[wp:Shanks%27s_square_forms_factorization|(SquFoF)]]. |
|||
Daniel Shanks's Square Form Factorization [[wp:Shanks%27s_square_forms_factorization|(SquFoF)]]. |
|||
Invented around 1975, ''‘On a 32-bit computer, SquFoF is the clear champion factoring algorithm'' |
Invented around 1975, ''‘On a 32-bit computer, SquFoF is the clear champion factoring algorithm'' |
||
''for numbers between 10<sup>10</sup> and 10<sup>18</sup>, and will likely remain so.’'' |
''for numbers between 10<sup>10</sup> and 10<sup>18</sup>, and will likely remain so.’'' |
||
An integral binary quadratic form [[wp:Binary_quadratic_form|(bqf)]] is a polynomial |
|||
An integral [[wp:Binary_quadratic_form|binary quadratic form]] is a polynomial |
|||
{{math|''f''(''x,y'') = ''ax<sup>2</sup>'' + ''bxy'' + ''cy<sup>2</sup>''}} |
{{math|''f''(''x,y'') = ''ax<sup>2</sup>'' + ''bxy'' + ''cy<sup>2</sup>''}} |
||
with integer coefficients and discriminant {{math|''D'' = ''b<sup>2</sup>'' – ''4ac''}}. |
with integer coefficients and discriminant {{math|''D'' = ''b<sup>2</sup>'' – ''4ac''}}. |
||
For each positive discriminant there are |
For each positive discriminant there are multiple forms {{math|(''a, b, c'')}}. |
||
The next form in a periodic sequence (cycle) of adjacent forms is found by applying a reduction operator |
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. |
|||
{{math|floor(''√N'')}}, rho constructs a ''principal form'' |
Using {{math|floor(''√N'')}}, rho constructs a ''principal form'' |
||
{{math|(''1, b, c'')}} with {{math|''D'' = ''4N''}}. |
{{math|(''1, b, c'')}} with {{math|''D'' = ''4N''}}. |
||
SquFoF |
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 {{math|('' |
They come in pairs of associated forms {{math|(''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 |
|||
{{math| |
{{math|(''ak'')''<sup>2</sup>'' – ''4ac'' = ''a''(''a·k<sup>2</sup>'' – ''4c'') = ''4N''}} |
||
and (if a is not equal to 1 or 2) N is split. |
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 |
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 |
is squared, that square form will always land in the principal cycle. Conversely, the square root of any |
||
form |
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 |
A square form is easy to find: the last coefficient ''c'' is a perfect square. This happens about once |
||
Line 33: | Line 34: | ||
Then ''a'' or ''a/2'' divides D and therefore N. |
Then ''a'' or ''a/2'' divides D and therefore N. |
||
To avoid trivial factorizations, Shanks created a list (queue) to |
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, |
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. |
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 proper square |
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. |
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. |
|||
and the program will duly report failure. |
|||
;Reference. |
;Reference. |
||
[https://homes.cerias.purdue.edu/~ssw/squfof.pdf] A detailed analysis of SquFoF (2007) |
|||
[https://www.ams.org/journals/mcom/2008-77-261/S0025-5718-07-02010-8/S0025-5718-07-02010-8.pdf] |
|||
A detailed analysis of SquFoF (2007) |
|||
__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.]] |
|||
<syntaxhighlight 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); |
|||
} |
|||
} |
|||
} |
|||
</syntaxhighlight>{{out}} |
|||
<pre> |
|||
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 |
|||
</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}}== |
=={{header|FreeBASIC}}== |
||
<syntaxhighlight lang="freebasic">' *********************************************** |
|||
<lang freebasic> |
|||
' *********************************************** |
|||
'subject: Shanks's square form factorization: |
'subject: Shanks's square form factorization: |
||
' ambiguous forms of discriminant 4N |
' ambiguous forms of discriminant 4N |
||
' give factors of N. |
' give factors of N. |
||
'tested : FreeBasic 1. |
'tested : FreeBasic 1.08.1 |
||
const qx = 50 |
|||
'queue size |
|||
'------------------------------------------------ |
'------------------------------------------------ |
||
const MxN = culngint(1) shl 62 |
const MxN = culngint(1) shl 62 |
||
' |
'input maximum |
||
const qx = (1 shl 5) - 1 |
|||
'queue size |
|||
type arg |
type arg |
||
'squfof arguments |
|||
as ulong m, f |
as ulong m, f |
||
as integer vb |
as integer vb |
||
Line 70: | Line 743: | ||
type bqf |
type bqf |
||
declare sub rho ( |
declare sub rho () |
||
'reduce indefinite form |
'reduce indefinite form |
||
declare function |
declare function issq (byref r as ulong) as integer |
||
'return -1 if c is square, set r:= sqrt(c) |
'return -1 if c is square, set r:= sqrt(c) |
||
declare sub qform (byref g as string, byval t as integer) |
declare sub qform (byref g as string, byval t as integer) |
||
'print binary quadratic form #t (a, 2b, c) |
'print binary quadratic form #t (a, 2b, c) |
||
as ulongint mN |
|||
as ulong rN, a, b, c |
as ulong rN, a, b, c |
||
as integer vb |
as integer vb |
||
Line 88: | Line 760: | ||
'return -1 if a proper square form is found |
'return -1 if a proper square form is found |
||
as |
as ulong a(qx), L, m |
||
as |
as integer k, t |
||
end type |
end type |
||
'global variables |
|||
'global variable |
|||
dim shared N as ulongint |
dim shared N as ulongint |
||
'the number to split |
|||
dim shared flag as integer |
dim shared flag as integer |
||
'signal to end all threads |
|||
dim shared as ubyte q1024(1023), q3465(3464) |
|||
'quadratic residue tables |
|||
'------------------------------------------------ |
'------------------------------------------------ |
||
sub bqf.rho ( |
sub bqf.rho () |
||
dim as ulong q, t |
dim as ulong q, t |
||
swap a, c |
swap a, c |
||
'residue |
'residue |
||
q = (rN + b) \ a |
q = culng(rN + b) \ a |
||
b = q * a - b |
t = b: b = q * a - b |
||
'pseudo-square |
|||
if sw then |
|||
c += q * (t - b) |
|||
c += q * (t - b) |
|||
else |
|||
'initialize |
|||
c = (mN - culngint(b) * b) \ a |
|||
end if |
|||
end sub |
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 |
|||
F.c = (mN - 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 1023) andalso q3465(c mod 3465) then |
|||
'98.6% non-squares filtered |
|||
if q55(c mod 55) then |
|||
r = culng(sqr(c)) |
|||
'>98% non-squares filtered |
|||
if r * r = c then return -1 |
|||
if c = r * r then return -1 |
|||
end if |
|||
end if |
|||
end if |
end if |
||
issq = 0 |
|||
end function |
end function |
||
sub bqf.qform (byref g as string, byval t as integer) |
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 u = a, v = b, w = c |
|||
'{D/4 = mN} |
|||
d = v * v + u * w |
|||
if mN - d then |
|||
print "fail:"; d: exit sub |
|||
end if |
|||
if t and 1 then |
if t and 1 then |
||
w = -w |
w = -w |
||
Line 163: | Line 816: | ||
#macro red(r, a) |
#macro red(r, a) |
||
r = iif(a and 1, a, a shr 1) |
r = iif(a and 1, a, a shr 1) |
||
if m > |
if m > 2 then |
||
r = iif(r mod m, r, r \ m) |
r = iif(r mod m, r, r \ m) |
||
end if |
end if |
||
Line 169: | Line 822: | ||
sub queue.enq (byref P as bqf) |
sub queue.enq (byref P as bqf) |
||
if t = qx then exit sub |
|||
dim s as ulong |
dim s as ulong |
||
red(s, P.c) |
red(s, P.c) |
||
if s < L then |
if s < L then |
||
'circular queue |
|||
k = (k + 2) and qx |
|||
if k > t then t = k |
|||
'enqueue P.b, P.c |
'enqueue P.b, P.c |
||
a( |
a(k) = P.b mod s |
||
a( |
a(k + 1) = s |
||
end if |
end if |
||
end sub |
end sub |
||
Line 196: | Line 850: | ||
dim as ulong L2, m, r, t, f = 1 |
dim as ulong L2, m, r, t, f = 1 |
||
dim as integer ix, i, j |
dim as integer ix, i, j |
||
dim as ulongint h |
dim as ulongint mN, h |
||
'principal and ambiguous cycles |
'principal and ambiguous cycles |
||
dim as bqf P, A |
dim as bqf P, A |
||
dim Q as queue |
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)) |
h = culngint(sqr(N)) |
||
if |
if h * h = N then |
||
'N is square |
'N is square |
||
rp->f = culng(h) |
rp->f = culng(h) |
||
Line 208: | Line 867: | ||
end if |
end if |
||
h = N |
|||
rp->f = 1 |
rp->f = 1 |
||
'multiplier |
|||
m = rp->m |
m = rp->m |
||
if m > 1 then |
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 |
'check overflow m * N |
||
if |
if N > (MxN \ m) then exit sub |
||
h *= m |
|||
end if |
end if |
||
mN = N * m |
|||
r = int(sqr(mN)) |
|||
P.mN = h |
|||
A.mN = h |
|||
r = int(sqr(h)) |
|||
'float64 fix |
'float64 fix |
||
if culngint(r) * r > |
if culngint(r) * r > mN then r -= 1 |
||
P.rN = r |
P.rN = r |
||
A.rN = r |
A.rN = r |
||
Line 230: | Line 892: | ||
if P.vb then print "r = "; r |
if P.vb then print "r = "; r |
||
Q. |
Q.k = -2: Q.t = -1: Q.m = m |
||
'Queue entry bounds |
'Queue entry bounds |
||
Q.L = int(sqr(r * 2)) |
Q.L = int(sqr(r * 2)) |
||
Line 237: | Line 899: | ||
'principal form |
'principal form |
||
P.b = r: P.c = 1 |
P.b = r: P.c = 1 |
||
rhoin(P) |
|||
P.qform("P", 1) |
P.qform("P", 1) |
||
ix = Q.L shl |
ix = Q.L shl 2 |
||
'higher if necessary |
|||
for i = 2 to ix |
for i = 2 to ix |
||
'search principal cycle |
'search principal cycle |
||
Line 247: | Line 908: | ||
if P.c < L2 then Q.enq(P) |
if P.c < L2 then Q.enq(P) |
||
P.rho |
P.rho |
||
if (i and 1) = 0 then |
if (i and 1) = 0 andalso P.issq(r) then |
||
'square form found |
|||
if P.issqr(r) then |
|||
'square form found |
|||
if Q.pro(P, r) then |
|||
'cycle too short |
|||
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 |
A.rho: j += 1 |
||
A.qform("A", j) |
|||
if A.b = t then |
|||
' |
'symmetry point |
||
A.qform("A", j) |
|||
A. |
red(f, A.a) |
||
if f = 1 then exit do |
|||
flag = -1 |
|||
'factor found |
|||
end if |
|||
loop until flag |
|||
if f = 1 then exit do |
|||
end if ' proper square |
|||
end if ' square form |
|||
'factor found |
|||
end if |
|||
loop until flag |
|||
end if ' proper square |
|||
end if ' square form |
|||
end if ' even indices |
|||
if flag then exit for |
if flag then exit for |
||
Line 300: | Line 955: | ||
data 7091569 |
data 7091569 |
||
data 13290059 |
data 13290059 |
||
data 23515517 |
|||
data 42854447 |
data 42854447 |
||
data 223553581 |
data 223553581 |
||
Line 319: | Line 975: | ||
data 1537228672809128917 |
data 1537228672809128917 |
||
data 4611686018427387877 |
data 4611686018427387877 |
||
data 0 |
|||
'main |
'main |
||
Line 324: | Line 981: | ||
const tx = 4 |
const tx = 4 |
||
dim as double tim = timer |
dim as double tim = timer |
||
dim as ulongint i, f |
|||
dim h(4) as any ptr |
dim h(4) as any ptr |
||
dim a(4) as arg |
dim a(4) as arg |
||
dim |
dim as ulongint f |
||
dim as integer s, t |
|||
width 64, 30 |
width 64, 30 |
||
cls |
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 |
a(0).vb = 0 |
||
Line 351: | Line 1,015: | ||
print "N = "; N |
print "N = "; N |
||
flag = 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 |
|||
for t = 1 to tx + 1 step 2 |
|||
if t < tx then |
|||
for t = 1 to tx |
|||
h(t) = threadcreate(@squfof, @a(t)) |
h(t) = threadcreate(@squfof, @a(t)) |
||
end if |
|||
squfof(@a( |
squfof(@a(t - 1)) |
||
f = a(t - 1).f |
|||
if t < tx then |
|||
for t = 1 to tx |
|||
threadwait(h(t)) |
threadwait(h(t)) |
||
if f = 1 then f = a(t).f |
if f = 1 then f = a(t).f |
||
end if |
|||
if f > 1 then exit for |
|||
next t |
|||
if N mod f then f = 1 |
|||
'assert |
|||
elseif f = N then |
|||
if N mod f then f = 1 |
|||
f = 1 |
|||
end if |
|||
if f = 1 then |
if f = 1 then |
||
Line 389: | Line 1,044: | ||
print "total time:"; csng(timer - tim); " s" |
print "total time:"; csng(timer - tim); " s" |
||
end</syntaxhighlight> |
|||
end |
|||
</lang> |
|||
{{out| |
{{out|Examples}} |
||
<pre> |
<pre> |
||
N = 2501 |
N = 2501 |
||
f = |
f = 61 N/f = 41 |
||
N = 12851 |
N = 12851 |
||
f = |
f = 71 N/f = 181 |
||
N = 13289 |
N = 13289 |
||
Line 420: | Line 1,074: | ||
N = 13290059 |
N = 13290059 |
||
f = 3119 N/f = 4261 |
f = 3119 N/f = 4261 |
||
N = 23515517 |
|||
f = 53 N/f = 443689 |
|||
N = 42854447 |
N = 42854447 |
||
Line 425: | Line 1,082: | ||
N = 223553581 |
N = 223553581 |
||
f = |
f = 11213 N/f = 19937 |
||
N = 2027651281 |
N = 2027651281 |
||
Line 434: | Line 1,091: | ||
N = 100895598169 |
N = 100895598169 |
||
f = |
f = 112303 N/f = 898423 |
||
N = 1002742628021 |
N = 1002742628021 |
||
Line 478: | Line 1,135: | ||
f = 343242169 N/f = 13435662733 |
f = 343242169 N/f = 13435662733 |
||
total time: 0. |
total time: 0.0170462 s |
||
</pre> |
|||
=={{header|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. |
|||
<syntaxhighlight 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) |
|||
} |
|||
}</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 0 fail |
|||
1152921505680588799 139001459 8294312261 |
|||
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]] |
|||
<syntaxhighlight 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 |
|||
</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 |
|||
0.039027 seconds (698 allocations: 38.312 KiB) |
|||
</pre> |
|||
=={{header|Nim}}== |
|||
{{trans|Phix}} |
|||
<syntaxhighlight 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}"</syntaxhighlight> |
|||
{{out}} |
|||
<pre>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</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> |
|||
Number Mult M/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 97 137 |
|||
75301 293 257 |
|||
120787 43 2809 |
|||
967009 1609 601 |
|||
997417 257 3881 |
|||
7091569 2663 2663 |
|||
5214317 73 71429 |
|||
20834839 3361 6199 |
|||
23515517 53 443689 |
|||
33409583 991 33713 |
|||
44524219 593 75083 |
|||
13290059 3119 4261 |
|||
223553581 11213 19937 |
|||
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|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> |
|||
=={{header|REXX}}== |
|||
This REXX version is a modified version of the '''C''' code in the Wikipedia's article |
|||
at [https://en.wikipedia.org/wiki/Shanks%27s_square_forms_factorization Shanks' square forms factorization]. |
|||
The only "failure" of this REXX version of Shanks' square forms factorization (when using the numbers supplied |
|||
within the |
|||
<br>REXX program) is because the number being tested is a prime (1,002,742,628,021). |
|||
<syntaxhighlight lang="rexx">/*REXX pgm factors an integer using Daniel Shanks' (1917-1996) square form factorization*/ |
|||
numeric digits 100 /*ensure enough decimal digits.*/ |
|||
call dMults 1,3,5,7,11,3*5,3*7,3*11,5*7,5*11,7*11, 3*5*7, 3*5*11, 3*7*11, 5*7*11, 3*5*7*11 |
|||
call dTests 2501, 12851, 13289, 75301, 120787, 967009, 997417, 7091569, 13290059, , |
|||
42854447, 223553581, 2027651281, 11111111111, 100895598169, 1002742628021, , |
|||
60012462237239, 287129523414791, 9007199254740931, 11111111111111111, , |
|||
314159265358979323, 384307168202281507, 419244183493398773, , |
|||
658812288346769681, 922337203685477563, 1000000000000000127, , |
|||
1152921505680588799, 1537228672809128917, 4611686018427387877 |
|||
w= length( commas(!.$) ) /*the max width of test numbers*/ |
|||
do tests=1 for !.0; n= !.tests; nc= commas(n) |
|||
f= ssff(n); fc= commas(f); wf= length(fc); if f\==0 then nf= commas(n%f) |
|||
if f\==0 then do; nfc= commas(n%f); wnfc= length(nfc); end |
|||
if f ==0 then _= " (Shank's square form factor failed.)" |
|||
else _= ' factors are: ' right( fc, max(w%2 , wf ) ) " and " , |
|||
right(nfc, max(w%2+4, wnfc) ) |
|||
say right(nc, w+5) _ |
|||
end /*tests*/ |
|||
exit 0 /*stick a fork in it, we're all done. */ |
|||
/*──────────────────────────────────────────────────────────────────────────────────────*/ |
|||
commas: parse arg ?; do jc=length(?)-3 to 1 by -3; ?=insert(',', ?, jc); end; return ? |
|||
dMults: @.$= 0; do j=1 for arg(); @.j= arg(j); @.$=max(@.$, @.j); end; @.0=j-1; return |
|||
dTests: !.$= 0; do j=1 for arg(); !.j= arg(j); !.$=max(!.$, !.j); end; !.0=j-1; return |
|||
gcd: procedure; parse arg x,y; do until _==0; _= x // y; x= y; y= _; end; return x |
|||
/*──────────────────────────────────────────────────────────────────────────────────────*/ |
|||
iSqrt: procedure; parse arg x; r=0; q=1; do while q<=x; q=q*4; end |
|||
do while q>1; q=q%4; _=x-r-q; r=r%2; if _>=0 then do;x=_;r=r+q; end; end |
|||
return r |
|||
/*──────────────────────────────────────────────────────────────────────────────────────*/ |
|||
ssff: procedure expose @.; parse arg n; n= abs(n); er= '***error***' |
|||
s= iSqrt(n); if s**2==n then return s; big= 2**digits() |
|||
do #=1 for @.0; k= @.# /*get a # from the list of low factors*/ |
|||
if n>big/k then do; say er 'number is too large: ' commas(k); exit 8; end |
|||
d= n*k; po= iSqrt(d); p= po |
|||
pprev= po; QQ= d - po*po |
|||
qprev= 1; BB= iSqrt(s+s)*6 |
|||
do i=2 while i<BB; b= (po+p)%QQ |
|||
p= b*QQ - p; q= QQ |
|||
QQ= qprev + b*(pprev-p); r= iSqrt(QQ) |
|||
if i//2==0 then if r*r==QQ then leave |
|||
qprev= q; pprev= p |
|||
end /*i*/ |
|||
if i>=BB then iterate |
|||
b= (po-p)%r; p= b*r + p |
|||
pprev= p; qprev= r |
|||
QQ= (d - pprev*pprev)%qprev |
|||
do until p==pprev; pprev= p |
|||
b= (po+p)%QQ; q= QQ; p= b*QQ - p |
|||
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</syntaxhighlight> |
|||
{{out|output|text= when using the internal default input:}} |
|||
<pre> |
|||
2,501 factors are: 61 and 41 |
|||
12,851 factors are: 71 and 181 |
|||
13,289 factors are: 137 and 97 |
|||
75,301 factors are: 293 and 257 |
|||
120,787 factors are: 43 and 2,809 |
|||
967,009 factors are: 1,609 and 601 |
|||
997,417 factors are: 257 and 3,881 |
|||
7,091,569 factors are: 2,663 and 2,663 |
|||
13,290,059 factors are: 3,119 and 4,261 |
|||
42,854,447 factors are: 9,689 and 4,423 |
|||
223,553,581 factors are: 11,213 and 19,937 |
|||
2,027,651,281 factors are: 46,061 and 44,021 |
|||
11,111,111,111 factors are: 21,649 and 513,239 |
|||
100,895,598,169 factors are: 112,303 and 898,423 |
|||
1,002,742,628,021 (Shank's square form factor failed.) |
|||
60,012,462,237,239 factors are: 6,862,753 and 8,744,663 |
|||
287,129,523,414,791 factors are: 6,059,887 and 47,381,993 |
|||
9,007,199,254,740,931 factors are: 10,624,181 and 847,801,751 |
|||
11,111,111,111,111,111 factors are: 2,071,723 and 5,363,222,357 |
|||
314,159,265,358,979,323 factors are: 317,213,509 and 990,371,647 |
|||
384,307,168,202,281,507 factors are: 415,718,707 and 924,440,401 |
|||
419,244,183,493,398,773 factors are: 48,009,977 and 8,732,438,749 |
|||
658,812,288,346,769,681 factors are: 62,222,119 and 10,588,072,199 |
|||
922,337,203,685,477,563 factors are: 110,075,821 and 8,379,108,103 |
|||
1,000,000,000,000,000,127 factors are: 111,756,107 and 8,948,056,861 |
|||
1,152,921,505,680,588,799 factors are: 139,001,459 and 8,294,312,261 |
|||
1,537,228,672,809,128,917 factors are: 26,675,843 and 57,626,245,319 |
|||
4,611,686,018,427,387,877 factors are: 343,242,169 and 13,435,662,733 |
|||
</pre> |
|||
=={{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> |
|||
=={{header|Wren}}== |
|||
{{libheader|Wren-long}} |
|||
{{libheader|Wren-big}} |
|||
{{libheader|Wren-fmt}} |
|||
This is based on the C code in the Wikipedia article and the FreeBASIC entry examples. |
|||
'0' is not actually an example here but is used by FreeBASIC to mark the end of the 'data' statement list so I've ignored that. |
|||
As Wren doesn't natively support unsigned 64-bit arithmetic, I've used the ULong class in the first above named module for this task. |
|||
Even so, there are two examples which fail (1000000000000000127 and 1537228672809128917) because the code is unable to process enough 'multipliers' before an overflow situation is encountered. To deal with this, the program automatically switches to BigInt to process such cases. |
|||
<syntaxhighlight lang="wren">import "./long" for ULong |
|||
import "./big" for BigInt |
|||
import "./fmt" for Fmt |
|||
var 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 |
|||
] |
|||
var squfof = Fn.new { |N| |
|||
var s = ULong.new((N.toNum.sqrt + 0.5).floor) |
|||
if (s*s == N) return s |
|||
for (multiplier in multipliers) { |
|||
var T = ULong |
|||
var n = N |
|||
if (n > ULong.largest/multiplier) { |
|||
T = BigInt |
|||
n = BigInt.new(n.toString) |
|||
} |
|||
var D = n * multiplier |
|||
var P = D.isqrt |
|||
var Pprev = P |
|||
var Po = Pprev |
|||
var Qprev = T.one |
|||
var Q = D - Po*Po |
|||
var L = (s * 8).isqrt.toSmall |
|||
var B = 3 * L |
|||
var i = 2 |
|||
var b = T.zero |
|||
var q = T.zero |
|||
var r = T.zero |
|||
while (i < B) { |
|||
b = (Po + P) / Q |
|||
P = b * Q - P |
|||
q = Q |
|||
Q = Qprev + b * (Pprev - P) |
|||
r = T.new((Q.toNum.sqrt + 0.5).floor) |
|||
if ((i & 1) == 0 && r*r == Q) break |
|||
Qprev = q |
|||
Pprev = P |
|||
i = i + 1 |
|||
} |
|||
if (i < B) { |
|||
b = (Po - P) / r |
|||
Pprev = P = b*r + P |
|||
Qprev = r |
|||
Q = (D - Pprev*Pprev) / Qprev |
|||
i = 0 |
|||
while (true) { |
|||
b = (Po + P) / Q |
|||
Pprev = P |
|||
P = b * Q - P |
|||
q = Q |
|||
Q = Qprev + b * (Pprev - P) |
|||
Qprev = q |
|||
i = i + 1 |
|||
if (P == Pprev) break |
|||
} |
|||
r = T.gcd(n, Qprev) |
|||
if (r != T.one && r != n) return (r is ULong) ? r : ULong.new(r.toString) |
|||
} |
|||
} |
|||
return ULong.zero |
|||
} |
|||
var examples = [ |
|||
"2501", |
|||
"12851", |
|||
"13289", |
|||
"75301", |
|||
"120787", |
|||
"967009", |
|||
"997417", |
|||
"7091569", |
|||
"13290059", |
|||
"42854447", |
|||
"223553581", |
|||
"2027651281", |
|||
"11111111111", |
|||
"100895598169", |
|||
"1002742628021", |
|||
"60012462237239", |
|||
"287129523414791", |
|||
"9007199254740931", |
|||
"11111111111111111", |
|||
"314159265358979323", |
|||
"384307168202281507", |
|||
"419244183493398773", |
|||
"658812288346769681", |
|||
"922337203685477563", |
|||
"1000000000000000127", |
|||
"1152921505680588799", |
|||
"1537228672809128917", |
|||
"4611686018427387877" |
|||
] |
|||
System.print("Integer Factor Quotient") |
|||
System.print("------------------------------------------") |
|||
for (example in examples) { |
|||
var N = ULong.new(example) |
|||
var fact = squfof.call(N) |
|||
var quot = (fact.isZero) ? "fail" : (N / fact).toString |
|||
Fmt.print("$-20s $-10s $s", N, fact, quot) |
|||
}</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> |
</pre> |
Latest revision as of 16:05, 9 February 2024
![Task](http://static.miraheze.org/rosettacodewiki/thumb/b/ba/Rcode-button-task-crushed.png/64px-Rcode-button-task-crushed.png)
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)
ALGOL 68
Assumes LONG INT is long enough - this is true in ALGOL 68G versioins 2 and 3.
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
- 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
C
Based on Wp
From the Wikipedia entry for Shanks's square forms factorization [[2]]
#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);
}
}
}
- 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.
//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);
}
}
- 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 ...
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;
}
}
- Output:
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
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
.
.
FreeBASIC
' ***********************************************
'subject: Shanks's square form factorization:
' ambiguous forms of discriminant 4N
' give factors of N.
'tested : FreeBasic 1.08.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 = culng(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
print
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
- 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.
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)
}
}
- 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
J
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:
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
}}
Task examples:
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
where
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
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;
}
- Output:
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
jq
Adapted from 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
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 ;
The Tasks
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
)
- 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] 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]
Julia
Modified from Wikipedia's article at [[3]]
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
- 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
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}"
- 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
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.
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.
- Output:
Number Mult M/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
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) }
}
- 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) if h