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

   Tests if the squares p^2 ,with primes p (p0 <= p <= p1),
   divide any integer of the form 4n! + 1.

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

   Usage:
   sqrdiv p0 p1
                        
*/

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


static int usage(void)
{
  printf("\nUsage: sqrdiv p0 p1\n\n");
  printf(" Tests if the squares p^2 ,with primes p (p0 <= p <= p1),\n");
  printf(" divide any integer of the form 4n! + 1.\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 modulus;
  mpz_init(modulus);
  uint64_t p, p0, p1, temp;
  mpz_t n;
  mpz_init(n);
  mpz_t test;
  mpz_init(test);

  if (argc < 3)
    return usage();

  p0 = strtoull(argv[1],NULL,0);
  p1 = strtoull(argv[2],NULL,0);

  // Print out what we are about to do..
  printf("\n Tests if  p^2 | 4n! + 1  for some integer n,\n");
  printf(" where %"PRIu64" <= p <= %"PRIu64", p is prime.\n\n", p0, p1);

  // Teh main loop.
  for (p = p0; p+1 < p1; p++)
  {
   if (isprime(p))
   {
    printf("Testing p = %"PRIu64"      \r",p); // Print out the prime we are testing.
    mpz_set_uint64(n,4);   // 4 -> n 
    mpz_set_uint64(modulus,p);
    mpz_mul_ui(modulus,modulus,p);
    for (temp = 2; temp < p; temp++)
    {
     mpz_mul_ui(n,n,temp);  // 4*temp! -> n
     mpz_mod(n,n,modulus); // (mod p*p)

     mpz_add_ui(test,n,1); // 4*temp! + 1 -> test
     mpz_mod(test,test,modulus); // (mod p*p)

     if (mpz_cmp_ui(test,0) == 0) printf(" %"PRIu64"^2 | 4*%"PRIu64"!+1\n",p,temp);
    }
   }
  }

  printf("                              \n Done.\n\n");

  return 0;
}