8M + 3 = (2A+1)^2 + (2B+1)^2 + (2C+1)^2 then M = trinum(A) + trinum(B) + trinum(C) #### # 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 factors. foreach my \$num (@screeners) { (\$M, \$powercount) = countpowers(\$M, \$num); if (\$powercount % 2 == 1) { print "Eliminating impossible case with odd-power of 4k-1 factor \$num\n"; return; } # In even cases, half the factors can be multiplied into the result. \$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); } }