1#!/usr/bin/env perl
2use strict;
3use warnings;
4use Math::Prime::Util qw/factor/;
5# Compare to Math::Factor::XS, which uses trial division.
6use Math::Factor::XS qw/prime_factors/;
7
8use Benchmark qw/:all/;
9use List::Util qw/min max reduce/;
10my $count = shift || -2;
11my $is64bit = (~0 > 4294967295);
12my $maxdigits = ($is64bit) ? 20 : 10;  # Noting the range is limited for max.
13my $semiprimes = 0;
14my $howmany = 1000;
15
16for my $d ( 3 .. $maxdigits ) {
17  print "Factor $howmany $d-digit numbers\n";
18  test_at_digits($d, $howmany);
19}
20
21sub test_at_digits {
22  my $digits = shift;
23  die "Digits has to be >= 1" unless $digits >= 1;
24  die "Digits has to be <= $maxdigits" if $digits > $maxdigits;
25  my $quantity = shift;
26
27  my @rnd = ndigit_rand($digits, $quantity);
28  my @smp = genrough($digits, $quantity);
29
30  # verify (can be _really_ slow for 18+ digits)
31  foreach my $p (@rnd, @smp) {
32    next if $p < 2;
33    verify_factor($p, [prime_factors($p)], [factor($p)], "Math::Prime::Util $Math::Prime::Util::VERSION");
34  }
35
36  #my $min_num = min @nums;
37  #my $max_num = max @nums;
38  #my $whatstr = "$digits-digit ", $semiprimes ? "semiprime" : "random";
39  #print "factoring 1000 $digits-digit ",
40  #      $semiprimes ? "semiprimes" : "random numbers",
41  #      " ($min_num - $max_num)\n";
42
43  my $lref = {
44    "MPU   random"    => sub { my@a=factor($_) for @rnd },
45    "MPU   nonsmooth" => sub { my@a=factor($_) for @smp },
46    "MFXS  random"    => sub { my@a=prime_factors($_) for @rnd },
47    "MFXS  nonsmooth" => sub { my@a=prime_factors($_) for @smp },
48  };
49  cmpthese($count, $lref);
50}
51
52sub verify_factor {
53  my ($n, $aref1, $aref2, $name) = @_;
54
55  return 1 if "@$aref1" eq "@$aref2";
56
57  my @master = @$aref1;
58  my @check  = @$aref2;
59  die "Factor $n master fail!" unless $n == reduce { $a * $b } @master;
60  die "Factor $n fail: $name" unless $#check == $#master;
61  die "Factor $n fail: $name" unless $n == reduce { $a * $b } @check;
62  for (0 .. $#master) {
63    die "Factor $n fail: $name" unless $master[$_] == $check[$_];
64  }
65  1;
66}
67
68sub genrough {
69  my ($digits, $num) = @_;
70
71  my @min_factors_by_digit = (2,2,3,5,7,13,23,47,97);
72  my $smallest_factor = $min_factors_by_digit[$digits];
73  $smallest_factor = $min_factors_by_digit[-1] unless defined $smallest_factor;
74
75  my @semiprimes;
76  foreach my $i (1 .. $num) {
77    my $n;
78    my @facn;
79    do {
80      $n = ndigit_rand($digits, 1);
81      @facn = Math::Prime::Util::trial_factor($n,$smallest_factor);
82    } while scalar(@facn) > 1;
83    push @semiprimes, $n;
84  }
85  return @semiprimes;
86}
87
88use Bytes::Random::Secure qw/random_string_from/;
89sub ndigit_rand {
90  my($digits, $howmany) = @_;
91  die "digits must be > 0" if $digits < 1;
92  $howmany = 1 unless defined $howmany;
93  # TODO: need to skip things larger than ~0 for this module
94  my @nums = map { random_string_from("123456789",1) . random_string_from("0123456789",$digits-1) } 1 .. $howmany;
95  if (10**$digits > ~0) {  @nums = map { Math::BigInt->new($_) } @nums;  }
96  else                  {  @nums = map { int($_) } @nums;                }
97  return wantarray ? @nums : $nums[0];
98}
99