1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
|
#!/usr/bin/env perl
use strict;
use warnings;
$| = 1; # fast pipes
use Math::Prime::Util qw/urandomm/;
use Math::Primality;
use Config;
my $nlinear = 10000;
my $nrandom = shift || 20000;
my $randmax = ~0;
my $rand_ndigit_gen = sub {
my $digits = shift;
die "Digits must be > 0" unless $digits > 0;
my $howmany = shift || 1;
my ($base, $max);
if ( 10**$digits < ~0) {
$base = ($digits == 1) ? 0 : int(10 ** ($digits-1));
$max = int(10 ** $digits);
$max = ~0 if $max > ~0;
} else {
$base = Math::BigInt->new(10)->bpow($digits-1);
$max = Math::BigInt->new(10)->bpow($digits) - 1;
}
my @nums = map { $base + urandomm($max-$base) } (1 .. $howmany);
return (wantarray) ? @nums : $nums[0];
};
if (1) {
print "OK for first 1";
my $dig = 1;
my $i = 9;
foreach my $n (2 .. $nlinear) {
die "MR(2) failure for $n" unless Math::Prime::Util::is_strong_pseudoprime($n,2) == Math::Primality::is_strong_pseudoprime($n,2);
die "SLPSP failure for $n" unless Math::Prime::Util::is_strong_lucas_pseudoprime($n) == Math::Primality::is_strong_lucas_pseudoprime($n);
die "Prime failure for $n" unless Math::Prime::Util::is_prime($n) == Math::Primality::is_prime($n);
if (--$i == 0) {
print "0";
$dig++;
$i = (10 ** $dig) - (10 ** ($dig-1));
}
}
print " numbers\n";
print "Testing random numbers from $nlinear to ", $randmax, "\n";
foreach my $r (1 .. $nrandom) {
my $n = $nlinear + 1 + int(rand($randmax - $nlinear));
my $rand_base = 2 + urandomm($n-4);
die "MR(2) failure for $n" unless Math::Prime::Util::is_strong_pseudoprime($n,2) == Math::Primality::is_strong_pseudoprime($n,2);
die "MR($rand_base) failure for $n" unless Math::Prime::Util::is_strong_pseudoprime($n,$rand_base) == Math::Primality::is_strong_pseudoprime($n,$rand_base);
die "SLPSP failure for $n" unless Math::Prime::Util::is_strong_lucas_pseudoprime($n) == Math::Primality::is_strong_lucas_pseudoprime($n);
my $ip1 = Math::Primality::is_prime($n);
my $ip2 = Math::Prime::Util::is_prime($n);
die "Prime failure for $n ($ip1,$ip2)" unless !!$ip1 == !!$ip2;
print "." if ($r % 256) == 0;
}
print "\n";
}
if (1) {
use bigint try => 'GMP,Pari';
my $big_base = 2**64 + 1;
my $range = 2**1024 - 1;
my $end_base = $big_base + $range;
print "Testing random numbers from $big_base to $end_base\n";
foreach my $r (1 .. int($nrandom/100)) {
my $n = $big_base + urandomm($range);
my $rand_base = 2 + urandomm($n-4);
die "MR(2) failure for $n" unless Math::Prime::Util::is_strong_pseudoprime($n,2) == Math::Primality::is_strong_pseudoprime("$n","2");
die "MR($rand_base) failure for $n" unless Math::Prime::Util::is_strong_pseudoprime($n,$rand_base) == Math::Primality::is_strong_pseudoprime($n,$rand_base);
die "SLPSP failure for $n" unless Math::Prime::Util::is_strong_lucas_pseudoprime($n) == Math::Primality::is_strong_lucas_pseudoprime("$n");
die "Prime failure for $n" unless (Math::Prime::Util::is_prime($n)?1:0) == Math::Primality::is_prime("$n");
#print "SUCCESS with $n\n";
print "." if ($r % 16) == 0;
}
print "\n";
}
print "\nBenchmarks\n";
my $num_rns = 100;
my $len_rns = 100;
my $count = -1;
use bigint try => 'GMP,Pari';
my @rns; # make the primality tests at least lift a finger.
while (@rns < $num_rns) {
my $n = $rand_ndigit_gen->($len_rns);
next unless $n%2 && $n%3 && $n%5 && $n%7 && $n%11 && $n%13;
push @rns, $n;
}
use Benchmark qw/:all/;
require Math::Prime::Util::PP;
print "Starting benchmarks, $num_rns $len_rns-digit random numbers...\n";
if (1) {
print "\nMiller-Rabin, one base:\n";
cmpthese($count, {
"MPU:PP" => sub { Math::Prime::Util::PP::is_strong_pseudoprime($_,2) for @rns; },
"MPU:GMP" => sub { Math::Prime::Util::GMP::is_strong_pseudoprime($_,2) for @rns; },
"MPU" => sub { Math::Prime::Util::is_strong_pseudoprime($_,2) for @rns; },
"MP" => sub { Math::Primality::is_strong_pseudoprime("$_","2") for @rns; },
});
}
if (1) {
print "\nStrong Lucas test:\n";
cmpthese($count, {
"MPU:PP" => sub { Math::Prime::Util::PP::is_strong_lucas_pseudoprime($_) for @rns;},
"MPU:GMP" => sub { Math::Prime::Util::GMP::is_strong_lucas_pseudoprime($_) for @rns;},
"MPU" => sub { Math::Prime::Util::is_strong_lucas_pseudoprime($_) for @rns;},
"MP" => sub { Math::Primality::is_strong_lucas_pseudoprime("$_") for @rns;},
});
}
if (1) {
print "\nBPSW test:\n";
cmpthese($count, {
"MPU:PP" => sub { my $sum = 0;
do { $sum += ( Math::Prime::Util::PP::is_strong_pseudoprime($_, 2) &&
Math::Prime::Util::PP::is_strong_lucas_pseudoprime($_) )
? 1 : 0 } for @rns; },
"MPU:GMP" => sub { Math::Prime::Util::GMP::is_prob_prime($_) for @rns; },
"MPU" => sub { Math::Prime::Util::is_prob_prime($_) for @rns;},
"MP" => sub { Math::Primality::is_prime("$_") for @rns;},
});
}
|
