/* Fred, Gadzooks! I've found another bug in the primality testing, so I'm back with another entry. (The bug would have only affected the outcome of the test for a few hundred numbers, but I wouldn't have wanted to see one of those numbers on my screen when I checked my mail tomorrow.) -- Jim */ /* POTM ENTRY: wookin_po_nubbas.c */ /* Jim Roche */ /* roche@ccr-p.ida.org */ /* */ /* cc -O wookin_po_nubbas.c -o wookin_po_nubbas */ /* Below are notes to myself (Jim) */ /** wookin8.c modified Sun., Jan. 15, 2:00 pm, to correct another * bug in the primality test (change "r-=n" to "r %= n" at * appropriate place). **/ /*** NOTE: Remember to comment out the definition of DEBUG for * the final version!!!! ***/ /** RECALL: The fancy division routine assumes that we divide by _odd_ numbers only; do NOT divide by 2 !!! **/ /** NOTE: I could improve the trial divison by unrolling loop * and making the code in-line, but I don't think that this would * have a huge effect on overall running time (probably less * than 10%, I imagine). **/ #include /*** NOTE: Remember to comment out (or remove) the definition of * DEBUG in final version!!! ***/ /* #define DEBUG 1 */ #define min(A,B) ((A) < (B) ? (A) : (B)) #define max(A,B) ((A) > (B) ? (A) : (B)) #define TESTBIT 0x1000000; #define ONES 0x1ffffff; /*** Most of the bit maps used in this program will be signaled * by the fact that their names begin with the letter "b." They * will generally be 25-bit maps right justified inside 32-bit * words. TESTBIT represents the leftmost bit in these maps, * corresponding to the upper left position in the grid. * Grid locations are generally represented as follows: * * 0 1 2 3 4 * 5 6 7 8 9 * 10 11 12 13 14 * 15 16 17 18 19 * 20 21 22 23 24 * * Occasionally it will also be useful to think of letting the * row and column indices each vary from 0 to 4, with the column * index varying more rapidly. ***/ /*** HIGHBIT will be used when computing r*r mod n as part of * a primality test adapted from Miller's test (which is itself * based on Fermat's Litle Theorem). ***/ #define HIGHBIT 0x80000000 /** BNBRS(bits) takes a bitmap representation of an arbitrary * subset of the 25 grid positions and returns a bitmap * representation of all grid positions that are adjacent to * at least one of the positions in the first subset. * Note that a grid location is NOT considered to be * adjacent to itself. (However, the two subsets can clearly * still share elements if the first subset has at least two * elements next to each other in the grid.) * The vector brow0[32] is an array of bitmaps. * Writing the index of brow0[] as a 5-bit number (most * significant bit first), we use the index as a bitmap * to select a subset of the 5 grid locations in Row 0. E.g., * the subscript 20 ( = 10100 in binary) represents the first * and third locations in the row (Locations 0 and 2, since we * begin indexing at 0). Then brow[20] (for example) is a bit * map representing the set of all grid locations (not just those * in Row 0) that are adjacent to either position (0,0 or (0,2) * (or to both). * The bit maps brow0[0], ... , brow0[31] are defined as above. * the arrays brow1[], ..., brow4[] are defined analogously. **/ #define BNBRS(bits) brow0[(bits) >> 20] | \ brow1[((bits) >> 15) & 0x1f] | \ brow2[((bits) >> 10) & 0x1f] | \ brow3[((bits) >> 5) & 0x1f] | \ brow4[(bits) & 0x1f] /*** All of the quantities with "SIZE" in their names are used * for sieving and creating a table with a modest number of primes. * When testing a number for primality, we'll first do trial * division by these tabulated primes, then jump to the Miller * test (if necessary). **/ #define TABSIZE 45 /* total number of primes needed*/ #define BOOTSIZE 7 /* number of primes to bootstrap*/ #define SIEVESIZE 220 /* size of sieve */ /** NOTE: Following increase in array size fixes bug in wookin6.c */ int P[TABSIZE + 10]; /* array of primes */ /* int S[SIEVESIZE]; */ int S[2*SIEVESIZE]; /* Changed in jim11.c **/ /**** The array of 64-bit floating-point numbers Pinv[] contains * the approximate reciprocals of the primes in P[]. These * are used in a turbocharged divisibility test that's * significantly faster than using the integer operator "%". The array dx[] is just a kludge (almost certainly * useless, too) introduced in an attempt to make sure that a * an instruction in the trial-division routine would be * executed correctly. ***/ double Pinv[TABSIZE]; double dx[50]; /* Used as trick to make fancy trial division work; * in reality, it's unnecessary, though, I think. */ int bpresent[10],bgeneral[11], b17, b39, b1379; int npart[11], digit[11]; int mindigit, maxdigit; int mod3[90]; /* Need >= 84 for buildmod3(). */ int digitsum[10]; /* Following added in jim10.c */ int bdist[11]; /* Array of bit masks. bdist[j] represents those grid positions that can reach a legal final digit (1,3,7, or 9) in EXACTLY j steps. */ /** Following variables added in jim12.c and jim13.c **/ int b0369, bnot0369, b147, bnot147; int is0369[10] = { 1,0,0,1,0,0,1,0,0,1 }; int is147[10] = { 0,1,0,0,1,0,0,1,0,0 }; int specmask; /* Special mask to do lookahead mod 3 */ int specialdist[10], bdist0369[10], bdist147[10]; int *ptr, *ptr1, *ptr2; int D[25]; /* D[i] is the digit in cell i */ /* bnbrs[i] is (octal) bit map of neighbors of cell i **/ int bnbrs[25] = { 043000000, 0123400000, 051600000, 024700000, 010300000, 0141060000, 0162470000, 071234000, 034516000, 014206000, 03021400, 03451600, 01624700, 0712340, 0304140, 060430, 071234, 034516, 016247, 06103, 01410, 01624, 0712, 0345, 0142 }; int brow1[32] = { 00,014206000,034516000,034716000, 071234000,075236000,075736000,075736000, 0162470000,0176676000,0176576000,0176776000, 0173674000,0177676000,0177776000,0177776000, 0141060000,0155266000,0175576000,0175776000, 0171274000,0175276000,0175776000,0175776000, 0163470000,0177676000,0177576000,0177776000, 0173674000,0177676000,0177776000,0177776000 }; int brow0[32], brow2[32], brow3[32],brow4[32]; /*** In bitinit(), we create a number of useful bit masks from * the digits that we read in. * For 0 <= i <= 9, the bit map bpresent[i] shows which of the * 25 grid locations contain the digit i. * The bit map b1379 indicates the set of grid locations which * contain 1,3,7, or 9. Other bit maps are defined analogously. **/ bitinit() { int onebit, i; onebit = TESTBIT; i = 0; ptr = D; /* Set pointer ptr to point to D[0] */ do { bpresent[*(ptr++)] += onebit; } while(onebit >>= 1); for(i = 9; i >= 0; i--) { if(bpresent[i]) { maxdigit = i; break; } } for(i = 0; i < 10; i++) { if(bpresent[i]) { mindigit = i; break; } } b17 = bpresent[1] | bpresent[7]; b39 = bpresent[3] | bpresent[9]; b1379 = b17 | b39; b0369 = bpresent[0] | bpresent[3]; b0369 |= bpresent[6] | bpresent[9]; bnot0369 = (~(b0369)) & ONES; b147 = bpresent[1] | bpresent[4] | bpresent[7]; bnot147 = (~(b147)) & ONES; } /*** The function buildrows() takes the array brow0[32] defined * previously and uses it to create the arrays brow1[] through * brow4[]. ***/ buildbrows() { #define BL *(ptr1++) = (*(ptr2++) << 5) & ONES #define BR *(ptr1++) = (*(ptr2++) >> 5) #define BL8 BL; BL; BL; BL; BL; BL; BL; BL #define BR8 BR; BR; BR; BR; BR; BR; BR; BR #define BL32 BL8; BL8; BL8; BL8 #define BR32 BR8; BR8; BR8; BR8 ptr1 = brow0; ptr2 = brow1; BL32; ptr1 = brow2; ptr2 = brow1; BR32; ptr1 = brow3; ptr2 = brow2; BR32; ptr1 = brow4; ptr2 = brow3; BR32; } /*** The function buildmod3() fills in the vector mod3[0..84]. * This vector is used to compute k%3 by table lookup for * 0 <= k <= 81. For an arbitrary positive integer n having at * most 9 (decimal) digits, we can compute n%3 by forming the * digit sum of n (which we can easily do as we gradually build * up the value of n by winding through the grid), then * looking up the mod-3 value of this digit sum. **/ buildmod3() { #define M3 *(ptr++) = 0; *(ptr++) = 1; *(ptr++) = 2 #define M21 M3; M3; M3; M3; M3; M3; M3 ptr = mod3; M21; M21; M21; M21; } getdigits(filename) char *filename; { int f, i; char buf[50]; f = open(filename, 0); read(f, buf, 50); close(f); for(i = 0; i < 25; i++) { D[i] = buf[2*i] - '0'; if((D[i] < 0) || (D[i] > 9)) { printf("Error in input file.\n"); exit(1); } } #ifdef DEBUG write(1, buf, 50); #endif } /** NOTE: suffixinit() should be called AFTER bitinit() and buildbrows() ***/ /*** The function suffixinit() is allows us to do "mod-2 * lookahead" and "mod-5 lookahead." Clearly, any prime * having two or more decimal digits cannot end with an even digit * or with a "5." Sometimes we can look ahead several digits * and see that we cannot reach a legal fina; digit (or "suffix") * in time. * For example, if we're in * the process of building up the integer n in the grid * (beginning with the most significant digit), and if after 7 * digits we are at least 3 moves away from the nearest 1,3,7, or 9 * the only possible final digits), * then we can prune our search tree at this point and try a * different branch. * This initialization function constructs the array of * bit maps bdist[0],..., bdist[10] indicating the grid * locations from which a legal final digit can be reached in * exactly j steps (not necessarily by the most direct route). * Since every grid is guaranteed * by the contest rules to contain a prime, we know that (except * in pathological cases handled elsewhere in the program), there * must be a 1,3,7, or 9 somewhere in the grid. Furthermore, it * is clear that from any given grid location, we can reach * any other in at most 4 steps (and in any number of steps >= 4, * if we meander). * **/ suffixinit() { int bnew, bold; int i,j; int onebit; bdist[0] = b1379; /* Bit mask for positions with legal suffix*/ bold = b1379; onebit = TESTBIT; /* 0x1000000 */ bnew = BNBRS(bold); /* Find all nbrs of posns in bold, put in bnew */ bdist[1] = bnew; bold = bnew; bnew = BNBRS(bold); bdist[2] = bnew; bold = bnew; bnew = BNBRS(bold); bdist[3] = bnew; bold = bnew; bnew = BNBRS(bold); bdist[4] = bnew; /** EVERYONE is within 5 steps of a legal suffix. **/ bdist[5] = ONES; bdist[6] = ONES; bdist[7] = ONES; bdist[8] = ONES; bdist[9] = ONES; bdist[10] = ONES; } /*** In the function specialinit(bsame, bdiff), the parameter bsame is * a right-justified bit vector of grid locations having the same particular * value mod 3, and bdiff is the 25-bit * right-justified complementary vector. * This may be the most important function in the program, since * it allows us to do mod-3 lookahead (as suffixinit() allows us to * do mod-2 and mod-5 lookahead), thereby bypassing tens of thousands * of divisibility tests that would otherwise be required in certain * grids. (These grids typically are filled mostly with 0's, 3's, * 6's, and 9's, with just a small number (e.g., 1) of other * digits. Similar, but less ferocious, grids can be formed by * filling up most of the grid with 1's, 4's, and 7's.) * This function computes specialdist[i] for 0 <= i <= 9. * If all digits in the number currently being built up are 0,3,6, or * 9, then specialdist[i] is a bitmap of positions from which it is * possible to complete the number (typically 9 digits long) * using exactly i more digits (not necessarily in the most * direct way) so as to have a chance of forming a prime. * In particular, from our current grid location, we must * be able to reach a digit NOT divisible by 3 (else the whole * number will be divisible by 3), AND we still must end with a * 1,3,7, or 9 as before. * If all digits in the number being built up are so far * in {1,4,7} (and if the length of the number is a multiple * of 3, as 9 is), then similar things happen. ***/ specialinit(bsame, bdiff) int bsame, bdiff; { /* extern int specialdist[10]; */ int loc, radius, onebit, i; int bold, bnew, btemp, mindist; for(i = 0; i < 10; i++) specialdist[i] = 0; onebit = TESTBIT; /* 0x1000000 */ onebit <<= 1; /* Cancel out shift inside loop */ for(loc = 0; loc < 25; loc++) { onebit >>= 1; /** The test below skips grid location if it doesn't hold digit * of type "same." **/ if((bsame & onebit) == 0) { continue; } /** mindist will end up holding the minimum number of steps we must * take to get from current location in grid to a digit of type * "different" and then to a legal final digit (1,3,7, or 9). * The second leg may have length 0, but not the first leg * (conditioned on the fact that we're about to execute the * following instruction). **/ mindist = 9; bold = onebit; /** radius is distance of other grid locations in a ring from the * location currently being tested. **/ for(radius = 1; radius < 5; radius++) { bnew = BNBRS(bold); btemp = bnew & bdiff; #ifdef DEBUG printf("loc = %d, radius = %d.\n", loc, radius); printf("bnew:\n"); displaybits(bnew); printf("btemp:\n"); displaybits(btemp); #endif bold = bnew; /* For next iteration */ if(btemp == 0) { continue; } /** If we get here, then there are digits of type "different" * within the current ring. **/ if(btemp & bdist[0]) mindist = min(mindist,radius); else if(btemp & bdist[1]) mindist = min(mindist,radius+1); else if(btemp & bdist[2]) mindist = min(mindist,radius+2); else if(btemp & bdist[3]) mindist = min(mindist,radius+3); else { mindist = min(mindist,radius+4); } } /* End for(radius) */ for(i = mindist; i < 10; i++) { specialdist[i] |= onebit; } } /* End for(loc) */ } /* End specialinit() */ /** The function displaybits(bits) is just for debugging. * Added in jim14.c **/ displaybits(bits) int bits; { int i, onebit, bitvec[25]; onebit = TESTBIT; /* 0x1000000 */ for(i = 0; i < 25; i++) { bitvec[i] = ((bits & onebit) != 0); onebit >>=1; } for(i = 0; i < 25; i++) { printf("%2d", bitvec[i]); if((i%5) == 4) printf("\n"); } printf("\n"); } main(ac, av) char **av; { int k, bits; int i; /* Added in jim13.c **/ getdigits(av[1]); bits = 0; for(k = 0; k < 25; k++) bits |= 1 << D[k]; /*** We begin by doing a few easy tests to see whether the largest * prime in the gris is either 2,3,5,7, or 11. It is not hard * to prove that the test is both necessary and sufficient in the * case of 11, but the situation isn't so clear for the single- * digit primes. ***/ if((bits|341) == 341) { printf("2\n"); } else if((bits|585) == 585) { printf("3\n"); } else if((bits|373) == 373) { printf("5\n"); } else if((bits|129) == 129) { printf("7\n"); } else if(bits == 2) { printf("11\n"); } else { pinit(); /* Initialize table of primes. */ bitinit(); buildbrows(); buildmod3(); /* Added in jim9.c */ suffixinit(); /* Added in jim10.c */ #ifdef DEBUG for(i = 0; i < 10; i++) { printf("bdist[%d]:\n", i); displaybits(bdist[i]); } #endif /** Now I'll compute the two arrays of bitmasks for mod-3 * lookahead, namely bdist0369[10] and bdist147[10] . **/ specialinit(b0369, bnot0369); for(i = 0; i < 10; i++) { bdist0369[i] = specialdist[i]; #ifdef DEBUG printf("bdist0369[%d]:\n", i); displaybits(bdist0369[i]); #endif } specialinit(b147, bnot147); for(i = 0; i < 10; i++) { bdist147[i] = specialdist[i]; #ifdef DEBUG printf("bdist147[%d]:\n", i); displaybits(bdist147[i]); #endif } /*** I handle length-1 case * separately. The largest prime in the grid cannot be 11 unless * the grid is all 1's (a case for which I test above), but * it's conceivable that the largest prime in the grid could * be 2,3,5, or 7 for some grid not covered by the tests above. * In such a case, my program might make a mistake without handling * the length-1 case separately, because I assume that anything * divisible by 2,3, or 5 is composite. ****/ for(k = (bits|146) == 146 ? 8 : 9; k >= 2; k--) search(k); /** If we get here, then largest prime is 2, 3, 5, or 7. **/ if(bpresent[7]) { printf("7\n"); exit(0); } else if(bpresent[5]) { printf("5\n"); exit(0); } else if(bpresent[3]) { printf("3\n"); exit(0); } else if(bpresent[2]) { printf("2\n"); exit(0); } } /*End else (the one indicating that we don't have a special grid) */ } /*** The following function, search(len), is the workhorse of * the program. It searches for the largest prime of length * len decimal digits in the grid. ***/ search(len) int len; { int i, n, level, ntemp, btemp, cur, mincur, bnow; int onebit; int down; int mod3temp; int bsuftemp; /** Following variables were added in jim12.c **/ int sofar0369[10], sofar147[10]; /** The sofar stuff implements mod-3 lookahead; added in jim12.c **/ for(i = 1; i < 10; i++) { sofar0369[i] = 0; sofar147[i] = 0; } sofar0369[0] = 1; if((len == 9) || (len == 6) || (len == 3)) { sofar147[0] = 1; } else { sofar147[0] = 0; } onebit = TESTBIT; bgeneral[0] = ONES; /* 0x1ffffff; */ npart[0] = 0; digitsum[0] = 0; /* Added in jim9.c */ for(i = 0; i < 10; i++) { digit[i] = maxdigit; } level = 1; while(level) { if(level == len) { btemp = bgeneral[level-1] & b1379; if(btemp == 0) { level--; digit[level]--; continue; } ntemp = 10*npart[level-1]; mod3temp = mod3[digitsum[level-1]]; #define PTEST if(n%7){if(primetest(n)) {printf("%d\n",n);exit(0);}} #define DOPTEST(A,B) if((btemp & bpresent[A])&&(mod3temp - B)) { \ n = ntemp + A; PTEST } DOPTEST(9,0) DOPTEST(7,2) DOPTEST(3,0) DOPTEST(1,2) /* If I get here, then all final digits yielded composite numbers. */ level--; digit[level]--; continue; /* Might be unnec, but seems safer * to do this in case I add * something later. */ } else { /* Here level != len */ mincur=((level==1)?max(1,mindigit):mindigit); btemp = bgeneral[level-1]; down = 0; /* Next statement does lookahead mod 2 and mod 5. */ bsuftemp = bdist[len - level]; for(cur = digit[level]; cur >= mincur; cur--) { bnow = btemp & bpresent[cur]; if(bnow == 0) continue; if((bnow &= bsuftemp) == 0) continue; /** The sofar... stuff below checks to see whether all the digits * so far (_including_ the current digit) have come from {0,3,6,9} * or from {1,4,7}. (In the latter * case, we only use the mod-3 lookahead if len is divisible by 3.) * If so, then we do mod-3 lookahead. **/ if(sofar0369[level-1] && is0369[cur]){ specmask = bdist0369[len-level]; if((bnow & specmask) == 0) continue; sofar0369[level] = 1; sofar147[level] = 0; } else if(sofar147[level-1] && is147[cur]){ specmask = bdist147[len-level]; if((bnow & specmask) == 0) continue; sofar147[level] = 1; sofar0369[level] = 0; } else { sofar0369[level] = 0; sofar147[level] = 0; } /** If bnow == 0 and cur = mincur, we'll deal with this after * for(cur) loop. **/ digit[level] = cur; npart[level]=10*npart[level-1]+cur; /** JIM: Maybe do mult above by hand, with shifts and adds. **/ digitsum[level]=digitsum[level-1]+cur; /*** * btemp = 0; * * btemp |= brow0[bnow >> 20]; * btemp |= brow1[(bnow >> 15) & 0x1f]; * btemp |= brow2[(bnow >> 10) & 0x1f]; * btemp |= brow3[(bnow >> 5) & 0x1f]; * btemp |= brow4[bnow & 0x1f]; ****/ /* In jim10.c, the 6 lines above were replaced by the line below.*/ btemp = BNBRS(bnow); bgeneral[level] = btemp; down = 1; level++; digit[level] = maxdigit; break; /* Out of for(cur) loop */ } /* End for(cur) */ if(down == 0) { level--; digit[level]--; } } /* End else [level != len] */ } /* End while(level) */ } /* End of function search(len) */ /* * BOOTSIZE primes are computed to start. We need all primes less than * the square root of prime number TABSIZE. Note that since primes greater * than 10 may not appear explicitly, we can only bootstrap primes up to * 113, prime number 30; so BOOTSIZE must be at most 30. Note also that * if BOOTSIZE is 16 or more then primes greater than 49 are being booted * and so the bootstrap test needs an additional q%7 in it. * * SIEVESIZE should be at least k larger than prime number TABSIZE, where * k is half as big as the last prime in the bootstrap table. */ /* * Here's a table relating BOOTSIZE to TABSIZE. The first column * contains legitimate values for BOOTSIZE and the last is the maximum * value of TABSIZE for that BOOTSIZE. (These values have not all * been double-checked, though.) * * k Pk P(k+1) ^2 prevP maxT * 5 11 13 169 167 39 * 7 17 19 361 359 72 * 9 23 29 841 839 146 */ /** In jim11.c, I moved the following defs to front of pgm. * * #define TABSIZE 45 * #define BOOTSIZE 7 * #define SIEVESIZE 220 * * int P[TABSIZE + 10]; * int S[2*SIEVESIZE]; **/ /* * Initialize the prime table. Note that there is a subtle point here * where the while loop begins: the value of q should be equal to -1 * (mod 6). */ pinit() { int n = 0, p, q = 1, i, e, f; P[n++] = ++q; P[n++] = ++q; q++; P[n++] = ++q; do { ++q; if(q%2 && q%3 && q%5) P[n++] = q; } while(n < BOOTSIZE); for(i = 2; i < BOOTSIZE; i++) { p = P[i]; e = p << 1; f = e << 1; p = -p; while(p < SIEVESIZE) { S[p+=e] = 1; S[p+=f] = 1; } } while(q < SIEVESIZE) { if(S[q+=2] == 0) P[n++] = q; if(S[q+=4] == 0) P[n++] = q; } /** The creation of the array Pinv[] is new in jim11.c **/ /* The array of 64-bit floating-point numbers below is used * to do fast trial division. A bit of elementary number theory * is needed to see why the test works. **/ for(i = 0; i < TABSIZE; i++) Pinv[i] = 1.00000000000001 / P[i]; } /** Moved following def to beginning of pgm. **/ /* #define HIGHBIT 0x80000000 */ /* * Test for primality. This will work for any number less than 10^9, * unless n is one of the first TABSIZE primes. (It would be easy to * test this, but would add a little time to the program, and it can * be shown that if any prime from 13 to 229 inclusive appears in the * grid, then we will find some larger prime in the grid before * ever getting to this smaller prime. It is also easy to show that * 11 is the largest prime in the grid only if the table consists * of all 1's, and we test for that condition at the beginning of * the program. Finally, we handle 2,3,5, and 7 as special cases. * First test directly * for divisibility by small primes (the first TABSIZE of them). If * still undecided, perform Miller's test for base 2 and, if * necessary, for base 3. */ primetest(n) { register int s, ss, i, ii, r, rr, a, b, e, r2, q; /** register int *ip; **/ register unsigned t, tt; double nn, ninv, r1, r3, d, p; /* volatile int *ip; */ /** My compiler at CCR didn't recognize the keyword "volatile", * so I'll cross my fingers and drop the "volatile". ***/ int *ip; double *dptr; /* My own attempt to simulate "volatile" */ /** extern double dx[]; **/ /** JIM: Maybe make the trial divisions in-line later. **/ nn = n; /* dptr = &d; */ dptr = dx; /* Try to make sure this compiles OK */ /* ip = 1 + (int *)&d; */ ip = 1 + (int *)dptr; /** NOTE: Do NOT start with i = 0 below, because this fancy * division routine does NOT work for even denominators * (and P[0] = 2). It works for all odd denominators in the * range from 3 to 2000, for all strictly positive n < 2^31. **/ /* Start dividing by P[3] = 7, since the calling * function already handles mod 2, 3, and 5. */ for(i = 3; i < TABSIZE; i++) { /* d = nn * Pinv[i]; */ dx[0] = nn * Pinv[i]; if((*(ip) & 0x000fff00) == 0) { return 0; } } /*** After doing trial division by "small" primes, do Miller's * test for primality. Allan Wilks and I showed after considerable * work that because of the restrictions on the integers that * can occur in the grids, we need only do Miller's test with * the two bases 2 and 3. The full-fledged Miller test sometimes * has to be done for the four bases 2, 3, 5, and 7. * (Remark: Allan and I worked together on the primality * testing but maintained the proverbial wall of separation * with regard to searching the grid.) ***/ t = n-1; s = -1; while((t&1) == 0) { t >>= 1; s++; } i = 32; while((t&HIGHBIT) == 0) { t <<= 1; i--; } ss = s; tt = t; ii = i; ninv = 1.0/nn; /* strongly pseudoprime to the base 2? */ if(i >= 4) { r = 1<<(t>>28); if(r >= n) r %= n; /* Fixed this in wookin8.c */ t <<= 4; i -= 4; } else r = 1; /** NOTE: The "else r = 1;" above is a bug fix in wookin7.c . **/ while(i--) { /** The following code computes r = (r*r)%n . **/ r1 = r>>8; r2 = r&0xff; r3 = (r<<1) - r2; d = r1*r1; q = (d+0.5)*ninv; d = (d-nn*q)*65536.0 + r3*r2; q = (d+0.5)*ninv; r = d-nn*q; if(t&HIGHBIT) if((r <<= 1) >= n) r -= n; t <<= 1; } if(r != 1 && r != n-1) while(s-- > 0) { r1 = r>>8; r2 = r&0xff; r3 = (r<<1) - r2; d = r1*r1; q = (d+0.5)*ninv; d = (d-nn*q)*65536.0 + r3*r2; q = (d+0.5)*ninv; r = d-nn*q; if(r == n-1) break; } if(s < 0) return 0; /* strongly pseudoprime to the base 3? */ s = ss; t = tt; i = ii; r = 1; while(i--) { r1 = r>>8; r2 = r&0xff; r3 = (r<<1) - r2; d = r1*r1; q = (d+0.5)*ninv; d = (d-nn*q)*65536.0 + r3*r2; q = (d+0.5)*ninv; r = d-nn*q; if(t&HIGHBIT) { rr = r; if((r <<= 1) >= n) r -= n; if((r += rr) >= n) r -= n; } t <<= 1; } if(r != 1 && r != n-1) while(s-- > 0) { r1 = r>>8; r2 = r&0xff; r3 = (r<<1) - r2; d = r1*r1; q = (d+0.5)*ninv; d = (d-nn*q)*65536.0 + r3*r2; q = (d+0.5)*ninv; r = d-nn*q; if(r == n-1) break; } if(s < 0) return 0; /* If we get here, n must be prime! */ return 1; }