Tri-48
January 30, 2018
Today’s exercise comes from Albert H. Beiler’s book Recreations in the Theory of Numbers Chapter 1, Problem 2):
One side of a right-angled triangle is 48. Find ten pairs of whole numbers which may represent the other two sides.
Your task is to find the solution to Beiler’s problem. When you are finished, you are welcome to read or run a suggested solution, or to post your own solution or discuss the exercise in the comments below.
Programming Praxis 
In Python. Look here here for solving the second part.
from ma.primee import divisor_gen # generate all divisors z = 48 z2 = z ** 2 # try x^2 + y^2 = z^2 x, y = z, 0 while x >= y > 0: test = x ** 2 + y ** 2 - z2 if test == 0: print(x, y) x -= 1 y += 1 elif test < 0: y += 1 else: x -= 1 # try x^2 - y^2 = z^2 for div in divisor_gen(z2): if div & 1: # as z = 4 * k continue if div >= z: continue u, v = div, z2 // div x = (u + v) // 2 y = (v - u) // 2 print(x, y)A Perl 1-liner…. we know the largest value of a side is (48^2/4 + 1) as {(n+1)^2-(n-1)^2 = 4n = 48^2} so we loop through all values – work out what the other side is… remove entries for which the other side is 0 or non-integer… finally display the results…
print map {"@{[sort{$a<=>$b}@{$_}]}\n"} grep {$_->[2]&&$_->[2]==int $_->[2]} map {[48,$_,sqrt abs$_*$_-2304]} 1..577;This gives:
14 48 50
20 48 52
36 48 60
48 55 73
48 64 80
48 90 102
48 140 148
48 189 195
48 286 290
48 575 577
Just realised there is a minor optimisation in mine to avoid duplicates… – can change the bottom of the loop to 34….
if you want to do for arbitrary N…
$N=shift; print map { "@{[sort{$a<=>$b}@{$_}]}\n"} grep{$_->[2]&&$_->[2]==int$_->[2]} map {[$N,$_,sqrt abs$_*$_-$N*$N]} ($N/sqrt 2)..($N*$N/2+1);This was a fun exercise, thanks! To keep it interesting, i avoided a brute force search. Instead, I relied on Euclid’s formula. https://en.wikipedia.org/wiki/Pythagorean_triple#Generating_a_triple
There is still optimization to be done, but this iterates the innermost loop 71 times.
def primes(): """Generate all prime numbers.""" yield 2 yield 3 yield 5 raise RuntimeError('How many primes are there, anyway?') def factors(n): if n == 1: yield 1 else: for p in primes(): count = 0 while n % p == 0: n //= p count += 1 if count: for f in factors(n): for i in range(count + 1): yield p**i * f break def pythagorean_triples(num): for f in factors(num): k = num // f if f % 2 == 0: for n in factors(f): m = f // 2 // n if m <= n: continue a = m**2 - n**2 b = 2 * m * n assert(b == f) c = m**2 + n**2 yield tuple((k * a, k * b, k * c)) for d in range(1, f*f): md2 = f + d**2 if md2 % (2 * d) == 0: m = md2 // (2 * d) n = m - d if n <= 0: break a = m**2 - n**2 assert(a == f) b = 2 * m * n c = m**2 + n**2 yield (k * a, k * b, k * c) for t in sorted(set(tuple(sorted(t)) for t in pythagorean_triples(48))): print(t) a, b, c = t assert(a**2 + b**2 == c**2)Let’s see: any k,a,b, with a and b coprime and of different parity (so a+b is odd) define a Pythagorean triple (k(a²+b²), 2kab, k(a²-b²)). If n = kj with j = a²+b², j must be an odd divisor of n, expressible as the sum of two squares; if n = 2kab, a and b are any coprime divisors of n/2 of different parity; if n = kj with j = a²-b² = (a+b)(a-b), j is product of two odd factors c,d of n, c > d, with a = (c+d)/2, b = (c-d)/2.
For n = 48, odd factors are just 1 and 3, even factors are 2,4,6,8,12,24,48.
Neither 1 or 3 are the sum of (non-zero) squares.
Even factors of 24 coprime to 1 are: 2,4,8,6,12,24 gives a,b,k = (2,1,12),(4,1,6),(8,1,3),(6,1,4),(12,1,2),(24,1,1)
Even factors of 24 coprime to 3 are: 2,4,8 gives a,b,k = (2,3,4),(4,3,3),(8,3,1)
Single pair of odd factors 1 and 3 give a,b,k = ((3+1)/2,(3-1)/2,12) = (2,1,16)
10 triples total.
Can argue the same way for any number of form 16p with p an odd prime to get 10 different legs. Might also get a hypotenuse if p is the sum of two squares (eg. if n = 80, p = 5, n = 16(1²+2²) – see https://en.wikipedia.org/wiki/Fermat%27s_theorem_on_sums_of_two_squares).
Here’s a brute force Python 2.7 solution that also works for other sized sides.
import itertools # https://en.wikipedia.org/wiki/Integer_square_root#Using_only_integer_division def isqrt(n): # Works for integers n >= 1 x = n while True: y = (x + (n / x)) / 2 if y - x in (0, 1): return x x = y def triples(x): """Returns the other sides of Pythagorean triples that include x.""" x2 = x ** 2 for a in itertools.count(1): a2 = a ** 2 b2 = x2 - a2 # treats x as the hypotenuse if b2 > a2: # avoids duplicates: b = isqrt(b2) # e.g., (3,4) and (4,3) for x == 5 if b * b == b2: yield (a, b) c2 = a2 + x2 # treats x as a leg c = isqrt(c2) if c * c == c2: yield (a, c) if a == c: return for x in triples(48): print xOutput: