#!/usr/bin/env perl
use strict;
use warnings;
$| = 1; # fast pipes
srand(377);
use Math::Prime::Util qw/factor/;
use Math::Factor::XS qw/prime_factors/;
use Math::Pari qw/factorint/;
use Benchmark qw/:all/;
use Data::Dumper;
use Config;
my $digits = shift || 15;
my $count = shift || -3;
my $rgen = sub {
my $range = shift;
return 0 if $range <= 0;
my $rbits = 0; { my $t = $range; while ($t) { $rbits++; $t >>= 1; } }
while (1) {
my $rbitsleft = $rbits;
my $U = 0;
while ($rbitsleft > 0) {
my $usebits = ($rbitsleft > $Config{randbits}) ? $Config{randbits} : $rbitsleft;
$U = ($U << $usebits) + int(rand(1 << $usebits));
$rbitsleft -= $usebits;
}
return $U if $U <= $range;
}
};
my @min_factors_by_digit = (2,2,3,3,5,11,17,47,97);
my $smallest_factor_allowed = $min_factors_by_digit[$digits];
$smallest_factor_allowed = $min_factors_by_digit[-1] unless defined $smallest_factor_allowed;
my $numprimes = 200;
die "Digits has to be >= 2" unless $digits >= 2;
die "Digits has to be <= 10" if (~0 == 4294967295) && ($digits > 10);
die "Digits has to be <= 19" if $digits > 19;
my $skip_mfxs = ($digits > 17);
# Construct some semiprimes of the appropriate number of digits
# There are much cleverer ways of doing this, using randomly selected
# nth_primes, and so on, but this works well until we get lots of digits.
print "Generating $numprimes random $digits-digit semiprimes (min factor $smallest_factor_allowed) ";
my @semiprimes;
foreach my $i ( 1 .. $numprimes ) {
my $base = int(10 ** ($digits-1));
my $add = int(10 ** ($digits)) - $base;
my @factors;
my $n;
while (1) {
$n = $base + $rgen->($add);
next if $n > (~0 - 4);
$n += (1,0,5,4,3,2,1,0,3,2,1,0,1,0,3,2,1,0,1,0,3,2,1,0,5,4,3,2,1,0)[$n%30];
@factors = factor($n);
next if scalar @factors != 2;
next if $factors[0] < $smallest_factor_allowed;
next if $factors[1] < $smallest_factor_allowed;
last if scalar @factors == 2;
}
die "ummm... $n != $factors[0] * $factors[1]\n" unless $n == $factors[0] * $factors[1];
#print "$n == $factors[0] * $factors[1]\n";
push @semiprimes, $n;
print "." if ($i % ($numprimes/10)) == 0;
}
print "done.\n";
print "Verifying Math::Prime::Util $Math::Prime::Util::VERSION ...";
foreach my $sp (@semiprimes) {
my @factors = factor($sp);
die "wrong for $sp\n" unless ($#factors == 1) && ($factors[0] * $factors[1]) == $sp;
}
print "OK\n";
if (!$skip_mfxs) {
print "Verifying Math::Factor::XS $Math::Factor::XS::VERSION ...";
foreach my $sp (@semiprimes) {
my @factors = prime_factors($sp);
die "wrong for $sp\n" unless ($#factors == 1) && ($factors[0] * $factors[1]) == $sp;
}
print "OK\n";
} else {
print "Math::Factor::XS is too slow for $digits digits. Skipping.\n";
}
print "Verifying Math::Pari $Math::Pari::VERSION ...";
foreach my $sp (@semiprimes) {
my @factors;
my ($pn,$pc) = @{factorint($sp)};
push @factors, (int($pn->[$_])) x $pc->[$_] for (0 .. $#{$pn});
die "wrong for $sp\n" unless ($#factors == 1) && ($factors[0] * $factors[1]) == $sp;
}
print "OK\n";
my %compare = (
'MPU' => sub { do { my @f = factor($_) } for @semiprimes; },
'MFXS' => sub { do { my @f = prime_factors($_) } for @semiprimes; },
'Pari' => sub { do { my ($pn,$pc) = @{factorint($_)}; my @f = map { int($pn->[$_]) x $pc->[$_] } 0 .. $#$pn; } for @semiprimes; },
);
delete $compare{'MFXS'} if $skip_mfxs;
cmpthese($count, \%compare);