# Find a triangle-number decomposition by converting to
# a Diophantine squares problem.
use strict;
# These have no quadratic residue for -1.
# If they appear as factors in M an odd number of times
# a solution is impossible.
our @screeners = (3, 7, 11, 19, 23, 31);
# The brute-force search for a quadratic solution.
# We have reduced the number as much as we can before this
# point.
sub quad2 {
my $M = shift;
my $top = int(sqrt($M));
my $last = int($top/2) + 1;
my ($i, $j, $rem);
for ($i = $top; $i >= $last; --$i) {
$rem = $M - $i*$i;
$j = int(sqrt($rem));
if ($j*$j == $rem) {
return ($i,$j);
}
}
}
# How many times does $n go into $M?
sub countpowers {
my ($M, $n) = @_;
my $powercount = 0;
while ($M % $n == 0) {
$M = $M/$n;
++$powercount;
}
return ($M, $powercount);
}
# Reduce even numbers before solving by brute force.
# Also, screen out some impossible cases.
# Don't try a full factorization for now.
sub quadreduce {
my $M = shift;
my $powercount;
my $multfactor = 1;
# Screen out odd impossible cases, and factor out pairs of 4k-1 factors.
foreach my $num (@screeners) {
($M, $powercount) = countpowers($M, $num);
if ($powercount % 2 == 1) {
print "Eliminating impossible case with odd-power of 4k-1 factor $num\n";
return;
}
# In even cases, half the factors can be multiplied into the result.
$powercount = $powercount/2;
for (my $i = 0; $i < $powercount; $i++) {
$multfactor *= $num;
}
}
# Count powers of 2. We can handle the case of odd powers of 2.
($M, $powercount) = countpowers($M, 2);
if ($powercount % 2 == 0) {
# Two goes in an even number of times.
my ($i, $j) = quad2($M);
if ($powercount > 0) {
$multfactor *= 1 << (($powercount)/2);
}
return ($i*$multfactor, $j*$multfactor);
}
# Two goes in an odd number of times. Use the i+j, i-j trick.
my ($i, $j) = quad2($M);
if ($powercount > 1) {
$multfactor *= 1 << (($powercount - 1)/2);
}
return (($i + $j)*$multfactor, ($i - $j)*$multfactor);
}
my $M;
$M = ($#ARGV >= 0)? $ARGV[0] : 987654321;
# Convert from a triangle-number problem to a square-number problem.
my $M2 = $M*8 + 3;
# The top-level loop will try odd squares by brute force.
my $top = int(sqrt($M2));
$top = $top - 1 if ($top % 2 == 0); # start with odd.
my $last = int($top/2) + 1;
my $k;
for ($k = $top; $k >= $last; $k -= 2) {
my $M3 = $M2 - $k*$k;
my ($i, $j) = quadreduce($M3);
if ($i) {
# Convert solution back to triangle-numbers.
my $new_i = ($i - 1)/2;
my $new_j = ($j - 1)/2;
my $new_k = ($k - 1)/2;
print "Solution: triangle numbers $new_i, $new_j, $new_k\n";
my $tri_i = ($new_i+$new_i*$new_i)/2;
my $tri_j = ($new_j+$new_j*$new_j)/2;
my $tri_k = ($new_k+$new_k*$new_k)/2;
my $sum = $tri_i + $tri_j + $tri_k;
print "Verifying $tri_i + $tri_j + $tri_k = $sum = $M\n";
exit(0);
}
}