#!/opt/bin/perl

use Math::BigInt try => "GMP";

$| = 1;

# http://miller-rabin.appspot.com/ (slightly related: http://www.mersenneforum.org/showthread.php?t=12209)
# https://arxiv.org/abs/1509.00864
my @sprp = (
   [                   "341531", "9345883071009581737"],
   [                "520924141", "15, 750068417525532"],
   [             "109134866497", "2, 45650740, 3722628058"],
   [           "47636622961201", "2, 2570940, 211991001, 3749873356"],
   [         "3770579582154547", "2, 2570940, 880937, 610386380, 4130785767"],
   [     "18446744073709551615", "2, 325, 9375, 28178, 450775, 9780504, 1795265022"],
   [ "318665857834031151167461", "2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37"],
   ["3317044064679887385961981", "2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41"],
);

$_ = [map +(new Math::BigInt $_), $_->[0], split /,\s*/, $_->[1]]
   for @sprp;

my @bases;

# requires @bases to be set properly before call
sub isprime($) {
   my ($n) = @_;

   $n % $_ || $n == $_ || return
      for 3, 5, 7; # 2 not tested

   return 1 if $n < 121;

   my ($d, $s) = $n - 1;
   $s++ until $d & (1<<$s);
   $d >>= $s;

   a:
   for my $a (@bases) {
      my $x = $a->copy->bmodpow ($d, $n);
      next if $x == 1 || $x == $n - 1;
      for (1 .. $s - 1) {
         $x = $x * $x % $n;
         return if $x == 1; # warn "$n == 0 mod ", Math::BigInt::bgcd ($a ** ($d * 2 ** $_ / 2) - 1, $n);
         next a if $x == $n - 1;
      }
      return;
   }

   1
}

my $n = (new Math::BigInt shift) | 1;

my $max = 0;

while () {
   while ($n > $max) {
      @sprp
         or die "range exceeded\n";
      ($max, @bases) = @{ shift @sprp };
   }
   print "$n\n" if isprime $n;
   $n += 2;
}