Tri-48
January 30, 2018
There are two possibilities. If 48 is the hypotenuse of the triangle, the problem is to find two integers with squares that sum to a target value, which we examined in a previous exercise; in this case, there are no solutions in triangles:
(define (squares n)
(let loop ((x (isqrt n)) (y 0) (zs '()))
(cond ((< x y) zs)
((< (+ (* x x) (* y y)) n) (loop x (+ y 1) zs))
((< n (+ (* x x) (* y y))) (loop (- x 1) y zs))
(else (loop (- x 1) (+ y 1) (cons (list x y) zs))))))
> (squares (* 48 48)) ((48 0))
The other possibility is that 48 is one of the other legs of the triangle, in which case we count up from 1, testing at each step if the hypotenuse is a square:
(define (tri n)
(let loop ((a 1))
(let ((c (+ (* a a) (* n n))))
(when (square? c)
(display (sort < (list a (isqrt c) n)))
(newline)))
(loop (+ a 1))))
> (tri 48) (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) CTRL-C break>
This is intentionally an infinite loop. Beiler shows that any number expressible as 16k (such as 48) always has exactly ten solutions, so the list is exhaustive.
You can run the program at https://ideone.com/yX11Fm:
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: