sub primes {
my ($upto, $func) = @_;
return large_primes
($upto, $func) if $upto > 25000; $upto -= 1;
my @flags = (1) x $upto;
foreach my $i (0..$upto - 1){
if ($flags[$i]){
my $prime = $i + 2;
my $k = $i + $prime;
while($k < $upto){
$flags[$k] = 0;
$k += $prime;
}
$func->($prime);
}
}
}
sub large_primes {
my ($upto, $func) = @_;
# The Algorithm is adapted from http://w...content-available-to-author-only...k.com/~jrm/printprimes.html
# This code is a translated version of Andreas Raab's #largePrimesUpTo:do: method
# which was written for the Squeak Smalltalk environment.
my $idx_limit = int(sqrt($upto)) + 1; my @flags = (0xFF
) x
(int(($upto + 2309) / 2310) * 60 + 60); my @primes_up_to_2310 = ();
primes
(2310, sub{ push(@primes_up_to_2310, $_[0]) }); my @mask_bit_idx;
$#mask_bit_idx = 2309;
$mask_bit_idx[0] = 0;
$mask_bit_idx[1] = 1;
my $bit_idx = 1;
foreach my $i (0..4){
$func->($primes_up_to_2310[$i]);
}
my $idx = 5;
foreach my $n (2..2309){
$idx += 1 while $primes_up_to_2310[$idx] < $n;
if( $n == $primes_up_to_2310[$idx] ){
$mask_bit_idx[$n] = $bit_idx += 1;
}
elsif( $n % 2 == 0 || $n % 3 == 0 || $n % 5 == 0 || $n % 7 == 0 || $n % 11 == 0 ){
$mask_bit_idx[$n] = 0;
}
else {
$mask_bit_idx[$n] = $bit_idx += 1;
}
}
for( my $n = 13; $n <= $upto; $n += 2 ){
if( my $mask_bit = $mask_bit_idx[$n % 2310] ){
my $byte_idx = int($n / 2310) * 60 + ($mask_bit - 1 >> 3); $bit_idx = 1 << ($mask_bit & 7);
if( $flags[$byte_idx] & $bit_idx ){
$func->($n);
if( $n < $idx_limit ){
$idx = $n * $n;
$idx += $n if !($idx & 1);
while( $idx <= $upto ){
if( $mask_bit = $mask_bit_idx[$idx % 2310] ){
$byte_idx = int($idx / 2310) * 60 + ($mask_bit-1 >> 3); $mask_bit = 255 - (1 << ($mask_bit & 7));
$flags[$byte_idx] = $flags[$byte_idx] & $mask_bit;
}
$idx += 2 * $n;
}
}
}
}
}
}
my $max = @ARGV == 1 ? $ARGV[0] : 100;
my $q = 0;
my $c = 0;
primes($max, sub{ $q = $_[0]; $c += 1; });
c3ViIHByaW1lcyB7CiAgICBteSAoJHVwdG8sICRmdW5jKSA9IEBfOwoKICAgIHJldHVybiBsYXJnZV9wcmltZXMoJHVwdG8sICRmdW5jKSBpZiAkdXB0byA+IDI1MDAwOwogICAgJHVwdG8gLT0gMTsKICAgIG15IEBmbGFncyA9ICgxKSB4ICR1cHRvOwogICAgZm9yZWFjaCBteSAkaSAoMC4uJHVwdG8gLSAxKXsKICAgICAgICBpZiAoJGZsYWdzWyRpXSl7CiAgICAgICAgICAgIG15ICRwcmltZSA9ICRpICsgMjsKICAgICAgICAgICAgbXkgJGsgPSAkaSArICRwcmltZTsKICAgICAgICAgICAgd2hpbGUoJGsgPCAkdXB0byl7CiAgICAgICAgICAgICAgICAkZmxhZ3NbJGtdID0gMDsKICAgICAgICAgICAgICAgICRrICs9ICRwcmltZTsKICAgICAgICAgICAgfQogICAgICAgICAgICAkZnVuYy0+KCRwcmltZSk7CiAgICAgICAgfQogICAgfQp9CgpzdWIgbGFyZ2VfcHJpbWVzIHsKICAgIG15ICgkdXB0bywgJGZ1bmMpID0gQF87CgogICAgIyBUaGUgQWxnb3JpdGhtIGlzIGFkYXB0ZWQgZnJvbSBodHRwOi8vdy4uLmNvbnRlbnQtYXZhaWxhYmxlLXRvLWF1dGhvci1vbmx5Li4uay5jb20vfmpybS9wcmludHByaW1lcy5odG1sCiAgICAjIFRoaXMgY29kZSBpcyBhIHRyYW5zbGF0ZWQgdmVyc2lvbiBvZiBBbmRyZWFzIFJhYWIncyAjbGFyZ2VQcmltZXNVcFRvOmRvOiBtZXRob2QKICAgICMgd2hpY2ggd2FzIHdyaXR0ZW4gZm9yIHRoZSBTcXVlYWsgU21hbGx0YWxrIGVudmlyb25tZW50LgoKICAgIG15ICRpZHhfbGltaXQgPSBpbnQoc3FydCgkdXB0bykpICsgMTsKICAgIG15IEBmbGFncyA9ICgweEZGKSB4IChpbnQoKCR1cHRvICsgMjMwOSkgLyAyMzEwKSAqIDYwICsgNjApOwogICAgbXkgQHByaW1lc191cF90b18yMzEwID0gKCk7CiAgICBwcmltZXMoMjMxMCwgc3VieyBwdXNoKEBwcmltZXNfdXBfdG9fMjMxMCwgJF9bMF0pIH0pOwogICAgbXkgQG1hc2tfYml0X2lkeDsKICAgICQjbWFza19iaXRfaWR4ID0gMjMwOTsKICAgICRtYXNrX2JpdF9pZHhbMF0gPSAwOwogICAgJG1hc2tfYml0X2lkeFsxXSA9IDE7CiAgICBteSAkYml0X2lkeCA9IDE7CiAgICBmb3JlYWNoIG15ICRpICgwLi40KXsKICAgICAgICAkZnVuYy0+KCRwcmltZXNfdXBfdG9fMjMxMFskaV0pOwogICAgfQogICAgbXkgJGlkeCA9IDU7CiAgICBmb3JlYWNoIG15ICRuICgyLi4yMzA5KXsKICAgICAgICAkaWR4ICs9IDEgd2hpbGUgJHByaW1lc191cF90b18yMzEwWyRpZHhdIDwgJG47CiAgICAgICAgaWYoICRuID09ICRwcmltZXNfdXBfdG9fMjMxMFskaWR4XSApewogICAgICAgICAgICAkbWFza19iaXRfaWR4WyRuXSA9ICRiaXRfaWR4ICs9IDE7CiAgICAgICAgfQogICAgICAgIGVsc2lmKCAkbiAlIDIgPT0gMCB8fCAkbiAlIDMgPT0gMCB8fCAkbiAlIDUgPT0gMCB8fCAkbiAlIDcgPT0gMCB8fCAkbiAlIDExID09IDAgKXsKICAgICAgICAgICAgJG1hc2tfYml0X2lkeFskbl0gPSAwOwogICAgICAgIH0KICAgICAgICBlbHNlIHsKICAgICAgICAgICAgJG1hc2tfYml0X2lkeFskbl0gPSAkYml0X2lkeCArPSAxOwogICAgICAgIH0KICAgIH0KICAgIGZvciggbXkgJG4gPSAxMzsgJG4gPD0gJHVwdG87ICRuICs9IDIgKXsKICAgICAgICBpZiggbXkgJG1hc2tfYml0ID0gJG1hc2tfYml0X2lkeFskbiAlIDIzMTBdICl7CiAgICAgICAgICAgIG15ICRieXRlX2lkeCA9IGludCgkbiAvIDIzMTApICogNjAgKyAoJG1hc2tfYml0IC0gMSA+PiAzKTsKICAgICAgICAgICAgJGJpdF9pZHggPSAxIDw8ICgkbWFza19iaXQgJiA3KTsKICAgICAgICAgICAgaWYoICRmbGFnc1skYnl0ZV9pZHhdICYgJGJpdF9pZHggKXsKICAgICAgICAgICAgICAgICRmdW5jLT4oJG4pOwogICAgICAgICAgICAgICAgaWYoICRuIDwgJGlkeF9saW1pdCApewogICAgICAgICAgICAgICAgICAgICRpZHggPSAkbiAqICRuOwogICAgICAgICAgICAgICAgICAgICRpZHggKz0gJG4gaWYgISgkaWR4ICYgMSk7CiAgICAgICAgICAgICAgICAgICAgd2hpbGUoICRpZHggPD0gJHVwdG8gKXsKICAgICAgICAgICAgICAgICAgICAgICAgaWYoICRtYXNrX2JpdCA9ICRtYXNrX2JpdF9pZHhbJGlkeCAlIDIzMTBdICl7CiAgICAgICAgICAgICAgICAgICAgICAgICAgICAkYnl0ZV9pZHggPSBpbnQoJGlkeCAvIDIzMTApICogNjAgKyAoJG1hc2tfYml0LTEgPj4gMyk7CiAgICAgICAgICAgICAgICAgICAgICAgICAgICAkbWFza19iaXQgPSAyNTUgLSAoMSA8PCAoJG1hc2tfYml0ICYgNykpOwogICAgICAgICAgICAgICAgICAgICAgICAgICAgJGZsYWdzWyRieXRlX2lkeF0gPSAkZmxhZ3NbJGJ5dGVfaWR4XSAmICRtYXNrX2JpdDsKICAgICAgICAgICAgICAgICAgICAgICAgfQogICAgICAgICAgICAgICAgICAgICAgICAkaWR4ICs9IDIgKiAkbjsKICAgICAgICAgICAgICAgICAgICB9CiAgICAgICAgICAgICAgICB9CiAgICAgICAgICAgIH0KICAgICAgICB9CiAgICB9Cn0KCm15ICRtYXggPSBAQVJHViA9PSAxID8gJEFSR1ZbMF0gOiAxMDA7Cm15ICRxID0gMDsKbXkgJGMgPSAwOwpwcmltZXMoJG1heCwgc3VieyAkcSA9ICRfWzBdOyAkYyArPSAxOyB9KTsKcHJpbnQgJHEsICIgIiwgJGMsICJcbiI7