#include <iostream>
#include <fstream>
#include <stdlib.h>
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <givaro/givintprime.h>
#include <givaro/givintnumtheo.h>
#include <givaro/givtimer.h>
#include <givaro/givinit.h>
#include <givaro/givzpzInt.h>
#include <givaro/givinteger.h>
#include <givaro/givrandom.h>
#include <givaro/givinteger.h>
#define GIVABSDIV(a,b) ((a)<(b)?((a)/(b)):((b)/(a)))
#ifndef GIVMIN
#define GIVMIN(a,b) ((a)<(b)?(a):(b))
#endif
#define ProbLucas_factor_first_primes(tmp,n) (tmp = IP.isZero(IP.mod(tmp,n,23))?23:( IP.isZero(IP.mod(tmp,n,19))?19:( IP.isZero(IP.mod(tmp,n,17))?17: (IP.isZero(IP.mod(tmp,n,2))?2:( IP.isZero(IP.mod(tmp,n,3))?3:( IP.isZero(IP.mod(tmp,n,5))?5:( IP.isZero(IP.mod(tmp,n,7))?7: ( IP.isZero(IP.mod(tmp,n,11))?11:13 ))))))))
#define ProbLucas_factor_second_primes(tmp,n) (tmp = IP.isZero(IP.mod(tmp,n,31))?31:( IP.isZero(IP.mod(tmp,n,29))?29: ( IP.isZero(IP.mod(tmp,n,37))?37: ( IP.isZero(IP.mod(tmp,n,41))?41:( IP.isZero(IP.mod(tmp,n,43))?43: ( IP.isZero(IP.mod(tmp,n,71))?71:( IP.isZero(IP.mod(tmp,n,67))?67:( IP.isZero(IP.mod(tmp,n,61))?61:( IP.isZero(IP.mod(tmp,n,59))?59: ( IP.isZero(IP.mod(tmp,n,53))?53:( IP.isZero(IP.mod(tmp,n,47))?47: ( IP.isZero(IP.mod(tmp,n,97))?97: ( IP.isZero(IP.mod(tmp,n,89))?89:( IP.isZero(IP.mod(tmp,n,83))?83:( IP.isZero(IP.mod(tmp,n,79))?79:73)))))))))))))))
{
g=1;
static Integer PROD_first_primes(223092870);
static Integer PROD_second_primes(
"10334565887047481278774629361");
if (
isOne(
gcd(y,n,PROD_first_primes))) {
if (
isOne(
gcd(y,n,PROD_second_primes))) {
unsigned long p(1);
for(
unsigned long c = 0;
isOne(g) && (++c < threshold); ) {
if( p == c ) {
x=y;
p <<= 1;
}
}
return g;
} else {
return ProbLucas_factor_second_primes(g,n);
}
} else {
return ProbLucas_factor_first_primes(g,n);
}
}
unsigned long Revert(
const Integer p,
const double epsilon,
double firstguess)
{
double t1, t4, t8, t34(firstguess), dL;
t1 = (double)p-1.0;
t4 = 2.0/t1+1.0;
t8 =
logtwo(p)*log(2.0)-log(2.0);
do {
double t3, t5, t7, t9, t11, t12, t18, t22, t23 ;
dL = t34;
t3 = 1.0/dL;
t5 = t3*t3;
t7 = 1.0-t5;
t9 = 2.0*log(dL);
t11 = t8/t9;
t18 = t9*t9;
t22 = log(t7);
t23 = (-t8/t18*t3*t22+t11*t5*t3/t7);
t34 = dL-( 1.0 - (1.0-epsilon)/(t4*t12))/(2.0*t23);
} while( (GIVABSDIV(t34,dL) < 0.95) && (dL<1048576.0) );
return (unsigned long)GIVMIN(dL,1048576.0);
}
unsigned long Revert(
const Integer p,
const double epsilon)
{
return Revert(p, epsilon, ::
sqrt(1.2/epsilon));
}
bool ProbLucas(
const Integer n,
const double orig_epsilon)
{
#ifdef __GMP_PLUSPLUS__
Integer::seeding( (unsigned long)BaseTimer::seed() );
#endif
double P = 1.0;
double epsilon = orig_epsilon;
unsigned long s=0;
for( ; !( (int)Q & 0x1) ; Q>>=1, ++s) { }
for(unsigned int i=0; ; ) {
IP.random(generator,a,n);
if (tmp!=1) {
if (tmp != nmu) {
for(i=(unsigned int) s; --i>0;) {
tmp *= tmp;
tmp %= n;
if ((tmp == 1) || (tmp == nmu))
break;
}
if (tmp != nmu)
return 0;
if (i == 0)
break;
}
else {
if (s == 1) break;
}
}
epsilon *= 4.0;
P /= 2.0;
}
unsigned long L = Revert(Q,epsilon);
q = 2;
expo = nmu;
MyPollard(generator, q, Q, L);
if (q == 1) {
q = Q;
expo = nmu; expo /= q;
IP.random(generator, a, n);
IP.powmod(tmp,a,expo,n);
if (tmp == 1) {
P /= (double)(q);
if (P < epsilon) {
std::cerr << "Composite with probability > 1-" << P << std::endl;
return 0;
}
}
Q = 1;
} else {
expo = nmu; expo /= q;
r=0;
Integer::divexact(b, Q, q);
Integer::divmod( b, r, Q, q );
}
L = Revert(Q,epsilon, (double)L);
}
while ( Q > trn ) {
IP.random(generator, a, n);
IP.powmod(tmp,a,expo,n);
if (tmp == 1) {
P /= (double)(q);
if (P < epsilon) {
std::cerr << "Composite with probability > 1-" << P << std::endl;
return 0;
}
} else {
if (IP.
gcd(q,(tmp-1UL),n) != 1) {
std::cerr << "Factor found : " << q << std::endl;
return 0;
}
MyPollard(generator, q, Q, L);
if (q == 1) {
q = Q;
expo = nmu; expo /= q;
IP.random(generator, a, n);
IP.powmod(tmp,a,expo,n);
if (tmp == 1) {
P /= (double)(q);
if (P < epsilon) {
std::cerr << "Composite with probability > 1-" << P << std::endl;
return 0;
}
}
break;
} else {
expo = nmu; expo /= q;
r=0;
Integer::divexact(b, Q, q);
Integer::divmod( b, r, Q, q );
}
L = Revert(Q,epsilon, (double)L);
}
}
}
std::cerr << "Prime with probability > 1-" << orig_epsilon << std::endl;
return 1;
}
int main (
int argc,
char * * argv)
{
if (argc > 1) P =
Integer(argv[1]);
else std::cin >> P;
double epsilon = argc > 2 ? atof(argv[2]) : 0.000001;
unsigned int NB = argc > 3 ? (unsigned int)atoi(argv[3]) : 1U;
bool a1(true);
for(unsigned int i = 0; i < NB; ++i) {
a1 = ProbLucas(P, epsilon);
tim += tt;
}
std::cout << (a1?"prime":"composite") << endl;
std::cerr << tim << std::endl;
return 0;
}