/* * Author: Dr Peter Dickman * Created: 20 December 2004 * by adaptation of an earlier program (see below). * adapted to allow for a single 1 in the square, and to * allow for repetition of prime values, in line with a * suggestion from: * Stan Schmidt, schmidtw@newpaltz.edu, USA * * Software this was built from was: * Author: Dr Peter Dickman * Created: 24th April 2001 * by modifying a less efficient version built on 5th February 2001 * * The original program was 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. * This version generates magic squares of the form * * a 1 c * d e f * g h i * * where the other eight values are (possibly non-distinct) primes. * We consider the nine positions of the square to be numbered for * indexing and commentary purposes: * * 1 2 3 * 4 5 6 * 7 8 9 * * (c) 2001,2004 University of Glasgow * All Rights Reserved * */ #include #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)) #define FILENAME_LENGTH 255 /* 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 #define DEFAULT_ROOTSIEVELEN 1000 /* * 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 */ /* * Once populated, the sieve contains TRUE at index i iff 2i+3 is prime */ /*========================================================================*/ 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 contains 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 /* The largest value we can get in a square of the structure we are looking for * is twice the centre value - 1. Since the largest value represented in our * sieve is 2*(SIEVELENGTH-1) + 3, this means that our centre value should * be limited to SIEVELENGTH+1 and this is represented by an index value * of SIEVELENGTH/2 - 1 */ #define SEARCHLENGTH ((SIEVELENGTH/2)-1) #define MAX(a,b) ((a)>=(b) ? (a) : (b)) void generate_magic_squares(SIEVEENTRYTYPE sieve[], unsigned long SIEVELENGTH, char * outputfile, char * timingfile, BOOLEAN * verbose, BOOLEAN * veryverbose) { long count = 0; long flushcount = 0; long timecount = 0; long p5; long delta, deltainc, tmp; time_t the_time = time(NULL); FILE * timefile; FILE * resultsfile; if (0 == strcmp(timingfile,"")) { timefile = fopen(TIME_FILENAME, "w"); printf("Timings are being written to %s\n", TIME_FILENAME); } else { timefile = fopen(timingfile, "w"); printf("Timings are being written to %s\n", timingfile); } if (0 == strcmp(outputfile,"")) { resultsfile = fopen(RESULTS_FILENAME, "w"); printf("Squares are being written to %s\n", RESULTS_FILENAME); } else { resultsfile = fopen(outputfile, "w"); printf("Squares are being written to %s\n", outputfile); } printf("\n\n"); fflush(stdout); fprintf(timefile, "Program started at %s", ctime(&the_time)); /* The key value for controlling progress is the centre value in the * square. The row/column/diagonal sum for the square will be three * times the centre value, so we work through all possible prime * centre values exploring as we go. * Incidentally, if we start with this lowest possible prime * in the centre the following square is almost acceptable: * * 5 1 3 * 1 3 5 * 3 5 1 * * Indeed, if the central value is a prime Q and 2Q-1 is prime * then the square 2Q-1 1 Q * 1 Q 2Q-1 * Q 2Q-1 1 * is almost a solution. But the repetition of the 1's isn't allowed. * This corresponds to the case epsilon=0 in the discussion below, * so we won't consider it further except to note that epsilon must * be non-zero. * * We could start with a carefully chosen initial central value, * hand-eliminating the early options, but the code is so fast it * is simpler to just start with the samllest odd prime in the * centre (i.e. 3, which is represented by an index of 0) and * let the code work its way up from there. * */ for (p5 = 0; p5 < SEARCHLENGTH; p5++) { /* Set the centre value for a given run */ /* then slowly specify more limited ranges for each of the others */ /* As we progress we ensure that we have a prime centre value and */ /* that the value opposite the 1 is also prime (detailed comment */ /* immediately after the test explains this if statement). */ /* Only in that situation do we do any more work for this centre. */ if (sieve[p5] && sieve[2*p5+1]) { /* we know we want a 1 in a mid-side position, so we force p2 */ /* and that immediately determines p8, so we calculate it */ /* and check whether it is prime.... */ /* Because the sieve contains TRUE in slot n to represent 2n+3 */ /* being prime, p5 represents the prime value 2*p5 + 3 and the */ /* square therefore has a row-sum of 6*p5 + 9. With 1 fixed in */ /* the middle of a side, we must therefore require that we */ /* have 4*p5+5 opposite it, which means that we want the slot */ /* 2*p5+1 in our sieve representation to contain TRUE (i.e. we */ /* want 2*(2*p5+1) + 3 = 4*p5 + 5 to be prime). */ /* Now we can simply explore the possible values for the */ /* remaining free parameter, which can be any corner. */ /* In the original version of this program I was looking */ /* for a triple of triples, with all nine indicated values */ /* being prime. This program works vaguely similarly, but */ /* we fix the smallest value (1), choose a central value, */ /* and calculate a max value (2*central - 1); we then need */ /* to find an integer epsilon such that: * 1 + epsilon and 1 + 2*epsilon are prime, * and centre - epsilon and centre + epsilon are prime * and max - epsilon and max - 2*epsilon are prime. * * Obviously, as discussed above, epsilon = 0 always gives * us an unacceptable magic square at this stage, as we've * already confirmed that the largest value in the square * is also prime but we'll end up with multiple 1's. * * The square we then build ends up as something like: * * centre+epsilon 1 2*centre-epsilon-1 * 2*centre-2*epsilon-1 centre 2*espilon+1 * 1+espilon 2*centre-1 centre-epsilon * * Now we can choose values of epsilon that are rather large, * so that 1 + 2*epsilon is bigger than centre, the limit is * that epsilon must fall in the range 0 < epsilon < centre-1 * In fact, because we can only use odd primes (2 isn't useful) * we can allow epsilon to progress in steps of 2, which is * the same as stepping through the sieve in steps of 1. * * In the code we use the variable delta, to represent an * index in the sieve corresponding to the value 2*delta+3 * We need to know if 1+epsilon and 1+2epsilon etc are all prime. * We allow delta to run through all possible indices from 0 * (corresponding to the prime 3, which is indicated by * 1+epsilon with epsilon=2) all the way to p5-1 * (corresponding to the possible prime 2*p5+1, which is the * largest allowed value of 1+epsilon as the next step causes * repetition of the centre value (and more importantly repetition * of 1 which isn't allowed) and larger steps would force us to * have some numbers below 1 in the square, which also isn't * allowed). */ for (delta = 2; delta < p5; delta++) { /* We use the following below, to optimise the sums */ deltainc = delta+1; /* check that 1+epsilon is OK, if not give up */ if (!sieve[delta]) continue; /* check that 1+2*epsilon is OK, if not give up */ /* index delta means 2delta+3, so epsilon is 2delta+2 * so 2*epsilon + 1 is 4delta+5, which matches an index * of 2*delta+1 (as we multiply the index by two then add three) */ if (!sieve[delta+deltainc]) continue; /* check that centre-epsilon is OK, if not give up */ /* centre is 2p5+3, epilson is 2delta+2, so centre-epsilon * is 2p5+1-2delta, which gives an index of p5-delta-1 */ if (!sieve[p5-deltainc]) continue; /* check that centre+epsilon is OK, if not give up */ /* centre is 2p5+3, epilson is 2delta+2, so centre+epsilon * is 2p5+5+2delta, which gives an index of p5+delta+1 */ if (!sieve[p5+deltainc]) continue; /* check that max-epsilon is OK, if not give up */ /* max is 2*centre-1, centre is 2p5+3, so max is 4p5+5, * epsilon is 2delta+2, so max-epsilon is 4p5+3-2delta * which matches an index in the sieve of 2p5-delta */ tmp = 2*p5 - delta; if (!sieve[tmp]) continue; /* check that max-2*epsilon is OK, if not give up */ /* from immediately above we see that max-2*epsilon * is 4p5+1-4delta, so index is 2p5-2delta-1 */ if (!sieve[tmp-deltainc]) continue; fprintf(resultsfile, "%ld %ld %ld\n", (long) 1, 2*delta+3, 2*p5+3); ++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); if (*verbose) { fprintf(stdout, "Generated %ld squares at %s\n", count, ctime(&the_time)); fflush(stdout); } } else if (*veryverbose) { fprintf(stdout, "Generated %ld squares\n", flushcount); fflush(stdout); } } } } printf("Found %ld prime magic squares with a 1, in all\n", count); fflush(stdout); } /*========================================================================*/ void analyse_args(int argc, char* argv[], unsigned long * ROOTSIEVELEN, char * outputfile, char * timingfile, BOOLEAN * verbose, BOOLEAN * veryverbose) { int i = 1; while (i < argc) { if ((0 == strcmp(argv[i], "-h")) || (0 == strcmp(argv[i], "-help"))) { printf("Usage: %s -h prints this message and exits\n\n", argv[0]); if (i > 1) printf(" NB: defaults given below are affected by preceding options\n"); printf("Other options are:\n"); printf(" -help also prints this message and exits\n"); printf(" -v produces a commentary on stdout as it works\n"); printf(" -vv produces a long commentary on stdout as it works\n"); printf(" -nv does not produce a commentary on stdout as it works\n"); printf(" the default is %s produce a %s commentary\n", *verbose ? "to" : "not to", *veryverbose ? "long" : "(brief)"); printf(" -s int adjusts the size of the underlying sieve of primes\n"); printf(" where sieve size is 2*int*(int+1) and this allows for\n"); printf(" prime numbers up to 4*int*(int+1) + 1\n"); printf(" the default is %ld\n", *ROOTSIEVELEN); printf(" -f file redirects the output into the indicated file\n"); printf(" the default is %s\n", RESULTS_FILENAME); printf(" -t file redirects the timing info into the indicated file\n"); printf(" the default is %s\n", TIME_FILENAME); printf("Filenames are at most %d chars and cannot includes spaces\n", FILENAME_LENGTH); printf("NB: Each option must be given separately.\n\n"); fflush(stdout); exit(EXIT_SUCCESS); } else if ((0 == strcmp(argv[i], "-s")) && (argc > i+1)) { if (1 != sscanf(argv[i+1], "%lu", ROOTSIEVELEN)) { fprintf(stderr, "Illegal argument to -s option: %s\n", argv[i+1]); fprintf(stderr, "Usage: %s -h will give usage info\n", argv[0]); fflush(stderr); exit (EXIT_FAILURE); } else if ((*ROOTSIEVELEN < 10) || (*ROOTSIEVELEN > 1000000)) { fprintf(stderr, "Nonsensical parameter to -s option: %s\n", argv[i+1]); fprintf(stderr, "Try a non-trivial, not too large, positive integer\n"); fprintf(stderr, "Usage: %s -h will give usage info\n", argv[0]); fflush(stderr); exit (EXIT_FAILURE); } i += 2; } else if ((0 == strcmp(argv[i], "-v"))) { *verbose = BOOL_TRUE; *veryverbose = BOOL_FALSE; i++; } else if ((0 == strcmp(argv[i], "-nv"))) { *verbose = BOOL_FALSE; *veryverbose = BOOL_FALSE; i++; } else if ((0 == strcmp(argv[i], "-vv"))) { *verbose = BOOL_TRUE; *veryverbose = BOOL_TRUE; i++; } else if ((0 == strcmp(argv[i], "-f")) && (argc > i+1)) { sscanf(argv[i+1], "%255[^ ]", outputfile); i += 2; } else if ((0 == strcmp(argv[i], "-t")) && (argc > i+1)) { sscanf(argv[i+1], "%255[^ ]", timingfile); i += 2; } else { fprintf(stderr, "Usage: %s -h will give usage info\n", argv[0]); fflush(stderr); exit (EXIT_FAILURE); } } } int main (int argc, char* argv[]) { unsigned long ROOTSIEVELEN = DEFAULT_ROOTSIEVELEN; unsigned long SIEVELENGTH; SIEVEENTRYTYPE *sieve; char outputfile[FILENAME_LENGTH+1] = ""; char timingfile[FILENAME_LENGTH+1] = ""; BOOLEAN verbose = BOOL_TRUE; BOOLEAN veryverbose = BOOL_FALSE; if (argc > 1) analyse_args(argc, argv, &ROOTSIEVELEN, outputfile, timingfile, &verbose, &veryverbose); 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 with a 1 Explorer\n"); printf("(c) 2001, 2004 Peter Dickman & University of Glasgow\n\n"); printf("Finds magic squares containing distinct primes and a 1\n"); printf("following a question from: Stan Schmidt of newpalzt.edu\n"); printf("\n"); printf("Generating magic squares (using primes up to %lu)\n\n\n", 2*(SIEVELENGTH-1) + 3); fflush(stdout); generate_magic_squares(sieve, SIEVELENGTH, outputfile, timingfile, &verbose, &veryverbose); printf("\n\n ...magic squares generated.\n\n"); fflush(stdout); exit (EXIT_SUCCESS); }