/* Look for primes of the form Rp = (b^p+-1)/(b+-1).  Such primes have p ones
   in their base-b representation.  If p has factors, so does Rp.

   We use GMP's native operations for pseudoprime tests (Miller-Rabin), as it
   seems difficult to adopt Lucas' test for these numbers.

   In order to avoid calling the expensive Miller-Rabin tests in some cases,
   we make a lame attempt to find small factors.  Such factors are of the form
   2kp+1, and must be congruent to 1, 3, 9, 13, 27, 31, 37 or 39 mod 40.  */

#include <sys/types.h>
#include <sys/wait.h>
#include <unistd.h>
#include <stdlib.h>
#include <string.h>
#include <stdio.h>
#include "gmp.h"

#ifndef BASE
#define BASE 10
#endif

#if defined (REVERSE_REPUNIT) && REVERSE_REPUNIT - 0 != 0
#undef REVERSE_REPUNIT
#define REVERSE_REPUNIT 1
#else
#define REVERSE_REPUNIT 0
#endif

#if REVERSE_REPUNIT
#define PLUSMINUS +
#else
#define PLUSMINUS -
#endif

#if __GNUC__
typedef unsigned int UQItype    __attribute__ ((mode (QI)));
typedef unsigned int UHItype    __attribute__ ((mode (QI)));
typedef          int SItype     __attribute__ ((mode (SI)));
typedef unsigned int USItype    __attribute__ ((mode (SI)));
typedef          int DItype     __attribute__ ((mode (DI)));
typedef unsigned int UDItype    __attribute__ ((mode (DI)));
#else
typedef unsigned char UQItype;
typedef	unsigned short int UHItype;
typedef	int SItype;
typedef unsigned int USItype;
typedef long long int DItype;
typedef long long unsigned int UDItype;
#endif

#if ! defined (ARITH)
#error "need to define ARITH to 32 or 64"
#endif

#if ARITH == 64
typedef UDItype UWtype;
typedef DItype SWtype;
typedef UHItype UHWtype;
#define W_TYPE_SIZE 64
#endif

#if ARITH == 32
typedef USItype UWtype;
typedef SItype SWtype;
typedef UHItype UHWtype;
#define W_TYPE_SIZE 32
#endif

#include "longlong.h"

#if defined (_MIPS_SIM)
#if _MIPS_SIM == _ABIN32
typedef unsigned long long word;
#define DEFINED_word
#endif
#endif

#ifndef DEFINED_word
typedef unsigned long word;
#endif

char tab40[] = {0,1,0,1,0,0,0,0,0,1,0,0,0,1,0,0,0,0,0,0,
		0,0,0,0,0,0,0,1,0,0,0,1,0,0,0,0,0,1,0,1};
char tab105[] = {0,1,1,0,1,0,0,0,1,0,0,1,0,1,0,0,1,1,0,1,0,
		 0,1,1,0,0,1,0,0,1,0,1,1,0,1,0,0,1,1,0,0,1,
		 0,1,1,0,1,1,0,0,0,0,1,1,0,0,0,0,1,1,0,1,1,
		 0,1,0,0,1,1,0,0,1,0,1,1,0,1,0,0,1,0,0,1,1,
		 0,0,1,0,1,1,0,0,1,0,1,0,0,1,0,0,0,1,0,1,1};


int powmod ();

struct proc_info
{
  pid_t pid;
  int fd;
  unsigned int p;
};

struct proc_info *proc_info;
int max_n_processes = 1;

void start_proc ();
unsigned int check_one ();
unsigned int mpz_ones_composite_p ();
unsigned int findfactor ();
int isprime ();

