#!/usr/bin/env perl
use strict;
use warnings;

# One of the things this test shows is the impact of the '-nobigint' option.
#    "MPU"   results are the timings without the '-nobigint' option
#    "MPUxs" results are the timings _with_  the '-nobigint' option
use Math::Prime::Util qw/factor/;

# Compare to Math::Factor::XS, which uses trial division.
use Math::Factor::XS qw/prime_factors/;

use Benchmark qw/:all/;
use List::Util qw/min max reduce/;
use Config;
my $count = shift || -2;
my $is64bit = (~0 > 4294967295);
my $maxdigits = ($is64bit) ? 20 : 10;  # Noting the range is limited for max.
my $semiprimes = 0;
my $howmany = 1000;

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;
  }
};
srand(29);

for my $d ( 3 .. $maxdigits ) {
  print "Factor $howmany $d-digit numbers\n";
  test_at_digits($d, $howmany);
}

sub test_at_digits {
  my $digits = shift;
  die "Digits has to be >= 1" unless $digits >= 1;
  die "Digits has to be <= $maxdigits" if $digits > $maxdigits;
  my $quantity = shift;

  my @rnd = genrand($digits, $quantity);
  my @smp = gensemi($digits, $quantity);

  # verify (can be _really_ slow for 18+ digits)
  foreach my $p (@rnd, @smp) {
    my @mpxs = prime_factors($p);  push @mpxs, $p if $p < 2;

    verify_factor($p, \@mpxs, [factor($p)], "Math::Prime::Util $Math::Prime::Util::VERSION");
  }


  #my $min_num = min @nums;
  #my $max_num = max @nums;
  #my $whatstr = "$digits-digit ", $semiprimes ? "semiprime" : "random";
  #print "factoring 1000 $digits-digit ",
  #      $semiprimes ? "semiprimes" : "random numbers",
  #      " ($min_num - $max_num)\n";

  my $lref = {
    "MPUxs rand"  => sub { Math::Prime::Util::_XS_factor($_) for @rnd },
    "MPUxs semi"  => sub { Math::Prime::Util::_XS_factor($_) for @smp },
    "MPU   rand"  => sub { factor($_) for @rnd },
    "MPU   semi"  => sub { factor($_) for @smp },
    "MFXS  rand"  => sub { prime_factors($_) for @rnd },
    "MFXS  semi"  => sub { prime_factors($_) for @smp },
  };
  cmpthese($count, $lref);
}

sub verify_factor {
  my $n = shift;
  my $aref_master = shift;
  my $aref_check = shift;
  my $name = shift;

  my @master = sort {$a<=>$b} @{$aref_master};
  my @check  = sort {$a<=>$b} @{$aref_check};
  die "Factor $n master fail!" unless $n == reduce { $a * $b } @master;
  die "Factor $n fail: $name" unless $#check == $#master;
  die "Factor $n fail: $name" unless $n == reduce { $a * $b } @check;
  for (0 .. $#master) {
    die "Factor $n fail: $name" unless $master[$_] == $check[$_];
  }
  1;
}

sub genrand {
  my $digits = shift;
  my $num = shift;

  my $base = ($digits == 1) ? 0 : int(10 ** ($digits-1));
  my $max = int(10 ** $digits);
  $max = ~0 if $max > ~0;
  my @nums = map { $base + $rgen->($max-$base) } (1 .. $num);
  return @nums;
}

sub gensemi {
  my $digits = shift;
  my $num = shift;

  my @min_factors_by_digit = (2,2,3,5,7,13,23,47,97);
  my $smallest_factor = $min_factors_by_digit[$digits];
  $smallest_factor = $min_factors_by_digit[-1] unless defined $smallest_factor;

  my $base = ($digits == 1) ? 0 : int(10 ** ($digits-1));
  my $max = int(10 ** $digits);
  $max = (~0-4) if $max > (~0-4);
  my @semiprimes;

  foreach my $i (1 .. $num) {
    my @factors;
    my $n;
    while (1) {
      $n = $base + $rgen->($max-$base);
      # Skip past all multiples of 2, 3, and 5.
      $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 = Math::Prime::Util::factor($n);
      next if scalar @factors != 2;
      next if $factors[0] < $smallest_factor;
      next if $factors[1] < $smallest_factor;
      last if scalar @factors == 2;
    }
    die "ummm... $n != $factors[0] * $factors[1]\n" unless $n == $factors[0] * $factors[1];
    push @semiprimes, $n;
  }
  return @semiprimes;
}