/* * Author: Dr Peter Dickman * Created: 24th April 2001 * by modifying a less efficient version built on 5th February * * 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 #include #include #include #include #include #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.timing" #define RESULTS_FILENAME "primesquares.results" #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"); FILE * resultsfile = fopen(RESULTS_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 - flushcount) > FLUSH_LIMIT) { flushcount = count; fflush(resultsfile); if ((count - timecount) > TIME_LIMIT) { timecount = count; the_time = time(NULL); fprintf(timefile, "Generated %ld squares at %s", count, 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); }