int
main (int argc, char **argv)
{
  unsigned int from, to;
  int n_processes;
  int i;
  unsigned int p;
  char out_file[20];

  from = strtoul (argv[1], 0, 0);
  to = strtoul (argv[2], 0, 0);
  if (argc == 4)
    max_n_processes = strtol (argv[3], 0, 0);
  sprintf (out_file, "repunit%d.out", BASE);
  unlink (out_file);

  n_processes = 0;

  proc_info = malloc (max_n_processes * sizeof *proc_info);

  for (i = 0; i < max_n_processes; i++)
    proc_info[i].pid = -1;

  p = from;
  while (1)
    {
      pid_t pid;
      int status;

      /* Start more processes if we don't already have enough.  */

      while (n_processes != max_n_processes)
	{
#ifdef NOT_JUST_PRIME_EXPONENTS
	  if ((p & 1) == 0)
	    p += 1;
#else
	  while (! isprime (p))
	    p += 1;
#endif
	  if (p > to)
	    break;

	  start_proc (p);
	  n_processes++;

	  p += 1;
	}

      /* We have the desired # of processes.  Wait for the first one to
	 terminate and store its result.  */

#ifdef NOFORK
      if (n_processes == 0)
	break;
#else
      pid = wait (&status);
      if (pid == -1)
	break;
      if (! WIFEXITED (status))
	{
	  fprintf (stderr, "subprocess was killed (signal %d)\n",
		   WTERMSIG (status));
	  kill (0, WTERMSIG (status));
	}
#endif

      for (i = 0; i < max_n_processes; i++)
	{
#ifndef NOFORK
	  if (proc_info[i].pid == pid)
#endif
	    {
	      char buf[11];
	      FILE *ofs;

	      read (proc_info[i].fd, buf, 10);
	      close (proc_info[i].fd);
	      ofs = fopen (out_file, "a+");
	      fprintf (ofs, "%u %s\n", proc_info[i].p, buf);
	      fclose (ofs);
	      proc_info[i].pid = -1; /* mark slot as free */
	      n_processes--;
	    }
	}
    }

  exit (0);
}

void
start_proc (unsigned int p)
{
  int channel[2];
  pid_t pid;
  int i;

  for (i = 0; i < max_n_processes; i++)
    {
      if (proc_info[i].pid == -1)
	break;
    }

  pipe (channel);
#ifdef NOFORK
  for (pid = 1; pid != (pid_t) -1; pid--)
#else
  pid = fork ();
#endif
  if (pid != 0)
    {
      /* Parent.  */
#if DEBUG
      fprintf (stderr, "starting: %d (%d) pid=%d\n", p, i, pid);
#endif
#ifndef NOFORK
      close (channel[1]);
#endif
      proc_info[i].fd = channel[0];
      proc_info[i].pid = pid;
      proc_info[i].p = p;
    }
  else
    {
      /* Child.  */
      unsigned int factor;
      char buf[11];

#ifndef NOFORK
      close (channel[0]);
#endif
      factor = check_one (p);
      if (factor == 0)
	sprintf (buf, "y");
      else if (factor == (unsigned int) -1)
	sprintf (buf, "n");
      else
	sprintf (buf, "%u", factor);
      write (channel[1], buf, strlen (buf) + 1);
      close (channel[1]);

#ifndef NOFORK
      _exit (0);
#endif
    }
}

