/* ex: set 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 -ffast-math -fomit-frame-pointer -fno-math-errno -o miller-rabin-distributed miller-rabin-distributed.c
 *
 *
 * pizza@debian:~/proj/primes-challenge$ gcc -W -Wall -O9 -ffast-math -fomit-frame-pointer -o miller-rabin-distributed miller-rabin-distributed.c
 * miller-rabin-distributed.c:28: warning: 'millertime2' defined but not used
 * pizza@debian:~/proj/primes-challenge$ time ./miller-rabin-distributed 1                                                         Using 1 workers...
 * [ 2693] 2012-03-23 13:54:35 start=2, stop=4294967295
 * [ 2693] 2012-03-23 13:54:35 launching 1 workers (4294967293 numbers each)...
 * [ 2694] 2012-03-23 13:54:35 searching 2-4294967295...
 * [ 2694] 2012-03-23 14:44:26 found 203280221 primes in 285564129 operations (0.07 ops/candidate, 1.40 ops/prime)
 * [ 2693] 2012-03-23 14:44:26 my work here is done.
 *
 * real    49m50.822s
 * user    49m31.686s
 * sys     0m10.516s
 * pizza@debian:~/proj/primes-challenge$
 *
 * if i use the first 5 primes instead of 7 we complete in 41 minutes but find 203280222 primes! one off!
 */

#include <assert.h>
#include <signal.h>
#include <stdio.h>
#include <stdlib.h>
#include <stdint.h>
#include <stdarg.h>
#include <sys/types.h>
#include <sys/wait.h>
#include <time.h>
#include <unistd.h>


static uint32_t Ops = 0; /* how many operations did i perform in my search? */

/* non-recursive, no faster */
static uint32_t millertime2(uint32_t a, uint32_t i, uint32_t n)
{
  uint32_t x[33], I[33];
  int k = 1;
  Ops++;
  x[0] = i;
  x[1] = i;
  I[0] = x[0];
  I[1] = x[1];
  while (x[k] > 1) {
     x[k+1] = x[k] / 2;
     I[k+1] = x[k+1];
    k++;
  }
  while (k--) {
    x[k] = ((uint64_t)x[k+1] * x[k+1]) % n;
    if ((x[k] == 1) & (x[k+1] != 1) & (x[k+1] != n - 1))
      return 0;
    if ((I[k+1] & 1))
      x[k] = ((uint64_t)a * x[k]) % n;
    if (0 == x[k])
      return 0;
  }
  return x[0];
}

/**
 * recursive, very simple version of the Miller-Rabin primality test.
 * @ref http://utilitymill.com/edit/Miller_Rabin_Primality_Test
 */
static uint32_t millertime(uint32_t a, uint32_t i, uint32_t n)
{
  uint32_t x, y;
  Ops++;
  if (i == 0)
    return 1;
  x = millertime(a, i / 2, n);
  if (x == 0)
    return 0;
  y = ((uint64_t)x * x) % n;
  if (y == 1 && x != 1 && x != n - 1)
    return 0;
  if ((i & 1))
    y = ((uint64_t)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(uint32_t n)
{
  static const uint32_t a[] = { 2, 3, 5, 7, 11, 13 }; /*, 17 };*/
  int i = 0;
  while (i < (sizeof a / sizeof a[0]) && a[i] < n && 1 == millertime(a[i++], n - 1, n))
    ;
  return (sizeof a / sizeof a[0]) == i || a[i] >= n;
}

static pid_t MyPid = 0;

/**
 * print with pid and date
 */
static void out(const char *s, ...)
{
  va_list args;
  time_t t = time(NULL);
  struct tm *l = localtime(&t);
  va_start(args, s);
  fprintf(stdout, "[%5u] %04u-%02u-%02u %02u:%02u:%02u ",
    (unsigned)MyPid, 1900 + l->tm_year, 1 + l->tm_mon, l->tm_mday, l->tm_hour, l->tm_min, l->tm_sec);
  vfprintf(stdout, s, args);
  va_end(args);
}

/**
 * search a range (inclusive) of integers for primes
 * @param min must be odd!
 */
static void search_range(uint32_t min, uint32_t max)
{
  uint32_t n = min,
           found = 0;
  out("searching %lu-%lu...\n",
    (unsigned long)n, (unsigned long)max);
  found += miller_rabin(n);
  n += 1 + (min & 1); /* go to next odd number */
  while (n <= max && n > 1) { /* prevent wraparound */
    found += miller_rabin(n);
    n += 2;
  }
  out("found %u primes in %lu operations (%.2f ops/candidate, %.2f ops/prime)\n",
    (unsigned long)found, (unsigned long)Ops, (double)Ops / (max - min), (double)Ops / found);
}

#define WORKERS_MAX 16
static int Worker_Cnt = 1;
static pid_t Workers[WORKERS_MAX];

/**
 * launch workers, divy up range between them
 */
static void run(void)
{
  uint32_t start = 2,
           stop = 0xFFFFFFFF,
           range;
  int w;
  MyPid = getpid();
  out("start=%lu, stop=%lu\n",
    (unsigned long)start, (unsigned long)stop);
  range = (stop - start) / Worker_Cnt;
  out("launching %d workers (%lu numbers each)...\n",
    Worker_Cnt, (unsigned long)range);
  for (w = 0; w < Worker_Cnt; w++) {
    pid_t pid = fork();
    if (0 == pid) { /* worker */
      MyPid = getpid();
      if (w == Worker_Cnt - 1 && stop - range > start) /* integer division can make us come up short... */
        range = stop - start;                          /* ensure we cover our range */
      search_range(start, start + range);
      exit(EXIT_SUCCESS);
    } else { /* boss */
      Workers[w] = pid;
      start += range;
    }
  }
  /* collect workers */
  for (w = 0; w < Worker_Cnt; w++)
    waitpid(Workers[w], NULL, 0);
  out("my work here is done.\n");
}

/**
 * make sure we work
 */
static void sanity_check(void)
{
  uint32_t notprime = (uint32_t)7 * 127 * 4830073; /* 4293934897, very high 32-bit composite */
  assert( miller_rabin(5));
  assert( miller_rabin(5));
  assert( miller_rabin(2));
  assert( miller_rabin(3));
  assert(!miller_rabin(6));
  assert( miller_rabin(7));
  assert(!miller_rabin(9));
  assert(!miller_rabin(10));
  assert(!miller_rabin(91));
  assert(!miller_rabin(121));
  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(notprime));
}

int main(int argc, char *argv[])
{
  if (2 == argc) {
    Worker_Cnt = atoi(argv[1]);
    assert(Worker_Cnt > 0 && Worker_Cnt <= WORKERS_MAX);
  }
  printf("Using %d workers...\n", Worker_Cnt);
  sanity_check();
  run();
  return 0;
}