/* brogen.c -- (C) Tapio Rajala, Jan 2011

   Calculates Legendre symbols to determine if the equation
                     k*n! + A = m^2
   can have integer solutions for n0 <= n <= n1.

   The Legendre symbols are of the form ( k*n!+A / p ),
   where the prime p > n1.

   Compile:
   gcc -Wall -O2 -o brogen brogen.c -lgmp

   Usage:
   brogen n0 n1 A

*/

#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <inttypes.h>
#include <gmp.h>

#define TESTS 40

static int usage(void)
{
  printf("\nUsage: brogen n0 n1 A k\n\n");
  printf(" Calculates Legendre symbols to determine if the equation\n");
  printf("                  k*n! + A = m^2\n");
  printf(" might have integer solutions for n0 <= n <= n1 < 2^64.\n\n");
  return 0;
}

static void mpz_set_uint64(mpz_t ROP, uint64_t a)
{
  if (sizeof(unsigned long) >= sizeof(uint64_t))
    mpz_set_ui(ROP,a);
  else
  {
    mpz_set_ui(ROP,a>>32);
    mpz_mul_2exp(ROP,ROP,32);
    mpz_add_ui(ROP,ROP,a);
  }
}

static int isprime(uint64_t p)
{
  uint64_t i, r;

  if (p <= 2 || p%2==0)
    return (p == 2);

  for (i = 3, r = sqrt(p); i <= r; i += 2)
    if (p % i == 0)
      return 0;

  return 1;
}

int main(int argc, char **argv)
{
  mpz_t prime[TESTS];
  mpz_t power[TESTS];
  mpz_t residue[TESTS];
  uint64_t n, n0, n1, temp, A, k;
  mpz_t test;
  mpz_init(test);

  if (argc < 5)
    return usage();

  n0 = strtoull(argv[1],NULL,0);
  n1 = strtoull(argv[2],NULL,0);
  A = strtoull(argv[3],NULL,0);
  k = strtoull(argv[4],NULL,0);

  // Print out what we are about to do..
  if (k == 1)
  {
   printf("\n Calculating Legendre symbols ( n!+%"PRIu64" / p ) to determine \n",A);  
   printf(" if the equation n! + %"PRIu64" = m^2 might have integer solutions\n for %"PRIu64" <= n <= %"PRIu64"\n\n",A,n0, n1); 
  } else {
   printf("\n Calculating Legendre symbols ( %"PRIu64"*n!+%"PRIu64" / p ) to determine \n",k,A);
   printf(" if the equation %"PRIu64"*n! + %"PRIu64" = m^2 might have integer solutions\n for %"PRIu64" <= n <= %"PRIu64"\n\n",k,A,n0, n1);
  }
  printf(" Using the following %i primes:\n",TESTS);

  // Find the needed primes.
  temp = 0;
  n = n1 + 1;
  while (temp < TESTS)
  {
   if (isprime(n))
   {
    printf(" %"PRIu64"",n); // Print out the prime.
    if (temp % 8 == 7) printf("\n");
    mpz_init(prime[temp]);
    mpz_set_uint64(prime[temp],n);
    mpz_init(power[temp]);
    mpz_set_uint64(power[temp],(n-1)/2);
    temp++;
   }
   n++;
  }

  printf("\n");

  // Initialize the residues.  
  for (temp = 0; temp < TESTS; temp++)
  {
    mpz_init(residue[temp]);
    mpz_set_uint64(residue[temp],k);
  }

  // The main loop.
  for (n = 2; n < n1; n++)
  {
   if (n % 100000 == 0) printf("Testing n = %"PRIu64"      \r",n);
   for (temp = 0; temp < TESTS; temp++) // Calculate n! (mod p[i])
   {
     mpz_mul_ui(residue[temp],residue[temp],n);
     mpz_mod(residue[temp],residue[temp],prime[temp]);
   }

   if (n+1 > n0) // Start calculating Legendre symbols if we are far enough
   {    
    temp = 0;
    while (temp < TESTS)
    {
      mpz_add_ui(test,residue[temp],A); // n! + A -> test
      mpz_mod(test,test,prime[temp]);
      if (mpz_cmp_ui(test,0) != 0) // Check if n!+A != 0 (mod p[i])
      {
       mpz_powm(test,test,power[temp],prime[temp]);
       if ((mpz_cmp_ui(test,1) != 0) && (mpz_cmp_ui(test,0) != 0)) // Check if n!+1^{(p[i]-1)/2} = -1 (mod p[i])
         temp = 1000; //OK. This should be more than TESTS...
      }
      temp++;
    }
    if (temp == TESTS)
     printf("Candidate n = %"PRIu64" passed the test!\n",n);
   }
  }
  printf("\nDone.\n\n");

  return 0;
}