#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));
}