use warnings;
use strict;

use Getopt::Long;
use Math::BigInt;

my %pars = ( verify => 0 );
GetOptions ( \%pars, 'verify', 'help', 'usage' )
    or die "Bad GetOptions";

die "primes number [ more numbers ] \nreturns the number's primality status"
    if $pars{help} || $pars{usage};

my @inputs = @ARGV 
    or die "no numbers given for primality test\n";

foreach ( @inputs ) {

    my $number = Math::BigInt->new ( $_ );
    my ( $primality_prob, $witness ) = miller_rabin_prime ( $number );

    print "$number is ",
          ( $primality_prob )
              ? "prime with prob $primality_prob\n"
              : "not prime, witness $witness\n";
}

sub miller_rabin_prime {

    my $n       = shift () or die "no n in MillerRabin";
    my $s       = shift () || 25;

    my $a = Math::BigInt->new ( 0 );

    for ( my $j = 1; $j <= $s; ++$j ) {

        # $a is a random int in [1,n-1]
        #
        my $ranlim = $n-1;
        do {
            $a = Math::BigInt->new ( sprintf "%d", rand ( $ranlim ) + 1 );
            $ranlim /= 2;
        } while ( 0 > $a );

        if ( 1 == Witness ( $a, $n ) ) {
            return 0, $a;
        }
    }
    my $prob = 1 - (3/4)**$s;
    return $prob, undef;
}

sub Witness {

    my $a = shift () or die "no a in Witness";
    my $n = shift () or die "no n in Witness";

    # n-1 = ( 2**t ) * u
    #
    my ( $t, $u ) = get_t_u ( $n - 1 );

    my @x;
    $x[0] = mod_exp ( $a, $u, $n );

    for ( my $i = 1; $i <= $t; ++$i ) {
        $x[$i] = mod_exp ( $x[$i-1], 2, $n );
        return 1
            if (1 == $x[$i]) && (1 != $x[$i-1]) && (( $n - 1 ) != $x[$i-1]);
    }
    return 1 
        if 1 != $x[$t];

    return 0;
}

sub get_t_u {

    my $m = shift () or die "no m in get_t_u";

    my ( $t, $u ) = ( 0, 0 );

    while ( 0 == $m % 2 ) {
        ++$t;
        $m /= 2;
    }
    $u = $m;
    return ( $t, $u );
}

sub mod_exp {

    # i**j mod n
    #
    # works with Math::BigInt in addition to regular scalars
    #
    my $i = shift () or die "no i in mod_exp";
    my $j = shift () or die "no j in mod_exp";
    my $n = shift () or die "no n in mod_exp";

    # return 1 for zero exponents
    #
    my $result = $i - $i + 1;
    return $result unless $j;

    my $pow2 = $i;

    while (1) {
        if ( $j % 2 ) {
            $result = ( $pow2 * $result ) % $n;
            return $result unless --$j;
        }
        $j /= 2;
        $pow2 = ( $pow2 * $pow2 ) % $n;
    }
}