#include <pari/pari.h>
/********************************************************************/
/**                                                                **/
/**       Is q = 2^p - 1 prime ? Straightforward Lucas test        **/
/**                                                                **/
/********************************************************************/
int lucas_isprime(long p);
int lucas_isprime2(long p);
int lucas_isprime3(long p);
int lucas_isprime4(long p);
int bench(long p);

int
main()
{
GEN p;

pari_init(100000000, 2); /* 100MB stack, no prime numbers */
printf("input a prime number: ");
p = lisGEN(stdin);
if (!BSW_psp(p))
err(talker, "sorry, %Z is not prime", p);

printf("... is %s\n", lucas_isprime(itos(p))? "prime": "not prime");
bench(itos(p));
}

/* basic, no garbage collecting */
int
lucas_isprime(long p)
{
GEN q = subis(shifti(gun, p), 1);
GEN u = stoi(4);
long k;
for (k = 3; k <= p; k++)
u = resii(subis(sqri(u), 2), q);
return gcmp0(u);
}

/* "inefficient" generic garbage collecting */
int
lucas_isprime2(long p)
{
pari_sp av0 = avma;
GEN q = subis(shifti(gun, p), 1);
pari_sp av = avma;
GEN u = stoi(4);
long k;
for (k = 3; k <= p; k++)
u = gerepilecopy(av, resii(subis(sqri(u), 2), q));
k = gcmp0(u); avma = av0; return k;
}

/* type-specific garbage collecting */
int
lucas_isprime3(long p)
{
GEN q = subis(shifti(gun, p), 1);
pari_sp av = avma;
GEN u = stoi(4);
long k;
for (k = 3; k <= p; k++)
u = gerepileuptoint(av, resii(subis(sqri(u), 2), q));
k = gcmp0(u); avma = av; return k;
}

/* "random" garbage collecting */
int
lucas_isprime4(long p)
{
GEN q = subis(shifti(gun, p), 1);
pari_sp av = avma, lim = stack_lim(av, 2);
GEN u = stoi(4);
long k;
for (k = 3; k <= p; k++)
{
u = resii(subis(sqri(u), 2), q);
#if 1
if ((k % 16) == 0) u = gerepileuptoint(av, u);
#else /* Variant */
if (low_stack(lim, stack_lim(av,2))) u = gerepileuptoint(av, u);
#endif
}
k = gcmp0(u); avma = av; return k;
}

int
bench(long p)
{
pari_timer T;
pari_sp av = avma;
long i, N = 10000000 / (p*p);

N = max(N, 1);
TIMERstart(&T); printf("N = %ld\n", N);
for (i=0; i < N; i++) { lucas_isprime(p); avma = av; }
avma = av; printf("Time lucas:  %ldms\n", TIMER(&T));
for (i=0; i < N; i++) { lucas_isprime2(p); avma = av; }
avma = av; printf("Time lucas2: %ldms\n", TIMER(&T));
for (i=0; i < N; i++) { lucas_isprime3(p); avma = av; }
avma = av; printf("Time lucas3: %ldms\n", TIMER(&T));
for (i=0; i < N; i++) { lucas_isprime4(p); avma = av; }
avma = av; printf("Time lucas4: %ldms\n", TIMER(&T));
}