#!/usr/bin/perl
use Math::BigInt lib => 'GMP';
use Math::PariInit qw/ nextprime /;

=pod

The algorithm can be written as follows:

    Inputs: n: a value to test for primality;
            k: a parameter that determines the accuracy of the test
    Output: composite if n is composite, otherwise probably prime
    write n - 1 as 2^s . d by factoring powers of 2 from n - 1
    repeat k times:

        pick a randomly in the range [1, n - 1]
        if a^d mod n != 1 and a^{2^r.d} mod n != -1 for all r
        in the range [0, s - 1] then return composite

    return probably prime 

=cut

$|++;

my $k = shift() || 25; # accuracy


for (my $n=91001; $n<93001; $n+=2) {
    
    my $d = $n - 1;
    my $s = 0;
    while( ! ($d & 1) ) { # using bit manipulation instead of division by 2 (this is a bit of a premature optimization...)
        $d >>=1;
        $s++;
    }

    my $comp = 1;
    my $witn = 0;

  TEST:
    for( 1 .. $k ) { # test in k bases
        my $a = int(rand($n-1))+1; # select a random base between 1 and n - 1
        my $b = Math::BigInt->new($a);

        if( $b->bmodpow($d, $n) == 1 ) { # a is not a witness for compositeness
            $comp = 0;  
            next TEST ;
        }

        for my $r( 0 .. $s-1 ) {
            if( $b->bmodpow(2**$r * $d, $n) == $n-1 ) { # a is not a witness
                $comp = 0;
                next TEST ;
            }
        }
        $witn++;
    }

    my $p = nextprime($n-1);
    print "$n is probably prime (with $witn witnesses out of $k)\n" unless $comp;
    print "$n is not a prime, however: next prime = $p\n" if !$comp and $n != $p;
    print "$n is a prime, but we missed it\n" if $comp and $n == $p;
}