File Coverage

blib/lib/Math/Primality/AKS.pm

Criterion	Covered	Total	%
statement	7	9	77.7
branch			n/a
condition			n/a
subroutine	3	3	100.0
pod			n/a
total	10	12	83.3

package Math::Primality::AKS;
{
  $Math::Primality::AKS::VERSION = '0.08';
}
use warnings;
use strict;

use Math::GMPz qw/:mpz/;
use Math::Primality::BigPolynomial;

use base 'Exporter';
use Carp qw/croak/;

use constant DEBUG => 0;

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

# ABSTRACT: Check for primality using the AKS (Agrawal-Kayal-Saxena) algorithm


our @EXPORT_OK = qw/is_aks_prime/;

our %EXPORT_TAGS = ( all => \@EXPORT_OK );


sub debug {
    if ( DEBUG or $ENV{DEBUG} ) {
      warn $_[0];
    }
}

sub is_aks_prime($) {
  # http://ece.gmu.edu/courses/ECE746/project/F06_Project_resources/Salembier_Southerington_AKS.pdf
  # http://islab.oregonstate.edu/koc/ece575/04Project2/Halim-Chanleudfa/Report.pdf
  # http://fatphil.org/maths/AKS/
  my $n = GMP->new($_[0]);
  # Step 1 - check if $n = m^d for some m, d
  # if $n is a power then return 0
  if (Rmpz_perfect_power_p($n)) {
    debug "fails step 1 of aks - $n is a perfect power";
    return 0;  # composite
  }
    my $r = Math::GMPz->new(2);
    my $logn = Rmpz_sizeinbase($n, 2);
    my $limit = Math::GMPz->new($logn * $logn);
    Rmpz_mul_ui($limit, $limit, 4);

    # Witness search

    OUTERLOOP: while (Rmpz_cmp($r, $n) == -1) {
    if(Rmpz_divisible_p($n, $r)) {
        debug "$n is divisible by $r\n";
        return 0;
    }

    if(Rmpz_probab_prime_p($n, 5)) {
        my $i = Math::GMPz->new(1);
        my $res = Math::GMPz->new(0);

        INNERLOOP: for ( ; Rmpz_cmp($n, $limit) <= 0; Rmpz_add_ui($i, $i, 1)) {
            Rmpz_powm($res, $n, $i, $r);
                if (Rmpz_cmp_ui($res, 1) == 0) {
                    last OUTERLOOP;
                }
        }

    }
    Rmpz_add_ui($r, $r, 1);
    }
    if (Rmpz_cmp($r, $n) == 0) {
        debug "Found $n is prime while checking for r\n";
        return 1;
    }

    # Polynomial check
    my $a;
    my $sqrtr = Math::GMPz->new($r);

    Rmpz_sqrt($sqrtr, $r);
    my $polylimit = Math::GMPz->new(0);
    Rmpz_add_ui($polylimit, $sqrtr, 1);
    Rmpz_mul_ui($polylimit, $polylimit, $logn);
    Rmpz_mul_ui($polylimit, $polylimit, 2);

    my $intr = Rmpz_get_ui($r);

    for($a = 1; Rmpz_cmp_ui($polylimit, $a) >= 0; $a++) {
        debug "Checking at $a\n";
        my $final_size = Math::GMPz->new(0);
        Rmpz_mod($final_size, $n, $r);
        my $compare = Math::Primality::BigPolynomial->new(Rmpz_get_ui($final_size));
        $compare->setCoef(Math::GMPz->new(1), $final_size);
        $compare->setCoef(Math::GMPz->new($a), 0);
        my $res = Math::Primality::BigPolynomial->new($intr);
        my $base = Math::Primality::BigPolynomial->new(1);
        $base->setCoef(Math::GMPz->new(0), $a);
        $base->setCoef(Math::GMPz->new(1), 1);

        Math::Primality::BigPolynomial::mpz_poly_mod_power($res, $base, $n, $n, $intr);


        if($res->isEqual($compare)) {
            debug "Found not prime at $a\n";
            return 0;
        }
    }
    return 1;
}


exp(0); # End of Math::Primality::AKS

__END__