/
Factoring.pm
223 lines (169 loc) · 5.14 KB
/
Factoring.pm
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
package Math::Factoring;
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
=head1 SYNOPSIS
use Math::Factoring;
my $n = 42;
my @factors = factor(42); # 2 3 7
=head1 AUTHOR
Jonathan "Duke" Leto, C<< <jonathan at leto.net> >>
=head1 BUGS
Please report any bugs or feature requests to C<bug-math-factoring at rt.cpan.org>, or through
the web interface at L<http://rt.cpan.org/NoAuth/ReportBug.html?Queue=Math::Factoring>. I will be notified, and then you'll
automatically be notified of progress on your bug as I make changes.
You can also submit patches or file a Github issue:
https://github.com/leto/Math-Factoring
=head1 SUPPORT
You can find documentation for this module with the perldoc command.
perldoc Math::Factoring
You can also look for information at:
=over 4
=item * RT: CPAN's request tracker
L<http://rt.cpan.org/NoAuth/Bugs.html?Dist=Math::Factoring>
=item * AnnoCPAN: Annotated CPAN documentation
L<http://annocpan.org/dist/Math::Factoring>
=item * CPAN Ratings
L<http://cpanratings.perl.org/d/Math::Factoring>
=item * Search CPAN
L<http://search.cpan.org/dist/Math::Factoring>
=back
=head1 ACKNOWLEDGEMENTS
=head1 COPYRIGHT & LICENSE
Copyright 2009-2012 Jonathan "Duke" Leto, all rights reserved.
This program is free software; you can redistribute it and/or modify it
under the same terms as Perl itself.
=cut
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;
}
=head2 factor($number)
Returns the prime factors of $number as an array. Currently this falls back to trial division.
The values in the array are plain Perl string variables, not Math::GMPz objects.
=cut
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;
}
=head2 factor_pollard_rho($number)
Return the factors of $number as an array using the Pollard-Rho algorithm.
=cut
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