File Coverage

blib/lib/Math/Factoring.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::Factoring;
{
  $Math::Factoring::VERSION = '0.02';
}

use warnings;
use strict;
use Math::GMPz qw/:mpz/;
use Math::Primality qw/is_prime/;
use base 'Exporter';
use constant GMP => 'Math::GMPz';
our @EXPORT_OK = qw/factor factor_trial/;
our @EXPORT = qw//;
use Data::Dumper;

# this gets rid of prototype warnings
sub factor_trial($);

my @small_primes = qw/
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
1p99     211     223     227     229     233     239   241     251     257
/;
my %small_primes = map { $_ => 1 } @small_primes;

# ABSTRACT: Math::Factoring - Advanced Factoring Algorithms


sub _random()
{
    my $n   = GMP->new(int rand(1e15) );
    my $state = rand_init($n);
    my $rand = GMP->new;
    Rmpz_urandomm($rand, $state, $n,1);
    #warn "return rand=$rand";
    return $rand;
}

# this is fast but only works for certain numbers
# which satisfy smoothness constraints that are currently not being
# checked for
sub _factor_pollard_rho($$$)
{
    my ($n,$a,$x0) = @_;
    warn "_factor_pollard($n,$a,$x0)\n";
    my ($x,$y,$q,$d) = map { GMP->new } ( 1 .. 4 );
    my ($i,$j) = (1,1);
    $q = 1; $x = $x0; $y = $x0;

    do {
        $x  = ($x*$x + $a ) % $n;
        $y  = ($y*$y + $a ) % $n;
        $y  = ($y*$y + $a ) % $n;
        $q *= ($x - $y);
        $q %= $n;

        $i++;
        $j = 1 if !$j;
        if( ($i % $j ) == 0 ) {
            $j++;
            Rmpz_gcd($d, $q, $n);
            if ($d != 1) {
                if (!is_prime($d)) {
                    no warnings 'prototype';
                    return _factor_pollard_rho( $d,
                            (_random() & 32) - 16,
                             _random() & 31 );
                } else {
                    return $d;
                }
            }
        }
    };
    return 0;

}


sub factor($)
{
    my ($n) = @_;
    if ($n <= 257 && $small_primes{$n} ){
        return ("$n");
    } else {
        return factor_trial($n);
    }
}

sub factor_trial($)
{
    my $n = shift;
    if ($n >= 0 and $n <= 3) {
        return ("$n");
    }
    $n   = GMP->new($n);
    my $sqrt = GMP->new;
    Rmpz_sqrt($sqrt, $n);

    # speed up factors of perfect squares

    if( Rmpz_perfect_square_p($n) ){
        my @root_factors = factor_trial($sqrt);
        return map { ("$_","$_") } @root_factors;
    }

    my @factors;
    my $cur = GMP->new(2);
    my ($mod,$square) = (GMP->new,GMP->new);
    Rmpz_mul($square,$cur,$cur);

    while( $square <= $n ) {
        Rmpz_mod($mod,$n,$cur);
        if( Rmpz_cmp_ui($mod,0) == 0 ) {
            push @factors,"$cur";
            Rmpz_tdiv_q($n,$n,$cur);  # $n = $n / $cur;
        } else {
            Rmpz_add_ui($cur,$cur,1); # $cur++
        }
        Rmpz_mul($square,$cur,$cur);
    }
    if (@factors == 0) {
        return ("$n");                # it was prime
    }
    if ( Rmpz_cmp_ui($n,1) ) {
        push @factors,"$n";           # add the last factor
    }
    return sort { $a <=> $b } @factors;
}


sub factor_pollard_rho($)
{
    my $n   = GMP->new($_[0]);
    my @factors;
    if ($n >= 0 and $n <= 3) {
        return "$n";
    } else {
        my ($a,$x0) = (1,3);
        my $t;
        while( !is_prime($n) ) {
            $t = _factor_pollard_rho($n,$a,$x0);
            warn "found t=$t,n=$n";
            last if $t == 0;
            push @factors, "$t";
            $n /= $t;
        }
        push @factors, "$n";
    }

    return sort { $a <=> $b } @factors;
}

1; # End of Math::Factoring

__END__