unsigned int
check_one (unsigned int p)
{
  int k;
  unsigned int klimit;
  unsigned long long int k2;
  mpz_t t;
  int result;

  mpz_init (t);

#if REVERSE_REPUNIT
  mpz_ui_pow_ui (t, (unsigned long int) (BASE), (unsigned long int) p);
  mpz_add_ui (t, t, 1L);
  mpz_tdiv_q_ui (t, t, (unsigned long int) ((BASE) + 1));
#else
  mpz_ui_pow_ui (t, (unsigned long int) (BASE), (unsigned long int) p);
  mpz_tdiv_q_ui (t, t, (unsigned long int) ((BASE) - 1));
#endif

  /* Set the factoring search limit approximately propotional to the time
     we will need for testing primality.  */
#ifdef FACTOR_AGGRESSIVELY
  k2 = p * p;
  if (k2 != (unsigned int) k2)
    klimit = (unsigned int) ~0;
  else
    {
      klimit = k2;
      if (klimit < 10000000)
	klimit = 10000000;
    }
#else
  k2 = p * p / 5;
  if (k2 != (unsigned int) k2)
    klimit = (unsigned int) ~0;
  else
    {
      klimit = k2;
      if (klimit < 100000)
	klimit = 100000;
    }
#endif

  k = findfactor (p, klimit);
  if (k)
    {
      mpz_t f;
      mpz_init (f);
      /* Check is the found factor is actually equal to the checked number.
	 In that case we really have a prime (not even a prp).  */
      mpz_set_ui (f, (unsigned long int) k);
      mpz_mul_ui (f, f, (unsigned long int) 2 * p);
      mpz_add_ui (f, f, 1L);
      if (mpz_cmp (t, f) == 0)
	k = 0;		/* prime/prp */
      mpz_clear (f);
      mpz_clear (t);
      return k;
    }

#if NOPRIMETEST
  mpz_clear (t);
  return (unsigned int) -1;
#endif

  if (mpz_cmp_ui (t, 10000) < 0)
    result = isprime ((unsigned int) mpz_get_ui (t));
  else
    {
      result = fermat_pseudoprime ((BASE) == 2 ? 3 : 2, t);
      if (result != 0)
	result = mpz_probab_prime_p (t, 3);
    }
  mpz_clear (t);
  return result != 0 ? 0 : (unsigned int) -1;
}

unsigned int
findfactor (unsigned int p, unsigned int klimit)
{
  static word pows[ARITH + 1];
  static int digits_per_word = 0;

  if (digits_per_word == 0)
    {
      word x = 1;
      word xx;
      int i;
      for (i = 0;; i++)
	{
	  pows[i] = x;
	  xx = x * BASE;
	  if (xx / BASE != x)
	    {
	      digits_per_word = i;
	      break;
	    }
	  x = xx;
	}
    }

  {
    unsigned int k;
    unsigned int pprev, pprem;
    int cnt;
    word powinit;

    pprev = 0;
    pprem = p;
    cnt = 0;
    while (pprem > digits_per_word)
      {
	pprev = (pprev << 1) | (pprem & 1);
	pprem >>= 1;
	cnt++;
      }

    powinit = pows[pprem];

    /* Look for k such that 2kp+1 | (b^p+1)/(b+1).  */
    for (k = 1; k < klimit; k++)
      {
#if ARITH == 64
	word d;
	int r;
	d = (word) 2 * p * k + 1;
#if BASE == 10 && ! REVERSE_REPUNIT
	if (! tab40[d % 40])
	  continue;
#endif
	if (! tab105[d % 105] && d > 7)
	  continue;
	r = powmod (pprev, powinit, cnt, d);
	/* If we found a factor, return it.  But watch out for factors that
	   just divide BASE+-1.  */
	if (r != 0 && (BASE PLUSMINUS 1) % d != 0)
	  return k;
#endif

#if ARITH == 32
	word dh, dl;
	int r;
	word i, dummy;
	umul_ppmm (dh, dl, k, 2 * p);
	add_ssaaaa (dh, dl, dh, dl, 0, 1);

#if BASE == 10 && ! REVERSE_REPUNIT
	if (UDIV_NEEDS_NORMALIZATION)
	  {
	    word xh, xl;
	    xh = dh % 40;  xh = (xh << 26) | (dl >> 6);  xl = dl << 26;
	    udiv_qrnnd (dummy, i, xh, xl, 40 << 26);
	    i >>= 26;
	  }
	else
	  udiv_qrnnd (dummy, i, dh % 40, dl, 40);
	if (! tab40[i])
	  continue;
#endif
	if (UDIV_NEEDS_NORMALIZATION)
	  {
	    word xh, xl;
	    xh = dh % 105;  xh = (xh << 25) | (dl >> 7);  xl = dl << 25;
	    udiv_qrnnd (dummy, i, xh, xl, 105 << 25);
	    i >>= 25;
	  }
	else
	  udiv_qrnnd (dummy, i, dh % 105, dl, 105);
	if (! tab105[i] && (dh != 0 || dl > 7))
	  continue;

	r = powmod (pprev, powinit, cnt, dh, dl);
	if (r != 0 && (dh != 0 || (BASE PLUSMINUS 1) % dl != 0))
	  return k;
#endif
      }
    return 0;
  }
}

