/////////////////////////////////////////////////////////////////////////////
// findprp - find a probable prime near a primorial.
//
// This command can be called as follows:
//
// findprp #  to find all probable primes above and near the primorial n#
// findprp #-1 to find all probable primes below and near the primorial n#
// findprp *# to find probable primes near the multiple of the primorial m*n#
//
// Primality checking is done using the Miller-Rabin algorithm.
//
// You can also start the search at a given offset by executing the command 
// witha positive or negative offset.  This allows you to start from a given
// point if the values have already been tested.
//
// For example: 'findprp 27509#-31253' will find all primes starting at that 
// and working down.
//
// This program requires the GMP library (http://www.gmplib.org) and can be 
// number compiled:
//
//	cc -o findprp findprp.c -lgmp
//
// Progress of the search will be sent to standard error.  When a probable
// prime is found, it will be send to standard output.  It will contain the
// number followed by a comma and the number of decimal digits in the number.
//
// Copyright, 2008, Paul Schmidt.  Can be used for research and personal use.
/////////////////////////////////////////////////////////////////////////////

#include 
#include 
#include 
#include 

void print_num(
FILE *f, 
unsigned long mult, 
unsigned long p, 
int sign, 
unsigned long c)
{
  if (mult != 1)
    fprintf(f, "%d*", mult);
  fprintf(f, "%d#", p);
  if (sign)
    fprintf(f, "+");
  else
    fprintf(f, "-");
  fprintf(f, "%d", c);
}

main (int argc, char *argv[])
{
  char *progname = argv[0];
  unsigned long p_ui;
  unsigned long pn;
  unsigned long pn1;
  unsigned long from = 1;
  unsigned long mult = 1;
  int above = 1;
  mpz_t primorial;
  mpz_init(primorial);
  mpz_t p;
  mpz_init(p);

  
  if (argc < 2)
  {
    fprintf(stderr, "Usage: %s [multiplier]*+/-[start]\n", progname);
    exit(1);
  }

  if (strchr(argv[1],'*'))
  {
    if (strchr(argv[1],'+'))
      sscanf(argv[1], "%U#*%U+%U", &mult, &pn, &from);
    else if (strchr(argv[1],'-'))
    {
      sscanf(argv[1], "%U#*%U-%U", &mult, &pn, &from);
      above = 0;
    }
    else
      sscanf(argv[1], "%U#*%U", &mult, &pn);
  }
  else
  {
    if (strchr(argv[1],'+'))
      sscanf(argv[1], "%U#+%U", &pn, &from);
    else if (strchr(argv[1],'-'))
    {
      sscanf(argv[1], "%U#-%U", &pn, &from);
      above = -1;
    }
    else
      sscanf(argv[1], "%U#", &pn);
  }

  mpz_set_ui(primorial, 2);
  for (mpz_set_ui(p, 3);  mpz_cmp_ui(p, pn) <= 0;  mpz_nextprime(p, p))
  {
    mpz_mul(primorial, primorial, p);
    p_ui = mpz_get_ui(p);
  }

  pn = p_ui;
  pn1 = mpz_get_ui(p);

  mpz_t number;
  mpz_init2(number, mpz_sizeinbase(primorial,2) + 32);
  mpz_mul_ui(number, primorial, mult);

  mpz_t tmp;
  mpz_init2(tmp, mpz_sizeinbase(number,2) + 1);

  mpz_t to;
  mpz_init_set_ui(to, pn1);
  mpz_mul_ui(to, to, pn1); 

  for (mpz_set_ui(p, from); mpz_cmp(p, to) < 0; mpz_nextprime(p, p))
  {
    fprintf(stderr, "\nTesting "); 
    print_num(stderr, mult, pn, above, mpz_get_ui(p)); 
    if (above)
      mpz_add(tmp, number, p);
    else
      mpz_sub(tmp, number, p);
    if (mpz_millerrabin(tmp, 5))
    {
      print_num(stdout, mult, pn, above, mpz_get_ui(p));
      fprintf(stdout, ",%U\n", mpz_sizeinbase(tmp, 10));
      fflush(stdout);
      fprintf(stderr, " is a prp!");
    }
    if (mpz_cmp_ui(p, 1) == 0)
      mpz_set_ui(p, pn1);
  }

  mpz_clear(primorial);
  mpz_clear(number);
  mpz_clear(p);
  mpz_clear(tmp);
}