Author: James Wanless (c) 2000-9
*/
#include
#include
#include
#include
#include
#include
gmp_randclass r (gmp_randinit_default);
mpz_class rand2(long bits)
/* returns random bits-bit integer */
{
return r.get_z_bits(bits);
}
int Rabin_Miller(mpz_class n)
/* given an integer n >= 3 returns 0 composite, 1 probably prime */
{
return mpz_probab_prime_p(n.get_mpz_t(), 10);
}
mpz_class gcd(mpz_class a, mpz_class b)
/* returns greatest common divisor of a and b */
{
mpz_t temp;
mpz_init(temp);
mpz_gcd(temp, a.get_mpz_t(), b.get_mpz_t());
mpz_class temp_class(temp);
mpz_clear(temp);
return temp_class;
}
mpz_class exp_mod(mpz_class x, mpz_class b, mpz_class n)
/* returns x ^ b mod n */
{
mpz_t temp;
mpz_init (temp);
mpz_powm(temp, x.get_mpz_t(), b.get_mpz_t(), n.get_mpz_t());
mpz_class temp_class(temp);
mpz_clear (temp);
return temp_class;
}
bool Wanless (mpz_class n)
{
mpz_class basestested = 0;
mpz_class T = 1;
mpz_class b = 1;
mpz_class A = 2;
mpz_class prodA = 2;
mpz_class wanless = 2;
long nbits = 1;
if (Rabin_Miller(n)) {
cout << n << "\n";
fflush(stdout);
return true;
}
if (n < 2)
return false;
if (n % 2 == 0) {
cout << "2\n";
fflush(stdout);
return Wanless(n / 2);
}
while (wanless < n) {
wanless *= 2;
nbits ++;
}
while (T == 1 || T == n) {
basestested++;
cout << "base#" << basestested << "\r";
fflush(stdout);
prodA = 1;
for (long i = 0; i < nbits; i++) {
A = rand2(nbits);
prodA *= A;
prodA = prodA % n;
}
b = exp_mod(prodA, wanless, n);
T = gcd(n, exp_mod(b, n, n) - b);
}
if (T > 1 && T < n) {
cout << "\n" << T << "\t[prodA=" << prodA << "]\n";
fflush(stdout);
Wanless(n / T);
}
return (false);
}
int main(void)
{
mpz_class N;
for (;;) {
cout << "number to be tested or 0 to quit:\n";
cin >> N;
if (N == 0) break;
Wanless(N);
}
return 0;
}
1 comment:
Once again (as expected) the #include's haven't survived the transfer to blog intact - but I imagine you can soon guess them to recreate the correct code...
J
Post a Comment