/*
  Author:  Pate Williams (c) 1997

  "Exercise 4.1. Strong pseudoprime test.
  Write the function abmodn(a, b, n) hinted at in
  the program above and incorporate it in a complete
  program for your computer. In order to save
  computing time, precede the psuedoprime test by
  trial divisions by the small primes <= G,
  thereby discarding many composites before the
  more time-consuming pseudoprime test is applied.
  Make computer experiments with various values
  of G in order to find its optimal value. Use the
  integers between 10 ^ 8 and 10 ^ 8 + 100 as test
  numbers in these experiments." -Hans Riesel-
  See "Prime Numbers and Computer Methods for
  Factorization" by Hans Riesel second edition
  page 94.
*/

#include <stdio.h>
#include <string.h>
#include "lip.h"

#define G 11

int test(long a, long d, long n, long n1, long s)
{
  int found;
  long r, t = 2;
  verylong za = 0, zb = 0, zn = 0;

  zintoz(a, &za);
  zintoz(n, &zn);
  zsexpmod(za, d, zn, &zb);
  found = zscompare(zb, 1) == 0 || zscompare(zb, n1) == 0;
  if (!found) {
    for (r = 1; !found && r < s; r++) {
      zsexpmod(za, d * t, zn, &zb);
      found = zscompare(zb, n1) == 0;
      t <<= 1;
    }
  }
  zfree(&za);
  zfree(&zb);
  zfree(&zn);
  return found;
}

void convert(char *s, verylong *zn)
{
  long i;
  verylong za = 0;

  zzero(zn);
  for (i = 0; i < strlen(s); i++) {
    zsmul(*zn, 10, &za);
    zsadd(za, s[i] - '0', zn);
  }
  zfree(&za);
}

int strong_pseudoprime_test(long n)
{
  int found;
  long d = n - 1, i, n1 = d, s = 0;
  long p[G] = {2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31};
  verylong za = 0, zn = 0;

  for (i = 0; i < G; i++)
    if (n % p[i] == 0) return 0;
  while (!(d & 1)) d >>= 1, s++;
  if (!test(2, d, n, n1, s)) return 0;
  if (!test(3, d, n, n1, s)) return 0;
  if (n < 1373653l) return 1;
  if (!test(5, d, n, n1, s)) return 0;
  if (n < 25326001l) return 1;
  if (!test(7, d, n, n1, s)) return 0;
  zintoz(n, &zn);
  convert("3215031751", &za);
  if (zcompare(zn, za) < 0) {
    zfree(&za);
    zfree(&zn);
    return 1;
  }
  if (!test(11, d, n, n1, s)) {
    zfree(&za);
    zfree(&zn);
    return 0;
  }
  convert("2152302898747", &za);
  if (zcompare(zn, za) < 0) {
    zfree(&za);
    zfree(&zn);
    return 1;
  }
  if (!test(13, d, n, n1, s)) {
    zfree(&za);
    zfree(&zn);
    return 0;
  }
  convert("3474749660383", &za);
  if (zcompare(zn, za) < 0) {
    zfree(&za);
    zfree(&zn);
    return 1;
  }
  if (!test(17, d, n, n1, s)) {
    zfree(&za);
    zfree(&zn);
    return 0;
  }
  convert("3415500717128321", &za);
  found = zcompare(za, zn) < 0;
  zfree(&za);
  zfree(&zn);
  return found;
}

int main(void)
{
  long a = 100000000l, b = 100000100l, c = 0, i;

  for (i = a; i <= b; i++) {
    if (strong_pseudoprime_test(i)) {
      c++;
      printf("%9ld ", i);
      if (c % 5 == 0) printf("\n");
    }
  }
  if (c % 5 != 0) printf("\n");
  printf("total number of primes = %ld\n", c);
  return 0;
}