Beefy Boxes and Bandwidth Generously Provided by pair Networks
Welcome to the Monastery
 
PerlMonks  

Triangle Numbers Revisited

by Limbic~Region (Chancellor)
on Oct 13, 2004 at 23:37 UTC ( [id://399054]=perlquestion: print w/replies, xml ) Need Help??

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 non-negative 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?

Cheers - L~R

* If anyone wants to know, feel free to ask. If anyone wants to offer their own explanation, feel free to post.

Replies are listed 'Best First'.
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

      All,
      In HTML comments, I have provided a spoiler explanation of how the code works.

      Update: 2013-10-24 When I first wrote this node, we didn't have spoiler tags. The HTML comments remain but I have also added spoiler for those who don't want to view the source.

      Cheers - L~R

      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.

      Cheers - L~R

        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

      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 non-triangular 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.

      Cheers - L~R

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]#
      tachyon,
      Nice, but since you decided to cheat, I am going to make you pay for it (compile time). I changed the gen_targets.pl to 10_000 random targets ranging from 1 .. 999_999_999. My solution completed in 25 seconds while yours took a whopping 36 seconds. I realize you only pay the compiling penalty once unless you change the code or run it from a different directory, but hey - I had to level the playing field some how. Incidently, your code ran in 10 seconds the second time around and would still win with the compiling penalty if the number of targets was large enough. Good Job.

      Cheers - L~R

      Update: If you want to go up to a billion, you are going to need to change the C code. I thought I was doing really well until I realized what was going on.

Re: Triangle Numbers Revisited
by hv (Prior) 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

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.
      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:

      Cheers - L~R

        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.
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 brute-force triangular number search.

    # 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 fac +tors. foreach my $num (@screeners) { ($M, $powercount) = countpowers($M, $num); if ($powercount % 2 == 1) { print "Eliminating impossible case with odd-power of 4k-1 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, 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); } }
      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?)
        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.

        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).
      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

      Cheers - L~R

Re: Triangle Numbers Revisited
by TedPride (Priest) on Oct 14, 2004 at 11:21 UTC
    EDIT: For some reason, a few numbers such as 1035 are causing errors (Can't take sqrt of -1, <STDIN> line 1.) Going to have to hunt down the bug and fix it...

    EDIT: Just needed a slight change to the first line after BING to prevent tget calls for 0:

    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-2); $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 3-100000, 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.

      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.

      Cheers - L~R

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.

    Cheers - L~R

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 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
Re: Triangle Numbers Revisited
by ambrus (Abbot) on Oct 24, 2013 at 15:30 UTC

Log In?
Username:
Password:

What's my password?
Create A New User
Domain Nodelet?
Node Status?
node history
Node Type: perlquestion [id://399054]
Approved by atcroft
Front-paged by FoxtrotUniform
help
Chatterbox?
and the web crawler heard nothing...

How do I use this?Last hourOther CB clients
Other Users?
Others meditating upon the Monastery: (5)
As of 2024-03-19 09:33 GMT
Sections?
Information?
Find Nodes?
Leftovers?
    Voting Booth?

    No recent polls found