#! /usr/bin/perl
###########################
#
# a little piece of perl code to find a prime
# number using Miller Rabin primality test
# and Math::BigInt library
#
# switcher - may 2007
###########################
# please install libmath-bigint-gmp-perl
###########################

use strict;
use Math::BigInt lib=>"GMP";
use Math::BigInt::Random qw/ random_bigint /;

use Time::HiRes;

print "\n\nPrime Number Tool - switcher, may 2007\n",
	"\t---------------------\n",
	"what do you want to do ?\n",
	"1. generate a prime number\n",
	"2. test the primality of a number\n\n",
	"choice > ";

my $choice=<STDIN>;

if ($choice==1) {
	my $keylength;
	do{
		print "size of the prime number (in bits) : ";
		$keylength=<STDIN>;
	}while($keylength !~ /\d+/);

	print "please wait.....\n";

	# benchmark clock
	my $globalTime = [Time::HiRes::gettimeofday()];
	my $bign = new Math::BigInt('0');
	my $state="caribou";
	do {
		# generate a random number
		$bign = random_bigint( length => $keylength, length_bin => 1);
		if(division($bign) eq "prime"){
			$state=miller($bign);
		}
	#repeat while we find a prime
	}while($state ne "prime");
	
	#benchmark clock
	my $global = Time::HiRes::tv_interval($globalTime); 
	print "prime number = $bign\nfound in $global\n";
}

if ($choice==2) {
	print "enter the number you want to test :\n> ";
	my $input=<STDIN>; chomp($input);
	my $bign = new Math::BigInt('1');
	$bign->bmul($input);

	print "please wait.....\n";
	my $state = miller($bign);
	if( $state eq "prime"){print "This number is prime ! \n";}
	else{print "This number seems not to be prime, $state\n";}
}


sub division{
	my $bign = shift;
	my @firstprime = qw(
	2 3 5 7 11 13 17 19 23 29 31 37 41 43 47 53 59 61 67 71 73 79 83 89 97 101
	103 107 109 113 127 131 137 139 149 151 157 163 167 173 179 181 191 193  197
	199 211 223 227 229 233 239 241 251 257 263 269 271 277 281 283 293 307  311
	);
		
	foreach my $prime (@firstprime){
		my $test = new Math::BigInt($bign);
		if ($test->bmod($prime) == 0){
			return "composed (division detection)";
		}
	}
	return "prime";

}

####################
# Rabin Miller primality check#
####################
# source : 
# http://en.wikipedia.org/wiki/Rabin-Miller_primality_test#Algorithm_and_running_time
####################
sub miller(){

	# get argument	
	my $bign = new Math::BigInt(shift);

	# define ($bign-1) in 2 variables
	my $d = new Math::BigInt($bign-1);

	# init "s" an integer
	my $s = 0;

	#Find d and a satisfying: n-1 = (2^s) * d {with d odd}
	while( !($d & 1) ){
		$d >>= 1;
		$s++;
	}
	
	# compute 100 tests
	for my $i (0..200){
		my $a = new Math::BigInt('0');
		# use a random number between 0 and bign-1 with pgcd a and bign == 1
		do{ $a = random_bigint(max => $bign-1);} 
		while ( Math::BigInt::bgcd($a,$bign-1) != 1 );

		# test 1 : if a^d mod n !=1, then do test2
		$a->bmodpow($d,$bign);
		if ($a != 1 and $a != $bign-1) {
			my $j=1;
			# test 2 : (quite explicit.... don't need any comment)
			while ($a<$s-1 and $a != $bign -1){
				$a->bmodpow(2,$bign); # compute y^2 mod bign
				if ($a == 1){ return "composed";} # return "composed", bign is not prime
				$j++;
			}
			if ( $a != $bign-1){ return "composed";} # return "composed", bign is not prime
		}
	}
	# all tests passed without a "composed" return, so $bign is strongly probably prime
	return "prime";
}