/*
 * Author:  Dr Peter Dickman
 * Created: 21st May 2001
 * Created: 24th April 2001
 * by modifying the equivalent program to avoid outputting squares
 *
 * 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);

}

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

#define TIME_FILENAME    "primesquares.puretiming"

#define FLUSH_LIMIT    100
#define TIME_LIMIT    1000

#define SEARCHLENGTH (SIEVELENGTH/3)

#define MAX(a,b) ((a)>=(b) ? (a) : (b))

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).
     */

    long count = 0;
    long flushcount = 0;
    long timecount = 0;

    long p1,p2,p5,p8;
    long delta, p1lim;

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

    time_t the_time = time(NULL);

    FILE * timefile    = fopen(TIME_FILENAME, "w");

    /* Find the first 9 valid numbers and note them as start positions */
    /* actually only need the ninth value, but calculate them all anyway */
    while (p_sctr < 9)
    {
	while (!sieve[p_sptr]) ++p_sptr;
        p_start[p_sctr++] = p_sptr++; 
    }

    fprintf(timefile, "Program started at %s", ctime(&the_time));

    for (p5 = p_start[4]; p5 < SEARCHLENGTH; p5++)
    {
	/* Set the centre value for a given run */
        /* then slowly specify more limited ranges for each of the others */
	if (sieve[p5]) 
	{
	    for (p2 = (p5-3); p2 > 0; p2--)
	    {
		if (sieve[p2])
		{
		    p1lim = MAX(0,(1 + 2*p2 - p5));

		    for (p1 = (p2-1); p1 >= p1lim; p1--)
		    {
			if (sieve[p1])
			{
			    delta = p2-p1;
			    if (2*delta == p5-p2) continue;

			    if (!sieve[p5+delta]) continue;
			    if (!sieve[p5-delta]) continue;
			    if (!sieve[p2+delta]) continue;

			    p8 = 2*p5 - p1 - delta;
			    if (!sieve[p8]) continue;
			    if (!sieve[p8+delta]) continue;
			    if (!sieve[p8-delta]) continue;

/* fprintf(resultsfile, "%ld %ld %ld\n", p1, p2, p5); */
			    ++count;
			}
		    }
		}
	    }


		if ((count - timecount) > TIME_LIMIT)
		{
		    timecount = count;
		    the_time = time(NULL);
		    fprintf(timefile, 
			    "Generated %ld squares (last centre %ld) at %s",
			    count,
			    2*p5 + 3,
			    ctime(&the_time));
		    fflush(timefile);
		}
	}
    }
    printf("Found %ld prime magic squares in all\n", count);
}


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

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

   }
}

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

    SIEVEENTRYTYPE *sieve;

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

    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");

    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);
/*    printf("Squares are being written to %s", RESULTS_FILENAME); */
    printf("Timings are being written to %s", TIME_FILENAME);

    generate_magic_squares(sieve, SIEVELENGTH);
    printf("\n\n                  ...magic squares generated.\n\n");

    exit (EXIT_SUCCESS);
}