/* Program to test that the smallest prime factor of 2^(8+48*k)+1 and
  2^(40+48*k)+1 are 257.  */

/* This program is put into the public domain.  Thomas König, 2017.  */

#include <gmp.h>
#include <stdio.h>

#define NMAX 100000000

/* Leave out two.  */

const unsigned int primes[] = {
  3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59,
  61, 67, 71, 73, 79, 83, 89, 97, 101, 103, 107, 109, 113, 127, 131,
  137, 139, 149, 151, 157, 163, 167, 173, 179, 181, 191, 193, 197,
  199, 211, 223, 227, 229, 233, 239, 241, 251, 257};

int main()
{
  mpz_t f, num, rem, two;
  unsigned int i, j, k;
  mpz_init(f);
  mpz_init(num);
  mpz_init(rem);
  mpz_init_set_ui (f, 1);
  mpz_init_set_ui (two, 2);
  int found_257;
  
  for (i=0; i<sizeof(primes)/sizeof(primes[0]); i++)
    mpz_mul_ui (f, f, primes[i]);

  for (k=0; k<NMAX; k+=48)
    {
      int found_257;

      found_257 = 0;
      j = 8 + k;
      mpz_powm_ui (num, two, j, f);
      mpz_add_ui (num, num, 1);
      mpz_mod (num, num, f);
      
      for (i=0; i<sizeof(primes)/sizeof(primes[0]); i++)
	if (mpz_divisible_ui_p (num, primes[i]))
	  {
	    found_257 = primes[i] == 257;
	    break;
	  }
      if (!found_257)
	printf("Conjecture 1 failed: %u %Zd\n", j, num);

      found_257 = 0;
      j = 40 + k;
      mpz_powm_ui (num, two, j, f);
      mpz_add_ui (num, num, 1);
      mpz_mod (num, num, f);
      
      for (i=0; i<sizeof(primes)/sizeof(primes[0]); i++)
	if (mpz_divisible_ui_p (num, primes[i]))
	  {
	    found_257 = primes[i] == 257;
	    break;
	  }
      if (!found_257)
	printf("Conjecture 2 failed: %u %Zd\n", j, num);
    }
  mpz_clear (f);
  mpz_clear (num);
  mpz_clear (rem);
  mpz_clear (two);
  return 0;
}