#!/usr/bin/env perl
use 5.010;
use warnings;
use strict;

use Carp qw/croak/;
use Time::HiRes qw/time/;
use Benchmark qw/cmpthese/;
use Math::Prime::Util qw/is_prime forprimes forcomposites/;
use Math::Primality;
use Math::Prime::XS;
use Math::Pari;

# Other options include:
#   Math::Prime::FastSieve
#   Math::GMPz (some undocumented functions)

# Toby's example
sub isPrime1 (_)
{
        my $num = shift;

        return if $num == 1;  # 1 is not prime

        croak "usage: isPrime(NATURAL NUMBER)"
                unless $num =~ /^[1-9][0-9]*$/;

        for my $div (2 .. sqrt $num)
        {
                return if $num % $div == 0;
        }
        return 1;
}

# Johngg's conversion of Toby's example to a mod-2 wheel
sub tobyinkPlus (_)
   {
       my $toTest = shift;

       return 0 if     $toTest == 1;
       return 1 if     $toTest == 2;
       return 0 unless $toTest %  2;

       my $sqrtLimit = sqrt $toTest;
       for ( my $div = 3; $div <= $sqrtLimit; $div += 2 )
       {
           return 0 unless $toTest % $div;
       }
       return 1;
   }

# Slightly faster version of the trivial method
sub isPrime1a (_)
{
  my($n) = @_;
  croak "isPrime must be a positive integer" unless $n !~ tr/0123456789//c;
  return (0,0,1,1,0,1,0,1,0,0)[$n] if $n < 10;

  for my $div (2..sqrt $n) {
    return 0 unless $n % $div;
  }
  1;
}

# A Mod-2 (i.e. odds-only) trial division
sub isPrime1b (_)
{
  my($n) = @_;
  croak "isPrime must be a positive integer" unless $n !~ tr/0123456789//c;
  return (0,0,1,1,0,1,0,1,0,0)[$n] if $n < 10;

  return unless $n & 1;
  my $limit = int(sqrt($n));
  for (my $div = 3; $div <= $limit; $div += 2) {
    return 0 unless $n % $div;
  }
  1;
}

# Trial division on a mod-6 wheel, so we skip multiples of 2 and 3
sub isPrime1c (_)
{
  my($n) = @_;
  croak "isPrime must be a positive integer" unless $n !~ tr/0123456789//c;
  return (0,0,1,1,0,1,0,1,0,0)[$n] if $n < 10;

  return 0 unless $n & 1;
  return 0 unless $n % 3;
  my $limit = int(sqrt($n));
  for (my $div = 5; $div <= $limit; $div += 6) {
    return 0 unless $n % $div;
    return 0 unless $n % ($div+2);
  }
  1;
}

# Johngg's SoE method
sub isPrime2
{
        my $toTest    = shift;
        my $sqrtLimit = sqrt $toTest;

        my $sieve = q{};
        vec( $sieve, 0      , 1 ) = 1;
        vec( $sieve, 1      , 1 ) = 1;
        vec( $sieve, $toTest, 1 ) = 0;

        my $marker = 1;
        while ( $marker < $sqrtLimit )
        {
                my $possPrime = $marker + 1;
                $possPrime ++ while vec( $sieve, $possPrime, 1 );

                my $fill = 2 * $possPrime;
                while ( $fill <= $toTest )
                {
                        vec( $sieve, $fill, 1 ) = 1;
                        $fill += $possPrime;
                }
                last if vec( $sieve, $toTest, 1 );

                $marker = $possPrime;
        }

        return not vec($sieve, $toTest, 1);
}

# Using better sieve.  Still a bad way to do primality.
sub isPrime2a
{
  my($n) = @_;
  croak "isPrime must be a positive integer" unless $n !~ tr/0123456789//c;
  return (0,0,1,1,0,1,0,1,0,0)[$n] if $n < 10;
  return 0 unless $n & 1;
  my ($sieve, $i, $limit, $s_end) = ( '', 3, int(sqrt($n)), $n >> 1 );
  while ( $i <= $limit ) {
    for (my $s = ($i*$i) >> 1; $s <= $s_end; $s += $i) {
      vec($sieve, $s, 1) = 1;
    }
    do { $i += 2 } while vec($sieve, $i >> 1, 1) != 0;
  }
  vec($sieve, $n >> 1, 1) ? 0 : 1;
}

