Square form factorization: Difference between revisions
C
m (timings) |
(C) |
||
Line 50:
__TOC__
=={{header|C}}==
From the Wikipedia entry for Shanks's square forms factorization [[//en.wikipedia.org/wiki/Shanks%27s_square_forms_factorization.]]
<lang c>#include <math.h>
#include <stdio.h>
#define nelems(x) (sizeof(x) / sizeof((x)[0]))
const unsigned long multiplier[] = {1, 3, 5, 7, 11, 3*5, 3*7, 3*11, 5*7, 5*11, 7*11, 3*5*7, 3*5*11, 3*7*11, 5*7*11, 3*5*7*11};
unsigned long long gcd(unsigned long long a, unsigned long long b)
{
while (b != 0)
{
a %= b;
a ^= b;
b ^= a;
a ^= b;
}
return a;
}
unsigned long long SQUFOF( unsigned long long N )
{
unsigned long long D, Po, P, Pprev, Q, Qprev, q, b, r, s;
unsigned long L, B, i;
s = (unsigned long long)(sqrtl(N)+0.5);
if (s*s == N) return s;
for (int k = 0; k < nelems(multiplier) && N <= 0xffffffffffffffff/multiplier[k]; k++) {
D = multiplier[k]*N;
Po = Pprev = P = sqrtl(D);
Qprev = 1;
Q = D - Po*Po;
L = 2 * sqrtl( 2*s );
B = 3 * L;
for (i = 2 ; i < B ; i++) {
b = (unsigned long long)((Po + P)/Q);
P = b*Q - P;
q = Q;
Q = Qprev + b*(Pprev - P);
r = (unsigned long long)(sqrtl(Q)+0.5);
if (!(i & 1) && r*r == Q) break;
Qprev = q;
Pprev = P;
};
if (i >= B) continue;
b = (unsigned long long)((Po - P)/r);
Pprev = P = b*r + P;
Qprev = r;
Q = (D - Pprev*Pprev)/Qprev;
i = 0;
do {
b = (unsigned long long)((Po + P)/Q);
Pprev = P;
P = b*Q - P;
q = Q;
Q = Qprev + b*(Pprev - P);
Qprev = q;
i++;
} while (P != Pprev);
r = gcd(N, Qprev);
if (r != 1 && r != N) return r;
}
return 0;
}
int main(int argc, char *argv[]) {
int i;
const unsigned long long data[] = {
2501,
12851,
13289,
75301,
120787,
967009,
997417,
7091569,
13290059,
42854447,
223553581,
2027651281,
11111111111,
100895598169,
1002742628021,
60012462237239,
287129523414791,
9007199254740931,
11111111111111111,
314159265358979323,
384307168202281507,
419244183493398773,
658812288346769681,
922337203685477563,
1000000000000000127,
1152921505680588799,
1537228672809128917,
4611686018427387877};
for(int i = 0; i < nelems(data); i++) {
unsigned long long example, factor, quotient;
example = data[i];
factor = SQUFOF(example);
if(factor == 0) {
printf("%llu was not factored.\n", example);
}
else {
quotient = example / factor;
printf("Integer %llu has factors %llu and %llu\n",
example, factor, quotient);
}
}
}
</lang>{{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>
| |||