Limbic~Region has asked for the
wisdom of the Perl Monks concerning the following question:
All,
YuckFoo once asked if there was a more efficient way to verify that any nonnegative whole integer could be represented as the sum of 3 triangular numbers. particle came up with a pretty cool solution using bit strings. While fast, it blew up when I asked it for 987654321.
My usual approach to these kind of problems is to translate how I would solve the problem with a pencil and paper in code, make any obvious performance changes, and see if it is fast enough. I usually don't start thinking about more efficient algorithms unless it is a "challenge" or my initial approach was just plain abysmal. Here was my go at it:
#!/usr/bin/perl
use strict;
use warnings;
my $target = defined $ARGV[0] ? $ARGV[0] : 987654321;
print join ', ' , get_three( $target );
sub get_three {
my $num = shift;
my $prev = p_tri( $num );
return ($num, 0, 0) if ! $prev;
while ( $prev ) {
my $p_tri = ($prev * $prev + $prev) / 2;
my $i = 0;
my $j = $num  $p_tri;
while ( $j >= $i ) {
return $p_tri, $i, $j if ! p_tri( $i ) && ! p_tri( $j );
++$i;
$j;
}
$prev;
}
die "Something went horribly wrong : $num\n";
}
sub p_tri {
my $num = shift;
my $x = ( sqrt( 8 * $num + 1 ) + 1 )/ 2;
my $t = int $x;
return $t == $x ? 0 : $t;
}
Now without explaining how this works^{*}, I will ask the same question as YuckFoo  Any suggestions to optimize the search method?
* If anyone wants to know, feel free to ask. If anyone wants to offer their own explanation, feel free to post.
Re: Triangle Numbers Revisited
by FoxtrotUniform (Prior) on Oct 13, 2004 at 23:54 UTC

Just off the top of my head:
In the inner while loop, it looks like you're testing each $i, $j to see if it's triangular (although I must confess, I don't really understand p_tri). Since we can enumerate triangle numbers (if Jobby is to be believed), why not iterate through those?
Edit: Are you sure you're getting correct solutions? It looks like $prev starts out as a triangular number, but then you just decrement it. (Aside: if this algorithm gives you consistently correct solutions, then the greatest triangular number lower than n is always present in n's triangular decomposition. I think.) Man, I'm an idiot. Since when is $prev a trinum?
Edit 2: Oh, I get it, p_tri returns the rank of the previous triangular number, not the number itself... although it returns zero when given a triangular number, which is weird but useful.
Edit 3: Here's the start of an implementation. It's not as fast as the code posted above (by about a factor of two, if Unix time is to be trusted), probably because it calculates a lot of trinums and makes a lot of function calls. I'm posting it more or less as a proof of concept.
Edit 4: Inlining the calls to &trinum results in slightly faster code than Limbic~Region's. Code updated, benchmarks added.

Yours in pedantry,
F
o
x
t
r
o
t
U
n
i
f
o
r
m
 [reply] [d/l] [select] 

 [reply] [d/l] 

FoxtrotUniform,
The trouble with selecting a single number to benchmark from is that it may favor one method over the other. With this in mind, I created another rudimentary benchmark that shows the two methods are fairly equal:
I haven't spent any time thinking about optimizations, but if I come up with any tomorrow I will post it along with a "real" benchmark.
 [reply] [d/l] [select] 

The trouble with selecting a single number to benchmark from is that it may favor one method over the other.
Good catch. Another problem I ran into was that the load on my test system was varying (lab machine, someone was logged in remotely doing a bunch of matlab foolery), so I ended up getting quite different "real" time results with the same input... which is probably what a casual reader would check.
It would be interesting to know what kinds of inputs favour whose method, and (ideally) why.

Yours in pedantry,
F
o
x
t
r
o
t
U
n
i
f
o
r
m
 [reply] 

FoxtrotUniform,
I couldn't sleep as I had several optimizations pop into my head I just had to try. I played around with caching known triangular numbers as well as results of known nontriangular numbers. I only got a marginal speed increase which I guessed would be the case earlier in the CB. Moving away from caching, I tried a combination of my method and your method with dramatic results:
From 1 minute to just over 3 seconds.
 [reply] [d/l] 
Re: Triangle Numbers Revisited
by tachyon (Chancellor) on Oct 14, 2004 at 06:53 UTC

