Beefy Boxes and Bandwidth Generously Provided by pair Networks
Perl-Sensitive Sunglasses
 
PerlMonks  

Comment on

( #3333=superdoc: print w/ replies, xml ) Need Help??
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

In reply to Re: Triangle Numbers Revisited by tall_man
in thread Triangle Numbers Revisited by Limbic~Region

Title:
Use:  <p> text here (a paragraph) </p>
and:  <code> code here </code>
to format your post; it's "PerlMonks-approved HTML":



  • Posts are HTML formatted. Put <p> </p> tags around your paragraphs. Put <code> </code> tags around your code and data!
  • Read Where should I post X? if you're not absolutely sure you're posting in the right place.
  • Please read these before you post! —
  • Posts may use any of the Perl Monks Approved HTML tags:
    a, abbr, b, big, blockquote, br, caption, center, col, colgroup, dd, del, div, dl, dt, em, font, h1, h2, h3, h4, h5, h6, hr, i, ins, li, ol, p, pre, readmore, small, span, spoiler, strike, strong, sub, sup, table, tbody, td, tfoot, th, thead, tr, tt, u, ul, wbr
  • Outside of code tags, you may need to use entities for some characters:
            For:     Use:
    & &amp;
    < &lt;
    > &gt;
    [ &#91;
    ] &#93;
  • Link using PerlMonks shortcuts! What shortcuts can I use for linking?
  • See Writeup Formatting Tips and other pages linked from there for more info.
  • Log In?
    Username:
    Password:

    What's my password?
    Create A New User
    Chatterbox?
    and the web crawler heard nothing...

    How do I use this? | Other CB clients
    Other Users?
    Others scrutinizing the Monastery: (6)
    As of 2014-11-28 07:37 GMT
    Sections?
    Information?
    Find Nodes?
    Leftovers?
      Voting Booth?

      My preferred Perl binaries come from:














      Results (194 votes), past polls