Skip to content

Commit

Permalink
Convert to use Math::GMPz
Browse files Browse the repository at this point in the history
  • Loading branch information
leto committed May 11, 2009
1 parent 62ca66a commit 8a752f8
Show file tree
Hide file tree
Showing 3 changed files with 40 additions and 31 deletions.
2 changes: 1 addition & 1 deletion Build.PL
Expand Up @@ -9,7 +9,7 @@ my $builder = Module::Build->new(
dist_version_from => 'lib/Math/Primality.pm',
build_requires => {
'Test::More' => 0,
'Math::BigInt::GMP' => 1.25,
'Math::GMPz' => 0.26,
},
add_to_cleanup => [ 'Math::Primality-*' ],
create_makefile_pl => 'traditional',
Expand Down
1 change: 0 additions & 1 deletion MANIFEST
Expand Up @@ -4,5 +4,4 @@ MANIFEST
README
lib/Math/Primality.pm
t/00-load.t
t/pod-coverage.t
t/pod.t
68 changes: 39 additions & 29 deletions lib/Math/Primality.pm
Expand Up @@ -2,13 +2,11 @@ package Math::Primality;
use warnings;
use strict;
use Data::Dumper;
use Math::BigInt;
use Math::BigInt qw/bgcd/;
use Math::BigInt::GMP;
use Math::GMPz qw/:mpz/;
use base 'Exporter';
our $DEBUG = 0;

use constant GMP => 'Math::BigInt::GMP';
use constant GMP => 'Math::GMPz';

=head1 NAME
Expand Down Expand Up @@ -61,19 +59,26 @@ sub debug {
sub is_pseudoprime
{
my ($n, $base) = @_;
# force to BigInts for now
return 0 unless $n;
$base ||= 2;
$base = Math::BigInt->new("$base");
$n = Math::BigInt->new("$n");
$base = GMP->new("$base");
$n = GMP->new("$n");

# if $n and $base are not coprime, than $base is a factor of $n
# $base > 2 && ( Math::BigInt::bgcd($n,$base) != 1 ) && return 0;

my $m = $n->copy->bdec; # m = n - 1
my $mod = $base->copy;
$mod->bmodpow($m,$n); # (base**exp) (mod n)
return ! $mod->bcmp(1); # pseudoprime if $mod = 1
my $m = _copy($n);
Rmpz_sub_ui($m, $m, 1); # m = n - 1

my $mod = _copy($base);
Rmpz_powm($mod, $base, $m, $n );
return ! Rmpz_cmp_ui($mod, 1); # pseudoprime if $mod = 1
}

sub _copy
{
my ($n) = @_;
return GMP->new("$n");
}

=head2 is_strong_pseudoprime($n,$b)
Expand All @@ -87,51 +92,57 @@ sub is_strong_pseudoprime
{
my ($n, $base) = @_;

# force to BigInts for now
$base ||= 2;
$base = GMP->_new("$base");
$n = GMP->_new("$n");
$base = GMP->new("$base");
$n = GMP->new("$n");

my $cmp = GMP->_cmp_ui($n, 2 );
my $cmp = Rmpz_cmp_ui($n, 2 );
return 1 if $cmp == 0;
return 0 if $cmp < 0;

# unnecessary but faster if $n is even
return 0 if GMP->_is_even($n);
return 0 if Rmpz_even_p($n);

my $m = GMP->_copy($n);
$m = GMP->_dec($m);
my $s = GMP->scan1($m,0);
my $d = $n->_new(0);
my $m = _copy($n);
Rmpz_sub_ui($m,$m,1);

$d = GMP->tdiv_q_2exp($d, $m,$s);
my $s = Rmpz_scan1($m,0);
my $d = GMP->new(0);

Rmpz_tdiv_q_2exp($d, $m,$s);
debug "m=$m, s=$s, d=$d";

my $residue = GMP->_modpow($base,$d, $n);
my $residue = GMP->new(0);
Rmpz_powm($residue, $base,$d, $n);
debug "$base^$d % $n = $residue";

# if $base^$d = +-1 (mod $n) , $n is a strong pseudoprime

if ( GMP->_cmp_ui($residue,1) == 0 ) {
$cmp = Rmpz_cmp_ui( $residue,1);

if ( $cmp == 0 ) {
debug "found $n as spsp since $base^$d % $n == $residue == 1\n";
return 1;
}
if ( GMP->_acmp($residue,$m) == 0 ) {
$cmp = Rmpz_cmp($residue,$m);

if ( $cmp == 0 ) {
debug "found $n as spsp since $base^$d % $n == $residue == $m\n";
return 1;
}

map {
# successively square $residue, $n is a strong psuedoprime
# if any of these are congruent to -1 (mod $n)
GMP->_mul($residue,$residue);
Rmpz_mul($residue,$residue,$residue);
debug "$_: r=$residue";

my $mod = GMP->_copy($residue);
GMP->_mod($mod,$n);
my $mod = _copy($residue);
Rmpz_mod($mod, $mod,$n);
debug "$_:$residue % $n = $mod ";
$mod = Rmpz_cmp($mod, $m);

if ( GMP->_acmp($mod, $m) == 0) {
if ($mod == 0) {
debug "$_:$mod == $m => spsp!";
return 1;
}
Expand All @@ -140,7 +151,6 @@ sub is_strong_pseudoprime
return 0;
}


=head1 AUTHOR
Jonathan Leto, C<< <jonathan at leto.net> >>
Expand Down

0 comments on commit 8a752f8

Please sign in to comment.