/*
 * Author:  Dr Peter Dickman
 * Created: 1st February 2001
 * by adding to a prime number generator that was
 * originally created: 20th March & 8th September 2000
 *
 * This program is designed for exploring prime magic squares,
 * i.e. 3x3 magic squares that contain 9 distinct prime numbers.
 * It includes a high-speed prime generator based
 * on the Sieve of Eratosthenes.
 *
 * (c) 2001  University of Glasgow
 * All Rights Reserved
 *
 */

#include <stdio.h>
#include <stdlib.h>
#include <limits.h>
#include <float.h>
#include <math.h>
#include <time.h>


#define BOOL_FALSE 0
#define BOOL_TRUE  1
typedef char BOOLEAN;

#define MIN(a,b) ((a)<(b) ? (a) : (b))

/* This code is built around a global array which represents the primes. */

/* various base types from which we build the array etc*/
#define SIEVEENTRYTYPE BOOLEAN
#define SIEVEINDEXTYPE unsigned long
#define PRIMECOUNTTYPE SIEVEINDEXTYPE

/* and specific values used in the construction */

#define ISNOTPRIME BOOL_FALSE
#define ISPRIME    BOOL_TRUE
#define MAYBEPRIME ISPRIME

/* 
 * To be confident the sieve is really tightly coded it's helpful
 * if we can be sure that:
 * (2*(ROOTSIEVELEN-1)+3)^2 = (2*(SIEVELENGTH-1)+3)
 * i.e. 4*(r-1)*(r-1) + 12*(r-1) + 9 = 2*s + 1
 * i.e. 2*(r*r - 2*r + 1) + 6*(r-1) + 4 = s
 * i.e. 2*r*r - 4*r + 2 + 6*r - 6 + 4 = s
 * i.e. 2*r*(r+1) = s
 */

/*========================================================================*/

void sieve_classical_alg(SIEVEENTRYTYPE sieve[], 
			 unsigned long SIEVELENGTH, 
			 unsigned long ROOTSIEVELEN)
{
    SIEVEINDEXTYPE ptr, step, current;

    PRIMECOUNTTYPE num_primes = 0;

    for (ptr = 0; ptr < SIEVELENGTH; sieve[ptr++] = MAYBEPRIME);

    for (ptr = 0; ptr < ROOTSIEVELEN; ptr++)
    {
	/* have we found another prime...?  */
	if (sieve[ptr])
	{
	    /* if so, which is it, that's how far apart our */
	    /* eliminations will be since we'll dispose of  */
	    /* every number that is 2*step apart, as this   */
	    /* sieve only conmtains odd numbers.....        */
	    step    = 2*ptr + 3;

	    /* the first value we eliminate is p*p, as smaller  */
	    /* multiples will have already been stomped on. Now */
	    /* p == step, so we want the index corresponding to */
	    /* step*step, i.e. we know that...                  */
	    /*      2*current + 3 = 4*ptr*ptr + 12*ptr + 9      */
	    /* so     current = 2*ptr*ptr + 6*ptr + 3           */
	    /* i.e.   current = step*ptr + step + ptr           */

	    current = step*ptr + step + ptr;

	    while (current < SIEVELENGTH)
	    {
		sieve[current] = ISNOTPRIME;
		current += step;
	    }

	    num_primes++;
	}
    }

    printf("Discovered %lu primes less than %lu\n", 
	   num_primes, 2*ROOTSIEVELEN + 3);

    for (ptr = ROOTSIEVELEN; ptr < SIEVELENGTH; ptr++)
	if (sieve[ptr]) num_primes++;

    printf("Found %lu primes less than %lu in all\n", 
	   num_primes, 2*SIEVELENGTH + 3);

}

/*========================================================================*/

void sieve_print1(SIEVEENTRYTYPE sieve[], unsigned long SIEVELENGTH)
{
    SIEVEINDEXTYPE ptr;
    for (ptr = 0; ptr < SIEVELENGTH; ptr++)
    {
	if (sieve[ptr]) printf("%lu\n", 2*ptr + 3);
    }
}

void sieve_print(SIEVEENTRYTYPE sieve[], unsigned long SIEVELENGTH) 
{
    sieve_print1(sieve, SIEVELENGTH);
}

/*========================================================================*/

static clock_t clockstart;

void print_time(void)
{
    time_t the_time;

    the_time = time(NULL);
    printf("The time is now: %s\n\n", ctime(&the_time));

#if 0
    double elapsed;
    clock_t clockend = clock();

    elapsed = (1e6*(clockend - clockstart))/CLOCKS_PER_SEC;
    printf("Elapsed time      = %12.3f secs\n\n", elapsed/1e6);
#endif
}

