Here is a new solution using
Math::Pari for big-number calculations and for factoring. If M is factored as (a^m * b^n * ...), a sum-of-two-squares for it can be found in O(log a + log b + ...) time (if possible: impossible cases can be identified and rejected instantly). It handles the 80-digit example I included in under a second.
# Find a triangle-number decomposition by converting to
# a Diophantine squares problem.
use strict;
use Math::Pari qw(:int factorint sqrtint PARI);
# Count the number of "guesses" to the a^2 + b^2 = M problem.
our $guesses = 0;
# The brute-force search for a quadratic solution.
# We have reduced the number as much as we can before this
# point. Use this method for small numbers only.
sub quad2 {
my $M = shift;
my $top = sqrtint($M);
my $last = int($top/2) + 1;
my ($i, $j, $rem);
for ($i = $top; $i >= $last; --$i) {
++$guesses;
$rem = $M - $i*$i;
$j = sqrtint($rem);
if ($j*$j == $rem) {
return ($i,$j);
}
}
}
# Use this algorithm for large numbers.
# This is an O(log M) algorithm like Euler's method for
# greatest common factor.
sub quad {
my $M = shift;
my ($a, $b, $k, $ap, $bp, $an, $bn, $atemp);
# Computation of sqrt(-1) mod M required for starting point.
# This exists for primes of the form 4k + 1.
# Start with: $a*$a + 1*1 = $k*$M
$a = PARI "component(sqrt(Mod(-1,$M)),2)";
$b = 1;
while (1) {
# Find a solution to $a*$a + $b*$b = $M
# At each step we will find $a and $b
# such that ($a*$a + $b*$b) = $k*$M,
# and then reduce mod $k so that there
# is a new $k less than the last one.
# Stop when $k is reduced to 1.
++$guesses;
$k = ($a*$a + $b*$b)/$M;
return ($a, $b) if ($k == 1);
$ap = $a % $k;
$an = ($k - $a) % $k;
$bp = $b % $k;
$bn = ($k - $b) % $k;
$ap = -$an if ($an < $ap);
$bp = -$bn if ($bn < $bp);
$atemp = ($a*$ap + $b*$bp)/$k;
$b = ($a*$bp - $b*$ap)/$k;
$a = $atemp;
$a = -$a if ($a < 0);
$b = -$b if ($b < 0);
}
}
# Reduce to factors before solving by brute force.
# Also, screen out some impossible cases.
# Now using full factorization with Math::Pari.
sub quadreduce {
my $M = shift;
my $powercount;
my $multfactor = 1;
my $doubler = 0;
my $a = 0;
my $b = 0;
my $factors = factorint($M);
my $rows = Math::Pari::matsize($factors)->[0];
# Screen out odd impossible cases, and factor out pairs of 4k-1 fac
+tors.
my (@goodfactors, @goodpowers);
for (my $i = 0; $i < $rows; ++$i) {
my $factor = $factors->[0][$i];
my $power = $factors->[1][$i];
# Special case for factor == 2. Reduce M to an odd number.
if ($factor == 2) {
$M = $M/(1 << $power);
if ($power % 2 == 1) {
$doubler = 1;
--$power;
}
$multfactor *= (1 << ($power/2)) if ($power > 0);
next;
}
# Otherwise, screen the factors that would make it impossible.
if ($factor % 4 == 3) {
if ($power % 2 == 1) {
print "Eliminating impossible case with odd-power $power o
+f 4k-1 factor $factor,\n";
return;
} else {
# In even cases, half the factors can be multiplied into t
+he result.
$M = $M/($factor ** $power);
$power = $power/2;
$multfactor *= ($factor ** $power);
}
} else {
push @goodfactors, $factor;
push @goodpowers, $power;
}
}
# The remaining factors are of the form 4k+1 and they
# must have a solution. We can build up an answer factor by factor
+.
my $pos = 0;
foreach (@goodfactors) {
my $fact = $_;
my $pow = $goodpowers[$pos];
if ($pow % 2 == 0) {
$multfactor *= ($fact ** ($pow/2));
} else {
$multfactor *= ($fact ** (($pow-1)/2)) if ($pow > 1);
my ($a1, $b1);
if ($fact > 1000) {
($a1, $b1) = quad($fact);
} else {
($a1, $b1) = quad2($fact);
}
if ($a == 0) {
# Capture the result from the first factor.
$a = $a1;
$b = $b1;
} else {
# Build the answer from the previous one. Uses the formula
+:
# ($a*$a + $b*$b)*($a1*$a1 + $b1*$b1) ==
# ($a*$a1 + $b*$b1)**2 + ($a*$b1 - $b*$a1)**2
my $atemp = ($a*$a1 + $b*$b1);
$b = ($a*$b1 - $b*$a1);
$a = $atemp;
}
}
$pos++;
}
# 2*($a*$a + $b*$b) == ($a + $b)**2 + ($a - $b)**2
# This is just a special case of the formula above, using 1*1 + 1*1
if ($doubler) {
return (($a+$b)*$multfactor, ($b - $a)*$multfactor) if ($a - $b)
+ < 0;
return (($a+$b)*$multfactor, ($a - $b)*$multfactor);
}
return ($a*$multfactor, $b*$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 = sqrtint($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) {
# UPDATE - bug fix. Sometimes get negatives.
$i = -$i if ($i < 0);
$j = -$j if ($j < 0);
# 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";
print "Found in $guesses guesses\n";
exit(0);
}
}
__END__
perl quad.pl 987654321000000000000000000000000000000000000000000000000
+00000000000000000000000
Solution: triangle numbers 104347715317115126631, 22509054715961300073
+, 14054567378613971410555040675654052625482
Verifying 5444222845950851406213673791756140268396 + 25332877210306982
+1564907070468155552701 + 98765432099999999999999999999999999999994302
+448381946078772221419137775704178903 = 987654321000000000000000000000
+00000000000000000000000000000000000000000000000000 = 9876543210000000
+0000000000000000000000000000000000000000000000000000000000000000
Found in 38 guesses
