Al Zimmermann sponsors programming contests about three times a year. There are many programming contests on the Internet, but Al's are special in many ways.
It's rude to gloat, but I can't help myself. The Al Zimmermann Programming Contests are real challenges, and winning one feels like a good accomplishment. Resumes are all about ``buzzwords'' and I suppose Simulated Annealing, Genetic Algorithm, and Primality testing might be the buzzwords for this contest, but I can claim no special expertise in any of those areas. Instead general skill and determination were far more relevant. Especially determination: I think everyone who's won one of Al's Contests devoted a lot of effort to it. I'm sure being semi-retired gave me an unfair advantage over contestants who were trying to squeeze in contest activity while doing a full-time job, but I'll take the win anyway!
I've never entered an Internet programming contest except for Al's. Some contests require, to win, that solutions be submitted within a few hours. (This would have worked for me when I was on ``Internet time'' but I'm semi-retired and taking time to raise my children and ``smell the flowers.'') Another ``contest'' requires that source code be submitted to solve a difficult problem with commercial value. One contest tests submissions for speed after compiling them without any compiler optimization! (Hunh? The Obfuscated C Contest makes more sense than that.)
This last point may be unclear to some. After all,
all submissions are subject to the same rule.
But techniques for optimizing code without compiler
assistance are quite different than normal techniques.
For example
#ifdef SANE_RULES
x = foo();
bar(x * 8);
#else
/* The following is uglier to read but better
* when we're not allowed to use the optimizer.
*/
bar((x = foo()) << 3);
#endif
(Gcc is an excellent compiler and would replace the multiply
with shift even without `-O', but it would still miss that
x was sitting in a register. And only with `-O'
does Gcc use the high-speed multiply by 5 available on most processors.)
Fill a 19-by-19 grid with the digits 0 through 9. Score k points for every distinct k-digit prime found in the grid, reading in either direction in any row, column or diagonal. Add a point for every even digit. Finally divide your score by the best score of any contestant. Repeat this for an 18-by-18 grid, 17-by-17, and so on all the way down to 3-by-3. (This is the ``Part B'' scoring, which I won. Mark Beyleveld won ``Part A'' with a slightly different scoring schedule.)
Sound silly? I was afraid to explain it to my drinking buddies until one of them mentioned that he used to hire airplanes to cross the South China Sea just to play in a Yahtzee (Balut) Tournament! Compared with that, entering an Al Zimmermann Programming Contest seems almost ... adult!
It is interesting that two of the best and most popular optimization heuristics are based on simulating nature: Genetic Algorithm (GA), Simulated Annealing. Each of these methods has an interesting variation: GA with 2-D genome, allean.
I probably wouldn't have participated in this contest at all, except that when I first saw it, what popped into my mind was ``I know how to do this! GA with the 2-D grid being the genome! Crossover is done by cutting at various angles, and mutations include shifts.'' That's how my earliest results were generated. I still think that might have been a good approach for the Part A contest but, for some reason, I decided to focus on Part B, where the 2-D genome made less sense. I still did use GA somewhat, e.g. in the pattbreed program.
I have a variation on Simulated Annealing which I call allean. It's simpler than Simulated Annealing, and the ``temperature'' is easier to tune. Maybe some day I'll have time to do experiments and see how it really compares for optimization speed.
It seems interesting that of 30-odd functions I needed to write for the main program (``allean''), over half of them were essentially tests for primality! (In some cases they were partial tests or tests for distinct primes.)
I had special processing for one-digit primes, a simple table-lookup routine for numbers between 10 and 100 million, and five different routines for numbers above 100 million, depending on which intermediate results required 32-, 64-, 96- or 128-bit arithmetic. The FSF Gcc compiler supports 64-bit (long long) arithmetic, but for larger numbers I wrote my own subroutines (though I guess Gcc's library has support for them too).
After checking divisibility by a few dozen small primes
(including 281 to dispose of the nasty pseudoprime 82929001)
I tested for primality with the simple Fermat equation
pow(2, p-1) % p == 1.
By a command line option (or when initially generating the
table for numbers less than 100 million),
the check could be extended to 3, 5, 7, etc.
but I found this to be a waste of time.
(I did waste the cycles
for final submissions and, at that point, never detected a composite
pseudo-prime.)
I maintained a (multi-compartment) hash table of all numbers encountered. This served two purposes: to avoid repeating a primality test, and to detect whether a prime was distinct or not.
Primes in the grid had to be maintained for the algorithm to function properly, but other entries in the hash table were ``soft'' and could be discarded to save space. (I don't recall encountering such a cache before, with both ``hard'' and ``soft'' entries.) For large N the dynamic hash table would eventually get too large, so I started discarding nonessential entries when it got to about 25 MB. This let me run two large-N alleans concurrently (though I needed to stop one whenever I wanted to run a hog like Acrobat in another window).
Preparing a high-scoring N-by-N grid was done in five stages:
Notes:
The `Stage 4' program (allean) was the most important part of the procedure. It operated by iteratively picking a grid cell randomly and computing the effect on the score of each of the nine possible changes at that cell. (Since the new grids differ at only one cell from the previous grid, it isn't necessary to score the entire grid, just to do (and undo) the numbers containing the given cell.) Improvements were incorporated into the grid under construction, along with changes to roughly mimic simulated annealing. Allean consisted of six C modules linked together:
Of these modules, the one with scoring functions is perhaps the most interesting, and the only one really related to the specific rules of this contest. I've included the source code to that module below. (If there's significant interest I may eventually provide source for other parts of the system as well.)
For this contest, coding and testing of software took much less time than directing and studying specific experiments. I used different techniques for different sizes. (For example, my submission for the 18x18 case was developed in a complicated way similar to the 13x13 grid described next, but the 19x19 grid was developed quickly by starting with the 18x18 grid.)
I will describe in detail how I derived the grid I submitted for the 13x13 case. I'll show the actual commands executed on my computer. Remember that a tree structure of related commands (or the same commands with slightly different parameters) was executed: the commands shown here are just the specific evolutionary path which led to the winning grid.
In the following detailed list of the specific steps which produced the winning grid, the following commands are used:
# Generate some nice patterns; this took a few hours.
# (Nevermind how Biglist.13 was generated.)
pattbreed 13 < Biglist.13 >> PB13.log2
# I extracted the 22 best-scoring patterns from the output
# of pattbreed 13.
vi Patt{A,B,C,D,E,F,G,H,I,J,K,L,M,N,O,P,Q,R,S,T,U,V}.13
# One of them looked particularly promising.
cp PattB.13 PattBz1.13
vi PattBz1.13
# Here are the contents of PattBz1.13.
# Note that six rows, seven columns, and both
# long diagonals contain 13-digit numbers drawn from
# "Biglist.13" -- a list of over 1000 high-scoring
# numbers. The '*'s were manually added by me to
# establish a set of 19 points to be optimized first.
# 3999173332933
# 9913396199933
# 9933739181939
# 7963**3**1*86
# 91237*3*1**97
# 3971*315***73
# 3391939373213
# 7752*789...83
# 63192*1.6..79
# 8611**9..7.91
# 1399**1...690
# 9933799319333
# 9333013961999
# First optimize the 19 points marked '*'
allean S 3 10 40 < PattBz1.13 >> PattBz1.13.log1
tail -13 PattBz1.13.log1 > PattBz2.13
# Next optimize all 169 points.
# Several such alleans were run; the best scored
# 8242 and was selected for further work.
allean -e 1 R 4 17 70 < PattBz2.13 >> PattBz2.13.log2
tail -13 PattBz2.13.log2 > B8242.13
# Compare the result of this optimization with
# the original pattern constructed by pattbreed.
# Often there would be very little overlap (though I
# guessed without proof that the "good" starting
# pattern might be helpful anyway).
# In this case however, seven of the original 13-digit
# numbers in PattB.13 are fully intact and some others
# are partially intact. I take this to be a very
# good sign.
cat B8242.13 PattB.13 | intersection 13 > PattBb1.13
# Treat this result as a new starting pattern and
# prepare to run 'allean S' again.
cp PattBb1.13 PattBb2.13
vi PattBb2.13
# Here are the contents of PattBb2.13.
# The '*' are set manually as in PattBz1.13 above,
# but now the digits are the result of the alleans
# preceding.
# 3999173332933
# 99*33*61**933
# 9933739181939
# .**3**3**1...
# 9**37*3*1..9.
# 3*71*315*...3
# 3391939373213
# 7**2*789***.3
# ***92*1*6*...
# ..*1**9**7.9.
# .3.9*.1...6..
# 9933799319333
# 9333.139..9.9
allean S 7 19 50 < PattBb2.13 >> PattBb2.13.log5
tail -13 PattBb2.13.log5 > PattBc2.13
# Since the digits in PattBb2.13 are *known* to be
# good, I do not allow them to change in the
# following.
allean T 5 19 50 < PattBc2.13 >> PattBc2.13.log3
tail -13 PattBc2.13.log3 > PattBc3.13
# Finally allean all 169 points.
allean R 7 19 50 < PattBc3.13 >> PattBc3.13.log5
tail -13 PattBc3.13.log5 > B8520.13
# This output was the winning submission.
# After the contest I noted that an ever better 13x13 grid
# (scoring 8579) is sometimes obtained by the last allean:
# 3339333673393
# 3933316493933
# 9933739181939
# 9793493911932
# 3993193513313
# 7341591589553
# 3391939373213
# 9772378937777
# 4939291139391
# 3931369991399
# 3399931639339
# 9333139679399
# 3333133181333
As seen, this module can be run with either the Part A or Part B scoring. (No other changes would be required for Part A, though as mentioned above, I'd intended to use a different approach if I'd gotten around to working on Part A.)
The module also supported a third heuristic scoring, useful for partial grids. For example, primes consisting mainly of the digits `3' and `9' (e.g., 397) will be very plentiful in a finished grid so there is little or need to look for them deliberately. This scoring was used in ``Stage 1'' (above), and some other experiments, but not in the final ``Stage 4'' program, so I've eliminated it in the version here.
``Do as I say, not as I do!'' The cntprimes() routine is complicated by excessive optimization (e.g., switching to 32-bit instead of 64-bit variables when I can), but, although I understand that excessive optimization is usually a flaw, I never met an optimization I didn't like!
(Your browser may corrupt the view of this C code, but you
can get a pristine version using ``Save as'' or ``View Source''
and deleting everything except the routine.)
/* This is the scoring routine
* for Primal Squares, the July 2005 Zimmerman contest.
* It assumes long and long long are 32- and 64-bits resp.
* Same module is used for Part_A and Part_B of the contest:
* define PART_B and recompile for the latter.
* The module contains 5 entry-points:
* cntprimes() -- score along one ray.
* pointscore() -- score contributed by one point
* gridscore() -- score entire grid
* singmask(), issing() -- helpers for 1-digit primes
* This module calls two external functions:
* dipri_i() deletes or adds a (32-bit) number to table
* and returns non-zero if it is/was a unique prime.
* dipri_l() is same but argument is 64-bit.
*
* Notes:
* (a) One-digit primes are handled separately.
* (b) allean.c calls gridscore() only once for initialization,
* so pointscore() is the real workhorse.
* (c) In addition to the ten digits 0-9, routines must
* handle an eleventh cell-value (Unset): this facilitates
* the annealing of partial grids.
* (d) Although the primality routines return quickly
* for multiples of small primes, multiples of 2,3,5,7
* are detected in cntprimes() for utmost speed.
* (This detection is trivial for multiples of 2 and 5,
* and is bypassed for 3 and 7 for numbers < 100 million.)
*/
typedef unsigned long long int Largenum;
typedef unsigned long int Integer;
/* Of the numerals
* {sp[0] sp[inc]},
* {sp[0] sp[inc] sp[2*inc]},
* {sp[0] sp[inc] sp[2*inc] ... sp[(mlen-1)*inc]}
* How many are prime?
* (Ignore numbers of fewer than (sklen+1) digits.)
* 1-digit numbers are handled by a separate routine.
*
* In this routine and next, `mode' is 1 for doing (adding primes),
* and -1 for undoing.
*/
#define PDIG(x) (1 << (x) & 0x0ac) /* 2,3,5,7 */
#define GDIG(x) (1 << (x) & 0x28a) /* 1,3,7,9 */
#define SSDG(x) (1 << (x) & 0x174) /* 2,4,5,6,8 */
#define NMU3(x) (1 << (x) & 0x6db6) /* 1,2,4,5,7,8,10,11,13,14 */
cntprimes(sp, sklen, mlen, inc, mode)
char *sp;
{
int dipri_i(), dipri_l();
int nump;
unsigned int dig, res3, res7;
Integer accum;
Largenum lacc;
nump = 0;
lacc = accum = *sp;
if (9 < accum || accum < 0)
goto done;
if (sklen >= mlen)
goto largsk;
if (sklen && accum == 0)
goto done;
#define NEXTDIG(REG) \
sp += inc; \
dig = *sp; \
if (9 < dig || dig < 0) goto done; \
REG = REG * 10 + dig
#ifdef PART_B
#define PSCORE(LEN) LEN
#else
#define PSCORE(LEN) 1
#endif
#define PTEST(LEN) (LEN < 10 ? dipri_i : dipri_l)
#define PROCDIG(REG, LEN) \
if (GDIG(dig) && (LEN > 5 ? NMU3(res3) && res7 : 1)) { \
if (mode > 0) { \
if (PTEST(LEN)(REG, 1)) \
nump += PSCORE(LEN); \
} else { \
if (PTEST(LEN)(REG, -1)) \
nump -= PSCORE(LEN); \
} \
} \
if (! --mlen) goto done
#define FIX37 \
if ((res3 += dig) >= 15) res3 -= 15; \
res7 = (res7 * 3 + dig) % 7
switch (sklen) {
default:
case 19:
largsk:
printf("Large skip %d\n", sklen);
return 0;
case 18: NEXTDIG(accum);
case 17: NEXTDIG(accum);
case 16: NEXTDIG(accum);
case 15: NEXTDIG(accum);
case 14: NEXTDIG(accum);
case 13: NEXTDIG(accum);
case 12: NEXTDIG(accum);
case 11: NEXTDIG(accum);
lacc = accum;
case 10: NEXTDIG(lacc);
case 9: NEXTDIG(lacc);
case 8: NEXTDIG(lacc);
case 7: NEXTDIG(lacc);
case 6: NEXTDIG(lacc);
case 5: NEXTDIG(lacc);
case 4: NEXTDIG(lacc);
case 3: NEXTDIG(lacc);
case 2: NEXTDIG(lacc);
case 1: NEXTDIG(lacc);
case 0:
break;
}
mlen -= sklen;
accum = lacc;
res3 = lacc % 3;
res7 = lacc % 7;
switch (sklen) {
case 0:
/* Do NOT Check 1-digit number for primality */
#ifdef PART_B
if (inc == 1) {
/* Do the even-digit scoring for Part B */
if (mode > 0)
nump += ! (accum & 1);
else
nump -= ! (accum & 1);
}
#endif
/* Disallow leading zero */
if (accum == 0) goto done;
if (! --mlen) goto done;
NEXTDIG(accum);
case 1:
PROCDIG(accum, 2);
NEXTDIG(accum);
case 2:
PROCDIG(accum, 3);
NEXTDIG(accum);
case 3:
PROCDIG(accum, 4);
NEXTDIG(accum);
case 4:
PROCDIG(accum, 5);
NEXTDIG(accum);
case 5:
/* Before here, table lookup makes us fast.
* After here, we'll skip the function call
* for multiples of 3 or 7.
*/
res3 = accum % 3;
res7 = accum % 7;
PROCDIG(accum, 6);
NEXTDIG(accum);
FIX37;
case 6:
PROCDIG(accum, 7);
NEXTDIG(accum);
FIX37;
case 7:
PROCDIG(accum, 8);
NEXTDIG(accum);
FIX37;
case 8:
PROCDIG(accum, 9);
lacc = accum;
NEXTDIG(lacc);
FIX37;
case 9:
PROCDIG(lacc, 10);
NEXTDIG(lacc);
FIX37;
case 10:
PROCDIG(lacc, 11);
NEXTDIG(lacc);
FIX37;
case 11:
PROCDIG(lacc, 12);
NEXTDIG(lacc);
FIX37;
case 12:
PROCDIG(lacc, 13);
NEXTDIG(lacc);
FIX37;
case 13:
PROCDIG(lacc, 14);
NEXTDIG(lacc);
FIX37;
case 14:
PROCDIG(lacc, 15);
NEXTDIG(lacc);
FIX37;
case 15:
PROCDIG(lacc, 16);
NEXTDIG(lacc);
FIX37;
case 16:
PROCDIG(lacc, 17);
NEXTDIG(lacc);
FIX37;
case 17:
PROCDIG(lacc, 18);
NEXTDIG(lacc);
FIX37;
case 18:
PROCDIG(lacc, 19);
done:
return nump;
default:
printf("Twenty digit number submitted ... aborting\n");
exit(1);
}
}
/* Return only the score contribution from cell xy.
* This routine is complicated enough that it
* should be written *after* coffee and *before* beer.
*/
pointscore(b, wid, hei, xy, mode)
char *b;
{
char *bb;
int x, y, score;
int x1, y1;
y1 = xy / wid;
x1 = xy % wid;
score = 0;
bb = b + wid * y1;
for (x = 0; x <= x1; x++, bb++)
score += cntprimes(bb, x1 - x, wid - x, 1, mode);
for (x--, bb--; x < wid; x++, bb++)
score += cntprimes(bb, x - x1, x + 1, -1, mode);
bb = b + x1;
for (y = 0; y <= y1; y++, bb += wid)
score += cntprimes(bb, y1 - y, hei - y, wid, mode);
for (y--, bb -= wid; y < hei; y++, bb += wid)
score += cntprimes(bb, y - y1, y + 1, -wid, mode);
if (x1 > y1) {
bb = b + x1 - y1;
for (y = 0; y <= y1; y++, bb += wid + 1)
score += cntprimes(bb, y1 - y,
wid + y1 - x1 - y, wid + 1, mode);
for (y--, bb -= wid + 1; y + x1 - y1 < hei; y++, bb += wid + 1)
score += cntprimes(bb, y - y1,
y + 1, -wid - 1, mode);
} else {
bb = b + (y1 - x1) * wid;
for (y = y1 - x1; y <= y1; y++, bb += wid + 1)
score += cntprimes(bb, y1 - y,
hei - y, wid + 1, mode);
for (y--, bb -= wid + 1; y < hei; y++, bb += wid + 1)
score += cntprimes(bb, y - y1,
y + x1 - y1 + 1, -wid - 1, mode);
}
if (x1 + y1 < hei) {
bb = b + x1 + y1;
for (y = 0; y <= y1; y++, bb += wid - 1)
score += cntprimes(bb, y1 - y,
x1 + y1 + 1 - y, wid - 1, mode);
for (y--, bb -= wid - 1; y <= x1 + y1; y++, bb += wid - 1)
score += cntprimes(bb, y - y1,
y + 1, -wid + 1, mode);
} else {
bb = b + (y1 + x1 - wid + 2) * wid - 1;
for (y = y1 + x1 - wid + 1; y <= y1; y++, bb += wid - 1)
score += cntprimes(bb, y1 - y,
hei - y, wid - 1, mode);
for (y--, bb -= wid - 1; y < hei; y++, bb += wid - 1)
score += cntprimes(bb, y - y1,
y - y1 - x1 + wid, -wid + 1, mode);
}
return score;
}
singmask(b, area)
char *b;
{
int x, y;
for (y = x = 0; x < area; x++) {
y |= (1 << b[x]) & 0xac;
if (y == 0xac)
break;
}
return y;
}
issing(v, msk)
{
int t;
t = (1 << v) & 0xac;
return t ? ! (t & msk) : 0;
}
/* Compute score for entire grid; assumes wid == hei */
gridscore(b, wid, hei)
char *b;
{
int x, y, score = 0;
char *bb;
for (bb = b, y = 0; y < hei; y++, bb += wid)
for (x = 0; x < wid; x++) {
score += cntprimes(bb + x, 0,
wid - x, 1, 1);
score += cntprimes(bb + x, 0,
x + 1, -1, 1);
score += cntprimes(bb + x, 0,
hei - y, wid, 1);
score += cntprimes(bb + x, 0,
y + 1, -wid, 1);
score += cntprimes(bb + x, 0,
x > y ? wid - x : hei - y, wid + 1, 1);
score += cntprimes(bb + x, 0,
x + 1 < hei - y ? x + 1 : hei - y, wid - 1, 1);
score += cntprimes(bb + x, 0,
wid - x < y + 1 ? wid - x : y + 1, -wid + 1, 1);
score += cntprimes(bb + x, 0,
x > y ? y + 1 : x + 1, -wid - 1, 1);
}
/* Locate single-digit primes */
y = singmask(b, wid * hei);
/* add 1 to score if grid contains 7, etc. */
score += !! (y & 0x80);
score += !! (y & 0x20);
score += !! (y & 0x08);
score += !! (y & 0x04);
return score;
}