# prime number generator
# Global associative array _primegen stores the state of all active prime
# number generators. It is a two-dimensional array with an identifier (a
# positive integer) in its first dimension and an element of state in its
# second dimension, intended only for internal use by the prime number
# generator; it should not be modified by client code. The possibilities are:
#
# _primegen[k, "init"] -- initial primes 2 and 3 to "prime" the pump (sorry!)
# _primegen[k, "pgen"] -- the k-number of the internal prime number generator
# _primegen[k, "c"] -- the next candidate prime number
# _primegen[k, "p"] -- the current largest sieving prime
# _primegen[k, p] -- the stride associated with sieving prime p where
# p in (3, 5, 7, 11, ...); notice no stride for 2
#
# Global variable _primegen_counter is the largest currently-unused k; it is
# initially zero.
#
# Function primegen_init() creates a new prime number generator that returns
# primes starting from zero. It returns a state number (the k of the _primegen
# array). Function primegen(k) returns the next prime number from prime number
# sequence k, and resets the _primegen[k] generator state for the next call to
# primegen(k).
function primegen_init( k) {
k = _primegen_counter++
_primegen[k,"init"] = 2
return k }
function primegen(k, pgen,c,p,q,s,m) {
# initial primes 2 and 3
if ((k,"init") in _primegen) {
if (_primegen[k,"init"] == 2) {
_primegen[k,"init"] = 3
return 2 }
else {
delete _primegen[k,"init"]
return 3 } }
# set up for remaining primes after 3
if (! (k,"pgen") in _primegen) {
pgen = primegen_init()
_primegen[k,"pgen"] = pgen
p = primegen(pgen)
p = primegen(pgen)
_primegen[k,"c"] = 3
_primegen[k,"p"] = p
_primegen[k,p] = p+p }
# main prime generation routine
pgen = _primegen[k,"pgen"]
c = _primegen[k,"c"]
p = _primegen[k,"p"]; q = p * p
while (1) {
c += 2
if ((k,c) in _primegen) {
s = _primegen[k,c]; m = c + s
delete _primegen[k,c]
while ((k,m) in _primegen) m += s
_primegen[k,m] = s }
else if (c < q) {
_primegen[k,"c"] = c
_primegen[k,"p"] = p
return c }
else {
s = p + p; m = c + s
while ((k,m) in _primegen) m += s
_primegen[k,m] = s
p = primegen(pgen); q = p * p } } }
BEGIN { ps = primegen_init()
for (p = primegen(ps); p < 100; p = primegen(ps)) {
for (x in _primegen) print x, _primegen[x] }
IyBwcmltZSBudW1iZXIgZ2VuZXJhdG9yCgojIEdsb2JhbCBhc3NvY2lhdGl2ZSBhcnJheSBfcHJpbWVnZW4gc3RvcmVzIHRoZSBzdGF0ZSBvZiBhbGwgYWN0aXZlIHByaW1lCiMgbnVtYmVyIGdlbmVyYXRvcnMuIEl0IGlzIGEgdHdvLWRpbWVuc2lvbmFsIGFycmF5IHdpdGggYW4gaWRlbnRpZmllciAoYQojIHBvc2l0aXZlIGludGVnZXIpIGluIGl0cyBmaXJzdCBkaW1lbnNpb24gYW5kIGFuIGVsZW1lbnQgb2Ygc3RhdGUgaW4gaXRzCiMgc2Vjb25kIGRpbWVuc2lvbiwgaW50ZW5kZWQgb25seSBmb3IgaW50ZXJuYWwgdXNlIGJ5IHRoZSBwcmltZSBudW1iZXIKIyBnZW5lcmF0b3I7IGl0IHNob3VsZCBub3QgYmUgbW9kaWZpZWQgYnkgY2xpZW50IGNvZGUuIFRoZSBwb3NzaWJpbGl0aWVzIGFyZToKIwojICAgICBfcHJpbWVnZW5baywgImluaXQiXSAtLSBpbml0aWFsIHByaW1lcyAyIGFuZCAzIHRvICJwcmltZSIgdGhlIHB1bXAgKHNvcnJ5ISkKIyAgICAgX3ByaW1lZ2VuW2ssICJwZ2VuIl0gLS0gdGhlIGstbnVtYmVyIG9mIHRoZSBpbnRlcm5hbCBwcmltZSBudW1iZXIgZ2VuZXJhdG9yCiMgICAgIF9wcmltZWdlbltrLCAiYyJdICAgIC0tIHRoZSBuZXh0IGNhbmRpZGF0ZSBwcmltZSBudW1iZXIKIyAgICAgX3ByaW1lZ2VuW2ssICJwIl0gICAgLS0gdGhlIGN1cnJlbnQgbGFyZ2VzdCBzaWV2aW5nIHByaW1lCiMgICAgIF9wcmltZWdlbltrLCBwXSAgICAgIC0tIHRoZSBzdHJpZGUgYXNzb2NpYXRlZCB3aXRoIHNpZXZpbmcgcHJpbWUgcCB3aGVyZQojICAgICAgICAgICAgICAgICAgICAgICAgICAgICBwIGluICgzLCA1LCA3LCAxMSwgLi4uKTsgbm90aWNlIG5vIHN0cmlkZSBmb3IgMgojCiMgR2xvYmFsIHZhcmlhYmxlIF9wcmltZWdlbl9jb3VudGVyIGlzIHRoZSBsYXJnZXN0IGN1cnJlbnRseS11bnVzZWQgazsgaXQgaXMKIyBpbml0aWFsbHkgemVyby4KIwojIEZ1bmN0aW9uIHByaW1lZ2VuX2luaXQoKSBjcmVhdGVzIGEgbmV3IHByaW1lIG51bWJlciBnZW5lcmF0b3IgdGhhdCByZXR1cm5zCiMgcHJpbWVzIHN0YXJ0aW5nIGZyb20gemVyby4gSXQgcmV0dXJucyBhIHN0YXRlIG51bWJlciAodGhlIGsgb2YgIHRoZSBfcHJpbWVnZW4KIyBhcnJheSkuIEZ1bmN0aW9uIHByaW1lZ2VuKGspIHJldHVybnMgdGhlIG5leHQgcHJpbWUgbnVtYmVyIGZyb20gcHJpbWUgbnVtYmVyCiMgc2VxdWVuY2UgaywgYW5kIHJlc2V0cyB0aGUgX3ByaW1lZ2VuW2tdIGdlbmVyYXRvciBzdGF0ZSBmb3IgdGhlIG5leHQgY2FsbCB0bwojIHByaW1lZ2VuKGspLgogICAgICAgIApmdW5jdGlvbiBwcmltZWdlbl9pbml0KCAgICBrKSB7CiAgICBrID0gX3ByaW1lZ2VuX2NvdW50ZXIrKwogICAgX3ByaW1lZ2VuW2ssImluaXQiXSA9IDIKICAgIHJldHVybiBrIH0KCmZ1bmN0aW9uIHByaW1lZ2VuKGssICAgIHBnZW4sYyxwLHEscyxtKSB7CgogICAgIyBpbml0aWFsIHByaW1lcyAyIGFuZCAzCiAgICBpZiAoKGssImluaXQiKSBpbiBfcHJpbWVnZW4pIHsKICAgICAgICBpZiAoX3ByaW1lZ2VuW2ssImluaXQiXSA9PSAyKSB7CiAgICAgICAgICAgIF9wcmltZWdlbltrLCJpbml0Il0gPSAzCiAgICAgICAgICAgIHJldHVybiAyIH0KICAgICAgICBlbHNlIHsKICAgICAgICAgICAgZGVsZXRlIF9wcmltZWdlbltrLCJpbml0Il0KICAgICAgICAgICAgcmV0dXJuIDMgfSB9CgogICAgIyBzZXQgdXAgZm9yIHJlbWFpbmluZyBwcmltZXMgYWZ0ZXIgMwogICAgaWYgKCEgKGssInBnZW4iKSBpbiBfcHJpbWVnZW4pIHsKICAgICAgICBwZ2VuID0gcHJpbWVnZW5faW5pdCgpCiAgICAgICAgX3ByaW1lZ2VuW2ssInBnZW4iXSA9IHBnZW4KICAgICAgICBwID0gcHJpbWVnZW4ocGdlbikKICAgICAgICBwID0gcHJpbWVnZW4ocGdlbikKICAgICAgICBfcHJpbWVnZW5baywiYyJdID0gMwogICAgICAgIF9wcmltZWdlbltrLCJwIl0gPSBwCiAgICAgICAgX3ByaW1lZ2VuW2sscF0gPSBwK3AgfQoKICAgICMgbWFpbiBwcmltZSBnZW5lcmF0aW9uIHJvdXRpbmUKICAgIHBnZW4gPSBfcHJpbWVnZW5baywicGdlbiJdCiAgICBjID0gX3ByaW1lZ2VuW2ssImMiXQogICAgcCA9IF9wcmltZWdlbltrLCJwIl07IHEgPSBwICogcAogICAgd2hpbGUgKDEpIHsKICAgICAgICBjICs9IDIKICAgICAgICBpZiAoKGssYykgaW4gX3ByaW1lZ2VuKSB7CiAgICAgICAgICAgIHMgPSBfcHJpbWVnZW5bayxjXTsgbSA9IGMgKyBzCiAgICAgICAgICAgIGRlbGV0ZSBfcHJpbWVnZW5bayxjXQogICAgICAgICAgICB3aGlsZSAoKGssbSkgaW4gX3ByaW1lZ2VuKSBtICs9IHMKICAgICAgICAgICAgX3ByaW1lZ2VuW2ssbV0gPSBzIH0KICAgICAgICBlbHNlIGlmIChjIDwgcSkgewogICAgICAgICAgICBfcHJpbWVnZW5baywiYyJdID0gYwogICAgICAgICAgICBfcHJpbWVnZW5baywicCJdID0gcAogICAgICAgICAgICByZXR1cm4gYyB9CiAgICAgICAgZWxzZSB7CiAgICAgICAgICAgIHMgPSBwICsgcDsgbSA9IGMgKyBzCiAgICAgICAgICAgIHdoaWxlICgoayxtKSBpbiBfcHJpbWVnZW4pIG0gKz0gcwogICAgICAgICAgICBfcHJpbWVnZW5bayxtXSA9IHMKICAgICAgICAgICAgcCA9IHByaW1lZ2VuKHBnZW4pOyBxID0gcCAqIHAgfSB9IH0KCkJFR0lOIHsgcHMgPSBwcmltZWdlbl9pbml0KCkKICAgICAgICBmb3IgKHAgPSBwcmltZWdlbihwcyk7IHAgPCAxMDA7IHAgPSBwcmltZWdlbihwcykpIHsKICAgICAgICAgICAgcGNvdW50Kys7IHByaW50ZiAiJWQgIiwgcCB9CiAgICAgICAgcHJpbnRmICIlZFxuIiwgcGNvdW50CiAgICAgICAgZm9yICh4IGluIF9wcmltZWdlbikgcHJpbnQgeCwgX3ByaW1lZ2VuW3hdIH0=