Codebase list libmath-prime-util-perl / 0fbfc57
Imported Upstream version 0.26 Salvatore Bonaccorso 11 years ago
15 changed file(s) with 2021 addition(s) and 294 deletion(s). Raw diff Collapse all Expand all
+19
-0
Changes less more
00 Revision history for Perl extension Math::Prime::Util.
1
2 0.26 21 April 2013
3
4 - Pure Perl factoring:
5 - real p-1 -- much faster and more effective
6 - Fermat (no better than HOLF)
7 - speedup for pbrent
8 - simple ECM
9 - redo factoring mix
10
11 - New functions:
12 prime_certificate produces a certificate of primality.
13 verify_prime checks a primality certificate.
14
15 - Pure perl primality proof now uses BLS75 instead of Lucas, so some
16 numbers will be much faster [n-1 only needs factoring to (n/2)^1/3].
17
18 - Math::Prime::Util::ECAffinePoint and ECProjectivePoint modules for
19 dealing with elliptic curves.
120
221 0.25 19 March 2013
322
+3
-0
MANIFEST less more
33 lib/Math/Prime/Util/PrimeArray.pm
44 lib/Math/Prime/Util/PP.pm
55 lib/Math/Prime/Util/ZetaBigFloat.pm
6 lib/Math/Prime/Util/ECAffinePoint.pm
7 lib/Math/Prime/Util/ECProjectivePoint.pm
68 LICENSE
79 Makefile.PL
810 MANIFEST
5557 examples/test-factor-gnufactor.pl
5658 examples/test-primes-script.pl
5759 examples/test-primes-script2.pl
60 examples/verify-gmp-eccp-cert.pl
5861 bin/primes.pl
5962 bin/factor.pl
6063 t/01-load.t
+1
-1
META.json less more
5757 "url" : "https://github.com/danaj/Math-Prime-Util"
5858 }
5959 },
60 "version" : "0.25"
60 "version" : "0.26"
6161 }
+1
-1
META.yml less more
3535 resources:
3636 homepage: https://github.com/danaj/Math-Prime-Util
3737 repository: https://github.com/danaj/Math-Prime-Util
38 version: 0.25
38 version: 0.26
+7
-3
README less more
0 Math::Prime::Util version 0.25
0 Math::Prime::Util version 0.26
11
22 A set of utilities related to prime numbers. These include multiple sieving
33 methods, is_prime, prime_count, nth_prime, approximations and bounds for
99 Math::Prime::XS
1010 Math::Prime::FastSieve
1111 Math::Factor::XS
12 Math::Big
1213 Math::Big::Factors
1314 Math::Factoring
1415 Math::Primality
1516 Math::Prime::TiedArray
17 Crypt::Primes
1618 For non-bignums, it is typically faster than Math::Pari (and doesn't
1719 require Pari to be installed).
1820
4850
4951 DEPENDENCIES
5052
51 Perl 5.6.2 or later. No modules outside of Core have been used.
53 Perl 5.6.2 or later (5.8 or later is preferred).
54
55 Bytes::Random::Secure 0.23 or later.
5256
5357
5458 COPYRIGHT AND LICENCE
5559
56 Copyright (C) 2011-2012 by Dana Jacobsen <dana@acm.org>
60 Copyright (C) 2011-2013 by Dana Jacobsen <dana@acm.org>
5761
5862 This library is free software; you can redistribute it and/or modify
5963 it under the same terms as Perl itself.
+5
-2
TODO less more
2020 - Test all functions return either native or bigints. Functions that return
2121 raw MPU::GMP results will return strings, which isn't right.
2222
23 - Make proper pminus1 in PP code, like factor.c.
24
2523 - An assembler version of mulmod for i386 would be _really_ helpful for
2624 all the non-x86-64 Intel machines.
2725
3331 - More efficient totient segment. Do we really need primes to n/2?
3432
3533 - Implement S2 calculation for LMO prime count.
34
35 - Add Sage output style verification. Looks like NZMATH has ECPP but doesn't
36 produce a certificate. Adding a Primo parser (see WraithX's code) would
37 be awesome, but may be a lot more work. It would still be nice to have yet
38 another independent codebase for this.
+27
-5
bin/factor.pl less more
55 use Math::Prime::Util qw/factor nth_prime/;
66 $| = 1;
77 no bigint;
8
9 # Allow execution of any of these functions in the command line
10 my @mpu_funcs = (qw/next_prime prev_prime prime_count nth_prime random_prime
11 random_ndigit_prime random_nbit_prime random_strong_prime
12 random_maurer_prime primorial pn_primorial moebius mertens
13 euler_phi jordan_totient exp_mangoldt divisor_sum
14 consecutive_integer_lcm/);
15 my %mpu_func_map;
816
917 my %opts;
1018 GetOptions(\%opts,
4149 sub eval_expr {
4250 my $expr = shift;
4351 die "$expr cannot be evaluated" if $expr =~ /:/; # Use : for escape
44 $expr =~ s/nth_prime\(/:1(/g;
45 $expr =~ s/log\(/:2(/g;
52 if (scalar(keys %mpu_func_map) == 0) {
53 my $n = 10;
54 foreach my $func (@mpu_funcs) {
55 $mpu_func_map{$func} = sprintf("%03d", $n++);
56 }
57 }
58 $expr =~ s/\blog\(/:001(/g;
59 foreach my $func (@mpu_funcs) {
60 $expr =~ s/\b$func\(/:$mpu_func_map{$func}(/g;
61 }
4662 die "$expr cannot be evaluated" if $expr =~ tr|-0123456789+*/() :||c;
47 $expr =~ s/:1/nth_prime/g;
48 $expr =~ s/:2/log/g;
49 $expr =~ s/(\d+)/ Math::BigInt->new($1) /g;
63 $expr =~ s/:001/log/g;
64 foreach my $func (@mpu_funcs) {
65 $expr =~ s/:$mpu_func_map{$func}\(/Math::Prime::Util::$func(/g;
66 }
67 $expr =~ s/(\d+)/ Math::BigInt->new("$1") /g;
5068 my $res = eval $expr; ## no critic
5169 die "Cannot eval: $expr\n" if !defined $res;
5270 $res = int($res->bstr) if ref($res) eq 'Math::BigInt' && $res <= ~0;
6179 Print the prime factors of each positive integer given on the command line,
6280 or reads numbers from standard input if called without arguments.
6381
82 Math expressions may be given as arguments, which will be evaluated before
83 factoring. This includes most Math::Prime::Util functions including things
84 like prime_count(#), nth_prime(#), primorial(#), random_nbit_prime(#), etc.
85
6486 --help displays this help message
6587 --version displays the version information
6688
+1
-1
examples/test-factor-gnufactor.pl less more
88 use autodie;
99 use Text::Diff;
1010 use Time::HiRes qw(gettimeofday tv_interval);
11 my $maxdigits = 50;
11 my $maxdigits = 100;
1212 $| = 1; # fast pipes
1313 srand(87431);
1414 my $num = 1000;
+61
-0
examples/verify-gmp-eccp-cert.pl less more
0 #!/usr/bin/env perl
1 use warnings;
2 use strict;
3 use Math::BigInt try=>"GMP,Pari";
4 use Math::Prime::Util qw/:all/;
5 use Data::Dump qw/dump/;
6
7 # Takes the output of GMP-ECPP, creates a certificate in the format used
8 # by MPU, and runs it through the verifier.
9 #
10 # Example:
11 #
12 # perl -MMath::Prime::Util -E 'say random_ndigit_prime(60)' | \
13 # gmp-ecpp -q | \
14 # perl examples/verify-gmp-eccp-cert.pl
15
16
17 my $early_check = 0;
18
19 my $N;
20 my ($n, $a, $b, $m, $q, $Px, $Py);
21 my @cert;
22
23 while (<>) {
24 if (/^N\[(\d+)\]\s*=\s*(\d+)/) {
25 $n = $2;
26 if ($1 == 0) {
27 die "Bad input" if defined $N;
28 $N = $n;
29 @cert = ($n, "AGKM");
30 }
31 }
32 elsif (/^a\s*=\s*(\d+)/) { $a = $1; }
33 elsif (/^b\s*=\s*(\d+)/) { $b = $1; }
34 elsif (/^m\s*=\s*(\d+)/) { $m = $1; }
35 elsif (/^q\s*=\s*(\d+)/) { $q = $1; }
36 elsif (/^P\s*=\s*\(\s*(\d+)\s*,\s*(\d+)\s*\)/) {
37 $Px = $1;
38 $Py = $2;
39 die "Bad input\n"
40 unless defined $N && defined $a && defined $b && defined $m
41 && defined $q && defined $Px && defined $Py;
42
43 # If for a given q value, is_prime returns 2, that indicates it can
44 # produce an n-1 primality proof very quickly, so we could stop now.
45 if ($early_check) {
46 my $bq = Math::BigInt->new("$q");
47 if (is_prime($bq) == 2) {
48 push @cert, [$n, $a, $b, $m, [prime_certificate($bq)], [$Px,$Py]];
49 last;
50 }
51 }
52 push @cert, [$n, $a, $b, $m, $q, [$Px,$Py]];
53 }
54 else {
55 last if /^proven prime/;
56 }
57 }
58
59 print dump(\@cert), "\n";
60 print verify_prime(@cert) ? "SUCCESS\n" : "FAILURE\n";
+3
-1
factor.c less more
77 #include "util.h"
88 #include "sieve.h"
99 #include "mulmod.h"
10 #include "cache.h"
1011
1112 /*
1213 * You need to remember to use UV for unsigned and IV for signed types that
622623 /* Pollard's P-1 */
623624 int pminus1_factor(UV n, UV *factors, UV B1, UV B2)
624625 {
625 UV q, f;
626 UV f;
627 UV q = 2;
626628 UV a = 2;
627629 UV savea = 2;
628630 UV saveq = 2;
+236
-0
lib/Math/Prime/Util/ECAffinePoint.pm less more
0 package Math::Prime::Util::ECAffinePoint;
1 use strict;
2 use warnings;
3 use Carp qw/carp croak confess/;
4
5 if (!defined $Math::BigInt::VERSION) {
6 eval { require Math::BigInt; Math::BigInt->import(try=>'GMP,Pari'); 1; }
7 or do { croak "Cannot load Math::BigInt"; };
8 }
9
10 BEGIN {
11 $Math::Prime::Util::ECAffinePoint::AUTHORITY = 'cpan:DANAJ';
12 $Math::Prime::Util::ECAffinePoint::VERSION = '0.26';
13 }
14
15 # Pure perl (with Math::BigInt) manipulation of Elliptic Curves
16 # in affine coordinates. Should be split into a point class.
17
18 sub new {
19 my ($class, $a, $b, $n, $x, $y) = @_;
20 $a = Math::BigInt->new("$a") unless ref($a) eq 'Math::BigInt';
21 $b = Math::BigInt->new("$b") unless ref($b) eq 'Math::BigInt';
22 $n = Math::BigInt->new("$n") unless ref($n) eq 'Math::BigInt';
23 $x = Math::BigInt->new("$x") unless ref($x) eq 'Math::BigInt';
24 $y = Math::BigInt->new("$y") unless ref($y) eq 'Math::BigInt';
25
26 croak "n must be >= 2" unless $n >= 2;
27 $a->bmod($n);
28 $b->bmod($n);
29
30 my $self = {
31 a => $a,
32 b => $b,
33 n => $n,
34 x => $x,
35 y => $y,
36 f => $n->copy->bone,
37 };
38
39 bless $self, $class;
40 return $self;
41 }
42
43 sub _add {
44 my ($self, $P1x, $P1y, $P2x, $P2y) = @_;
45 my $n = $self->{'n'};
46
47 if ($P1x == $P2x) {
48 my $t = ($P1y + $P2y) % $n;
49 return (Math::BigInt->bzero,Math::BigInt->bone) if $t == 0;
50 }
51 my $deltax = ($P2x - $P1x) % $n;
52 $deltax->bmodinv($n);
53 return (Math::BigInt->bzero,Math::BigInt->bone) if $deltax eq "NaN";
54
55 my $deltay = ($P2y - $P1y) % $n;
56 my $m = ($deltay * $deltax) % $n; # m = deltay / deltax
57
58 my $x = ($m*$m - $P1x - $P2x) % $n;
59 my $y = ($m*($P1x - $x) - $P1y) % $n;
60 return ($x,$y);
61 }
62
63 sub _double {
64 my ($self, $P1x, $P1y) = @_;
65 my $n = $self->{'n'};
66
67 my $m = 2*$P1y;
68 $m->bmodinv($n);
69 return (Math::BigInt->bzero,Math::BigInt->bone) if $m eq "NaN";
70
71 $m = ((3*$P1x*$P1x + $self->{'a'}) * $m) % $n;
72
73 my $x = ($m*$m - 2*$P1x) % $n;
74 my $y = ($m*($P1x - $x) - $P1y) % $n;
75 return ($x,$y);
76 }
77
78 sub mul {
79 my ($self, $k) = @_;
80 my $x = $self->{'x'};
81 my $y = $self->{'y'};
82 my $n = $self->{'n'};
83 my $f = $self->{'f'};
84
85 my $Bx = $n->copy->bzero;
86 my $By = $n->copy->bone;
87 my $v = 1;
88
89 while ($k > 0) {
90 if ( ($k % 2) != 0) {
91 $k--;
92 $f->bmul($Bx-$x)->bmod($n);
93 if ($x == 0 && $y == 1) { }
94 elsif ($Bx == 0 && $By == 1) { ($Bx,$By) = ($x,$y); }
95 else { ($Bx,$By) = $self->_add($x,$y,$Bx,$By); }
96 } else {
97 $k >>= 1;
98 $f->bmul(2*$y)->bmod($n);
99 ($x,$y) = $self->_double($x,$y);
100 }
101 }
102 $f = Math::BigInt::bgcd( $f, $n);
103 $self->{'x'} = $Bx;
104 $self->{'y'} = $By;
105 $self->{'f'} = $f;
106 return $self;
107 }
108
109 sub add {
110 my ($self, $other) = @_;
111 croak "add takes a EC point"
112 unless ref($other) eq 'Math::Prime::Util::ECAffinePoint';
113 croak "second point is not on the same curve"
114 unless $self->{'a'} == $other->{'a'} &&
115 $self->{'b'} == $other->{'b'} &&
116 $self->{'n'} == $other->{'n'};
117
118 ($self->{'x'}, $self->{'y'}) = $self->_add($self->{'x'}, $self->{'y'},
119 $other->{'x'}, $other->{'y'});
120 return $self;
121 }
122
123
124 sub a { return shift->{'a'}; }
125 sub b { return shift->{'b'}; }
126 sub n { return shift->{'n'}; }
127 sub x { return shift->{'x'}; }
128 sub y { return shift->{'y'}; }
129 sub f { return shift->{'f'}; }
130
131 sub is_infinity {
132 my $self = shift;
133 return ($self->{'x'}->is_zero() && $self->{'y'}->is_one());
134 }
135
136 1;
137
138 __END__
139
140
141 # ABSTRACT: Elliptic curve operations for affine points
142
143 =pod
144
145 =encoding utf8
146
147
148 =head1 NAME
149
150 Math::Prime::Util::ECAffinePoint - Elliptic curve operations for affine points
151
152
153 =head1 VERSION
154
155 Version 0.26
156
157
158 =head1 SYNOPSIS
159
160 # Create a point on a curve (a,b,n) with coordinates 0,1
161 my $ECP = Math::Prime::Util::ECAffinePoint->new($a, $b, $n, 0, 1);
162
163 # scalar multiplication by k.
164 $ECP->mul($k)
165
166 # add two points on the same curve
167 $ECP->add($ECP2)
168
169 say "P = O" if $ECP->is_infinity();
170
171 =head1 DESCRIPTION
172
173 This really should just be in Math::EllipticCurve.
174
175 To write.
176
177
178 =head1 FUNCTIONS
179
180 =head2 new
181
182 $point = Math::Prime::Util::ECAffinePoint->new(a, b, n, x, y);
183
184 Returns a new point at C<(x,y,1)> on the curve C<(a,b,n)>.
185
186 =head2 a
187
188 =head2 b
189
190 =head2 n
191
192 Returns the C<a>, C<b>, or C<n> values that describe the curve.
193
194 =head2 x
195
196 =head2 y
197
198 Returns the C<x> or C<y> values that define the point on the curve.
199
200 =head2 f
201
202 Returns a possible factor found during EC multiplication.
203
204 =head2 add
205
206 Takes another point on the same curve as an argument and adds it this point.
207
208 =head2 mul
209
210 Takes an integer and performs scalar multiplication of the point.
211
212 =head2 is_infinity
213
214 Returns true if the point is (0,1), which is the point at infinity for
215 the affine coordinates.
216
217
218 =head1 SEE ALSO
219
220 L<Math::EllipticCurve::Prime>
221
222 This really should just be in a L<Math::EllipticCurve> module.
223
224 =head1 AUTHORS
225
226 Dana Jacobsen E<lt>dana@acm.orgE<gt>
227
228
229 =head1 COPYRIGHT
230
231 Copyright 2012-2013 by Dana Jacobsen E<lt>dana@acm.orgE<gt>
232
233 This program is free software; you can redistribute it and/or modify it under the same terms as Perl itself.
234
235 =cut
+297
-0
lib/Math/Prime/Util/ECProjectivePoint.pm less more
0 package Math::Prime::Util::ECProjectivePoint;
1 use strict;
2 use warnings;
3 use Carp qw/carp croak confess/;
4
5 if (!defined $Math::BigInt::VERSION) {
6 eval { require Math::BigInt; Math::BigInt->import(try=>'GMP,Pari'); 1; }
7 or do { croak "Cannot load Math::BigInt"; };
8 }
9
10 BEGIN {
11 $Math::Prime::Util::ECProjectivePoint::AUTHORITY = 'cpan:DANAJ';
12 $Math::Prime::Util::ECProjectivePoint::VERSION = '0.26';
13 }
14
15 # Pure perl (with Math::BigInt) manipulation of Elliptic Curves
16 # in projective coordinates.
17
18 sub new {
19 my ($class, $c, $n, $x, $z) = @_;
20 $c = Math::BigInt->new("$c") unless ref($c) eq 'Math::BigInt';
21 $n = Math::BigInt->new("$n") unless ref($n) eq 'Math::BigInt';
22 $x = Math::BigInt->new("$x") unless ref($x) eq 'Math::BigInt';
23 $z = Math::BigInt->new("$z") unless ref($z) eq 'Math::BigInt';
24
25 croak "n must be >= 2" unless $n >= 2;
26 $c->bmod($n);
27
28 my $self = {
29 c => $c,
30 d => ($c + 2) >> 2,
31 n => $n,
32 x => $x,
33 z => $z,
34 f => $n-$n+1,
35 };
36
37 bless $self, $class;
38 return $self;
39 }
40
41 sub _addx {
42 my ($x1, $x2, $xin, $n) = @_;
43
44 my $u = ($x2 - 1) * ($x1 + 1);
45 my $v = ($x2 + 1) * ($x1 - 1);
46
47 my $upv2 = ($u + $v) ** 2;
48 my $umv2 = ($u - $v) ** 2;
49
50 return ( $upv2 % $n, ($umv2*$xin) % $n );
51 }
52
53 sub _add3 {
54 my ($x1, $z1, $x2, $z2, $xin, $zin, $n) = @_;
55
56 my $u = ($x2 - $z2) * ($x1 + $z1);
57 my $v = ($x2 + $z2) * ($x1 - $z1);
58
59 my $upv2 = ($u + $v) ** 2;
60 my $umv2 = ($u - $v) ** 2;
61
62 return ( ($upv2*$zin) % $n, ($umv2*$xin) % $n );
63 }
64
65 sub _double {
66 my ($x, $z, $n, $d) = @_;
67
68 my $u = ($x + $z) ** 2;
69 my $v = ($x - $z) ** 2;
70 my $w = $u - $v;
71 my $t = $d * $w + $v;
72 return ( ($u * $v) % $n , ($w * $t) % $n );
73 }
74
75 sub mul {
76 my ($self, $k) = @_;
77 my $x = $self->{'x'};
78 my $z = $self->{'z'};
79 my $n = $self->{'n'};
80 my $d = $self->{'d'};
81
82 my ($x1, $x2, $z1, $z2);
83
84 my $r = --$k;
85 my $l = -1;
86 while ($r != 1) { $r >>= 1; $l++ }
87 if ($k & (1 << $l)) {
88 ($x2, $z2) = _double($x, $z, $n, $d);
89 ($x1, $z1) = _add3($x2, $z2, $x, $z, $x, $z, $n);
90 ($x2, $z2) = _double($x2, $z2, $n, $d);
91 } else {
92 ($x1, $z1) = _double($x, $z, $n, $d);
93 ($x2, $z2) = _add3($x, $z, $x1, $z1, $x, $z, $n);
94 }
95 $l--;
96 while ($l >= 1) {
97 if ($k & (1 << $l)) {
98 ($x1, $z1) = _add3($x1, $z1, $x2, $z2, $x, $z, $n);
99 ($x2, $z2) = _double($x2, $z2, $n, $d);
100 } else {
101 ($x2, $z2) = _add3($x2, $z2, $x1, $z1, $x, $z, $n);
102 ($x1, $z1) = _double($x1, $z1, $n, $d);
103 }
104 $l--;
105 }
106 if ($k & 1) {
107 ($x, $z) = _double($x2, $z2, $n, $d);
108 } else {
109 ($x, $z) = _add3($x2, $z2, $x1, $z1, $x, $z, $n);
110 }
111
112 $self->{'x'} = $x;
113 $self->{'z'} = $z;
114 return $self;
115 }
116
117 sub add {
118 my ($self, $other) = @_;
119 croak "add takes a EC point"
120 unless ref($other) eq 'Math::Prime::Util::ECProjectivePoint';
121 croak "second point is not on the same curve"
122 unless $self->{'c'} == $other->{'c'} &&
123 $self->{'n'} == $other->{'n'};
124
125 ($self->{'x'}, $self->{'z'}) = _add3($self->{'x'}, $self->{'z'},
126 $other->{'x'}, $other->{'z'},
127 $self->{'x'}, $self->{'z'},
128 $self->{'n'});
129 return $self;
130 }
131
132 sub double {
133 my ($self) = @_;
134 ($self->{'x'}, $self->{'z'}) = _double($self->{'x'}, $self->{'z'}, $self->{'n'}, $self->{'d'});
135 return $self;
136 }
137
138 sub _extended_gcd {
139 my ($a, $b) = @_;
140 my $zero = $a-$a;
141 my ($x, $lastx, $y, $lasty) = ($zero, $zero+1, $zero+1, $zero);
142 while ($b != 0) {
143 my $q = int($a/$b);
144 ($a, $b) = ($b, $a % $b);
145 ($x, $lastx) = ($lastx - $q*$x, $x);
146 ($y, $lasty) = ($lasty - $q*$y, $y);
147 }
148 return ($a, $lastx, $lasty);
149 }
150
151 sub normalize {
152 my ($self) = @_;
153 my $n = $self->{'n'};
154 my $z = $self->{'z'};
155 #my ($f, $u, undef) = _extended_gcd( $z, $n );
156 my $f = Math::BigInt::bgcd( $z, $n );
157 my $u = $z->copy->bmodinv($n);
158 $self->{'x'} = ( $self->{'x'} * $u ) % $n;
159 $self->{'z'} = $n-$n+1;
160 $self->{'f'} = ($f * $self->{'f'}) % $n;
161 return $self;
162 }
163
164 sub c { return shift->{'c'}; }
165 sub d { return shift->{'d'}; }
166 sub n { return shift->{'n'}; }
167 sub x { return shift->{'x'}; }
168 sub z { return shift->{'z'}; }
169 sub f { return shift->{'f'}; }
170
171 sub is_infinity {
172 my $self = shift;
173 return ($self->{'x'}->is_zero() && $self->{'z'}->is_one());
174 }
175
176 sub copy {
177 my $self = shift;
178 return Math::Prime::Util::ECProjectivePoint->new(
179 $self->{'c'}, $self->{'n'}, $self->{'x'}, $self->{'z'});
180 }
181
182 1;
183
184 __END__
185
186
187 # ABSTRACT: Elliptic curve operations for projective points
188
189 =pod
190
191 =encoding utf8
192
193
194 =head1 NAME
195
196 Math::Prime::Util::ECProjectivePoint - Elliptic curve operations for projective points
197
198
199 =head1 VERSION
200
201 Version 0.26
202
203
204 =head1 SYNOPSIS
205
206 # Create a point on a curve (a,b,n) with coordinates 0,1
207 my $ECP = Math::Prime::Util::ECProjectivePoint->new($c, $n, 0, 1);
208
209 # scalar multiplication by k.
210 $ECP->mul($k)
211
212 # add two points on the same curve
213 $ECP->add($ECP2)
214
215 say "P = O" if $ECP->is_infinity();
216
217 =head1 DESCRIPTION
218
219 This really should just be in Math::EllipticCurve.
220
221 To write.
222
223
224 =head1 FUNCTIONS
225
226 =head2 new
227
228 $point = Math::Prime::Util::ECProjectivePoint->new(c, n, x, z);
229
230 Returns a new point on the curve defined by the Montgomery parameter c.
231
232 =head2 c
233
234 =head2 n
235
236 Returns the C<c>, C<d>, or C<n> values that describe the curve.
237
238 =head2 d
239
240 Returns the precalculated value of C<int( (c + 2) / 4 )>.
241
242 =head2 x
243
244 =head2 z
245
246 Returns the C<x> or C<z> values that define the point on the curve.
247
248 =head2 f
249
250 Returns a possible factor found after L</normalize>.
251
252 =head2 add
253
254 Takes another point on the same curve as an argument and adds it this point.
255
256 =head2 double
257
258 Double the current point on the curve.
259
260 =head2 mul
261
262 Takes an integer and performs scalar multiplication of the point.
263
264 =head2 is_infinity
265
266 Returns true if the point is (0,1), which is the point at infinity for
267 the affine coordinates.
268
269 =head2 copy
270
271 Returns a copy of the point.
272
273 =head2 normalize
274
275 Performs an extended gcd operation to make C<z=1>. If a factor of C<n> is
276 found it is put in C<f>.
277
278
279 =head1 SEE ALSO
280
281 L<Math::EllipticCurve::Prime>
282
283 This really should just be in a L<Math::EllipticCurve> module.
284
285 =head1 AUTHORS
286
287 Dana Jacobsen E<lt>dana@acm.orgE<gt>
288
289
290 =head1 COPYRIGHT
291
292 Copyright 2012-2013 by Dana Jacobsen E<lt>dana@acm.orgE<gt>
293
294 This program is free software; you can redistribute it and/or modify it under the same terms as Perl itself.
295
296 =cut
+700
-142
lib/Math/Prime/Util/PP.pm less more
44
55 BEGIN {
66 $Math::Prime::Util::PP::AUTHORITY = 'cpan:DANAJ';
7 $Math::Prime::Util::PP::VERSION = '0.24';
7 $Math::Prime::Util::PP::VERSION = '0.26';
88 }
99
1010 # The Pure Perl versions of all the Math::Prime::Util routines.
881881 # It's now time to perform the Lucas pseudoprimality test using $D.
882882
883883 if (ref($n) ne 'Math::BigInt') {
884 if (!defined $MATH::BigInt::VERSION) {
884 if (!defined $Math::BigInt::VERSION) {
885885 eval { require Math::BigInt; Math::BigInt->import(try=>'GMP,Pari'); 1; }
886886 or do { croak "Cannot load Math::BigInt "; }
887887 }
10301030 sub is_aks_prime {
10311031 my $n = shift;
10321032
1033 if (!defined $MATH::BigInt::VERSION) {
1033 if (!defined $Math::BigInt::VERSION) {
10341034 eval { require Math::BigInt; Math::BigInt->import(try=>'GMP,Pari'); 1; }
10351035 or do { croak "Cannot load Math::BigInt "; }
10361036 }
1037 if (!defined $MATH::BigFloat::VERSION) {
1037 if (!defined $Math::BigFloat::VERSION) {
10381038 eval { require Math::BigFloat; Math::BigFloat->import(); 1; }
10391039 or do { croak "Cannot load Math::BigFloat "; }
10401040 }
11361136 @factors;
11371137 }
11381138
1139 my $_holf_r;
1140 my @_fsublist = (
1141 sub { my $n = shift; return ($n) unless $n < $_half_word;
1142 holf_factor ($n, 64*1024, $_holf_r); $_holf_r += 64*1024; },
1143 sub { prho_factor (shift, 8*1024, 3) },
1144 sub { pbrent_factor (shift, 32*1024, 1) },
1145 sub { pminus1_factor(shift, 10_000); },
1146 sub { pminus1_factor(shift, 600_000); },
1147 sub { pbrent_factor (shift, 512*1024, 7) },
1148 sub { ecm_factor (shift, 1_000, 5_000, 10) },
1149 sub { pminus1_factor(shift, 4_000_000); },
1150 sub { pbrent_factor (shift, 512*1024, 11) },
1151 sub { ecm_factor (shift, 10_000, 50_000, 10) },
1152 sub { holf_factor (shift, 256*1024, $_holf_r); $_holf_r += 256*1024; },
1153 sub { pminus1_factor(shift,20_000_000); },
1154 sub { ecm_factor (shift, 100_000, 800_000, 10) },
1155 sub { holf_factor (shift, 512*1024, $_holf_r); $_holf_r += 512*1024; },
1156 sub { pbrent_factor (shift, 2048*1024, 13) },
1157 sub { holf_factor (shift, 2048*1024, $_holf_r); $_holf_r += 2048*1024; },
1158 sub { ecm_factor (shift, 1_000_000, 1_000_000, 10) },
1159 sub { pminus1_factor(shift, 100_000_000, 500_000_000); },
1160 );
1161
11391162 sub factor {
11401163 my($n) = @_;
11411164 _validate_positive_integer($n);
11671190 #print "Looking at $n with stack ", join(",",@nstack), "\n";
11681191 while ( ($n >= (31*31)) && !_is_prime7($n) ) {
11691192 my @ftry;
1170 my $holf_rounds = 0;
1171 if ($n < $_half_word) {
1172 $holf_rounds = 64*1024;
1173 #warn "trying holf 64k on $n\n";
1174 @ftry = holf_factor($n, $holf_rounds);
1175 }
1176 if (scalar @ftry < 2) {
1177 foreach my $add (3, 5, 7, 11, 13) {
1178 #warn "trying prho 64k {$add} on $n\n" if scalar @ftry < 2;
1179 @ftry = prho_factor($n, 64*1024, $add) if scalar @ftry < 2;
1180 }
1181 }
1182 if (scalar @ftry < 2) {
1183 #warn "trying holf 128k on $n\n";
1184 @ftry = holf_factor($n, 128*1024, $holf_rounds);
1185 $holf_rounds += 128*1024;
1186 }
1187 if (scalar @ftry < 2) {
1188 #warn "trying prho 128k {17} on $n\n";
1189 @ftry = prho_factor($n, 128*1024, 17);
1193 $_holf_r = 1;
1194 foreach my $sub (@_fsublist) {
1195 last if scalar @ftry >= 2;
1196 @ftry = $sub->($n);
11901197 }
11911198 if (scalar @ftry > 1) {
11921199 #print " split into ", join(",",@ftry), "\n";
12051212 sort {$a<=>$b} @factors;
12061213 }
12071214
1215 sub _found_factor {
1216 my($f, $n, $what, @factors) = @_;
1217 if ($f == 1 || $f == $n) {
1218 push @factors, $n;
1219 } else {
1220 # Perl 5.6.2 needs things spelled out for it.
1221 my $f2 = (ref($n) eq 'Math::BigInt') ? $n->copy->bdiv($f)->as_int
1222 : int($n/$f);
1223 push @factors, $f;
1224 push @factors, $f2;
1225 croak "internal error in $what" unless $f * $f2 == $n;
1226 # MPU::GMP prints this type of message if verbose, so do the same.
1227 print "$what found factor $f\n" if Math::Prime::Util::prime_get_config()->{'verbose'} > 0;
1228 }
1229 @factors;
1230 }
1231
12081232 # TODO:
1209 sub fermat_factor { trial_factor(@_) }
12101233 sub squfof_factor { trial_factor(@_) }
12111234
12121235 sub prho_factor {
12281251 $V = $n->copy->bzero->badd($V);
12291252 for my $i (1 .. $rounds) {
12301253 # Would use bmuladd here, but old Math::BigInt's barf with scalar $a.
1231 #$U->bmuladd($U, $a); $U->bmod($n);
1232 #$V->bmuladd($V, $a); $V->bmod($n);
1233 #$V->bmuladd($V, $a); $V->bmod($n);
1234 $U->bmul($U); $U->badd($a); $U->bmod($n);
1235 $V->bmul($V); $V->badd($a); $V->bmod($n);
1236 $V->bmul($V); $V->badd($a); $V->bmod($n);
1254 $U->bmul($U)->badd($a)->bmod($n);
1255 $V->bmul($V)->badd($a);
1256 $V->bmul($V)->badd($a)->bmod($n);
12371257 my $f = Math::BigInt::bgcd( ($U > $V) ? $U-$V : $V-$U, $n);
12381258 if ($f == $n) {
12391259 last if $inloop++; # We've been here before
12401260 } elsif ($f != 1) {
1241 my $f2 = $n->copy->bdiv($f)->as_int;
1242 push @factors, $f;
1243 push @factors, $f2;
1244 croak "internal error in prho" unless ($f * $f2) == $n;
1245 return @factors;
1261 return _found_factor($f, $n, "prho", @factors);
12461262 }
12471263 }
12481264
12561272 if ($f == $n) {
12571273 last if $inloop++; # We've been here before
12581274 } elsif ($f != 1) {
1259 push @factors, $f;
1260 push @factors, int($n/$f);
1261 croak "internal error in prho" unless ($f * int($n/$f)) == $n;
1262 return @factors;
1275 return _found_factor($f, $n, "prho", @factors);
12631276 }
12641277 }
12651278
12661279 } else {
12671280
12681281 for my $i (1 .. $rounds) {
1269 # U^2+a % n
1270 $U = _mulmod($U, $U, $n);
1271 $U = (($n-$U) > $a) ? $U+$a : $U+$a-$n;
1272 # V^2+a % n
1273 $V = _mulmod($V, $V, $n);
1274 $V = (($n-$V) > $a) ? $V+$a : $V+$a-$n;
1275 # V^2+a % n
1276 $V = _mulmod($V, $V, $n);
1277 $V = (($n-$V) > $a) ? $V+$a : $V+$a-$n;
1282 $U = _mulmod($U, $U, $n); $U = (($n-$U) > $a) ? $U+$a : $U-$n+$a;
1283 $V = _mulmod($V, $V, $n); $V = (($n-$V) > $a) ? $V+$a : $V-$n+$a;
1284 $V = _mulmod($V, $V, $n); $V = (($n-$V) > $a) ? $V+$a : $V-$n+$a;
12781285 my $f = _gcd_ui( ($U > $V) ? $U-$V : $V-$U, $n );
12791286 if ($f == $n) {
12801287 last if $inloop++; # We've been here before
12811288 } elsif ($f != 1) {
1282 push @factors, $f;
1283 push @factors, int($n/$f);
1284 croak "internal error in prho" unless ($f * int($n/$f)) == $n;
1285 return @factors;
1289 return _found_factor($f, $n, "prho", @factors);
12861290 }
12871291 }
12881292
12921296 }
12931297
12941298 sub pbrent_factor {
1295 my($n, $rounds) = @_;
1299 my($n, $rounds, $a) = @_;
12961300 _validate_positive_integer($n);
12971301 $rounds = 4*1024*1024 unless defined $rounds;
1302 $a = 3 unless defined $a;
12981303
12991304 my @factors = _basic_factor($n);
13001305 return @factors if $n < 4;
13011306
1302 my $a = 1;
13031307 my $Xi = 2;
13041308 my $Xm = 2;
13051309
13061310 if ( ref($n) eq 'Math::BigInt' ) {
13071311
1308 $Xi = $n->copy->bzero->badd($Xi);
1309 $Xm = $n->copy->bzero->badd($Xm);
1310 for my $i (1 .. $rounds) {
1311 $Xi->bmul($Xi); $Xi->badd($a); $Xi->bmod($n);
1312 my $f = Math::BigInt::bgcd( ($Xi > $Xm) ? $Xi-$Xm : $Xm-$Xi, $n);
1313 if ( ($f != 1) && ($f != $n) ) {
1314 my $f2 = $n->copy->bdiv($f)->as_int;
1315 push @factors, $f;
1316 push @factors, $f2;
1317 croak "internal error in pbrent" unless ($f * $f2) == $n;
1318 return @factors;
1319 }
1320 $Xm = $Xi->copy if ($i & ($i-1)) == 0; # i is a power of 2
1312 # Same code as the GMP version, but runs *much* slower. Even with
1313 # Math::BigInt::GMP it's >200x slower. With the default Calc backend
1314 # it's thousands of times slower.
1315 my $inner = 256;
1316 my $zero = $n->copy->bzero;
1317 my $saveXi;
1318 my $f;
1319 $Xi = $zero->copy->badd($Xi);
1320 $Xm = $zero->copy->badd($Xm);
1321 my $r = 1;
1322 while ($rounds > 0) {
1323 my $rleft = ($r > $rounds) ? $rounds : $r;
1324 while ($rleft > 0) {
1325 my $dorounds = ($rleft > $inner) ? $inner : $rleft;
1326 my $m = $zero->copy->bone;
1327 $saveXi = $Xi->copy;
1328 foreach my $i (1 .. $dorounds) {
1329 $Xi->bmul($Xi)->badd($a)->bmod($n);
1330 $m->bmul($Xi - $Xm);
1331 }
1332 $rleft -= $dorounds;
1333 $rounds -= $dorounds;
1334 $m->bmod($n);
1335 $f = Math::BigInt::bgcd( $m, $n);
1336 last if $f != 1;
1337 }
1338 if ($f == 1) {
1339 $r *= 2;
1340 $Xm = $Xi->copy;
1341 next;
1342 }
1343 if ($f == $n) { # back up to determine the factor
1344 $Xi = $saveXi->copy;
1345 do {
1346 $Xi->bmul($Xi)->badd($a)->bmod($n);
1347 $f = Math::BigInt::bgcd( ($Xi > $Xm) ? $Xi-$Xm : $Xm-$Xi, $n);
1348 } while ($f != 1 && $r-- != 0);
1349 last if $f == 1 || $f == $n;
1350 }
1351 return _found_factor($f, $n, "pbrent", @factors);
13211352 }
13221353
13231354 } elsif ($n < $_half_word) {
13251356 for my $i (1 .. $rounds) {
13261357 $Xi = ($Xi * $Xi + $a) % $n;
13271358 my $f = _gcd_ui( ($Xi > $Xm) ? $Xi-$Xm : $Xm-$Xi, $n );
1328 if ( ($f != 1) && ($f != $n) ) {
1329 push @factors, $f;
1330 push @factors, int($n/$f);
1331 croak "internal error in pbrent" unless ($f * int($n/$f)) == $n;
1332 return @factors;
1333 }
1359 return _found_factor($f, $n, "pbrent", @factors) if $f != 1 && $f != $n;
13341360 $Xm = $Xi if ($i & ($i-1)) == 0; # i is a power of 2
13351361 }
13361362
13411367 $Xi = _mulmod($Xi, $Xi, $n);
13421368 $Xi = (($n-$Xi) > $a) ? $Xi+$a : $Xi+$a-$n;
13431369 my $f = _gcd_ui( ($Xi > $Xm) ? $Xi-$Xm : $Xm-$Xi, $n );
1344 if ( ($f != 1) && ($f != $n) ) {
1345 push @factors, $f;
1346 push @factors, int($n/$f);
1347 croak "internal error in pbrent" unless ($f * int($n/$f)) == $n;
1348 return @factors;
1349 }
1370 return _found_factor($f, $n, "pbrent", @factors) if $f != 1 && $f != $n;
13501371 $Xm = $Xi if ($i & ($i-1)) == 0; # i is a power of 2
13511372 }
13521373
13551376 @factors;
13561377 }
13571378
1358 # This code is bollocks. See a proper implementation in factor.c
13591379 sub pminus1_factor {
1360 my($n, $rounds) = @_;
1380 my($n, $B1, $B2) = @_;
13611381 _validate_positive_integer($n);
1362 $rounds = 4*1024*1024 unless defined $rounds;
13631382
13641383 my @factors = _basic_factor($n);
13651384 return @factors if $n < 4;
13661385
1367 if ( ref($n) eq 'Math::BigInt' ) {
1368 my $kf = $n->copy->bzero->badd(13);
1369 for my $i (1 .. $rounds) {
1370 $kf->bmodpow($i,$n);
1371 $kf = $n if $kf == 0;
1372 my $f = Math::BigInt::bgcd( $kf-1, $n );
1373 if ( ($f != 1) && ($f != $n) ) {
1374 my $f2 = $n->copy->bdiv($f)->as_int;
1375 push @factors, $f;
1376 push @factors, $f2;
1377 croak "internal error in pminus1" unless ($f * $f2) == $n;
1386 if ( ref($n) ne 'Math::BigInt' ) {
1387 # Stage 1 only
1388 $B1 = 10_000_000 unless defined $B1;
1389 my $a = 2;
1390 my $f = 1;
1391 my($pc_beg, $pc_end, @bprimes);
1392 $pc_beg = 2;
1393 $pc_end = $pc_beg + 100_000;
1394 my $sqrtb1 = int(sqrt($B1));
1395 while (1) {
1396 $pc_end = $B1 if $pc_end > $B1;
1397 @bprimes = @{ primes($pc_beg, $pc_end) };
1398 while (@bprimes) {
1399 my $q = shift @bprimes;
1400 my $k = $q;
1401 if ($q <= $sqrtb1) {
1402 my $kmin = int($B1 / $q);
1403 while ($k <= $kmin) { $k *= $q; }
1404 }
1405 $a = _powmod($a, $k, $n);
1406 if ($a == 0) { push @factors, $n; return @factors; }
1407 my $f = _gcd_ui( $a-1, $n );
1408 return _found_factor($f, $n, "pminus1", @factors) if $f != 1;
1409 }
1410 last if $pc_end >= $B1;
1411 $pc_beg = $pc_end+1;
1412 $pc_end += 500_000;
1413 }
1414 push @factors, $n;
1415 return @factors;
1416 }
1417
1418 # Stage 2 isn't really any faster than stage 1 for the examples I've tried.
1419
1420 if (!defined $B1) {
1421 for my $mul (1, 100, 1000, 10_000, 100_000, 1_000_000) {
1422 $B1 = 1000 * $mul;
1423 $B2 = 1*$B1;
1424 #warn "Trying p-1 with $B1 / $B2\n";
1425 my @nf = pminus1_factor($n, $B1, $B2);
1426 if (scalar @nf > 1) {
1427 push @factors, @nf;
13781428 return @factors;
13791429 }
13801430 }
1381 } else {
1382 my $kf = 13;
1383 for my $i (1 .. $rounds) {
1384 $kf = _powmod($kf, $i, $n);
1385 $kf = $n if $kf == 0;
1386 my $f = _gcd_ui( $kf-1, $n );
1387 if ( ($f != 1) && ($f != $n) ) {
1388 push @factors, $f;
1389 push @factors, int($n/$f);
1390 croak "internal error in pminus1" unless ($f * int($n/$f)) == $n;
1391 return @factors;
1392 }
1393 }
1394 }
1395 push @factors, $n;
1396 @factors;
1431 push @factors, $n;
1432 return @factors;
1433 }
1434 $B2 = 1*$B1 unless defined $B2;
1435
1436 my $one = $n->copy->bone;
1437 my ($j, $q, $saveq) = (1, 2, 2);
1438 my $t = $one->copy;
1439 my $a = $one->copy->badd(1);
1440 my $savea = $a->copy;
1441 my $f = 1;
1442 my($pc_beg, $pc_end, @bprimes);
1443
1444 $pc_beg = 2;
1445 $pc_end = $pc_beg + 100_000;
1446 while (1) {
1447 $pc_end = $B1 if $pc_end > $B1;
1448 @bprimes = @{ primes($pc_beg, $pc_end) };
1449 while (@bprimes) {
1450 $q = shift @bprimes;
1451 my $k = $q;
1452 my $kmin = int($B1 / $q);
1453 while ($k <= $kmin) { $k *= $q; }
1454 $t *= $k; # accumulate powers for a
1455 if ( ($j++ % 32) == 0) {
1456 $a->bmodpow($t, $n);
1457 $t = $one->copy;
1458 if ($a == 0) { push @factors, $n; return @factors; }
1459 $f = Math::BigInt::bgcd( $a-1, $n );
1460 last if $f == $n;
1461 if ($f != 1) {
1462 push @factors, $f, $n/$f;
1463 return @factors;
1464 }
1465 $saveq = $q;
1466 $savea = $a->copy;
1467 }
1468 }
1469 last if $f != 1 || $pc_end >= $B1;
1470 $pc_beg = $pc_end+1;
1471 $pc_end += 500_000;
1472 }
1473 undef @bprimes;
1474 $a->bmodpow($t, $n);
1475 if ($a == 0) { push @factors, $n; return @factors; }
1476 $f = Math::BigInt::bgcd( $a-1, $n );
1477 if ($f == $n) {
1478 $q = $saveq;
1479 $a = $savea->copy;
1480 while ($q <= $B1) {
1481 my $k = $q;
1482 my $kmin = int($B1 / $q);
1483 while ($k <= $kmin) { $k *= $q; }
1484 $a->bmodpow($k, $n);
1485 my $f = Math::BigInt::bgcd( $a-1, $n );
1486 if ($f == $n) { push @factors, $n; return @factors; }
1487 last if $f != 1;
1488 $q = next_prime($q);
1489 }
1490 }
1491 # STAGE 2
1492 if ($f == 1 && $B2 > $B1) {
1493 my $bm = $a->copy;
1494 my $b = $one->copy;
1495 my @precomp_bm;
1496 $precomp_bm[0] = ($bm * $bm) % $n;
1497 foreach my $j (1..19) {
1498 $precomp_bm[$j] = ($precomp_bm[$j-1] * $bm * $bm) % $n;
1499 }
1500 $a->bmodpow($q, $n);
1501 my $j = 1;
1502 $pc_beg = $q+1;
1503 $pc_end = $pc_beg + 100_000;
1504 while (1) {
1505 $pc_end = $B2 if $pc_end > $B2;
1506 @bprimes = @{ primes($pc_beg, $pc_end) };
1507 while (@bprimes) {
1508 my $lastq = $q;
1509 $q = shift @bprimes;
1510 my $qdiff = ($q - $lastq) / 2 - 1;
1511 if (!defined $precomp_bm[$qdiff]) {
1512 $precomp_bm[$qdiff] = $bm->copy->bmodpow($q-$lastq, $n);
1513 }
1514 $a->bmul($precomp_bm[$qdiff])->bmod($n);
1515 if ($a == 0) { push @factors, $n; return @factors; }
1516 $b->bmul($a-1);
1517 if (($j++ % 64) == 0) {
1518 $b->bmod($n);
1519 $f = Math::BigInt::bgcd( $b, $n );
1520 last if $f != 1;
1521 }
1522 }
1523 last if $f != 1 || $pc_end >= $B2;
1524 $pc_beg = $pc_end+1;
1525 $pc_end += 500_000;
1526 }
1527 $f = Math::BigInt::bgcd( $b, $n );
1528 }
1529 return _found_factor($f, $n, "pminus1", @factors);
13971530 }
13981531
13991532 sub holf_factor {
14011534 _validate_positive_integer($n);
14021535 $rounds = 64*1024*1024 unless defined $rounds;
14031536 $startrounds = 1 unless defined $startrounds;
1537 $startrounds = 1 if $startrounds < 1;
14041538
14051539 my @factors = _basic_factor($n);
14061540 return @factors if $n < 4;
14091543 for my $i ($startrounds .. $rounds) {
14101544 my $ni = $n->copy->bmul($i);
14111545 my $s = $ni->copy->bsqrt->bfloor->as_int;
1412 $s->binc if ($s * $s) != $ni;
1413 my $m = $s->copy->bmul($s)->bmod($n);
1546 if ($s * $s == $ni) {
1547 # s^2 = n*i, so m = s^2 mod n = 0. Hence f = GCD(n, s) = GCD(n, n*i)
1548 my $f = Math::BigInt::bgcd($ni, $n);
1549 return _found_factor($f, $n, "HOLF", @factors);
1550 }
1551 $s->binc;
1552 my $m = ($s * $s) - $ni;
14141553 # Check for perfect square
14151554 my $mc = int(($m & 31)->bstr);
14161555 next unless $mc==0||$mc==1||$mc==4||$mc==9||$mc==16||$mc==17||$mc==25;
14171556 my $f = $m->copy->bsqrt->bfloor->as_int;
14181557 next unless ($f*$f) == $m;
14191558 $f = Math::BigInt::bgcd( ($s > $f) ? $s-$f : $f-$s, $n);
1420 last if $f == 1 || $f == $n; # Should never happen
1421 my $f2 = $n->copy->bdiv($f)->as_int;
1422 push @factors, $f;
1423 push @factors, $f2;
1424 croak "internal error in HOLF" unless ($f * $f2) == $n;
1425 # print "HOLF found factors in $i rounds\n";
1426 return @factors;
1559 return _found_factor($f, $n, "HOLF ($i rounds)", @factors);
14271560 }
14281561 } else {
14291562 for my $i ($startrounds .. $rounds) {
14361569 my $f = int(sqrt($m));
14371570 next unless $f*$f == $m;
14381571 $f = _gcd_ui( ($s > $f) ? $s - $f : $f - $s, $n);
1439 last if $f == 1 || $f == $n; # Should never happen
1440 push @factors, $f;
1441 push @factors, int($n/$f);
1442 croak "internal error in HOLF" unless ($f * int($n/$f)) == $n;
1443 # print "HOLF found factors in $i rounds\n";
1444 return @factors;
1572 return _found_factor($f, $n, "HOLF ($i rounds)", @factors);
14451573 }
14461574 }
14471575 push @factors, $n;
14481576 @factors;
14491577 }
14501578
1451
1579 sub fermat_factor {
1580 my($n, $rounds) = @_;
1581 _validate_positive_integer($n);
1582 $rounds = 64*1024*1024 unless defined $rounds;
1583
1584 my @factors = _basic_factor($n);
1585 return @factors if $n < 4;
1586
1587 if ( ref($n) eq 'Math::BigInt' ) {
1588 my $a = $n->copy->bsqrt->bfloor->as_int;
1589 return _found_factor($a, $n, "Fermat", @factors) if $a*$a == $n;
1590 $a++;
1591 my $b2 = $a*$a - $n;
1592 my $lasta = $a + $rounds;
1593 while ($a <= $lasta) {
1594 my $mc = int(($b2 & 31)->bstr);
1595 if ($mc==0||$mc==1||$mc==4||$mc==9||$mc==16||$mc==17||$mc==25) {
1596 my $s = $b2->copy->bsqrt->bfloor->as_int;
1597 if ($s*$s == $b2) {
1598 my $i = $a-($lasta-$rounds)+1;
1599 return _found_factor($a - $s, $n, "Fermat ($i rounds)", @factors);
1600 }
1601 }
1602 $a++;
1603 $b2 = $a*$a-$n;
1604 }
1605 } else {
1606 my $a = int(sqrt($n));
1607 return _found_factor($a, $n, "Fermat", @factors) if $a*$a == $n;
1608 $a++;
1609 my $b2 = $a*$a - $n;
1610 my $lasta = $a + $rounds;
1611 while ($a <= $lasta) {
1612 my $mc = $b2 & 31;
1613 if ($mc==0||$mc==1||$mc==4||$mc==9||$mc==16||$mc==17||$mc==25) {
1614 my $s = int(sqrt($b2));
1615 if ($s*$s == $b2) {
1616 my $i = $a-($lasta-$rounds)+1;
1617 return _found_factor($a - $s, $n, "Fermat ($i rounds)", @factors);
1618 }
1619 }
1620 $a++;
1621 $b2 = $a*$a-$n;
1622 }
1623 }
1624 push @factors, $n;
1625 @factors;
1626 }
1627
1628
1629 sub ecm_factor {
1630 my($n, $B1, $B2, $ncurves) = @_;
1631 _validate_positive_integer($n);
1632
1633 my @factors = _basic_factor($n);
1634 return @factors if $n < 4;
1635
1636 $ncurves = 10 unless defined $ncurves;
1637
1638 if (!defined $B1) {
1639 for my $mul (1, 10, 100, 1000, 10_000, 100_000, 1_000_000) {
1640 $B1 = 100 * $mul;
1641 $B2 = 10*$B1;
1642 #warn "Trying ecm with $B1 / $B2\n";
1643 my @nf = ecm_factor($n, $B1, $B2, $ncurves);
1644 if (scalar @nf > 1) {
1645 push @factors, @nf;
1646 return @factors;
1647 }
1648 }
1649 push @factors, $n;
1650 return @factors;
1651 }
1652
1653 $B2 = 10*$B1 unless defined $B2;
1654 my $sqrt_b1 = int(sqrt($B1)+1);
1655
1656 # Affine code. About 3x slower than the projective, and no stage 2.
1657 #
1658 #if (!defined $Math::Prime::Util::ECAffinePoint::VERSION) {
1659 # eval { require Math::Prime::Util::ECAffinePoint; 1; }
1660 # or do { croak "Cannot load Math::Prime::Util::ECAffinePoint"; };
1661 #}
1662 #my @bprimes = @{ primes(2, $B1) };
1663 #my $irandf = Math::Prime::Util::_get_rand_func();
1664 #foreach my $curve (1 .. $ncurves) {
1665 # my $a = $irandf->($n-1);
1666 # my $b = 1;
1667 # my $ECP = Math::Prime::Util::ECAffinePoint->new($a, $b, $n, 0, 1);
1668 # foreach my $q (@bprimes) {
1669 # my $k = $q;
1670 # if ($k < $sqrt_b1) {
1671 # my $kmin = int($B1 / $q);
1672 # while ($k <= $kmin) { $k *= $q; }
1673 # }
1674 # $ECP->mul($k);
1675 # my $f = $ECP->f;
1676 # if ($f != 1) {
1677 # last if $f == $n;
1678 # warn "ECM found factors with B1 = $B1 in curve $curve\n";
1679 # return _found_factor($f, $n, "ECM B1=$B1 curve $curve", @factors);
1680 # }
1681 # last if $ECP->is_infinity;
1682 # }
1683 #}
1684
1685 if (!defined $Math::Prime::Util::ECProjectivePoint::VERSION) {
1686 eval { require Math::Prime::Util::ECProjectivePoint; 1; }
1687 or do { croak "Cannot load Math::Prime::Util::ECProjectivePoint"; };
1688 }
1689
1690 # With multiple curves, it's better to get all the primes at once.
1691 # The downside is this can kill memory with a very large B1.
1692 my @bprimes = @{ primes(3, $B1) };
1693 my @b2primes = ($B2 > $B1) ? @{primes($B1+1, $B2)} : ();
1694 my $irandf = Math::Prime::Util::_get_rand_func();
1695
1696 foreach my $curve (1 .. $ncurves) {
1697 my $sigma = $irandf->($n-1-6) + 6;
1698 my ($u, $v) = ( ($sigma*$sigma - 5) % $n, (4 * $sigma) % $n );
1699 my ($x, $z) = ( ($u*$u*$u) % $n, ($v*$v*$v) % $n );
1700 my $b = (4 * $x * $v) % $n;
1701 my $a = ( (($v-$u)**3) * (3*$u + $v) ) % $n;
1702 my $f = Math::BigInt::bgcd( $b, $n );
1703 $f = Math::BigInt::bgcd( $z, $n ) if $f == 1;
1704 next if $f == $n;
1705 return _found_factor($f,$n, "ECM B1=$B1 curve $curve", @factors) if $f != 1;
1706 $u = $b->copy->bmodinv($n);
1707 $a = (($a*$u) - 2) % $n;
1708
1709 my $ECP = Math::Prime::Util::ECProjectivePoint->new($a, $n, $x, $z);
1710 my $fm = $n-$n+1;
1711 my $i = 15;
1712
1713 for (my $q = 2; $q < $B1; $q *= 2) { $ECP->double(); }
1714 foreach my $q (@bprimes) {
1715 my $k = $q;
1716 if ($k < $sqrt_b1) {
1717 my $kmin = int($B1 / $q);
1718 while ($k <= $kmin) { $k *= $q; }
1719 }
1720 $ECP->mul($k);
1721 $fm = ($fm * $ECP->x() ) % $n;
1722 if ($i++ % 32 == 0) {
1723 $f = Math::BigInt::bgcd($fm, $n);
1724 last if $f != 1;
1725 }
1726 }
1727 $f = Math::BigInt::bgcd($fm, $n);
1728 next if $f == $n;
1729
1730 if ($f == 1 && $B2 > $B1) { # BEGIN STAGE 2
1731 my $D = int(sqrt($B2/2)); $D++ if $D % 2;
1732 my $one = $n - $n + 1;
1733 my $g = $one;
1734
1735 my $S2P = $ECP->copy->normalize;
1736 $f = $S2P->f;
1737 if ($f != 1) {
1738 next if $f == $n;
1739 #warn "ECM S2 normalize f=$f\n" if $f != 1;
1740 return _found_factor($f, $n, "ECM S2 B1=$B1 curve $curve");
1741 }
1742 my $S2x = $S2P->x;
1743 my $S2d = $S2P->d;
1744 my @nqx = ($n-$n, $S2x);
1745
1746 foreach my $i (2 .. 2*$D) {
1747 my($x2, $z2);
1748 if ($i % 2) {
1749 ($x2, $z2) = Math::Prime::Util::ECProjectivePoint::_addx($nqx[($i-1)/2], $nqx[($i+1)/2], $S2x, $n);
1750 } else {
1751 ($x2, $z2) = Math::Prime::Util::ECProjectivePoint::_double($nqx[$i/2], $one, $n, $S2d);
1752 }
1753 $nqx[$i] = $x2;
1754 #($f, $u, undef) = _extended_gcd($z2, $n);
1755 $f = Math::BigInt::bgcd( $z2, $n );
1756 last if $f != 1;
1757 $u = $z2->copy->bmodinv($n);
1758 $nqx[$i] = ($x2 * $u) % $n;
1759 }
1760 if ($f != 1) {
1761 next if $f == $n;
1762 #warn "ECM S2 1: B1 $B1 B2 $B2 curve $curve f=$f\n";
1763 return _found_factor($f, $n, "ECM S2 B1=$B1 curve $curve", @factors);
1764 }
1765
1766 $x = $nqx[2*$D-1];
1767 my $m = 1;
1768 while ($m < ($B2+$D)) {
1769 if ($m != 1) {
1770 my $oldx = $S2x;
1771 my ($x1, $z1) = Math::Prime::Util::ECProjectivePoint::_addx($nqx[2*$D], $S2x, $x, $n);
1772 $f = Math::BigInt::bgcd( $z1, $n );
1773 last if $f != 1;
1774 $u = $z1->copy->bmodinv($n);
1775 $S2x = ($x1 * $u) % $n;
1776 $x = $oldx;
1777 last if $f != 1;
1778 }
1779 if ($m+$D > $B1) {
1780 my @p = grep { $_ >= $m-$D && $_ <= $m+$D } @b2primes;
1781 foreach my $i (@p) {
1782 last if $i >= $m;
1783 $g = ($g * ($S2x - $nqx[$m+$D-$i])) % $n;
1784 }
1785 foreach my $i (@p) {
1786 next unless $i > $m;
1787 next if is_prime($m+$m-$i);
1788 $g = ($g * ($S2x - $nqx[$i-$m])) % $n;
1789 }
1790 $f = Math::BigInt::bgcd($g, $n);
1791 #warn "ECM S2 3: found $f in stage 2\n" if $f != 1;
1792 last if $f != 1;
1793 }
1794 $m += 2*$D;
1795 }
1796 } # END STAGE 2
1797
1798 next if $f == $n;
1799 if ($f != 1) {
1800 #warn "ECM found factors with B1 = $B1 in curve $curve\n";
1801 return _found_factor($f, $n, "ECM B1=$B1 curve $curve", @factors);
1802 }
1803 # end of curve loop
1804 }
1805 push @factors, $n;
1806 @factors;
1807 }
1808
1809
1810 sub primality_proof_lucas {
1811 my ($n) = shift;
1812 my @composite = (0, []);
1813
1814 # Since this can take a very long time with a composite, try some easy cuts
1815 return @composite if !defined $n || $n < 2;
1816 return (2, [$n]) if $n < 4;
1817 return @composite if miller_rabin($n,2,15,325) == 0;
1818
1819 if (!defined $Math::BigInt::VERSION) {
1820 eval { require Math::BigInt; Math::BigInt->import(try=>'GMP,Pari'); 1; }
1821 or do { croak "Cannot load Math::BigInt"; };
1822 }
1823
1824 my $nm1 = $n-1;
1825 my @factors = factor($nm1);
1826 { # remove duplicate factors and make a sorted array of bigints
1827 my %uf;
1828 undef @uf{@factors};
1829 @factors = sort {$a<=>$b} map { Math::BigInt->new("$_") } keys %uf;
1830 }
1831 for (my $a = 2; $a < $nm1; $a++) {
1832 my $ap = Math::BigInt->new("$a");
1833 # 1. a must be coprime to n
1834 next unless Math::BigInt::bgcd($ap, $n) == 1;
1835 # 2. a^(n-1) = 1 mod n.
1836 next unless $ap->copy->bmodpow($nm1, $n) == 1;
1837 # 3. a^((n-1)/f) != 1 mod n for all f.
1838 next if (scalar grep { $_ == 1 }
1839 map { $ap->copy->bmodpow(int($nm1/$_),$n); }
1840 @factors) > 0;
1841 # Verify each factor and add to proof
1842 my @fac_proofs;
1843 foreach my $f (@factors) {
1844 my ($isp, $fproof) = Math::Prime::Util::is_provable_prime_with_cert($f);
1845 if ($isp != 2) {
1846 carp "could not prove primality of $n.\n";
1847 return (1, []);
1848 }
1849 push @fac_proofs, (scalar @$fproof == 1) ? $fproof->[0] : $fproof;
1850 }
1851 my @proof = ("$n", "Pratt", [@fac_proofs], $a);
1852 return (2, [@proof]);
1853 }
1854 return @composite;
1855 }
1856
1857 sub primality_proof_bls75 {
1858 my ($n) = shift;
1859 my @composite = (0, []);
1860
1861 # Since this can take a very long time with a composite, try some easy cuts
1862 return @composite if !defined $n || $n < 2;
1863 return (2, [$n]) if $n < 4;
1864 return @composite if ($n & 1) == 0;
1865 return @composite if miller_rabin($n,2,15,325) == 0;
1866
1867 if (!defined $Math::BigInt::VERSION) {
1868 eval { require Math::BigInt; Math::BigInt->import(try=>'GMP,Pari'); 1; }
1869 or do { croak "Cannot load Math::BigInt"; };
1870 }
1871
1872 $n = Math::BigInt->new("$n") unless ref($n) eq 'Math::BigInt';
1873 my $nm1 = $n-1;
1874 my $A = $nm1->copy->bone; # factored part
1875 my $B = $nm1->copy; # unfactored part
1876 my @factors = (2);
1877 croak "BLS75 error: n-1 not even" unless $nm1->is_even();
1878 while ($B->is_even) { $B /= 2; $A *= 2; }
1879 foreach my $f (3,5,7,11,13,17,19,23,29,31,37,41,43,47,53,59,61,67,71,73,79) {
1880 last if $f*$f > $B;
1881 if (($B % $f) == 0) {
1882 push @factors, $f;
1883 do { $B /= $f; $A *= $f; } while (($B % $f) == 0);
1884 }
1885 }
1886 my @nstack;
1887 # nstack should only hold composites
1888 if ($B == 1) {
1889 # Completely factored. Nothing.
1890 } elsif (Math::Prime::Util::is_prob_prime($B)) {
1891 push @factors, $B;
1892 $A *= $B; $B /= $B; # completely factored already
1893 } else {
1894 push @nstack, $B;
1895 }
1896 while (@nstack) {
1897 my ($s,$r) = $B->copy->bdiv($A->copy->bmul(2));
1898 my $fpart = ($A+1) * (2*$A*$A + ($r-1) * $A + 1);
1899 last if $n < $fpart;
1900
1901 my $m = pop @nstack;
1902 # Don't use bignum if it has gotten small enough.
1903 $m = int($m->bstr) if ref($m) eq 'Math::BigInt' && $m <= ~0;
1904 # Try to find factors of m, using the default set of factor subs.
1905 my @ftry;
1906 $_holf_r = 1;
1907 foreach my $sub (@_fsublist) {
1908 last if scalar @ftry >= 2;
1909 @ftry = $sub->($m);
1910 }
1911 # If we couldn't find a factor, skip it.
1912 next unless scalar @ftry > 1;
1913 # Process each factor
1914 foreach my $f (@ftry) {
1915 croak "Invalid factoring" if $f == 1 || $f == $m || ($B%$f) != 0;
1916 if (Math::Prime::Util::is_prob_prime($f)) {
1917 push @factors, $f;
1918 do { $B /= $f; $A *= $f; } while (($B % $f) == 0);
1919 } else {
1920 push @nstack, $f;
1921 }
1922 }
1923 }
1924 # Just in case:
1925 foreach my $f (@factors) {
1926 while (($B % $f) == 0) {
1927 $B /= $f; $A *= $f;
1928 }
1929 }
1930 { # remove duplicate factors and make a sorted array of bigints
1931 my %uf = map { $_ => 1 } @factors;
1932 @factors = sort {$a<=>$b} map { Math::BigInt->new("$_") } keys %uf;
1933 }
1934 # Did we factor enough?
1935 my ($s,$r) = $B->copy->bdiv($A->copy->bmul(2));
1936 my $fpart = ($A+1) * (2*$A*$A + ($r-1) * $A + 1);
1937 return (1,[]) if $n >= $fpart;
1938 # Check we didn't mess up
1939 croak "BLS75 error: $A * $B != $nm1" unless $A*$B == $nm1;
1940 croak "BLS75 error: $A not even" unless $A->is_even();
1941 croak "BLS75 error: A and B not coprime" unless Math::BigInt::bgcd($A, $B)==1;
1942
1943 my $rtest = $r*$r - 8*$s;
1944 my $rtestroot = $rtest->copy->bsqrt;
1945 return @composite if $s != 0 && ($rtestroot*$rtestroot) == $rtest;
1946
1947 my @fac_proofs;
1948 my @as;
1949 foreach my $f (@factors) {
1950 my $success = 0;
1951 foreach my $a (2 .. 10000) {
1952 my $ap = Math::BigInt->new("$a");
1953 next unless $ap->copy->bmodpow($nm1, $n) == 1;
1954 next unless Math::BigInt::bgcd($ap->copy->bmodpow($nm1/$f, $n)->bsub(1), $n) == 1;
1955 push @as, $a;
1956 $success = 1;
1957 last;
1958 }
1959 return @composite unless $success;
1960 my ($isp, $fproof) = Math::Prime::Util::is_provable_prime_with_cert($f);
1961 if ($isp != 2) {
1962 carp "could not prove primality of $n.\n";
1963 return (1, []);
1964 }
1965 push @fac_proofs, (scalar @$fproof == 1) ? $fproof->[0] : $fproof;
1966 }
1967 my @proof = ($n, "n-1", [@fac_proofs], [@as]);
1968 return (2, [@proof]);
1969 }
14521970
14531971
14541972 my $_const_euler = 0.57721566490153286060651209008240243104215933593992;
14711989 my $wantbf = 0;
14721990 my $xdigits = 17;
14731991 if (defined $bignum::VERSION || ref($x) =~ /^Math::Big/) {
1474 if (!defined $MATH::BigFloat::VERSION) {
1992 if (!defined $Math::BigFloat::VERSION) {
14751993 eval { require Math::BigFloat; Math::BigFloat->import(); 1; }
14761994 or do { croak "Cannot load Math::BigFloat "; }
14771995 }
15812099 my $wantbf = 0;
15822100 my $xdigits = 17;
15832101 if (defined $bignum::VERSION || ref($x) =~ /^Math::Big/) {
1584 if (!defined $MATH::BigFloat::VERSION) {
2102 if (!defined $Math::BigFloat::VERSION) {
15852103 eval { require Math::BigFloat; Math::BigFloat->import(); 1; }
15862104 or do { croak "Cannot load Math::BigFloat "; }
15872105 }
17132231 my $wantbf = 0;
17142232 my $xdigits = 17;
17152233 if (defined $bignum::VERSION || ref($x) =~ /^Math::Big/) {
1716 if (!defined $MATH::BigFloat::VERSION) {
2234 if (!defined $Math::BigFloat::VERSION) {
17172235 eval { require Math::BigFloat; Math::BigFloat->import(); 1; }
17182236 or do { croak "Cannot load Math::BigFloat "; }
17192237 }
18042322 my $wantbf = 0;
18052323 my $xdigits = 17;
18062324 if (defined $bignum::VERSION || ref($x) =~ /^Math::Big/) {
1807 if (!defined $MATH::BigFloat::VERSION) {
2325 if (!defined $Math::BigFloat::VERSION) {
18082326 eval { require Math::BigFloat; Math::BigFloat->import(); 1; }
18092327 or do { croak "Cannot load Math::BigFloat "; }
18102328 }
19002418
19012419 =head1 VERSION
19022420
1903 Version 0.14
2421 Version 0.26
19042422
19052423
19062424 =head1 SYNOPSIS
20602578 count of primes between the ranges (e.g. C<(13,17)> returns 2, C<14,17>
20612579 and C<13,16> return 1, and C<14,16> returns 0).
20622580
2063 The Lehmer method is used for large values, which speeds up results greatly.
2581 The Lehmer method is used for large values, which greatly speeds up results.
20642582
20652583
20662584 =head2 nth_prime
21092627 prime using the AKS primality test. This code is included for completeness
21102628 and as an example, but is impractically slow.
21112629
2630 =head2 primality_proof_lucas
2631
2632 Given a positive number C<n> as input, performs a full factorization of C<n-1>,
2633 then attempts a Lucas test on the result. A Pratt-style certificate is
2634 returned. Note that if the input is composite, this will take a B<very> long
2635 time to return.
2636
2637 =head2 primality_proof_bls75
2638
2639 Given a positive number C<n> as input, performs a partial factorization of
2640 C<n-1>, then attempts a proof using theorem 5 of Brillhart, Lehmer, and
2641 Selfridge's 1975 paper. This can take a long time to return if given a
2642 composite, though it should not be anywhere near as long as the Lucas test.
2643
21122644
21132645 =head1 UTILITY FUNCTIONS
21142646
21602692
21612693 my @factors = fermat_factor($n);
21622694
2163 Currently unimplementated in Perl.
2695 Produces factors, not necessarily prime, of the positive number input. This
2696 uses Fermat's classic algorithm, without any special improvements (e.g.
2697 Lehman or McKee). This is typically slow and usually L</holf_factor> will
2698 be faster.
21642699
21652700 =head2 holf_factor
21662701
21782713
21792714 my @factors = squfof_factor($n);
21802715
2181 Currently unimplementated in Perl.
2716 Currently not implemented in Perl.
21822717
21832718 =head2 prho_factor
21842719
21852720 =head2 pbrent_factor
2186
2187 =head2 pminus1_factor
21882721
21892722 my @factors = prho_factor($n);
21902723
21922725 my @factors = prho_factor($n, 1000);
21932726
21942727 Produces factors, not necessarily prime, of the positive number input. An
2195 optional number of rounds can be given as a second parameter. These attempt
2196 to find a single factor using one of the probabilistic algorigthms of
2197 Pollard Rho, Brent's modification of Pollard Rho, or Pollard's C<p - 1>.
2198 These are more specialized algorithms usually used for pre-factoring very
2199 large inputs, or checking very large inputs for naive mistakes. If the
2728 optional number of rounds can be given as a second parameter. An optional
2729 additive factor may be given as a third parameter (default 3). These methods
2730 attempt to find a single factor using Pollard's Rho algorithm (or Brent's
2731 alternative which is typically faster).
2732 These methods can quickly find small factors. If the
22002733 input is prime or they run out of rounds, they will return the single
22012734 input value. On some inputs they will take a very long time, while on
22022735 others they succeed in a remarkably short time.
2736
2737 =head2 pminus1_factor
2738
2739 my @factors = pminus1_factor($n);
2740 my @factors = pminus1_factor($n, 1_000); # set B1 smoothness
2741 my @factors = pminus1_factor($n, 1_000, 50_000); # set B1 and B2
2742
2743 Produces factors, not necessarily prime, of the positive number input. This
2744 is Pollard's C<p-1> method, using two stages. The default with no B1 or B2
2745 given is to run with increasing B1 factors.
2746 This method can rapidly find a factor C<p> of C<n> where C<p-1> is smooth
2747 (i.e. C<p-1> has no large factors).
2748
2749 =head2 ecm_factor
2750
2751 my @factors = ecm_factor($n);
2752 my @factors = ecm_factor($n, 1000); # B1 = B2 = 1000, curves = 10
2753 my @factors = ecm_factor($n, 1000, 5000, 10); # Set B1, B2, and ncurves
2754
2755 Given a positive number input, tries to discover a factor using ECM. The
2756 resulting array will contain either two factors (it succeeded) or the original
2757 number (no factor was found). In either case, multiplying @factors yields the
2758 original input. The B1 and B2 smoothness factors for stage 1 and stage 2
2759 respectively may be supplied, as can be number of curves to try.
2760
22032761
22042762
22052763
23502908
23512909 =head1 COPYRIGHT
23522910
2353 Copyright 2012 by Dana Jacobsen E<lt>dana@acm.orgE<gt>
2911 Copyright 2012-2013 by Dana Jacobsen E<lt>dana@acm.orgE<gt>
23542912
23552913 This program is free software; you can redistribute it and/or modify it under the same terms as Perl itself.
23562914
+607
-138
lib/Math/Prime/Util.pm less more
55
66 BEGIN {
77 $Math::Prime::Util::AUTHORITY = 'cpan:DANAJ';
8 $Math::Prime::Util::VERSION = '0.25';
8 $Math::Prime::Util::VERSION = '0.26';
99 }
1010
1111 # parent is cleaner, and in the Perl 5.10.1 / 5.12.0 core, but not earlier.
1414 our @EXPORT_OK =
1515 qw( prime_get_config prime_set_config
1616 prime_precalc prime_memfree
17 is_prime is_prob_prime is_provable_prime
17 is_prime is_prob_prime is_provable_prime is_provable_prime_with_cert
18 prime_certificate verify_prime
1819 is_strong_pseudoprime is_strong_lucas_pseudoprime
1920 is_aks_prime
2021 miller_rabin
2425 prime_count_lower prime_count_upper prime_count_approx
2526 nth_prime nth_prime_lower nth_prime_upper nth_prime_approx
2627 random_prime random_ndigit_prime random_nbit_prime
27 random_strong_prime random_maurer_prime
28 random_strong_prime random_maurer_prime random_maurer_prime_with_cert
2829 primorial pn_primorial consecutive_integer_lcm
2930 factor all_factors
3031 moebius mertens euler_phi jordan_totient exp_mangoldt
832833 }
833834
834835 sub random_maurer_prime {
836 my ($n, $cert) = random_maurer_prime_with_cert(@_);
837 croak "maurer prime $n failed certificate verification!"
838 unless verify_prime(@$cert);
839 return $n;
840 }
841
842 sub random_maurer_prime_with_cert {
835843 my($k) = @_;
836844 _validate_positive_integer($k, 2);
845 my @cert;
837846 if ($] < 5.008 && $_Config{'maxbits'} > 32) {
838 return random_nbit_prime($k) if $k <= 49;
847 if ($k <= 49) {
848 my $n = random_nbit_prime($k);
849 return ($n, [$n]);
850 }
839851 croak "Random Maurer not supported on old Perl";
840852 }
841853
844856 # returning 2. This should be reasonably fast to ~128 bits with MPU::GMP.
845857 my $p0 = $_Config{'maxbits'};
846858
847 return random_nbit_prime($k) if $k <= $p0;
859 if ($k <= $p0) {
860 my $n = random_nbit_prime($k);
861 my ($isp, $cert) = is_provable_prime_with_cert($n);
862 croak "small nbit prime could not be proven" if $isp != 2;
863 return ($n, $cert);
864 }
848865
849866 if (!defined $Math::BigInt::VERSION) {
850867 eval { require Math::BigInt; Math::BigInt->import(try=>'GMP,Pari'); 1; }
875892 }
876893
877894 # I've seen +0, +1, and +2 here. Maurer uses +0. Menezes uses +1.
878 my $q = random_maurer_prime( ($r * $k)->bfloor + 1 );
895 my ($q, $certref) = random_maurer_prime_with_cert( ($r * $k)->bfloor + 1 );
879896 $q = Math::BigInt->new("$q") unless ref($q) eq 'Math::BigInt';
880897 my $I = Math::BigInt->new(2)->bpow($k-2)->bdiv($q)->bfloor->as_int();
881898 print "B = $B r = $r k = $k q = $q I = $I\n" if $verbose && $verbose != 3;
953970 # 2) make history by finding the first known BPSW pseudo-prime
954971 croak "Maurer prime $n=2*$R*$q+1 failed BPSW" unless is_prob_prime($n);
955972
956 return $n;
973 # Build up cert, knowing n-1 = 2*q*R, q > sqrt(n).
974 # We'll need to find the right a value for the factor 2.
975 foreach my $f2a (2 .. 200) {
976 $a = Math::BigInt->new($f2a);
977 next unless $a->copy->bmodpow($n-1, $n) == 1;
978 next unless Math::BigInt::bgcd($a->copy->bmodpow(($n-1)/2, $n)->bsub(1), $n) == 1;
979 @cert = ("$n", "n-1", [2, [@$certref]], [$f2a, $try_a]);
980 return ($n, \@cert);
981 }
957982 }
958983 # Didn't pass the selected a values. Try another R.
959984 }
15981623 return ($n <= 18446744073709551615) ? 2 : 1;
15991624 }
16001625
1626 # Return just the non-cert portion.
16011627 sub is_provable_prime {
16021628 my($n) = @_;
16031629 return 0 if defined $n && $n < 2;
16041630 _validate_positive_integer($n);
16051631
1606 # Shortcut some of the calls.
16071632 return _XS_is_prime($n) if ref($n) ne 'Math::BigInt' && $n <= $_XS_MAXVAL;
1608 return Math::Prime::Util::GMP::is_provable_prime($n) if $_HAVE_GMP
1609 && defined &Math::Prime::Util::GMP::is_provable_prime;
1610
1611 my $is_prob_prime = is_prob_prime($n);
1612 return $is_prob_prime unless $is_prob_prime == 1;
1613
1614 # At this point we know it is almost certainly a prime, but we need to
1615 # prove it. We should do ECPP or APR-CL now, or failing that, do the
1616 # Brillhart-Lehmer-Selfridge test, or Pocklington-Lehmer. Until those
1617 # are written here, we'll do a Lucas test, which is super simple but may
1618 # be very slow. We have AKS code, but it's insanely slow.
1633 return Math::Prime::Util::GMP::is_provable_prime($n)
1634 if $_HAVE_GMP && defined &Math::Prime::Util::GMP::is_provable_prime;
1635
1636 my ($is_prime, $cert) = is_provable_prime_with_cert($n);
1637 return $is_prime;
1638 }
1639
1640 # Return just the cert portion.
1641 sub prime_certificate {
1642 my($n) = @_;
1643 my ($is_prime, $cert) = is_provable_prime_with_cert($n);
1644 return @$cert;
1645 }
1646
1647
1648 sub is_provable_prime_with_cert {
1649 my($n) = @_;
1650 return 0 if defined $n && $n < 2;
1651 _validate_positive_integer($n);
1652
1653 # Set to 0 if you want the proof to go down to 11.
1654 if (1) {
1655 if (ref($n) ne 'Math::BigInt' && $n <= $_XS_MAXVAL) {
1656 my $isp = _XS_is_prime($n);
1657 return ($isp == 2) ? ($isp, [$n]) : ($isp, []);
1658 }
1659 if ($_HAVE_GMP && defined &Math::Prime::Util::GMP::is_provable_prime_with_cert) {
1660 return Math::Prime::Util::GMP::is_provable_prime_with_cert($n);
1661 }
1662
1663 my $isp = is_prob_prime($n);
1664 if ($isp != 1) {
1665 return ($isp == 2) ? ($isp, [$n]) : ($isp, []);
1666 }
1667 } else {
1668 if ($n <= 10) {
1669 if ($n==2||$n==3||$n==5||$n==7) {
1670 return (2, [$n]);
1671 }
1672 return 0;
1673 }
1674 }
1675
1676 # Choice of methods for proof:
1677 # ECPP needs a fair bit of programming work
1678 # APRCL needs a lot of programming work
1679 # BLS75 combo Corollary 11 of BLS75. Trial factor n-1 and n+1 to B, find
1680 # factors F1 of n-1 and F2 of n+1. Quit when:
1681 # B > (N/(F1*F1*(F2/2)))^1/3 or B > (N/((F1/2)*F2*F2))^1/3
1682 # BLS75 n+1 Requires factoring n+1 to (n/2)^1/3 (theorem 19)
1683 # BLS75 n-1 Requires factoring n-1 to (n/2)^1/3 (theorem 5 or 7)
1684 # Pocklington Requires factoring n-1 to n^1/2 (BLS75 theorem 4)
1685 # Lucas Easy, requires factoring of n-1 (BLS75 theorem 1)
1686 # AKS horribly slow
16191687 # See http://primes.utm.edu/prove/merged.html or other sources.
16201688
1621 # It shouldn't be possible to get here without BigInt already loaded.
1689 #my ($isp, $pref) = Math::Prime::Util::PP::primality_proof_lucas($n);
1690 my ($isp, $pref) = Math::Prime::Util::PP::primality_proof_bls75($n);
1691 carp "proved $n is not prime\n" if !$isp;
1692 return ($isp, $pref);
1693 }
1694
1695
1696 sub verify_prime {
1697 my @pdata = @_;
1698 my $verbose = $_Config{'verbose'};
1699
1700 # Empty input = no certificate = not prime
1701 return 0 if scalar @pdata == 0;
1702
1703 # Handle case of being handed a reference to the certificate.
1704 @pdata = @{$pdata[0]} if scalar @pdata == 1 && ref($pdata[0]) eq 'ARRAY';
1705
1706 my $n = shift @pdata;
1707 if (length($n) == 1) {
1708 return 1 if $n =~ /^[2357]$/;
1709 print "primality fail: $n tiny and not prime\n" if $verbose;
1710 return 0;
1711 }
1712
16221713 if (!defined $Math::BigInt::VERSION) {
16231714 eval { require Math::BigInt; Math::BigInt->import(try=>'GMP,Pari'); 1; }
16241715 or do { croak "Cannot load Math::BigInt"; };
16251716 }
1626 my $nm1 = $n-1;
1627 my @factors = factor($nm1);
1628 # Remember that we have to prove the primality of every factor.
1629 if ( (scalar grep { is_provable_prime($_) != 2 } @factors) > 0) {
1630 carp "could not prove primality of $n.\n";
1717 $n = Math::BigInt->new("$n") if ref($n) ne 'Math::BigInt';
1718 if ($n->is_even) {
1719 print "primality fail: $n even\n" if $verbose;
1720 return 0;
1721 }
1722
1723 my $method = (scalar @pdata > 0) ? shift @pdata : 'BPSW';
1724
1725 if ($method eq 'BPSW') {
1726 if ($n > Math::BigInt->new("18446744073709551615")) {
1727 print "primality fail: $n too large for BPSW proof\n" if $verbose;
1728 return 0;
1729 }
1730 my $bpsw = 0;
1731 if ($_HAVE_GMP) {
1732 $bpsw = Math::Prime::Util::GMP::is_prob_prime($n);
1733 } else {
1734 $bpsw = Math::Prime::Util::PP::miller_rabin($n, 2)
1735 && Math::Prime::Util::PP::is_strong_lucas_pseudoprime($n);
1736 }
1737 if (!$bpsw) {
1738 print "primality fail: BPSW indicated $n is composite\n" if $verbose;
1739 return 0;
1740 }
1741 print "primality success: $n by BPSW\n" if $verbose > 1;
16311742 return 1;
16321743 }
16331744
1634 for (my $a = 2; $a < $nm1; $a++) {
1635 my $ap = Math::BigInt->new("$a");
1636 # 1. a^(n-1) = 1 mod n.
1637 next if $ap->copy->bmodpow($nm1, $n) != 1;
1638 # 2. a^((n-1)/f) != 1 mod n for all f.
1639 next if (scalar grep { $_ == 1 }
1640 map { $ap->copy->bmodpow(int($nm1/$_),$n); }
1641 @factors) > 0;
1642 return 2;
1643 }
1644 carp "proved $n is not prime\n";
1745 if ($method eq 'Pratt' || $method eq 'Lucas') {
1746 # Based on Lucas primality test, which requires full n-1 factorization.
1747 if (scalar @pdata != 2 || (ref($pdata[0]) ne 'ARRAY') || (ref($pdata[1]) eq 'ARRAY')) {
1748 warn "verify_prime: incorrect Pratt format, must have factors and a value\n";
1749 return 0;
1750 }
1751 my @factors = @{shift @pdata};
1752 my $a = Math::BigInt->new(shift @pdata);
1753 my $nm1 = $n - 1;
1754 my $B = $nm1; # Unfactored part
1755
1756 my @prime_factors;
1757 my %factors_seen;
1758 foreach my $farray (@factors) {
1759 my $f;
1760 if (ref($farray) eq 'ARRAY') {
1761 $f = Math::BigInt->new("$farray->[0]");
1762 return 0 unless verify_prime(@$farray);
1763 } else {
1764 $f = $farray;
1765 return 0 unless verify_prime($f);
1766 }
1767 next if defined $factors_seen{"$f"}; # repeated factors
1768 if (($B % $f) != 0) {
1769 print "primality fail: given factor $f does not divide $nm1\n" if $verbose;
1770 return 0;
1771 }
1772 while (($B % $f) == 0) {
1773 $B /= $f;
1774 }
1775 push @prime_factors, $f;
1776 $factors_seen{"$f"} = 1;
1777 }
1778 croak "Pratt error: n-1 not completely factored" unless $B == 1;
1779
1780 # 1. a must be co-prime to n.
1781 if (Math::BigInt::bgcd($a, $n) != 1) {
1782 print "primality fail: a and n not coprime\n" if $verbose;
1783 return 0;
1784 }
1785 # 2. n is a psp base a
1786 if ($a->copy->bmodpow($nm1, $n) != 1) {
1787 print "primality fail: n is not a psp base a\n" if $verbose;
1788 return 0;
1789 }
1790 # 3. For each factor f of n-1, a^((n-1)/f) != 1 mod n
1791 foreach my $f (@prime_factors) {
1792 if ($a->copy->bmodpow(int($nm1/$f),$n) == 1) {
1793 print "primality fail: factor f fails a^((n-1)/f) != 1 mod n\n" if $verbose;
1794 return 0;
1795 }
1796 }
1797 print "primality success: $n by Lucas test\n" if $verbose > 1;
1798 return 1;
1799 }
1800
1801 if ($method eq 'n-1') {
1802 # BLS75 or generalized Pocklington
1803 # http://www.ams.org/journals/mcom/1975-29-130/S0025-5718-1975-0384673-1/S0025-5718-1975-0384673-1.pdf
1804 if (scalar @pdata != 2 || (ref($pdata[0]) ne 'ARRAY') || (ref($pdata[1]) ne 'ARRAY')) {
1805 warn "verify_prime: incorrect n-1 format, must have factors and a values\n";
1806 return 0;
1807 }
1808 my @factors = @{shift @pdata};
1809 my @as = @{shift @pdata};
1810 if ($#factors != $#as) {
1811 warn "verify_prime: incorrect n-1 format, must have a value for each factor\n";
1812 return 0;
1813 }
1814
1815 my $nm1 = $n - 1;
1816 my $A = $n-$n+1; # Factored part (F_1 in BLS paper)
1817 my $B = $nm1; # Unfactored part (R_1 in BLS paper)
1818
1819 my @prime_factors;
1820 my @pfas;
1821 my %factors_seen;
1822 foreach my $farray (@factors) {
1823 my $f;
1824 my $a = shift @as;
1825 if (ref($farray) eq 'ARRAY') {
1826 $f = Math::BigInt->new("$farray->[0]");
1827 return 0 unless verify_prime(@$farray);
1828 } else {
1829 $f = Math::BigInt->new("$farray");
1830 return 0 unless verify_prime($f);
1831 }
1832 next if defined $factors_seen{"$f"}; # repeated factors
1833 if (($B % $f) != 0) {
1834 print "primality fail: given factor $f does not divide $nm1\n" if $verbose;
1835 return 0;
1836 }
1837 while (($B % $f) == 0) {
1838 $B /= $f;
1839 $A *= $f;
1840 }
1841 push @prime_factors, $f;
1842 push @pfas, $a;
1843 $factors_seen{"$f"} = 1;
1844 }
1845 croak "BLS75 error: $A * $B != $nm1" unless $A*$B == $nm1;
1846
1847 # The theorems state that A is the even portion, so we are requiring 2 be
1848 # listed as a factor.
1849 if ($A->is_odd) {
1850 print "primality fail: 2 must be included as a factor" if $verbose;
1851 return 0;
1852 }
1853
1854 # TODO: consider: if B=1 and a single a is given, then Lucas test.
1855
1856 if (Math::BigInt::bgcd($A, $B) != 1) {
1857 print "primality fail: A and B not coprime\n" if $verbose;
1858 return 0;
1859 }
1860 # Theorem 5, m = 1, page 624
1861 {
1862 my ($s,$r) = $B->copy->bdiv($A->copy->bmul(2));
1863 my $fpart = ($A+1) * (2*$A*$A + ($r-1) * $A + 1);
1864 if ($n >= $fpart) {
1865 print "primality fail: not enough factors\n" if $verbose;
1866 return 0;
1867 }
1868 my $rtest = $r*$r - 8*$s;
1869 my $rtestroot = $rtest->copy->bsqrt;
1870 if ($s != 0 && ($rtestroot*$rtestroot) == $rtest) {
1871 print "primality fail: BLS75 theorem 5: s=$s, r=$r indicates composite\n" if $verbose;
1872 return 0;
1873 }
1874 }
1875 # Now verify (I), page 623
1876 foreach my $i (0 .. $#prime_factors) {
1877 my $f = $prime_factors[$i];
1878 my $a = Math::BigInt->new("$pfas[$i]");
1879 if ($a->copy->bmodpow($nm1, $n) != 1 ||
1880 Math::BigInt::bgcd($a->copy->bmodpow($nm1/$f, $n)->bsub(1), $n) != 1) {
1881 print "primality fail: BLS75 factor=$f, a=$a failed.\n" if $verbose;
1882 return 0;
1883 }
1884 }
1885 print "primality success: $n by BLS75 theorem 5\n" if $verbose > 1;
1886 return 1;
1887 }
1888
1889 if ($method eq 'ECPP' || $method eq 'AGKM') {
1890 # EC cert: Atkin-Morain etc.
1891 # Normally we'd have the q values set up recursively, but to follow the
1892 # standard trend, we have this set up as a list:
1893 # n, "AGKM", [n,a,b,m,q,P], [n1,a,b,m,q,P], ...
1894 #
1895 # Examples:
1896 # (100000000000000000039, "AGKM", [100000000000000000039, 31484432173069852672, 39553474583282556928, 100000000014867206541, 539348143913549, [39164891430400385024,86449249723524901718]])
1897 # (677826928624294778921, "AGKM", [677826928624294778921, 404277700094248015180, 599134911995823048257, 677826928656744857936, 104088901820753203, [2293544533, 356794037129589115041]])
1898 # Ux,Uy should be 600992528322000913770, 206075883056439332684
1899 # Vx,Vy should be 0, 1
1900 if (scalar @pdata < 1) {
1901 warn "verify_prime: incorrect AGKM format\n";
1902 return 0;
1903 }
1904 my ($ni, $a, $b, $m, $P);
1905 my ($qval, $q) = ($n, $n);
1906 foreach my $block (@pdata) {
1907 if (ref($block) ne 'ARRAY' || scalar @$block != 6) {
1908 warn "verify_prime: incorrect AGKM block format\n";
1909 return 0;
1910 }
1911 if ($block->[0] != $q) {
1912 warn "verify_prime: incorrect AGKM block format: block n != q\n";
1913 return 0;
1914 }
1915 ($ni, $a, $b, $m, $qval, $P) = @$block;
1916 $q = ref($qval) eq 'ARRAY' ? $qval->[0] : $qval;
1917 if (ref($P) ne 'ARRAY' || scalar @$P != 2) {
1918 warn "verify_prime: incorrect AGKM block point format\n";
1919 return 0;
1920 }
1921 my($Px, $Py) = @$P;
1922 $ni = $n->copy->bzero->badd("$ni") unless ref($ni) eq 'Math::BigInt';
1923 $a = $n->copy->bzero->badd("$a") unless ref($a) eq 'Math::BigInt';
1924 $b = $n->copy->bzero->badd("$b") unless ref($b) eq 'Math::BigInt';
1925 $m = $n->copy->bzero->badd("$m") unless ref($m) eq 'Math::BigInt';
1926 $q = $n->copy->bzero->badd("$q") unless ref($q) eq 'Math::BigInt';
1927 $Px = $n->copy->bzero->badd("$Px") unless ref($Px) eq 'Math::BigInt';
1928 $Py = $n->copy->bzero->badd("$Py") unless ref($Py) eq 'Math::BigInt';
1929 if (Math::BigInt::bgcd($ni, 6) != 1) {
1930 warn "verify_prime: AGKM block n '$ni' is divisible by 2 or 3\n";
1931 return 0;
1932 }
1933 my $c = $a*$a*$a * 4 + $b*$b * 27;
1934 if (Math::BigInt::bgcd($c, $ni) != 1) {
1935 warn "verify_prime: AGKM block gcd 4a^3+27b^2,n incorrect\n";
1936 return 0;
1937 }
1938 if ($q <= $ni->copy->broot(4)->badd(1)->bpow(2)) {
1939 warn "verify_prime: AGKM block q is too small\n";
1940 return 0;
1941 }
1942 # Final check, check that we've got a bound above and below (Hasse)
1943 if (!defined $Math::Prime::Util::ECAffinePoint::VERSION) {
1944 eval { require Math::Prime::Util::ECAffinePoint; 1; }
1945 or do { croak "Cannot load Math::Prime::Util::ECAffinePoint"; };
1946 }
1947 my $ECP = Math::Prime::Util::ECAffinePoint->new($a, $b, $ni, $Px, $Py);
1948 # Compute U = (m/q)P, check U != point at infinity
1949 $ECP->mul( $m->copy->bdiv($q)->as_int );
1950 if ($ECP->is_infinity) {
1951 warn "verify_prime: AGKM point does not multiply correctly.\n";
1952 return 0;
1953 }
1954 # Compute V = qU, check V = point at infinity
1955 $ECP->mul( $q );
1956 if (! $ECP->is_infinity) {
1957 warn "verify_prime: AGKM point does not multiply correctly.\n";
1958 return 0;
1959 }
1960 }
1961 # Check primality of last q
1962 return 0 unless verify_prime($qval);
1963
1964 print "primality success: $n by A-K-G-M elliptic curve\n" if $verbose > 1;
1965 return 1;
1966 }
1967
1968 warn "verify_prime: Unknown method: '$method'.\n";
16451969 return 0;
16461970 }
16471971
19912315
19922316 =head1 VERSION
19932317
1994 Version 0.25
2318 Version 0.26
19952319
19962320
19972321 =head1 SYNOPSIS
21352459
21362460 =over 4
21372461
2138 =item primes.pl displays primes between start and end values, with many
2139 options for filtering (e.g. twin, safe, circular, good, lucky, etc.).
2140
2141 =item factor.pl operates similar to the GNU C<factor> program. It supports
2462 =item *
2463
2464 primes.pl displays primes between start and end values or expressions,
2465 with many options for filtering (e.g. twin, safe, circular, good, lucky,
2466 etc.). Use C<--help> to see all the options.
2467
2468 =item *
2469
2470 factor.pl operates similar to the GNU C<factor> program. It supports
21422471 bigint and expression inputs.
21432472
21442473 =back
22282557 test. If L<Math::Prime::Util::GMP> is installed, some quick primality proofs
22292558 are run on larger numbers, so will return 2 for many of those also.
22302559
2231 Also see the L</"is_prob_prime"> function, which will never do additional
2232 tests, and the L</"is_provable_prime"> function which will try very hard to
2560 Also see the L</is_prob_prime> function, which will never do additional
2561 tests, and the L</is_provable_prime> function which will try very hard to
22332562 return only 0 or 2 for any input.
22342563
22352564 For native precision numbers (anything smaller than C<2^64>, all three
22362565 functions are identical and use a deterministic set of Miller-Rabin tests.
2237 While L</"is_prob_prime"> and L</"is_prime"> return probable prime results
2566 While L</is_prob_prime> and L</is_prime> return probable prime results
22382567 for larger numbers, they use the strong Baillie-PSW test, which has had
22392568 no counterexample found since it was published in 1980 (though certainly they
22402569 exist).
23072636 In contrast, even primesieve using 12 cores would take over a week on this
23082637 same computer to determine C<Pi(10^16)>.
23092638
2310 Also see the function L</"prime_count_approx"> which gives a very good
2311 approximation to the prime count, and L</"prime_count_lower"> and
2312 L</"prime_count_upper"> which give tight bounds to the actual prime count.
2639 Also see the function L</prime_count_approx> which gives a very good
2640 approximation to the prime count, and L</prime_count_lower> and
2641 L</prime_count_upper> which give tight bounds to the actual prime count.
23132642 These functions return quickly for any input, including bigints.
23142643
23152644
24822811
24832812 Takes a positive number as input and returns back either 0 (composite),
24842813 2 (definitely prime), or 1 (probably prime). This gives it the same return
2485 values as C<is_prime> and C<is_prob_prime>.
2486
2487 The current implementation uses a Lucas test requiring a complete factorization
2488 of C<n-1>, which may not be possible in a reasonable amount of time. The GMP
2489 version uses the BLS (Brillhart-Lehmer-Selfridge) method, requiring C<n-1> to
2490 be factored to the cube root of C<n>, which is more likely to succeed and will
2491 usually take less time, but can still fail. Hence you should always test that
2492 the result is C<2> to ensure the prime is proven.
2814 values as L</is_prime> and L</is_prob_prime>.
2815
2816 The current implementation of both the Perl and GMP proofs is using theorem 5
2817 of BLS75 (Brillhart-Lehmer-Selfridge), requiring C<n-1> to be factored to
2818 C<(n/2)^(1/3))>. This takes less time than factoring to C<n^0.5> as required
2819 by the generalized Pocklington test or C<n-1> for the Lucas test. However it
2820 is possible a factor cannot be found in a reasonable amount of time, so you
2821 should always test that the result in C<2> to ensure it was proven.
2822
2823 A later implementation will use an ECPP test for larger inputs.
2824
2825
2826 =head2 prime_certificate
2827
2828 my @cert = prime_certificate($n);
2829 say verify_prime(@cert) ? "proven prime" : "not prime";
2830
2831 Given a positive integer C<n> as input, returns either an empty array (we could
2832 not prove C<n> prime) or an array representing a certificate of primality.
2833 This may be examined or given to L</verify_prime> for verification. The latter
2834 function contains the description of the format.
2835
2836
2837 =head2 is_provable_prime_with_cert
2838
2839 Given a positive integer as input, returns a two element array containing the
2840 result of L</is_provable_prime> and an array reference containing the primality
2841 certificate like L</prime_certificate>. The certificate will be an empty
2842 array reference if the result is not 2 (definitely prime).
2843
2844
2845 =head2 verify_prime
2846
2847 my @cert = prime_certificate($n);
2848 say verify_prime(@cert) ? "proven prime" : "not prime";
2849
2850 Given an array representing a certificate of primality, returns either 0 (not
2851 verified), or 1 (verified). The computations are all done using pure Perl
2852 Math::BigInt and should not be time consuming (the Pari or GMP backends will
2853 help with large inputs).
2854
2855 A certificate is an array holding an C<n-cert>. An C<n-cert> is one of:
2856
2857 n
2858 implies n,"BPSW"
2859
2860 n,"BPSW"
2861 the number n is small enough to be proven with BPSW. This
2862 currently means smaller than 2^64.
2863
2864 n,"Pratt",[n-cert, ...],a
2865 A Pratt certificate. We are given n, the method "Pratt" or
2866 "Lucas", a list of n-certs that indicate all the unique factors
2867 of n-1, and an 'a' value to be used in the Lucas primality test.
2868 The certificate passes if:
2869 1 all factor n-certs can be verified
2870 2 all n-certs are factors of n-1 and none are missing
2871 3 a is coprime to n
2872 4 a^(n-1) = 1 mod n
2873 5 a^((n-1)/f) != 1 mod n for each factor
2874
2875 n,"n-1",[n-cert, ...],[a,...]
2876 An n-1 certificate suitable for the generalized Pocklington or the
2877 BLS75 (Brillhart-Lehmer-Selfridge 1975, theorem 5) test. The
2878 proof is performed using BLS75 theorem 5 which requires n-1 to be
2879 factored up to (n/2)^1/3. If n-1 is factored to more than
2880 sqrt(n), then the conditions are identical to the generalized
2881 Pocklington test.
2882 The certificate passes if:
2883 1 all factor n-certs can be verified
2884 2 all factor n-certs are factors of n-1
2885 3 there must be a corresponding 'a' for each factor n-cert
2886 4 given A (the factored part of n-1), B = (n-1)/A (the
2887 unfactored part), s = int(B/(2A)), r = B-s*2A:
2888 - n < (A+1)(2*A*A+(r-a)A+a) [ n-1 factored to (n/2)^1/3 ]
2889 - s = 0 or r*r-8s not a perfect square
2890 - A and B are coprime
2891 5 for each pair (f,a) representing a factor n-cert and its 'a':
2892 - a^(n-1) = 1 mod n
2893 - gcd( a^((n-1)/f)-1, n ) = 1
2894
2895 n,"AGKM",[ec-block],[ec-block],...
2896 An Elliptic Curve certificate. We are given n, the method "AGKM"
2897 or "ECPP", and one or more 6-element blocks representing a
2898 standard ECPP or Atkin-Goldwasser-Kilian-Morain certificate.
2899 In its traditional form, it is non-recursive, with each q value
2900 being proved by successive blocks (this makes it easy to use for
2901 programs like Sage and GMP-ECPP). A q value is also allowed to
2902 be an n-cert, which allows an alternative proof for the last q.
2903 Every ec-block has 6 elements:
2904 N the N value this block proves prime if q is prime
2905 a value describing the elliptic curve to be used
2906 b value describing the elliptic curve to be used
2907 m order of the curve
2908 q a probable prime > (N^1/4+1)^2 (may be an n-cert)
2909 P a point [x,y] on the curve (affine coordinates)
2910 The certificate passes if:
2911 - the final q can be proved with BPSW.
2912 - for each block:
2913 - N is the same as the preceeding block's q
2914 - N is not divisible by 2 or 3
2915 - gcd( 4a^3 + 27b^2, N ) == 1;
2916 - q > (N^1/4+1)^2
2917 - U = (m/q)P is not the point at infinity
2918 - V = qU is the point at infinity
24932919
24942920
24952921 =head2 is_aks_prime
25312957 say "Mertens(10M) = ", mertens(10_000_000); # = 1037
25322958
25332959 Returns M(n), the Mertens function for a non-negative integer input. This
2534 function is defined as C<sum(moebius(1..n))>. This is a much more efficient
2535 solution for larger inputs. For example, computing Mertens(100M) takes:
2960 function is defined as C<sum(moebius(1..n))>, but calculated more efficiently
2961 for large inputs. For example, computing Mertens(100M) takes:
25362962
25372963 time approx mem
25382964 0.4s 0.1MB mertens(100_000_000)
25472973 which involves calculating all moebius values to C<n^1/2>, which in turn will
25482974 require prime sieving to C<n^1/4>.
25492975
2550 Various methods exist for this, using differing quantities of μ(n). The
2976 Various algorithms exist for this, using differing quantities of μ(n). The
25512977 simplest way is to efficiently sum all C<n> values. Benito and Varona (2008)
25522978 show a clever and simple method that only requires C<n/3> values. Deléglise
25532979 and Rivat (1996) describe a segmented method using only C<n^1/3> values. The
2554 current implementation does a simple non-segmented C<n^1/2> version of this.
2555 Kuznetsov (2011) gives an alternate method that he indicates is even faster.
2556 Lastly, one of the advanced prime count algorithms could be theoretically used
2557 to create a faster solution.
2980 current implementation does a simple non-segmented C<n^1/2> version of their
2981 method. Kuznetsov (2011) gives an alternate method that he indicates is even
2982 faster. Lastly, one of the advanced prime count algorithms could be
2983 theoretically used to create a faster solution.
25582984
25592985
25602986 =head2 euler_phi
25903016
25913017 say "exp(lambda($_)) = ", exp_mangoldt($_) for 1 .. 100;
25923018
2593 Returns Λ(n), the Mangoldt function (also known as von Mangoldt's function) for
2594 an integer value. It is equal to log p if n is prime or a power of a prime,
3019 Returns exp(Λ(n)), the exponential of the Mangoldt function (also known
3020 as von Mangoldt's function) for an integer value.
3021 It is equal to log p if n is prime or a power of a prime,
25953022 and 0 otherwise. We return the exponential so all results are integers.
2596 Hence the return value for C<exp_mangoldt> is:
2597 C<p> if C<n = p^m> for some prime C<p> and integer C<m E<gt>= 1>
3023 Hence the return value for C<exp_mangoldt> is:
3024
3025 p if n = p^m for some prime p and integer m >= 1
25983026 1 otherwise.
25993027
26003028
26043032
26053033 Returns θ(n), the first Chebyshev function for a non-negative integer input.
26063034 This is the sum of the logarithm of each prime where C<p E<lt>= n>. An
2607 alternate computation is as the logarithm of n primorial, hence:
2608
2609 use List::Util qw/sum/;
3035 alternate computation is as the logarithm of n primorial.
3036 Hence these functions:
3037
3038 use List::Util qw/sum/; use Math::BigFloat;
3039
26103040 sub c1a { 0+sum( map { log($_) } @{primes(shift)} ) }
2611 use Math::BigFloat;
26123041 sub c1b { Math::BigFloat->new(primorial(shift))->blog }
26133042
26143043 yield similar results, albeit slower and using more memory.
26243053 Another alternate computation is as the logarithm of lcm(1,2,...,n).
26253054 Hence these functions:
26263055
2627 use List::Util qw/sum/;
3056 use List::Util qw/sum/; use Math::BigFloat;
3057
26283058 sub c2a { 0+sum( map { log(exp_mangoldt($_)) } 1 .. shift ) }
2629 use Math::BigFloat;
26303059 sub c2b { Math::BigFloat->new(consecutive_integer_lcm(shift))->blog }
26313060
26323061 yield similar results, albeit slower and using more memory.
26583087
26593088 divisor_sum( $n, sub { my $d=shift; $d**5 * moebius($n/$d); } );
26603089
2661 though in the last case we have a function L</"jordan_totient"> to compute
3090 though in the last case we have a function L</jordan_totient> to compute
26623091 it more efficiently.
26633092
26643093 This function is useful for calculating things like aliquot sums, abundant
27423171 uniformity but results in many fewer bits of randomness being consumed as
27433172 well as being much faster.
27443173
2745 If an C<irand> function has been set via L</"prime_set_config">, it will be
3174 If an C<irand> function has been set via L</prime_set_config>, it will be
27463175 used to construct any ranged random numbers needed. The function should
27473176 return a uniformly random 32-bit integer, which is how the irand functions
27483177 exported by L<Math::Random::Secure>, L<Math::Random::MT>,
28433272 advances in factoring and attacks, for practical purposes, using large random
28443273 primes offer security equivalent to using strong primes.
28453274
2846 Similar to L</"random_nbit_prime">, the result will be a BigInt if the
3275 Similar to L</random_nbit_prime>, the result will be a BigInt if the
28473276 number of bits is greater than the native bit size. For better performance
28483277 with large bit sizes, install L<Math::Prime::Util::GMP>.
28493278
3279
28503280 =head2 random_maurer_prime
28513281
28523282 my $bigprime = random_maurer_prime(512);
28533283
28543284 Construct an n-bit provable prime, using the FastPrime algorithm of
28553285 Ueli Maurer (1995). This is the same algorithm used by L<Crypt::Primes>.
2856 Similar to L</"random_nbit_prime">, the result will be a BigInt if the
3286 Similar to L</random_nbit_prime>, the result will be a BigInt if the
28573287 number of bits is greater than the native bit size. For better performance
28583288 with large bit sizes, install L<Math::Prime::Util::GMP>.
28593289
28603290 The differences between this function and that in L<Crypt::Primes> are
28613291 described in the L</"SEE ALSO"> section.
3292
3293 Internally this additionally runs the BPSW probable prime test on every
3294 partial result, and constructs a primality certificate for the final
3295 result, which is verified. These add additional checks that the resulting
3296 value has been properly constructed.
3297
3298
3299 =head2 random_maurer_prime_with_cert
3300
3301 my($n, $cert_ref) = random_maurer_prime_with_cert(512)
3302
3303 As with L</random_maurer_prime>, but returns a two-element array containing
3304 the n-bit provable prime along with a primality certificate. The certificate
3305 is the same as produced by L</prime_certificate> or
3306 L</is_provable_prime_with_cert>, and can be parsed by L</verify_prime> or
3307 any other software that can parse the certificate (the "n-1" form is described
3308 in detail in L</verify_prime>).
28623309
28633310
28643311
29193366
29203367 Allows setting of some parameters. Currently the only parameters are:
29213368
2922 xs Allows turning off the XS code, forcing the Pure Perl code
2923 to be used. Set to 0 to disable XS, set to 1 to re-enable.
3369 xs Allows turning off the XS code, forcing the Pure Perl
3370 code to be used. Set to 0 to disable XS, set to 1 to
3371 re-enable. You probably will never want to do this.
3372
3373 gmp Allows turning off the use of L<Math::Prime::Util::GMP>,
3374 which means using Pure Perl code for big numbers. Set
3375 to 0 to disable GMP, set to 1 to re-enable.
29243376 You probably will never want to do this.
29253377
2926 gmp Allows turning off the use of L<Math::Prime::Util::GMP>,
2927 which means using Pure Perl code for big numbers. Set to
2928 0 to disable GMP, set to 1 to re-enable.
2929 You probably will never want to do this.
2930
2931 assume_rh Allows functions to assume the Riemann hypothesis is true
2932 if set to 1. This defaults to 0. Currently this setting
2933 only impacts prime count lower and upper bounds, but could
2934 later be applied to other areas such as primality testing.
2935 A later version may also have a way to indicate whether
2936 no RH, RH, GRH, or ERH is to be assumed.
2937
2938 irand Takes a code ref to an irand function returning a uniform
2939 number between 0 and 2**32-1. This will be used for all
2940 random number generation in the module.
3378 assume_rh Allows functions to assume the Riemann hypothesis is
3379 true if set to 1. This defaults to 0. Currently this
3380 setting only impacts prime count lower and upper
3381 bounds, but could later be applied to other areas such
3382 as primality testing. A later version may also have a
3383 way to indicate whether no RH, RH, GRH, or ERH is to
3384 be assumed.
3385
3386 irand Takes a code ref to an irand function returning a
3387 uniform number between 0 and 2**32-1. This will be
3388 used for all random number generation in the module.
29413389
29423390
29433391 =head1 FACTORING FUNCTIONS
29573405 a few rounds of Pollard's Rho, SQUFOF, Pollard's p-1, Hart's OLF, a long
29583406 run of Pollard's Rho, and finally trial division if anything survives. This
29593407 process is repeated for each non-prime factor. In practice, it is very rare
2960 to require more than the first Rho + SQUFOF to find a factor.
3408 to require more than the first Rho + SQUFOF to find a factor, and I have not
3409 seen anything go to the last step.
29613410
29623411 Factoring bigints works with pure Perl, and can be very handy on 32-bit
29633412 machines for numbers just over the 32-bit limit, but it can be B<very> slow
33083757
33093758 =over 4
33103759
3311 =item L<Math::Prime::FastSieve>
3312 This is the alternative module I use for basic functionality with small
3313 integers. It's fast and simple, and has a good set of features.
3314
3315 =item L<Math::Primality>
3316 This is the alternative module I use for primality testing on bigints.
3317
3318 =item L<Math::Pari>
3319 If you want the kitchen sink and can install it and handle using it. There
3320 are still some functions it doesn't do well (e.g. prime count and nth_prime).
3760 =item * L<Math::Prime::FastSieve> is the alternative module I use for basic
3761 functionality with small integers. It's fast and simple, and has a good
3762 set of features.
3763
3764 =item * L<Math::Primality> is the alternative module I use for primality
3765 testing on bigints.
3766
3767 =item * L<Math::Pari> if you want the kitchen sink and can install it and
3768 handle using it. There are still some functions it doesn't do well
3769 (e.g. prime count and nth_prime).
33213770
33223771 =back
33233772
33243773
33253774 L<Math::Prime::XS> has C<is_prime> and C<primes> functionality. There is no
3326 bigint support. The C<is_prime> function is a well-written trial division,
3775 bigint support. The C<is_prime> function uses well-written trial division,
33273776 meaning it is very fast for small numbers, but terribly slow for large
33283777 64-bit numbers. The prime sieve is an unoptimized non-segmented SoE which
33293778 which returns an array. It works well for 32-bit values, but speed and
33583807 What Crypt::Primes has that MPU does not is support for returning a generator.
33593808
33603809 L<Math::Factor::XS> calculates prime factors and factors, which correspond to
3361 the L</"factor"> and L<"all_factors"> functions of MPU. These functions do
3810 the L</factor> and L</all_factors> functions of MPU. These functions do
33623811 not support bigints. Both are implemented with trial division, meaning they
33633812 are very fast for really small values, but quickly become unusably slow
33643813 (factoring 19 digit semiprimes is over 700 times slower). It has additional
33753824 trial division). It supports bigints. Unfortunately it is extremely slow on
33763825 any input that isn't comprised entirely of small factors. Even 7 digit inputs
33773826 can take hundreds or thousands of times longer to factor than MPU or
3378 L<Math::Factor::XS>. 19-digit semiprimes will take hours vs. MPU's single
3827 L<Math::Factor::XS>. 19-digit semiprimes will take I<hours> vs. MPU's single
33793828 milliseconds.
33803829
33813830 L<Math::Factoring> is a placeholder module for bigint factoring. Version 0.02
33993848 while MPU was originally targeted to native integers, but both have added
34003849 better support for the other. The main differences are extra functionality
34013850 (MPU has more functions) and performance. With native integer inputs, MPU
3402 is generally much faster, especially with L</"prime_count">. For bigints,
3851 is generally much faster, especially with L</prime_count>. For bigints,
34033852 MPU is slower unless the L<Math::Prime::Util::GMP> module is installed, in
34043853 which case they have somewhat similar speeds. L<Math::Primality> also installs
34053854 a C<primes.pl> program, but it has much less functionality than the one
34113860 count, factor count, Euler totient, primorials, etc. Math::NumSeq is mainly
34123861 set up for accessing these values in order, rather than for arbitrary values,
34133862 though some sequences support that. The primary advantage I see is the
3414 uniform access mechanism for a <I>lot</I> of sequences. For those methods that
3415 overlap, MPU is usually much faster. Importantly, most all the sequences in
3863 uniform access mechanism for a I<lot> of sequences. For those methods that
3864 overlap, MPU is usually much faster. Importantly, most of the sequences in
34163865 Math::NumSeq are limited to 32-bit indices.
34173866
34183867 L<Math::Pari> supports a lot of features, with a great deal of overlap. In
34233872
34243873 =over 4
34253874
3426 =item isprime Similar to MPU's L<is_prob_prime> or L<is_prime> functions.
3875 =item isprime
3876
3877 Similar to MPU's L<is_prob_prime> or L<is_prime> functions.
34273878 MPU is deterministic for native integers, and uses a strong
34283879 BPSW test for bigints (with a quick primality proof tried as well). The
34293880 default version of Pari used by L<Math::Pari> (2.1.7) uses 10 random M-R
34333884 can take longer (though will be much faster than the BLS75 proof used in
34343885 MPU's L<is_provable_prime> routine).
34353886
3436 =item primepi Similar to MPU's L<prime_count> function. Pari uses a naive
3887 =item primepi
3888
3889 Similar to MPU's L<prime_count> function. Pari uses a naive
34373890 counting algorithm with its precalculated primes, so this is not very useful.
34383891
3439 =item primes Doesn't support ranges, requires bumping up the precalculated
3892 =item primes
3893
3894 Doesn't support ranges, requires bumping up the precalculated
34403895 primes for larger numbers, which means knowing in advance the upper limit
34413896 for primes. Support for numbers larger than 400M requires making with Pari
34423897 version 2.3.5. If that is used, sieving is about 2x faster than MPU, but
34433898 doesn't support segmenting.
34443899
3445 =item factorint Similar to MPU's L<factor> though with a different return (I
3900 =item factorint
3901
3902 Similar to MPU's L<factor> though with a different return (I
34463903 find the result value quite inconvenient to work with, but others may like
34473904 its vector of factor/exponent format). Slower than MPU for all 64-bit inputs
34483905 on an x86_64 platform, it may be faster for large values on other platforms.
34493906 Bigints are slightly faster with Math::Pari for "small" values, and MPU slows
34503907 down rapidly (the difference is noticeable with 30+ digit numbers).
34513908
3452 =item eulerphi Similar to MPU's L<euler_phi>. MPU is 2-5x faster for native
3909 =item eulerphi
3910
3911 Similar to MPU's L<euler_phi>. MPU is 2-5x faster for native
34533912 integers. There is also support for a range, which can be much more efficient.
34543913 Without L<Math::Prime::Util::GMP> installed, MPU is very slow with bigints.
34553914 With it installed, it is about 2x slower than Math::Pari.
34563915
3457 =item moebius Similar to MPU's L<moebius>. Comparisons are similar to
3458 C<eulerphi>.
3459
3460 =item sumdiv Similar to MPU's L<divisor_sum>. The standard sum (sigma_1) is
3916 =item moebius
3917
3918 Similar to MPU's L<moebius>. Comparisons are similar to C<eulerphi>.
3919
3920 =item sumdiv
3921
3922 Similar to MPU's L<divisor_sum>. The standard sum (sigma_1) is
34613923 very fast in MPU. Giving it a sub makes it much slower, and for numbers
3462 with very many factors, Pari is <I>much</I> faster.
3463
3464 =item eint1 Similar to MPU's L<ExponentialIntegral>.
3465
3466 =item zeta A more substantial version of MPU's L<RiemannZeta> function.
3924 with very many factors, Pari is I<much> faster.
3925
3926 =item eint1
3927
3928 Similar to MPU's L<ExponentialIntegral>.
3929
3930 =item zeta
3931
3932 A more feature-rich version MPU's L<RiemannZeta> function (supports negative
3933 and complex inputs).
34673934
34683935 =back
34693936
36154082 faster. Without the L<Math::Prime::Util::GMP> module, almost all actions on
36164083 numbers greater than native scalars will be much faster in Pari.
36174084
3618 The presentation here:
3619 L<http://math.boisestate.edu/~liljanab/BOISECRYPTFall09/Jacobsen.pdf>
4085 L<This slide presentation|http://math.boisestate.edu/~liljanab/BOISECRYPTFall09/Jacobsen.pdf>
36204086 has a lot of data on 64-bit and GMP factoring performance I collected in 2009.
36214087 Assuming you do not know anything about the inputs, trial division and
36224088 optimized Fermat or Lehman work very well for small numbers (<= 10 digits),
36294095 L<msieve|http://sourceforge.net/projects/msieve/>,
36304096 L<gmp-ecm|http://ecm.gforge.inria.fr/>,
36314097 L<GGNFS|http://sourceforge.net/projects/ggnfs/>,
3632 and L<Pari|http://pari.math.u-bordeaux.fr/>.
4098 and L<Pari|http://pari.math.u-bordeaux.fr/>. The latest yafu should cover most
4099 uses, with GGNFS likely only providing a benefit for numbers large enough to
4100 warrant distributed processing.
36334101
36344102
36354103 The primality proving algorithms leave much to be desired. If you have
3636 numbers larger than C<2^128>, I recommend Pari's C<isprime(n, 2)> which
3637 will run a fast APRCL test, or
3638 L<GMP-ECPP|http://http://sourceforge.net/projects/gmp-ecpp/>. Either one
4104 numbers larger than C<2^128>, I recommend C<isprime(n, 2)> from Pari 2.3+
4105 which will run a fast APRCL test,
4106 L<GMP-ECPP|http://http://sourceforge.net/projects/gmp-ecpp/>, or
4107 L<Primo|http://www.ellipsa.eu/>. Any of those
36394108 will be much faster than the Lucas or BLS algorithms used in MPU for large
3640 inputs.
4109 inputs. For very large numbers, Primo is the one to use.
36414110
36424111
36434112 =head1 AUTHORS
36484117 =head1 ACKNOWLEDGEMENTS
36494118
36504119 Eratosthenes of Cyrene provided the elegant and simple algorithm for finding
3651 the primes.
4120 primes.
36524121
36534122 Terje Mathisen, A.R. Quesada, and B. Van Pelt all had useful ideas which I
36544123 used in my wheel sieve.
36704139
36714140 =item *
36724141
3673 Pierre Dusart, "Estimates of Some Functions Over Primes without R.H.", preprint, 2010. L<http://arxiv.org/abs/1002.0442/>
4142 Pierre Dusart, "Estimates of Some Functions Over Primes without R.H.", preprint, 2010. Updates to the best non-RH bounds for prime count and nth prime. L<http://arxiv.org/abs/1002.0442/>
36744143
36754144 =item *
36764145
36784147
36794148 =item *
36804149
3681 Gabriel Mincu, "An Asymptotic Expansion", <Journal of Inequalities in Pure and Applied Mathematics>, v4, n2, 2003. A very readable account of Cipolla's 1902 nth prime approximation. L<http://www.emis.de/journals/JIPAM/images/153_02_JIPAM/153_02.pdf>
4150 Gabriel Mincu, "An Asymptotic Expansion", I<Journal of Inequalities in Pure and Applied Mathematics>, v4, n2, 2003. A very readable account of Cipolla's 1902 nth prime approximation. L<http://www.emis.de/journals/JIPAM/images/153_02_JIPAM/153_02.pdf>
36824151
36834152 =item *
36844153
37064175
37074176 =item *
37084177
3709 Ueli M. Maurer, "Fast Generation of Prime Numbers and Secure Public-Key Cryptographic Parameters", 1995. L<http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.26.2151>
4178 Ueli M. Maurer, "Fast Generation of Prime Numbers and Secure Public-Key Cryptographic Parameters", 1995. Generating random provable primes by building up the prime. L<http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.26.2151>
37104179
37114180 =item *
37124181
3713 Pierre-Alain Fouque and Mehdi Tibouchi, "Close to Uniform Prime Number Generation With Fewer Random Bits", pre-print, 2011. L<http://eprint.iacr.org/2011/481>
4182 Pierre-Alain Fouque and Mehdi Tibouchi, "Close to Uniform Prime Number Generation With Fewer Random Bits", pre-print, 2011. Describes random prime distributions, their algorithm for creating random primes using few random bits, and comparisons to other methods. Definitely worth reading for the discussions of uniformity. L<http://eprint.iacr.org/2011/481>
37144183
37154184 =item *
37164185
37184187
37194188 =item *
37204189
3721 L<OEIS: Primorial|http://oeis.org/wiki/Primorial>.
4190 L<OEIS: Primorial|http://oeis.org/wiki/Primorial>
37224191
37234192 =item *
37244193
37304199
37314200 =item *
37324201
3733 Manuel Benito and Juan L. Varona, "Recursive formulas related to the summation of the Möbius function", I<The Open Mathematics Journal>, v1, pp 25-34, 2007. Presents a nice algorithm to compute the Mertens functions with only n/3 Möbius values. L<http://www.unirioja.es/cu/jvarona/downloads/Benito-Varona-TOMATJ-Mertens.pdf>
4202 Manuel Benito and Juan L. Varona, "Recursive formulas related to the summation of the Möbius function", I<The Open Mathematics Journal>, v1, pp 25-34, 2007. Among many other things, shows a simple formula for computing the Mertens functions with only n/3 Möbius values (not as fast as Deléglise and Rivat, but really simple). L<http://www.unirioja.es/cu/jvarona/downloads/Benito-Varona-TOMATJ-Mertens.pdf>
37344203
37354204 =back
37364205
+53
-0
t/81-bignum.t less more
7575 + scalar(keys %allfactors)
7676 + 5 # moebius, euler_phi, jordan totient
7777 + 15 # random primes
78 + 3 # verify_prime
79 + 2*$extra # "big" prime verification
7880 + 0;
7981
8082 # Using GMP makes these tests run about 2x faster on some machines
112114 random_nbit_prime
113115 random_strong_prime
114116 random_maurer_prime
117 verify_prime
115118 /;
116119 # TODO: is_strong_lucas_pseudoprime
117120 # ExponentialIntegral
285288 cmp_ok( ($pcapf - $pclo)/$pcapf, '<=', $percent , "PC lower");
286289 cmp_ok( ($pcup - $pcapf)/$pcapf, '<=', $percent , "PC upper");
287290 }
291
292 ###############################################################################
293
294 # Provable primes
295 {
296 my @proof;
297
298 @proof = (20907001, "Pratt", [ 2,
299 3,
300 5,
301 [23,"Pratt",[2,[11,"Pratt",[2,5],2]],5],
302 [101, "Pratt", [2, 5], 2],
303 ], 14 );
304 ok( verify_prime(@proof), "simple Lucas/Pratt proof verified" );
305
306 @proof = ('3364125245431456304736426076174232972735419017865223025179282077503701', 'n-1',
307 [2,5,127, ['28432789963853652887491983185920687231739655787', 'n-1',
308 [ 2,3,163,650933, [ '44662634059309451871488121651101494489', 'n-1',
309 [ 2,3,23,4021,2321273 ],
310 [ 11, 2, 2, 2, 2 ]
311 ] ],
312 [ 2, 2, 2, 2, 2 ]
313 ],
314 '9316417838190714313' ],
315 [ 2, 2, 2, 2, 2 ]);
316 ok( verify_prime(@proof), "simple n-1 proof verified" );
317
318 @proof = ('677826928624294778921',"AGKM", ['677826928624294778921', '404277700094248015180', '599134911995823048257', '677826928656744857936', '104088901820753203', ['2293544533', '356794037129589115041']], ['104088901820753203', '0', '73704321689372825', '104088902465395836', '1112795797', ['3380482019', '53320146243107032']], ['1112795797', '0', '638297481', '1112860899', '39019', ['166385704', '356512285']]);
319 ok( verify_prime(@proof), "ECPP primality proof of 677826928624294778921 verified" );
320 }
321
322 if ($extra) {
323 my @proof;
324 #skip "Skipping bigger proofs without fast BigInt library", 2
325 # if $bigintlib !~ /^(GMP|Pari)/;
326
327 @proof = ('6864797660130609714981900799081393217269435300143305409394463459185543183397656052122559640661454554977296311391480858037121987999716643812574028291115057151', 'n-1',
328 [ 2,3,5,11,17,31,41,53,131,157,521,1613,61681,8191,42641,858001,51481, '7623851', '308761441' ],
329 [ 3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3 ] );
330 ok( verify_prime(@proof), "2**521-1 primality proof verified" );
331
332 @proof = ('531137992816767098689588206552468627329593117727031923199444138200403559860852242739162502265229285668889329486246501015346579337652707239409519978766587351943831270835393219031728127', 'n-1',
333 [ 2,3,7,607,'112102729', '341117531003194129', '7432339208719',
334 ['845100400152152934331135470251', 'n-1',
335 [2,5,11,31,41,101,251,601,1801],
336 [2,3,3,3,3,2,3,3,3]
337 ] ],
338 [ 3,5,3,2,3,3,3,3 ] );
339 ok( verify_prime(@proof), "2**607-1 primality proof verified" );
340 }