Here is one way to make it 30x 150x faster using your 5000 random number bench. Kinda cheating of course and still room for more optimisation.
[root@devel3 root]# time ./trinum.pl >out
real 0m41.310s
user 0m41.310s
sys 0m0.000s
[root@devel3 root]# time ./trinum_c2.pl >out
real 0m0.276s
user 0m0.280s
sys 0m0.000s
[root@devel3 root]# more trinum_c2.pl
#!/usr/bin/perl
use Inline 'C';
my $file = $ARGV[0]  'targets.txt';
open (INPUT, '<', $file) or die "Unable to open $file for reading : $!
+";
while ( <INPUT> ) {
chomp;
get_three( $_ );
}
__END__
__C__
#include <math.h>
#define TRI(i) (i*(i+1)/2)
int get_three ( int num ) {
register max,target,i,j;
max = (int)sqrt((double)(2*num));
while ( max != 0 ) {
target = num  TRI(max);
for( i=0;TRI(i)<=target;i++ ) {
for( j=0;j<=i;j++ ) {
if ( target == (int)(TRI(i)+TRI(j)) ) {
printf( "%d, %d, %d\n", TRI(max), TRI(i), TRI(j) )
+;
return 1;
}
}
}
max;
}
printf( "Something went horribly wrong : %d\n", num );
return 0;
}
[root@devel3 root]#
 [reply] [d/l] [select] 

 [reply] 
Re: Triangle Numbers Revisited
by hv (Parson) on Oct 14, 2004 at 12:47 UTC

I'm not sure quite how you might take advantage of it for optimisation, but the possible decompositions are restricted by the mod 3 equivalences. That is:
T_n == 1 (mod 3) when n == 1 (mod 3)
T_n == 0 (mod 3) otherwise
so if we split the T_n sequence into A_n (== 0) and B_n (== 1) the possible decompositions are restricted such that: if n == 0 (mod 3), require A A A or B B B
if n == 1 (mod 3), require A A B
if n == 2 (mod 3), require A B B
You can get similar restrictions by considering other prime moduli, but 3 is likely to be the most beneficial because it has a shorter than possible cycle in T_n.
Also, a word of warning on your p_tri() routine  the final result relies on comparing a floating point number for equality, normally considered a bad idea. It would be preferable to do the test instead by calculating from $t back up to the last known integer value, something like: sub p_tri {
my $num = shift;
my $x = 8 * $num + 1;
my $t = int((sqrt($x) + 1)/2);
return +((2 * $t  1) * (2 * $t  1) == $x) ? 0 : $t;
}
Hugo  [reply] [d/l] [select] 
Re: Triangle Numbers Revisited
by TedPride (Priest) on Oct 14, 2004 at 10:25 UTC

use strict;
use warnings;
my ($i, $j, $inum, $jnum, $left);
$ = 1; print "Enter the number to check : ";
my $inp = <STDIN>;
BING:
for ($i = tget($inp); $i; $i) { # BUG ON THIS LINE  SEE BELOW FOR F
+IX
$inum = tmake($i);
$left = tget($inp  $inum  1) + 1;
for ($j = 1; $j < $left; $j++) {
$jnum = tmake($j);
last BING if (tis($inp  $inum  $jnum));
}
}
print " Your triangles are : $inum $jnum "
. ($inp  $inum  $jnum) . "\n";
sub tget { return int(sqrt(1 + $_[0] * 2)  .5); }
sub tmake { return (.5 * $_[0] * ($_[0] + 1)); }
sub tis { my $n = shift; return ($n == tmake(tget($n))); }
This is blazingly fast with numbers even several orders of magnitude larger than 987654321, and requires almost no memory to run.  [reply] [d/l] 

TedPride,
This is blazingly fast...
I wanted to see just how fast, so I modified your code to fit into the benchmark FoxtrotUniform and I were playing with. It turns out your code is only about half as fast as my second version and it has a bug I didn't bother to track down. To see the bug, try entering 3133756 with your code.
I say your code is only about half as fast because it got less than halfway through (line 2235 of 5000) when it blew up and it took approximately the same amount of time (3.619 seconds) to get there.
In case you want to run some more tests yourself, take a look at gen_targets.pl and then use this modified code:
 [reply] [d/l] 

I found and fixed the bug. See my second version. As to it being blazingly fast, it is blazingly fast  just not as blazingly fast as it could be if I took the time to optimize the algorithm. Congrats if yours is faster.
 [reply] 
Re: Triangle Numbers Revisited
by tall_man (Parson) on Oct 14, 2004 at 18:09 UTC

There is a theorem by Gauss which shows that finding M as the sum of three triangle numbers is equivalent to finding 8M+3 as the sum of three odd squares. That is, if
8M + 3 = (2A+1)^2 + (2B+1)^2 + (2C+1)^2
then
M = trinum(A) + trinum(B) + trinum(C)
The computation starts with a bigger number, but the calculations might be easier this way. I'll try it when I get a chance.
Update: Here is code that takes advantage of a few things that are known about Diophantine quadratic problems. For example, it can convert multiples of two into smaller problems. It also screens out a few factors which would make a solution impossible. In many cases, it should be faster than the bruteforce triangular number search.
# Find a trianglenumber 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 bruteforce 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 4k1 fac
+tors.
foreach my $num (@screeners) {
($M, $powercount) = countpowers($M, $num);
if ($powercount % 2 == 1) {
print "Eliminating impossible case with oddpower of 4k1 fac
+tor $num\n";
return;
}
# In even cases, half the factors can be multiplied into the res
+ult.
$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, ij 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 trianglenumber problem to a squarenumber problem.
my $M2 = $M*8 + 3;
# The toplevel 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 trianglenumbers.
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);
}
}
 [reply] [d/l] [select] 

