#!/usr/bin/perl
#
# Use Bit::vector->primes?

use Bit::Vector::Overload;

my @primes = ();
{ my @factors = (undef, Bit::Vector->new(0));
  my %is_prime;
  sub factor {
    my $n = shift;
    return $factors[$n] if defined $factors[$n];
    for my $pi (0 .. $#primes) {
      my $p = $primes[$pi];
      next if $n % $p;          # does not divide;
      last if $p*$p > $n;      # $n is prime
      return $factors[$n] = mjd_xor(factor($n/$p), $is_prime{$p});
    }
    # New prime
    push @primes, $n unless $is_prime{$n};
    my $v = $factors[$n] = $is_prime{$n} = Bit::Vector->new(scalar(@primes));
    $v->Bit_On($#primes);
    $v;
  }
}

# Extend shorter vector to length of longer
# then set $r to exclusive or of result
sub mjd_xor {
  my ($a, $b) = @_;
  my ($as, $bs) = ($a->Size, $b->Size);
  if    ($as > $bs) { $b->Resize($as) }
  elsif ($as < $bs) { $a->Resize($bs); $as = $bs }
  my $r = Bit::Vector->new($as);
  $r->ExclusiveOr($a, $b);
  $r;
}

sub Bit::Vector::show {
  print "Vector: ";
  for ($_[0]->Index_List_Read) {
    print $primes[$_], " ";
  }
  print "\n";
}


my $INFINITY = 999999999;
sub Graham {
  my $n = shift;
  my $MAX = $INFINITY;
  search($n, factor($n), \$MAX);
}

# Return smallest max
sub search {
  my ($lo, $t, $hiR) = @_;
  return if $lo > $$hiR;
  return $lo unless $t;
  # Largest prime still present in factorization
  my $maxprime_bit = $t->Max();
  my $maxprime = $primes[$maxprime_bit];
  # Smallest multiple of $maxprime that is > $lo
  my $lopm = (int($lo/$maxprime)+1);

  my $best = $$hiR + 1;  # Best so far is not good enough
  for ($lopm .. int($$hiR/$maxprime)) {
    my $try = $_ * $maxprime;
    last if $try > $$hiR;
    # Try adding $_ * $maxprime to the set
    my $f = factor($try);
    next unless $f->contains($maxprime_bit); # Not useful
    my $n = mjd_xor($t, $f);
#    $n->Bit_Off($maxprime_bit);
    my $res = search($try, $n, $hiR);
    if ($res < $best) {
      $$hiR = $best = $res;
    }
  }
  return $best > $$hiR ? () : $best;
}

my $z = Graham(2);