void init_time(void)
{
    clockstart = clock();
    print_time();
}


/*========================================================================*/

#define PVAL(x)          (2*x + 3)
#define SUMPRIMES(x,y,z) (PVAL(x) + PVAL(y) + PVAL(z))

void printAsquare(SIEVEENTRYTYPE sieve[], long target,
		  long p1, long p2, long p3, long p4, 
		  long p5, long p6, long p7, long p8, long p9) 
{
    printf("These primes form a magic square (type A):\n");
    printf("%ld %ld %ld %ld %ld %ld %ld %ld %ld\n",
	   PVAL(p1), PVAL(p2), PVAL(p3), PVAL(p4), PVAL(p5), 
           PVAL(p6), PVAL(p7), PVAL(p8), PVAL(p9));
    printf("The square's value is %ld\n", target);

    printf("\n  %6ld  %6ld  %6ld\n  %6ld  %6ld  %6ld\n  %6ld  %6ld  %6ld\n\n",
           PVAL(p6), PVAL(p7), PVAL(p2), PVAL(p1), PVAL(p5), 
           PVAL(p9), PVAL(p8), PVAL(p3), PVAL(p4));
}


void printBsquare(SIEVEENTRYTYPE sieve[], long target,
		  long p1, long p2, long p3, long p4, 
		  long p5, long p6, long p7, long p8, long p9) 
{
    printf("These primes form a magic square (type B):\n");
    printf("%ld %ld %ld %ld %ld %ld %ld %ld %ld\n",
	   PVAL(p1), PVAL(p2), PVAL(p3), PVAL(p4), PVAL(p5), 
           PVAL(p6), PVAL(p7), PVAL(p8), PVAL(p9));
    printf("The square's value is %ld\n", target);

    printf("\n  %6ld  %6ld  %6ld\n  %6ld  %6ld  %6ld\n  %6ld  %6ld  %6ld\n\n",
           PVAL(p7), PVAL(p6), PVAL(p2), PVAL(p1), PVAL(p5), 
           PVAL(p9), PVAL(p8), PVAL(p4), PVAL(p3));
}


int validate_potential_square(SIEVEENTRYTYPE sieve[], 
                              long p1, long p2, long p3, long p4, 
                              long p5, long p6, long p7, long p8, long p9) 
{
/* Given indexes for nine distinct primes, in ascending sequence, 
 * determine whether they form a magic square....
 * Any such collection of numbers will have to be laid out as follows
 * if it is to form a magic square (modulo symmetries of the square):
 *
 *      6 7 2       7 6 2
 *      1 5 9       1 5 9
 *      8 3 4       8 4 3
 *
 * So the following triples must have the same sum:
 * p1+p5+p9  p6+p5+p4  p8+p5+p2  p7+p5+p3  p2+p6+p7  p3+p4+p8
 * as well as one or other of the following pairs of triples:
 * p2+p3+p9  p1+p7+p8        or        p2+p4+p9  pp1+p6+p8
 *
 * returns 0 if no square is formed
 * returns 1 if a side-lowest square is formed
 * returns 2 if a corner-lowest square if formed
 */

  long target = SUMPRIMES(p1,p5,p9);

  if (target != SUMPRIMES(p4,p5,p6)) return 0;
  if (target != SUMPRIMES(p2,p5,p8)) return 0;
  if (target != SUMPRIMES(p3,p5,p7)) return 0;
  if (target != SUMPRIMES(p2,p6,p7)) return 0;
  if (target != SUMPRIMES(p3,p4,p8)) return 0;
  
  /* So far so good, now test for final triples */

  if ((target == SUMPRIMES(p2,p3,p9)) && (target == SUMPRIMES(p1,p7,p8)))
  {
      printBsquare(sieve,target,p1,p2,p3,p4,p5,p6,p7,p8,p9);
      return 1;
  }
  else if ((target == SUMPRIMES(p2,p4,p9)) && (target == SUMPRIMES(p1,p6,p8)))
  {
      printAsquare(sieve,target,p1,p2,p3,p4,p5,p6,p7,p8,p9);
      return 2;
  }
  else return 0;
}