Quick question regarding screening out "a few factors which would make a solution impossible": Didn't Gauss prove Fermat's polygonal theorem for triangle numbers, showing that every positive integer can be represented as a sum of at most three triangle numbers? If so, how can there be cases for which a solution is impossible? (Or are you saying that the method you're using only works for some cases?)
 [reply] 

I'm talking about an intermediate step. I'm looking for three odd squares that add to the number 8*M+3. I pick the first number, k, by brute force working down from the square root. So then I have to solve:
N = 8*M + 3  k^2
i^2 + j^2 = N
There are choices for k that don't work. I want to eliminate them quickly and move on to the next k in the loop instead of spending time trying all combinations of i and j. Eventually I will find an answer, but that choice of k won't be part of it.
For example, if N is a multiple of an odd power of 3, the quadratic problem can't be solved in integers. So I can eliminate about 1/3 of the possible choices for k.  [reply] [d/l] 

I hadn't known it, but you're right, Gauss proved that in his diary (which wasn't discovered until he had been dead 50 years).
 [reply] 

tall_man,
Wow! It took a bit to make your code with my benchmark, but the results were quite impressive (25_000 random targets between 1 .. 9_999_999
$ time ./lr.pl
real 0m21.139s
user 0m20.279s
sys 0m0.040s
$ time ./tm.pl
real 0m3.935s
user 0m3.585s
sys 0m0.300s
 [reply] [d/l] 
Re: Triangle Numbers Revisited
by TedPride (Priest) on Oct 14, 2004 at 11:21 UTC

use strict;
use warnings;
my ($i, $j, $inum, $jnum, $left);
$ = 1; print "Enter the number to check : ";
my $inp = <STDIN>;
BING:
for ($i = tget($inp2); $i; $i) {
$inum = tmake($i);
$left = tget($inp  $inum  1) + 1;
for ($j = 1; $j < $left; $j++) {
$jnum = tmake($j);
last BING if (tis($inp  $inum  $jnum));
}
}
print " Your triangles are : $inum $jnum "
. ($inp  $inum  $jnum) . "\n";
sub tget { return int(sqrt(1 + $_[0] * 2)  .5); }
sub tmake { return (.5 * $_[0] * ($_[0] + 1)); }
sub tis { my $n = shift; return ($n == tmake(tget($n))); }
Incidently, I ran a looped version of this for all numbers 3100000, and only the following could not be made with sums of three triangles:
4, 6, 11, 20, 29
As far as I can tell, all integers after 29 can be made with at least one sum of three triangles. The number of sums increases as you go along.  [reply] [d/l] 

TedPride,
Good catch! I didn't have time to track it down myself this morning as I was supposed to be getting ready for work. For the purposes of the original thread, it appears that 0 was being allowed as a triangular number. With that allowance, all numbers can indeed be the sum of 3 triangles.
As far as speed, you are right that this is quite fast. If it was not fast enough in a real problem you would likely use C as tachyon, the cheater, did ;)
I have no idea about how to figure out the big O of either solution, but I did want to mention that my solution is more than a constant factor faster in case you decided to spend the time optimizing. At 5_000 random targets between 1 and 5_000_000 it was twice as fast, but for 10_000 random targets between 1 and 987_654_321 it was three times as fast. I didn't want to come off as offensive  I just like silly challenges like this.
 [reply] [d/l] 
Re: Triangle Numbers Revisited
by Limbic~Region (Chancellor) on Oct 14, 2004 at 20:49 UTC

All,
I was having a conversation with blokhead in the CB concerning the most efficient solution to this problem. Not having any idea how to determine the big O notation, I ran a few tests to see how well my solution scaled.
The first thing I noticed was the average number of guesses doubled every time I increased the search group by a factor of 10. Then I noticed something really cool  it also corresponded to a power of 2. This allowed me to calculate the average case scenario for any number
So for 12_345 you can expect 16.42, and for 123_456_789 you can expect 262.67 As you can see this scales quite well. While this is the average case, the worst case scenario should still scale reasonably well.
 [reply] [d/l] 
Re: Triangle Numbers Revisited
by tall_man (Parson) on Oct 19, 2004 at 18:33 UTC

Here is a new solution using Math::Pari for bignumber calculations and for factoring. If M is factored as (a^m * b^n * ...), a sumoftwosquares 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 80digit example I included in under a second.
# Find a trianglenumber 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 bruteforce 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 4k1 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 oddpower $power o
+f 4k1 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 ** (($pow1)/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 trianglenumber problem to a squarenumber problem.
my $M2 = $M*8 + 3;
# The toplevel 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 trianglenumbers.
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
 [reply] [d/l] 
Re: Triangle Numbers Revisited
by ambrus (Abbot) on Oct 24, 2013 at 15:30 UTC

 [reply] 