warn "Validating tests on small primes\n";
forprimes {
  die "isPrime1 fail $_" unless isPrime1($_);
  die "isPrime1a fail $_" unless isPrime1a($_);
  die "isPrime1b fail $_" unless isPrime1b($_);
  die "isPrime1c fail $_" unless isPrime1c($_);
  die "is_prime fail $_" unless is_prime($_);
  die "isPrime2 fail $_" unless $_ > 1000 || isPrime2($_);
  die "isPrime2a fail $_" unless $_ > 1000 || isPrime2a($_);
} 100_000;
warn "Validating tests on small composites\n";
forcomposites {
  die "isPrime1 fail $_" if isPrime1($_);
  die "isPrime1a fail $_" if isPrime1a($_);
  die "isPrime1b fail $_" if isPrime1b($_);
  die "isPrime1c fail $_" if isPrime1c($_);
  die "is_prime fail $_" if $_ < 1000 && is_prime($_);
  die "isPrime2a fail $_" if $_ < 1000 && isPrime2a($_);
} 100_000;

my @small = qw/75 169 229 367 369 8794 9227 10807 11939 14803 19937 19938 39783 47083/;
my @medium = qw/526619 737159 872353 3022457 3932021 10180609 46051073 286914127 305669779 334473319 471623353 806431723 1615845311 1817131787/;
my @large = qw/15242882507 15442725479 33130600963 83197255529 86460333013 136073779057 626171233423 891666642109 1669672223953 1754764041599 3493560177931 3794553266983 43249947787433 112325626292951/;
my @ttest = qw/75 169 229 367 369 8794 9227 10807 11939 14803 19937 19938 39783 47083 199933 75640351 760149189769/;

say q[                   TOBYINK     DIV2        DIV6        MPU];
say q[       NUMBER   :  RES TIME    RES TIME    RES TIME    RES  TIME];
for my $num (qw(75 169 229 367 369 8794 9227 10807 11939 14803 19937 19938 39783 47083 199933 75640351 760149189769 635921898906263))
{
  printf('%15d :', $num);
  timeit($num, $_) for \&isPrime1, \&isPrime1b, \&isPrime1c, \&is_prime;
  print "\n";
}

sub timeit {
  my ($n, $f) = @_;
  my $start = time();
  print($f->($n) ? "  Y " : "  N ");
  my $end = time();
  printf '%.06f', ($end - $start);
}

say "\nSmall inputs (2-5 digits)\n";
cmpthese(-1, {
  "Toby Div1"   => sub { isPrime1($_) for @small },
  "JohnGG Div2" => sub { tobyinkPlus($_) for @small },
  "DJ   Div1"   => sub { isPrime1a($_) for @small },
  "DJ   Div2"   => sub { isPrime1b($_) for @small },
  "DJ   Div6"   => sub { isPrime1c($_) for @small },
  "JohnGG SoE"  => sub { isPrime2($_) for @small },
  "DJ SoE"      => sub { isPrime2a($_) for @small },
  "Module:MPU"  => sub { is_prime($_) for @small },
  "Module:MP"   => sub { Math::Primality::is_prime($_) for @small },
  "Module:MPXS" => sub { Math::Prime::XS::is_prime($_) for @small },
  "Module:Pari" => sub { Math::Pari::isprime($_) for @small },
});

say "\nMedium inputs (6-10 digits)\n";
cmpthese(-2, {
  "Toby Div1"   => sub { isPrime1($_) for @medium },
  "JohnGG Div2" => sub { tobyinkPlus($_) for @medium },
  "DJ   Div1"   => sub { isPrime1a($_) for @medium },
  "DJ   Div2"   => sub { isPrime1b($_) for @medium },
  "DJ   Div6"   => sub { isPrime1c($_) for @medium },
  "JohnGG SoE"  => sub { isPrime2($_) for @small },
  "DJ SoE"      => sub { isPrime2a($_) for @small },
  "Module:MPU"  => sub { is_prime($_) for @medium },
  "Module:MP"   => sub { Math::Primality::is_prime($_) for @medium },
  "Module:MPXS" => sub { Math::Prime::XS::is_prime($_) for @medium },
  "Module:Pari" => sub { Math::Pari::isprime($_) for @medium },
});

say "\nMedium inputs (10-15 digits)\n";
cmpthese(-8, {
  "Toby Div1"   => sub { isPrime1($_) for @large },
  "JohnGG Div2" => sub { tobyinkPlus($_) for @large },
  "DJ   Div1"   => sub { isPrime1a($_) for @large },
  "DJ   Div2"   => sub { isPrime1b($_) for @large },
  "DJ   Div6"   => sub { isPrime1c($_) for @large },
  "Module:MPU"  => sub { is_prime($_) for @large },
  "Module:MP"   => sub { Math::Primality::is_prime($_) for @large },
  "Module:MPXS" => sub { Math::Prime::XS::is_prime($_) for @large },
  "Module:Pari" => sub { Math::Pari::isprime($_) for @large },
});