Beefy Boxes and Bandwidth Generously Provided by pair Networks
Think about Loose Coupling
 
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. 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


Comment on Re: Triangle Numbers Revisited
Download Code

Log In?
Username:
Password:

What's my password?
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 surveying the Monastery: (6)
As of 2015-07-05 08:01 GMT
Sections?
Information?
Find Nodes?
Leftovers?
    Voting Booth?

    The top three priorities of my open tasks are (in descending order of likelihood to be worked on) ...









    Results (61 votes), past polls