#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include <time.h>
#ifdef _WIN32
#include <process.h>
#endif
#include <string.h>
#include "giants.h"
#define D 100
#define NUM_PRIMES 6542
#ifdef _WIN32
#pragma warning( disable : 4127 4706 )
#endif
extern giant scratch2;
int pr[NUM_PRIMES];
giant xr = NULL, xs = NULL, zs = NULL, zr = NULL, xorg = NULL,
zorg = NULL, t1 = NULL, t2 = NULL, t3 = NULL, N = NULL,
gg = NULL, An = NULL, Ad = NULL;
giant xb[D+2], zb[D+2], xzb[D+2];
int modmode = 0, Q, modcount = 0;
void ell_even(giant x1, giant z1, giant x2, giant z2, giant An,
giant Ad, giant N);
void ell_odd(giant x1, giant z1, giant x2, giant z2, giant xor,
giant zor, giant N);
void ell_mul(giant xx, giant zz, int n, giant An, giant Ad, giant N);
int least_2(int n);
void dot(void);
int psi_rand();
int
psi_rand(
void
)
{
unsigned short hi;
unsigned short low;
time_t tp;
int result;
time(&tp);
low = (unsigned short)rand();
hi = (unsigned short)rand();
result = ((hi << 16) | low) ^ ((int)tp);
return (result & 0x7fffffff);
}
void
set_random_seed(
void
)
{
time_t tp;
time(&tp);
srand((int)tp + (int)getpid());
}
int
isprime(
int odd
)
{
int j;
int p;
for (j=1; ; j++)
{
p = pr[j];
if (p*p > odd)
return(1);
if (odd % p == 0)
return(0);
}
}
int
primeq(
int p
)
{
register int j=3;
if ((p&1)==0)
return ((p==2)?1:0);
if (j>=p)
return (1);
while ((p%j)!=0)
{
j+=2;
if (j*j>p)
return(1);
}
return(0);
}
void
s_modg(
giant N,
giant t
)
{
++modcount;
switch (modmode)
{
case 0:
modg(N, t);
break;
case -1:
mersennemod(Q, t);
break;
case 1:
fermatmod(Q, t);
break;
}
}
void
reset_mod(
giant x,
giant N
)
{
divg(x, N);
modmode = 0;
}
void
ell_even(
giant x1,
giant z1,
giant x2,
giant z2,
giant An,
giant Ad,
giant N
)
{
gtog(x1, t1);
addg(z1, t1);
squareg(t1);
s_modg(N, t1);
gtog(x1, t2);
subg(z1, t2);
squareg(t2);
s_modg(N, t2);
gtog(t1, t3);
subg(t2, t3);
gtog(t2, x2);
mulg(t1, x2);
gshiftleft(2, x2);
s_modg(N, x2);
mulg(Ad, x2);
s_modg(N, x2);
mulg(Ad, t2);
gshiftleft(2, t2);
s_modg(N, t2);
gtog(t3, t1);
mulg(An, t1);
s_modg(N, t1);
addg(t1, t2);
mulg(t3, t2);
s_modg(N, t2);
gtog(t2,z2);
}
void
ell_odd(
giant x1,
giant z1,
giant x2,
giant z2,
giant xor,
giant zor,
giant N
)
{
gtog(x1, t1);
subg(z1, t1);
gtog(x2, t2);
addg(z2, t2);
mulg(t1, t2);
s_modg(N, t2);
gtog(x1, t1);
addg(z1, t1);
gtog(x2, t3);
subg(z2, t3);
mulg(t3, t1);
s_modg(N, t1);
gtog(t2, x2);
addg(t1, x2);
squareg(x2);
s_modg(N, x2);
gtog(t2, z2);
subg(t1, z2);
squareg(z2);
s_modg(N, z2);
mulg(zor, x2);
s_modg(N, x2);
mulg(xor, z2);
s_modg(N, z2);
}
void
ell_mul(
giant xx,
giant zz,
int n,
giant An,
giant Ad,
giant N
)
{
unsigned int c = (unsigned int)0x80000000;
if (n==1)
return;
if (n==2)
{
ell_even(xx, zz, xx, zz, An, Ad, N);
return;
}
gtog(xx, xorg);
gtog(zz, zorg);
ell_even(xx, zz, xs, zs, An, Ad, N);
while((c&n) == 0)
{
c >>= 1;
}
c>>=1;
do
{
if (c&n)
{
ell_odd(xs, zs, xx, zz, xorg, zorg, N);
ell_even(xs, zs, xs, zs, An, Ad, N);
}
else
{
ell_odd(xx, zz, xs, zs, xorg, zorg, N);
ell_even(xx, zz, xx, zz, An, Ad, N);
}
c >>= 1;
} while(c);
}
void
choose12(
giant x,
giant z,
int k,
giant An,
giant Ad,
giant N
)
{
itog(k, zs);
gtog(zs, xs);
squareg(xs);
itog(5, t2);
subg(t2, xs);
s_modg(N, xs);
addg(zs, zs);
addg(zs, zs);
s_modg(N, zs);
gtog(xs, x);
squareg(x);
s_modg(N, x);
mulg(xs, x);
s_modg(N, x);
gtog(zs, z);
squareg(z);
s_modg(N, z);
mulg(zs, z);
s_modg(N, z);
gtog(zs, t2);
subg(xs, t2);
gtog(t2, t3);
squareg(t2);
s_modg(N, t2);
mulg(t3, t2);
s_modg(N, t2);
gtog(xs, t3);
addg(t3, t3);
addg(xs, t3);
addg(zs, t3);
s_modg(N, t3);
mulg(t3, t2);
s_modg(N, t2);
gtog(zs, t3);
mulg(xs, t3);
s_modg(N, t3);
squareg(xs);
s_modg(N, xs);
mulg(xs, t3);
s_modg(N, t3);
addg(t3, t3);
addg(t3, t3);
s_modg(N, t3);
gtog(t3, Ad);
gtog(t2, An);
}
void
ensure(
int q
)
{
int nsh, j;
N = newgiant(INFINITY);
if(!q)
{
gread(N,stdin);
q = bitlen(N) + 1;
}
nsh = q/4;
if (!xr)
xr = newgiant(nsh);
if (!zr)
zr = newgiant(nsh);
if (!xs)
xs = newgiant(nsh);
if (!zs)
zs = newgiant(nsh);
if (!xorg)
xorg = newgiant(nsh);
if (!zorg)
zorg = newgiant(nsh);
if (!t1)
t1 = newgiant(nsh);
if (!t2)
t2 = newgiant(nsh);
if (!t3)
t3 = newgiant(nsh);
if (!gg)
gg = newgiant(nsh);
if (!An)
An = newgiant(nsh);
if (!Ad)
Ad = newgiant(nsh);
for (j=0;j<D+2;j++)
{
xb[j] = newgiant(nsh);
zb[j] = newgiant(nsh);
xzb[j] = newgiant(nsh);
}
}
int
bigprimeq(
giant x
)
{
itog(1, t1);
gtog(x, t2);
subg(t1, t2);
itog(5, t1);
powermodg(t1, t2, x);
if (isone(t1))
return(1);
return(0);
}
void
dot(
void
)
{
printf(".");
fflush(stdout);
}
main(
int argc,
char *argv[]
)
{
int j, k, C, nshorts, cnt, count,
limitbits = 0, pass, npr, rem;
long B;
int randmode = 0;
if (!strcmp(argv[argc-1], "-r"))
{
randmode = 1;
if (argc > 4)
limitbits = atoi(argv[argc-2]);
}
else
{
randmode = 0;
}
modmode = 0;
if (argc > 2)
{
modmode = atoi(argv[1]);
Q = atoi(argv[2]);
}
if (modmode==0)
Q = 0;
ensure(Q);
if (modmode)
{
itog(1, N);
gshiftleft(Q, N);
itog(modmode, t1);
addg(t1, N);
}
pr[0] = 2;
for (k=0, npr=1;; k++)
{
if (primeq(3+2*k))
{
pr[npr++] = 3+2*k;
if (npr >= NUM_PRIMES)
break;
}
}
if (randmode == 0)
{
printf("Sieving...\n");
fflush(stdout);
for (j=0; j < NUM_PRIMES; j++)
{
gtog(N, t1);
rem = idivg(pr[j], t1);
if (rem == 0)
{
printf("%d ", pr[j]);
gtog(t1, N);
if (isone(N))
{
printf("\n");
exit(0);
}
else
{
printf("* ");
fflush(stdout);
}
--j;
}
}
if (bigprimeq(N))
{
gout(N);
exit(0);
}
printf("\n");
printf("Commencing Pollard rho...\n");
fflush(stdout);
itog(1, gg);
itog(3, t1); itog(3, t2);
for (j=0; j < 15000; j++)
{
if((j%100) == 0)
{
dot();
gcdg(N, gg);
if (!isone(gg))
break;
}
squareg(t1);
iaddg(2, t1);
s_modg(N, t1);
squareg(t2);
iaddg(2, t2);
s_modg(N, t2);
squareg(t2);
iaddg(2, t2);
s_modg(N, t2);
gtog(t2, t3);
subg(t1, t3);
t3->sign = abs(t3->sign);
mulg(t3, gg);
s_modg(N, gg);
}
gcdg(N, gg);
if ((gcompg(N,gg) != 0) && (!isone(gg)))
{
fprintf(stdout,"\n");
gout(gg);
reset_mod(gg, N);
if (isone(N))
{
printf("\n");
exit(0);
}
else
{
printf("* ");
fflush(stdout);
}
if (bigprimeq(N))
{
gout(N);
exit(0);
}
}
printf("\n");
printf("Commencing Pollard (p-1)...\n");
fflush(stdout);
itog(1, gg);
itog(3, t1);
for (j=0; j< NUM_PRIMES; j++)
{
cnt = (int)(8*log(2.0)/log(1.0*pr[j]));
if (cnt < 2)
cnt = 1;
for (k=0; k< cnt; k++)
{
powermod(t1, pr[j], N);
}
itog(1, t2);
subg(t1, t2);
mulg(t2, gg);
s_modg(N, gg);
if (j % 100 == 0)
{
dot();
gcdg(N, gg);
if (!isone(gg))
break;
}
}
gcdg(N, gg);
if ((gcompg(N,gg) != 0) && (!isone(gg)))
{
fprintf(stdout,"\n");
gout(gg);
reset_mod(gg, N);
if (isone(N))
{
printf("\n");
exit(0);
}
else
{
printf("* ");
fflush(stdout);
}
if (bigprimeq(N))
{
gout(N);
exit(0);
}
}
}
printf("\n");
printf("Commencing ECM...\n");
fflush(stdout);
if (randmode)
set_random_seed();
pass = 0;
while (++pass)
{
if (randmode == 0)
{
if (pass <= 3)
{
B = 1000;
}
else if (pass <= 10)
{
B = 10000;
}
else if (pass <= 100)
{
B = 100000L;
} else
{
B = 1000000L;
}
}
else
{
B = 2000000L;
}
C = 50*((int)B);
cnt = psi_rand();
choose12(xr, zr, cnt, An, Ad, N);
printf("Choosing curve %d, with s = %d, B = %d, C = %d:\n", pass,cnt, B, C); fflush(stdout);
cnt = 0;
nshorts = 1;
count = 0;
for (j=0;j<nshorts;j++)
{
ell_mul(xr, zr, 1<<16, An, Ad, N);
ell_mul(xr, zr, 3*3*3*3*3*3*3*3*3*3*3, An, Ad, N);
ell_mul(xr, zr, 5*5*5*5*5*5*5, An, Ad, N);
ell_mul(xr, zr, 7*7*7*7*7*7, An, Ad, N);
ell_mul(xr, zr, 11*11*11*11, An, Ad, N);
ell_mul(xr, zr, 13*13*13*13, An, Ad, N);
ell_mul(xr, zr, 17*17*17, An, Ad, N);
}
k = 19;
while (k<B)
{
if (isprime(k))
{
ell_mul(xr, zr, k, An, Ad, N);
if (k<100)
ell_mul(xr, zr, k, An, Ad, N);
if (cnt++ %100==0)
dot();
}
k += 2;
}
count = 0;
gtog(zr, gg);
gcdg(N, gg);
if ((!isone(gg))&&(bitlen(gg)>limitbits))
{
fprintf(stdout,"\n");
gwriteln(gg, stdout);
fflush(stdout);
reset_mod(gg, N);
if (isone(N))
{
printf("\n");
exit(0);
}
else
{
printf("* ");
fflush(stdout);
}
if (bigprimeq(N))
{
gout(N);
exit(0);
}
continue;
}
else
{
printf("\n");
fflush(stdout);
}
k = ((int)B)/D;
gtog(xr, xb[0]);
gtog(zr, zb[0]);
ell_mul(xb[0], zb[0], k*D+1 , An, Ad, N);
gtog(xr, xb[D+1]);
gtog(zr, zb[D+1]);
ell_mul(xb[D+1], zb[D+1], (k+2)*D+1 , An, Ad, N);
for (j=1; j <= D; j++)
{
gtog(xr, xb[j]);
gtog(zr, zb[j]);
ell_mul(xb[j], zb[j], 2*j , An, Ad, N);
gtog(zb[j], xzb[j]);
mulg(xb[j], xzb[j]);
s_modg(N, xzb[j]);
}
modcount = 0;
printf("\nCommencing second stage, curve %d...\n",pass); fflush(stdout);
count = 0;
itog(1, gg);
while (1)
{
gtog(zb[0], xzb[0]);
mulg(xb[0], xzb[0]);
s_modg(N, xzb[0]);
mulg(zb[0], gg);
s_modg(N,gg);
for (j = 1; j < D; j++)
{
if (!isprime(k*D+1+ 2*j))
continue;
gtog(xb[0], t1);
subg(xb[j], t1);
gtog(zb[0], t2);
addg(zb[j], t2);
mulg(t1, t2);
s_modg(N, t1);
subg(xzb[0], t2);
addg(xzb[j], t2);
s_modg(N, t2);
--modcount;
mulg(t2, gg);
s_modg(N, gg);
if((++count)%1000==0)
dot();
}
k += 2;
if(k*D > C)
break;
gtog(xb[D+1], xs);
gtog(zb[D+1], zs);
ell_odd(xb[D], zb[D], xb[D+1], zb[D+1], xb[0], zb[0], N);
gtog(xs, xb[0]);
gtog(zs, zb[0]);
}
gcdg(N, gg);
if((!isone(gg))&&(bitlen(gg)>limitbits))
{
fprintf(stdout,"\n");
gwriteln(gg, stdout);
fflush(stdout);
reset_mod(gg, N);
if (isone(N))
{
printf("\n");
exit(0);
}
else
{
printf("* ");
fflush(stdout);
}
if (bigprimeq(N))
{
gout(N);
exit(0);
}
continue;
}
printf("\n");
fflush(stdout);
}
return 0;
}