/* ex: set ff=dos ts=2 et: */
/* Copyright 2008 Ryan Flynn */
/* <##c Skapare> vixey: write a program to calculate all primes that fit in a
 * 32 bit unsigned int without using sieve or division/modulus
 */
/*
 * gcc -W -Wall -O3 -march=prescott -ffast-math -fomit-frame-pointer -o miller-rabin miller-rabin.c
 */

#include <assert.h>
#include <stdio.h>
#include <stdlib.h>
#include <stdint.h>
#include <inttypes.h>
#include <time.h>

typedef uint32_t u32;
typedef uint64_t u64;

/**
 * recursive, very simple version of the Miller-Rabin primality test.
 */
static u32 millertime(u32 a, u32 i, u32 n)
{
  u32 x, y;
  if (i == 0)
    return 1;
  x = millertime(a, i / 2, n);
  if (x == 0)
    return 0;
  y = ((u64)x * x) % n;
  if (y == 1 && x != 1 && x != n - 1)
    return 0;
  if ((i & 1))
    y = ((u64)a * y) % n;
  return y;
}

/**
 * "multiple Rabin tests using the first 7 primes (using 8 gives no improvement)
 * are valid for every number up to 3.4e10(sic)"
 * @ref http://mathworld.wolfram.com/Rabin-MillerStrongPseudoprimeTest.html
 */
static int miller_rabin(u32 n)
{
  static const u32 a[] = { 2, 3, 5, 7, 11, 13 }; /*, 17 };*/
# define N_A (sizeof a / sizeof a[0])
  unsigned i = 0;
  while (
    i < N_A &&
    a[i] < n &&
    1 == millertime(a[i++], n - 1, n)
  );
  return N_A == i || a[i] >= n;
}

static void print(unsigned nth, u32 n)
{
  time_t t = time(NULL);
  struct tm *l = localtime(&t);
  fprintf(stderr, "%04u-%02u-%02u %02u:%02u:%02u %u %" PRIu32 "\n",
    1900 + l->tm_year, 1 + l->tm_mon, l->tm_mday, l->tm_hour, l->tm_min, l->tm_sec, nth, n);
}

static void short_output(void)
{
  u32 i = 5,
      primes_found = 2,
      last_print = 2;
  printf("1 2\n");
  printf("2 3\n");
  while (i > 1) {
    if (miller_rabin(i)) {
      primes_found++;
      if (primes_found == last_print * 2) {
        print(primes_found, i);
        last_print = primes_found;
      }
    }
    i += 2;
  }
}

int main(void)
{
  srand(time(NULL));
  printf("RAND_MAX=%u\n", RAND_MAX);
  assert(miller_rabin(2));
  assert(miller_rabin(3));
  assert( miller_rabin(5));
  assert(!miller_rabin(6));
  assert( miller_rabin(7));
  assert(!miller_rabin(10));
  assert(!miller_rabin(232));
  assert( miller_rabin(257));
  assert( miller_rabin(45821));
  assert( miller_rabin(45821));
  assert( miller_rabin(51539));
  assert( miller_rabin(63499));
  assert( miller_rabin(65521));
  assert( miller_rabin(65543));
  assert( miller_rabin(65537));
  assert( miller_rabin(70853));
  assert( miller_rabin(92431));
  assert( miller_rabin(104729));
  assert( miller_rabin(2146435153));
  assert(!miller_rabin((u32)7 * 127 * 4830073)); /* 4293934897; highest u32 non-prime */
  short_output();
	return 0;
}