void generate_magic_squares(SIEVEENTRYTYPE sieve[], unsigned long SIEVELENGTH) 
{
    /* The first nine odd primes are:      3,5,7,11,13,17,19,23,29
     * so their indexes in the sieve are:  0,1,2, 4, 5, 7, 8,10,13
     * These are therefore the start values for the loop counters
     * that run up through the primes finding all the combinations
     * of nine distinct primes that are then offered for validation.
     *
     * This could be done programmatically, based on the sieve, and
     * should be for generality (to allow the magic squares to be
     * generated from any given collection of integers).
     */

    int count = 0;
    long p0,p1,p2,p3,p4,p5,p6,p7,p8;

    long p_start[9];
    long p_sctr = 0;
    long p_sptr = 0;

    /* Find the first 9 valid numbers and note them as start positions */
    while (p_sctr < 9)
    {
	while (!sieve[p_sptr]) ++p_sptr;
	p_start[p_sctr++] = p_sptr++;
    }

    init_time();

    for (p8 = p_start[8]; p8 < SIEVELENGTH; p8++)
    {
	/* Set the maximum allowable prime index for a given run */
        /* then slowly specify more limited ranges for each of the others */
	if (sieve[p8]) for (p7 =  p_start[7]; p7 < p8; p7++)
	if (sieve[p7]) for (p6 =  p_start[6]; p6 < p7; p6++)
	if (sieve[p6]) for (p5 =  p_start[5]; p5 < p6; p5++)
	if (sieve[p5]) for (p4 =  p_start[4]; p4 < p5; p4++)
	if (sieve[p4]) for (p3 =  p_start[3]; p3 < p4; p3++)
	if (sieve[p3]) for (p2 =  p_start[2]; p2 < p3; p2++)
	if (sieve[p2]) for (p1 =  p_start[1]; p1 < p2; p1++)
	if (sieve[p1]) for (p0 =  p_start[0]; p0 < p1; p0++)
        if (sieve[p0]) 
	    if (validate_potential_square(sieve,p0,p1,p2,p3,p4,p5,p6,p7,p8))
	    { count++; print_time(); fflush(stdout); }
    }

    printf("Found %d prime magic squares in all\n", count);

}


/*========================================================================*/

void analyse_args(int argc, char* argv[], 
		  unsigned long * ROOTSIEVELEN, BOOLEAN * doprint)
{
   int i = 1;
   while (i < argc)
   {
       if (0 == strcmp(argv[i], "-p"))
       {
	   *doprint = BOOL_TRUE;
	   i++;
       }
       else if ((0 == strcmp(argv[i], "-s")) && (argc > i+1))
       {
	   sscanf(argv[i+1], "%lu", ROOTSIEVELEN);
	   i += 2;
       }
       else
       {
	   fprintf(stderr, "Usage: %s [-p] [-s size]\n", argv[0]);
	   exit (EXIT_FAILURE);
       }

   }
}

int main (int argc, char* argv[])
{
    unsigned long ROOTSIEVELEN = 1000;
    unsigned long SIEVELENGTH;
    BOOLEAN doprint = BOOL_FALSE;

    SIEVEENTRYTYPE *sieve;

    if (argc > 1) analyse_args(argc, argv, &ROOTSIEVELEN, &doprint);

    SIEVELENGTH = 2*ROOTSIEVELEN*(ROOTSIEVELEN+1);

    sieve = (SIEVEENTRYTYPE *) malloc(SIEVELENGTH*sizeof(SIEVEENTRYTYPE));

    /* calculate the primes */

    printf("\n\nThe Prime Exploratorium\n");
    printf("(c) 2000  Peter Dickman & University of Glasgow\n\n");
    printf("Working with a sieve of length %lu\n", SIEVELENGTH);
    printf("So maximum potential prime is %lu\n", 2*(SIEVELENGTH-1) + 3);
    printf("Generating primes... (up to %lu)\n", 2*(SIEVELENGTH-1) + 3);
    sieve_classical_alg(sieve, SIEVELENGTH, ROOTSIEVELEN);
    printf("                 ...primes generated.\n\n");


    if (doprint) sieve_print(sieve, SIEVELENGTH);   


    printf("\n\nThe Prime Magic Square Explorer\n");
    printf("(c) 2001 Peter Dickman & University of Glasgow\n\n");
    printf("Finds magic squares containing distinct primes");
    printf("\n");

    printf("Generating magic squares (using primes up to %lu)\n\n\n", 
	   2*(SIEVELENGTH-1) + 3);
    generate_magic_squares(sieve, SIEVELENGTH);
    printf("\n\n                  ...magic squares generated.\n\n");

    exit (EXIT_SUCCESS);
}