# 
# Hello,
# 
#     My solution is a combination of brute-force searching and prudent
# shortcuts.  I used an object class to represent the sequence so far;
# the object encapsulates the lower and upper bound of the sequence, the
# numbers in the list, and the prime factors of the product of the
# sequence.  Prime factors are represented as a hash of prime =>
# exponent.
# 
#     Factoring is done by a simple brute-force division by all primes
# up to 100, which means that numbers up to 10,000 can be factored.
# Factors are memoized, which greatly speeds up searches.
# 
#     The function graham_seq initializes the sequence to contain the
# input number.  It then finds the largest unpaired prime factor in the
# input, and sets initial upper bound to $n + $p (the input number plus
# that largest prime). For example:
# 
#     input:   28        ( 2^2 * 7^1)
#     largest unpaired prime:  7
#     start looking at: 28 + 7 = 35
# 
# since the series must have another multiple of 7 in it.
# 
#     Then graham_seq calls delve() to recursively delve into the
# sequence in search of a sequence of numbers that work.  At each
# iteration, delve() finds the largest unpaired prime number in the
# sequence, and loops over multiples of that prime number in the range.
# If at any point, delve() reaches a solution, it returns true and that
# propagates back up.  Otherwise, it keeps searching or returns 0 if it
# has exhausted its possibilities.
# 
#     If delve() returns false, graham_seq increments the upper bound
# and tries again.  This seems to me to be the weak point of the
# program.  Intuitively, it seems to me that one could do better than to
# step by 1, but I couldn't think of a way that was mathematically
# sound.
# 
#     Because the program skips ranges of numbers that can't possibly
# work, it's sorta smart; because it's otherwise a brute-force solution,
# it's kinda dumb.  Some people had reported long run times for prime
# input numbers, such as 19 or 337; my program blasts over those. On 
# the other hand, it takes a very long time for some other numbers.
# 918 took over 7 seconds; 4433 took over an hour before I stopped it.
# 
# 
#     The program is organized as follows: First there's a couple of
# utility functions, prime_factors and is_square.  These are memoized
# because they are called repeatedly on the same input, and their return
# values don't depend on anything else.  Then comes the Sequence object
# class, which I think is pretty straightforward.  When a number is 
# added to the class, it is factored, and the exponents of its prime
# factors are added to those of the product of the sequence.
#     Then the main package follows, with the delve and graham_seq
# functions as described above.
# 
# Here's the code:
# 
#!/usr/bin/perl

use strict;

package Utility;

use Memoize;
memoize $_ for qw(prime_factors is_square);

my @primes = (2, 3, 5, 7, 9, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 73, 79, 83, 89, 97);
sub prime_factors
	{
	my $n = shift;
	die "Can't factor numbers (like $n) larger than 10,000.\n" if $n > 10_000;
	my %exponent;
	Prime: foreach my $p (@primes)
		{
		while ($n % $p == 0)
			{
			$n /= $p;
			$exponent{$p}++;
			last Prime unless $n > 1;
			}
		}
	$exponent{$n} = 1  if $n > 1;
	return \%exponent;
	}

sub is_square
	{
	my $n = shift;
	my $nf = prime_factors $n;

	# Check whether they're all even
	my ($p,$e);
	while ( ($p,$e) = each %$nf)
		{
		return 0 if $e % 2;
		}
	return 1;
	}


package Sequence;

sub new
	{
	my $class = shift;
	my $start = shift;
	my $factors = Utility::prime_factors($start);
	$factors = { %$factors };    # Make deep copy because of caching
	bless { lower   => $start,
			upper   => undef,
			numbers => {$start => 1},
			factors => $factors }, $class;
	}
sub lower
	{
	my $self = shift;
	return $self->{lower};
	}
sub upper
	{
	my $self = shift;
	return $self->{upper};
	}
sub set_upper
	{
	my $self = shift;
	my $up   = shift;    # new upper bound
	$self->remove($self->{upper});    # Remove old upper bound, if any.
	$self->{upper} = $up;
	$self->add($up);
	}
sub contains
	{
	my $self = shift;
	my $num  = shift;
	return exists $self->{numbers}{$num};
	}
sub list
	{
	my $self = shift;
	return sort {$a <=> $b} keys %{$self->{numbers}};
	}
sub add
	{
	my $self = shift;
	my $new  = shift;
	return if exists $self->{numbers}{$new};    # already there; don't re-add
	return if $new > $self->{upper};   # bounding condition

	my $newf = Utility::prime_factors($new);

	$self->{numbers}{$new} = 1;    # Add it to list
	for my $prime (keys %$newf)
		{
		$self->{factors}{$prime} += $newf->{$prime};
		}
	}
sub remove
	{
	my $self = shift;
	my $gone = shift;
	return unless exists $self->{numbers}{$gone};  # Can't remove if not there.

	delete $self->{numbers}{$gone};    # Remove it
	my $oldf = Utility::prime_factors($gone);
	for my $prime (keys %$oldf)
		{
		$self->{factors}{$prime} -= $oldf->{$prime};
		}
	}
sub is_square
	{
	my $self = shift;
	for my $val (values %{$self->{factors}})
		{
		return 0 if $val % 2;
		}
	return 1;
	}
sub max_odd_prime
	{
	my $self = shift;
	my $max;
	my ($prime,$exp);
	while ( ($prime,$exp) = each %{$self->{factors}} )
		{
		next unless $exp % 2;
		$max = $prime if $max < $prime;
		}
	return $max;
	}
		

package main;
use Data::Dumper;

# Returns true if a square is found.
sub delve
	{
	my $seq = shift;

	# Find the largest odd-powerd prime factor in the sequence.
	# If there is none, the sequence represents a square number and we're done.
	my $p = $seq->max_odd_prime($seq) or return 1;
	my $q;
	my ($min, $max) = ($seq->lower, $seq->upper);

	# If $p is greater than min,max range, there cannot possibly be another.
	return 0 if $p >= $max - $min;

	# Loop over multiples of $p within the input range.
	for ($q = $min+$p - $min%$p; $q < $max; $q += $p)
		{
		next if $seq->contains($q);        # Skip if it's already there.
		next if Utility::is_square($q);    # Don't bother with squares.

		# Add it to the sequence and see if we can make the seq work with it.
		$seq->add($q);
		delve($seq) and return 1;
		$seq->remove($q);     # Nope, didn't fine a working set.
		}
	return 0;
	}

# $seq = graham_seq ($number);
sub graham_seq
	{
	my $n = shift;
	my $seq = new Sequence($n);

	# Find largest odd-powered prime factor.  If none, input is square.
	my $p = $seq->max_odd_prime or return $seq;
	my $m;

	# Start at $n+$p; no sequence is possible with a smaller upper bound.
	for ($m = $n+$p; ; $m++)
		{
		next if Utility::is_square($m);    # Don't bother with squares.
		$seq->set_upper($m);
		delve($seq) and return $seq;
		}
	}

# $num = graham_number($number);
sub graham_number
	{
	my $seq = graham_seq(@_);
	return $seq->upper;
	}

# ----------------------------------------------------------------
# Main program
BEGIN { *Graham = \&graham_number }
#my $n = shift;
#my $g = graham_number $n;
#print "$g\n";

# To print the sequence instead of the max number:
# my $s = graham_seq $n;
# print join(', ', $s->list), "\n";

__END__

-- 
Eric J. Roode                                              eric@myxa.com
Senior Software Engineer, Myxa Corporation
tr j, j ,j for @japh = (qw b lre h, uJ p, ekca tona, ts reh b, $/.r);
print scalar reverse sort @japh;
