#include <iostream>
#include <fstream>
#include <stdlib.h>
#include <stdio.h>
using namespace std;
#include <stdlib.h>
#include <cmath>
#include <givaro/givintprime.h>
#include <givaro/givintnumtheo.h>
#include <givaro/givtimer.h>
#include <givaro/givinit.h>
#include <givaro/modular-integer.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)))))))))))))))
Integer&
MyPollard(GivRandom& gen, Integer& g,
const Integer& n,
const unsigned long threshold)
{
g=1;
Integer m(0U), x, y, t;
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))) {
for(
unsigned long c = 0;
isOne(g) && (++
c < threshold); ) {
x=y;
}
}
return g;
} else {
}
} else {
}
}
unsigned long Revert(
const Integer
p,
const double epsilon,
double firstguess)
{
double t1, t4, t8, t34(firstguess), dL;
t4 = 2.0/t1+1.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)
{
}
bool ProbLucas(
const Integer n,
const double orig_epsilon)
{
#ifdef __GMP_PLUSPLUS__
Integer::seeding( (unsigned long)BaseTimer::seed() );
#endif
GivRandom generator;
Integer Q=n-1,
a,q,tmp(1);
Integer nmu=Q;
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);
Integer trn, r,
b, expo;
root( trn, n, 3); trn *= trn;
q = 2;
expo = nmu;
if (q == 1) {
q = Q;
expo = nmu; expo /= q;
IP.random(generator,
a, 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);
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-1U),n) != 1) {
std::cerr << "Factor found : " << q << std::endl;
return 0;
}
if (q == 1) {
q = Q;
expo = nmu; expo /= q;
IP.random(generator,
a, 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)
{
Integer P;
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);
Timer tim; tim.clear();
for(
unsigned int i = 0; i <
NB; ++i) {
Timer tt; tt.clear(); tt.start();
tt.stop();
tim += tt;
}
std::cout << (a1?"prime":"composite") << endl;
std::cerr << tim << std::endl;
return 0;
}
#define GIVABSDIV(a, b)
Definition: ProbLucas.C:44
Integer & MyPollard(GivRandom &gen, Integer &g, const Integer &n, const unsigned long threshold)
Definition: ProbLucas.C:54
unsigned long Revert(const Integer p, const double epsilon, double firstguess)
Definition: ProbLucas.C:92
#define ProbLucas_factor_second_primes(tmp, n)
Definition: ProbLucas.C:50
#define ProbLucas_factor_first_primes(tmp, n)
Definition: ProbLucas.C:48
IntFactorDom IP
Definition: ProbLucas.C:43
bool ProbLucas(const Integer n, const double orig_epsilon)
Definition: ProbLucas.C:126
int main(int argc, char **argv)
Definition: benchmark-recint_exp.C:20
#define NB
Definition: gfq_atomic.C:32
#define GIVMIN(a, b)
Definition: gfq_atomic.C:24
Namespace in which the whole Givaro library resides.
Definition: all_field.C:23
bool root(Integer &q, const Integer &, uint32_t n)
Definition: gmp++_int_misc.C:71
int32_t isOne(const Integer &a)
Definition: gmp++_int_compare.C:523
Integer gcd(const Integer &a, const Integer &b)
Definition: gmp++_int_gcd.C:43
Integer & pow(Integer &Res, const Integer &n, const int64_t l)
Definition: gmp++_int_pow.C:61
double logtwo(const Integer &a)
Definition: gmp++_int_misc.C:127
Integer sqrt(const Integer &p)
Definition: gmp++_int_misc.C:59
int32_t isZero(const Integer &a)
Definition: gmp++_int_compare.C:538
MG & a
Definition: rmadd.h:143
MG const rmint< K, MG > const T & c
Definition: rmadd.h:143
MG const rmint< K, MG > & b
Definition: rmadd.h:143
a p
Definition: rmadd.h:147