use v6-alpha;
# This function multiplies the squaring factors of $n and $m to receive
# the squaring factors of ($n*$m)
sub multiply_squaring_factors($n,$m)
{
return sort { $^a <=> $^b }, one(@$n,@$m).values;
}
my %gsf_cache = (1 => []);
# This function gets the squaring factors of $n.
# The squaring factors are those prime numbers that need to be multiplied
# by $n to reach a perfect square. They are the minimal such number.
sub get_squaring_factors($n,$start_from=2)
{
if (%gsf_cache.exists($n))
{
return %gsf_cache{$n};
}
my $p = [//] grep { $n % $_ == 0 }, ($start_from .. $n);
# This function is recursive to make better use of the Memoization
# feature.
return
(%gsf_cache{$n} =
multiply_squaring_factors(
[$p],
get_squaring_factors(int($n / $p), $p)
)
);
}
sub Graham($n)
{
my $n_sq_factors = get_squaring_factors($n);
# The graham number of a perfect square is itself.
if ($n_sq_factors.elems() == 0)
{
return ($n, [$n]);
}
# Cheating:
# Check if between n and n+largest_factor we can fit
# a square of SqFact{n*(n+largest_factor)}. If so, return
# n+largest_factor.
#
# So, for instance, if n = p than n+largest_factor = 2p
# and so SqFact{p*(2p)} = 2 and it is possible to see if
# there's a 2*i^2 between p and 2p. That way, p*2*i^2*2p is
# a square number.
{
my $largest_factor = $n_sq_factors.[-1];
my ($lower_bound, $lb_sq_factors);
$lower_bound = $n + $largest_factor;
loop
{
$lb_sq_factors = get_squaring_factors($lower_bound);
if (grep { $_ == $largest_factor }, @$lb_sq_factors)
{
last;
}
$lower_bound += $largest_factor;
}
my $n_times_lb =
multiply_squaring_factors($n_sq_factors, $lb_sq_factors);
my $rest_of_factors_product = [*] @$n_times_lb;
my $low_square_val = int(sqrt($n/$rest_of_factors_product));
my $high_square_val = int(sqrt($lower_bound/$rest_of_factors_product));
if ($low_square_val != $high_square_val)
{
return ($lower_bound, [$n, ($low_square_val+1)*($low_square_val+1)*$rest_of_factors_product,$lower_bound]);
}
}
# %primes_to_ids_map maps each prime number to its ID. IDs are consective.
my (%primes_to_ids_map);
my $next_id = 0;
# @base is an array that for each ID of a prime number holds the
# controlling vector for this number.
#
# This is in fact a matrix that is kept stair-shaped and canonized.
my (@base);
# The integers whose multiplication compose this vector
my (@base_composition);
# Register all the primes in the squaring factors of $n
for @$n_sq_factors -> $p
{
%primes_to_ids_map{$p} = ($next_id++);
}
# $n_vec is used to determine if $n can be composed out of the base's
# vectors.
my $n_vec = $n_sq_factors;
my $n_composition = [ $n ];
my $i;
loop ($i = $n+1 ; ; $i++)
{
my $i_sq_factors = get_squaring_factors($i);
# Skip perfect squares - they do not add to the solution
if ($i_sq_factors.elems() == 0)
{
next;
}
# Check if $i is a prime number
# We need n > 2 because for n == 2 it does include a prime number.
#
# Prime numbers cannot be included because 2*n is an upper bound
# to G(n) and so if there is a prime p > n than its next mulitple
# will be greater than G(n).
if (($n > 2) && ($i_sq_factors[0] == $i))
{
next;
}
# $final_vec is the new vector to add after it was
# stair-shaped by all the controlling vectors in the base.
my $final_vec = $i_sq_factors;
my $final_composition = [ $i ];
for @$i_sq_factors -> $p
{
if (!%primes_to_ids_map.exists($p))
{
# Register a new prime number
%primes_to_ids_map{$p} = ($next_id++);
}
else
{
my $id = %primes_to_ids_map{$p};
if (defined(@base[$id]))
{
$final_vec = multiply_squaring_factors($final_vec, @base[$id]);
$final_composition = multiply_squaring_factors($final_composition, @base_composition[$id]);
}
}
}
# Get the minimal ID and its corresponding prime number
# in $final_vec.
my $min_id = -1;
my $min_p = 0;
for @$final_vec -> $p
{
my $id = %primes_to_ids_map{$p};
if (($min_id < 0) || ($min_id > $id))
{
$min_id = $id;
$min_p = $p;
}
}
if ($min_id >= 0)
{
# Assign $final_vec as the controlling vector for this prime
# number
@base[$min_id] = $final_vec;
@base_composition[$min_id] = $final_composition;
# Canonize the rest of the vectors with the new vector.
my $j;
loop ($j=0 ; $j < @base.elems() ; $j++)
{
if (($j == $min_id) || (!defined(@base[$j])))
{
next;
}
if (grep { $_ == $min_p }, @base[$j][])
{
@base[$j] = multiply_squaring_factors(@base[$j], $final_vec);
@base_composition[$j] = multiply_squaring_factors(@base_composition[$j], $final_composition);
}
}
}
# A closure to print the base. It is not used but can prove useful.
my $print_base = sub {
print "Base=\n\n";
my $j;
loop ($j=0 ; $j < @base.elems() ; $j++)
{
next if (! defined(@base[$j]));
print "base[$j] (" ~ join(" * ", @base[$j][]) ~ ")\n";
}
print "\n\n";
};
# Check if we can form $n
while ($n_vec.elems())
{
# Assing $id as the minimal ID of the squaring factors of $p
my @ids_vec = (sort { $^a <=> $^b }, %primes_to_ids_map{@$n_vec});
my $id = @ids_vec[0];
# Mulitply by the controlling vector of this ID if such one exists
# or terminate if there isn't such.
last if (!defined(@base[$id]));
$n_vec = multiply_squaring_factors($n_vec, @base[$id]);
$n_composition = multiply_squaring_factors($n_composition, @base_composition[$id]);
}
if ($n_vec.elems() == 0)
{
return ($i, $n_composition);
}
}
}
my $start = int(@*ARGS.shift);
my $end = int(@*ARGS.shift || $start);
my $n;
for $start..$end -> $n
{
my ($g_val, $composition) = Graham($n);
say "G($n) = $g_val [{join(",", @$composition)}]";
}