more useful options PerlMonks

### Re: Triangle Numbers Revisited

by tall_man (Parson)
 on Oct 19, 2004 at 18:33 UTC ( #400619=note: print w/ replies, xml ) Need Help??

in reply to Triangle Numbers Revisited

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.
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.
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.
\$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.
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) {
} else {
}
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;
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__

+00000000000000000000000
Solution: triangle numbers 104347715317115126631, 22509054715961300073
+, 14054567378613971410555040675654052625482
Verifying 5444222845950851406213673791756140268396 + 25332877210306982
+1564907070468155552701 + 98765432099999999999999999999999999999994302
+448381946078772221419137775704178903 = 987654321000000000000000000000
+00000000000000000000000000000000000000000000000000 = 9876543210000000
+0000000000000000000000000000000000000000000000000000000000000000
Found in 38 guesses

Comment on Re: Triangle Numbers Revisited

Create A New User
Node Status?
node history
Node Type: note [id://400619]
help
Chatterbox?
and the web crawler heard nothing...

How do I use this? | Other CB clients
Other Users?
Others cooling their heels in the Monastery: (10)
As of 2015-11-25 18:29 GMT
Sections?
Information?
Find Nodes?
Leftovers?
Voting Booth?