/* * 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 #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); } /*========================================================================*/ 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); }