#!perl

use warnings;
use strict;


use Inline "C", <<__ENDC;

/* mulmod(a, b, m) gives (a*b)%c */
long mulmod (long a, long b, long m) {
 long long t = (long long)a * b;
 return (long)(t % m);
}

__ENDC


{
# Miller--Rabin test

my @smallprimes = (
 2, 3, 5, 7, 11, 13, 17, 19, 
);


##</code><code>##

# primep($x) returns whether $x is a prime
sub primep_rm {
 my($cand, $base, $p, $iter, $exp, $odd, $t, $x, $xprev);
 $cand = 0 + $_[0];
 1 < $cand && $cand == int($cand) or
  die "error: invalid number to primep_rm: $cand";
 for $p (@smallprimes) {
  0 != $cand % $p or
   return $p == $cand;
 }
 ($exp, $odd) = (0, $cand - 1);
 while (0 == $odd % 2) {
  ($exp, $odd) = ($exp + 1, $odd >> 1);
 }
 for $iter (1 .. 52) {
  $base = 1 + int(rand($cand - 1));
  $x = powmod($base, $odd, $cand);
  SQUARING: for $t (1 .. $exp) {
   ($xprev, $x) = ($x, mulmod($x, $x, $cand));
   $x == 1 and do {
    $xprev == 1 || $xprev == $cand - 1 or
     return;
    last;
   }
  }
  # KEY STEP: $x now contains $base**($cand - 1) % $cand
  # which should be 1 because of the Fermat theorem.
  $x == 1 or
   return;
 }
 return 1;
}

##</code><code>##


# powmod($a, $f, $m) calculates $a**$f % $m,
# precondition: 1 < $m, 0 < $f, $a are integers representable as a perl integer
sub powmod {
 my($a, $f, $m) = @_;
 my $r = 1;
 for (;;) {
  0 != ($f & 1) and
   $r = mulmod($r, $a, $m);
  $f >>= 1;
  0 == $f and
   last;
  $a = mulmod($a, $a, $m);
 }
 $r;
}

}


{
# simple factoring

my @primes;
{
 my($n, $p, $sqrp);
 NLOOP: for $n (2 .. 100_100) {
  $sqrp = sqrt($n) + 0.5;
  for $p (@primes) {
   $sqrp < $p and last; 
   0 == $n % $p and next NLOOP;
  }
  push @primes, $n
 } 
}
my $greatest_prime_squared = $primes[-1]**2;


##</code><code>##

sub primep_div {
 my $cand = 0 + $_[0];
 my($p, $sqr);
 1 < $cand && $cand == int($cand) or
  die "error: invalid number to primep_div: $cand";
  $cand < $greatest_prime_squared or
  die "error: number too large in primep_div: $cand";
 $sqr = sqrt($cand) + 0.5;
 for $p (@primes) {
  $p <= $sqr or
   return 1;
  0 == $cand % $p and
   return $p == $cand;
 }
 die "internal error: shouldn't fall through here";
}

##</code><code>##


}


if (1) {
 my $test = sub {
  my($n, $p_rm, $p_div);
  $n = $_[0];
  $p_rm = !!primep_rm($n);
  $p_div = !!primep_div($n);
  $p_rm == $p_div or
   die "error: primality tests don't match: $n";
  
 };
 my($n, $r, $l);
 warn "begin serial tests";
 for $n (2 .. 100000) {
  &$test($n);
 }
 warn "begin random tests";
 for $l (3 .. 9) {
  for $r (1 .. 10000) {
   $n = 2 + int(rand(10**$l));
   &$test($n);
  }
 }
 warn "all ok";
}


use Benchmark ":all";

if (1) {
 warn "start benchmarks";
 for my $len (3 .. 9) {
  my($r_rm, $r_div);
  my $ntries = $len < 8 ? 50_000 : 20_000;
  $ntries /= 5;
  my @tries = map {
   10**($len - 1) + 1 + 2 * int(rand((10**$len - 10**($len - 1))/2));
   } 0 .. $ntries;
  warn "benchmarking $len digit odd numbers, $ntries in total";
  cmpthese(1, {
   "primep_rm", sub { 
    my $r; 
    $r += primep_rm($_) ? 1 : 0 for
     @tries;
    $r_rm = $r;
   },
   "primep_div", sub { 
    my $r; 
    $r += primep_div($_) ? 1 : 0 for
     @tries;
    $r_div = $r;
   },
  });
  $r_rm == $r_div or
   die "error: primality tests don't agree in total";
 }
}


__END__

##</code><code>##

begin serial tests at p line 125.
begin random tests at p line 129.
all ok at p line 136.
start benchmarks at p line 143.
benchmarking 3 digit odd numbers, 10000 in total at p line 151.
            (warning: too few iterations for a reliable count)
            (warning: too few iterations for a reliable count)
benchmarking 4 digit odd numbers, 10000 in total at p line 151.
           s/iter  primep_rm primep_div
primep_rm    11.2         --       -97%
primep_div  0.300      3637%         --
            (warning: too few iterations for a reliable count)
            (warning: too few iterations for a reliable count)
benchmarking 5 digit odd numbers, 10000 in total at p line 151.
           s/iter  primep_rm primep_div
primep_rm    10.0         --       -96%
primep_div  0.370      2614%         --
            (warning: too few iterations for a reliable count)
            (warning: too few iterations for a reliable count)
benchmarking 6 digit odd numbers, 10000 in total at p line 151.
           s/iter  primep_rm primep_div
primep_rm    9.63         --       -95%
primep_div  0.490      1865%         --
            (warning: too few iterations for a reliable count)
            (warning: too few iterations for a reliable count)
benchmarking 7 digit odd numbers, 10000 in total at p line 151.
           s/iter  primep_rm primep_div
primep_rm    9.04         --       -92%
primep_div  0.740      1122%         --
            (warning: too few iterations for a reliable count)
            (warning: too few iterations for a reliable count)
benchmarking 8 digit odd numbers, 4000 in total at p line 151.
           s/iter  primep_rm primep_div
primep_rm    8.87         --       -85%
primep_div   1.34       562%         --
            (warning: too few iterations for a reliable count)
            (warning: too few iterations for a reliable count)
benchmarking 9 digit odd numbers, 4000 in total at p line 151.
           s/iter  primep_rm primep_div
primep_rm    3.36         --       -68%
primep_div   1.08       211%         --
            (warning: too few iterations for a reliable count)
            (warning: too few iterations for a reliable count)
           s/iter  primep_rm primep_div
primep_rm    3.41         --       -30%
primep_div   2.37        44%         --