#if ARITH == 64
#ifndef __alpha
word
__mpn_invert_normalized_limb (word d)
{
  word di, dummy;
  /* Special case for DIVISOR_LIMB == 100...000.  */
  if (d << 1 == 0)
    di = ~(word) 0;
  else
    udiv_qrnnd (di, dummy, -d, 0, d);
  return di;
}
#endif

int
powmod (unsigned int pprev, word x, int cnt, word d)
{
  int i;
  word s, sh, sl;
  word di;
  word dummy;
  int nc;

  count_leading_zeros (nc, d);
  di = __mpn_invert_normalized_limb (d << nc);

  s = x % d;

  for (i = cnt; i != 0; i--)
    {
      umul_ppmm (sh, sl, s, s << nc);
      udiv_qrnnd_preinv (dummy, s, sh, sl, d << nc, di);
      s >>= nc;

      if ((pprev & 1) != 0)
	{
	  s *= BASE;
	  if (BASE >= 256)
	    fprintf (stderr, "error cannot handle BASE >= 256\n");
	  if (BASE >= 128 && s >= 128*d)  s -= 128*d;
	  if (BASE >= 64 && s >= 64*d)    s -= 64*d;
	  if (BASE >= 32 && s >= 32*d)    s -= 32*d;
	  if (BASE >= 16 && s >= 16*d)    s -= 16*d;
	  if (BASE >= 8 && s >= 8*d)	  s -= 8*d;
	  if (BASE >= 4 && s >= 4*d)	  s -= 4*d;
	  if (s >= 2*d)			  s -= 2*d;
	  if (s >= d)			  s -= d;
	}
      pprev >>= 1;
    }

  if (REVERSE_REPUNIT)
    s = d - s;
  return s == 1;
}
#endif

#if ARITH == 32
int
powmod (unsigned int pprev, word x, int cnt, word dh, word dl)
{
  int i;
  word sh, sl;
  word rh, rl;

  sl = x;
  if (dh == 0)
    sl %= dl;
  sh = 0;

  for (i = cnt; i != 0; i--)
    {
      mulmod (sh, sl, sh, sl, dh, dl, &rh, &rl);
      sh = rh;
      sl = rl;

      if ((pprev & 1) != 0)
	{
	  mulmod (sh, sl, 0, BASE, dh, dl, &rh, &rl);
	  sh = rh;
	  sl = rl;
	}
      pprev >>= 1;
    }

  if (REVERSE_REPUNIT)
    sub_ddmmss (sh, sl, dh, dl, sh, sl);
  return sh == 0 && sl == 1;
}
#endif

int
isprime (unsigned int t)
{
  unsigned int q, r, d;

  if (t < 3 || (t & 1) == 0)
    return t == 2;

  for (d = 3, r = 1; r != 0; d += 2)
    {
      q = t / d;
      r = t - q * d;
      if (q < d)
	return 1;
    }
  return 0;
}

int
fermat_pseudoprime (unsigned int a, mpz_t t)
{
  mpz_t tmp, az, res;
  int ret;

  mpz_init (tmp);
  mpz_init (res);
  mpz_init_set_ui (az, (unsigned long int) a);
  mpz_sub_ui (tmp, t, 1L);
  mpz_powm (res, az, tmp, t);
  ret = mpz_cmp_ui (res, 1L);
  mpz_clear (az);
  mpz_clear (res);
  mpz_clear (tmp);
  return ret == 0;
}