Project Euler Problem #9:

Problem #9 features what is possibly the most famous theorem in all of mathematics: The Pythagorean Theorem. However, unlike typical discussion of the Pythagorean Theorem which involves right triangles in Euclidean geometry, this problem is more concerned with investigating the theorem as a Diophantine Equation. The problem reads:

Project Euler Problem 9: Special Pythagorean triplet 
A Pythagorean triplet is a set of three natural numbers, a < b < c, for which

a^2 + b^2 = c^2 

For example, 3^2+4^2 = 9 + 16 = 25 = 5^2. 

There exists exactly one Pythagorean triplet for which a+b+c = 1000. 
Find the product abc.

Luckily, there are some useful facts about Pythagorean triplets which we can use to our advantage while interpreting the theorem as a Diophantine equation. I chose to generalize this problem by finding the maximum product of all Pythagorean triplets with a given sum or returning -1 if no such triplets exist. Here is the solution I found:

Solution #1: General Pythagorean Triplet

As it turns out, prior investigation of Pythagorean triples has resulted in a useful method of generating Pythagorean triples. One can easily verify that for all positive integers m and n such that m>n, the triple (m^2-n^2, 2mn, m^2+n^2) is a Pythagorean triple. According to Wikipedia, this generator is one of the many facts that is called Euclid’s Formula.

It is a bit more difficult to show that this generator can generate all primitive Pythagorean triples. A primitive Pythagorean triple is a triple (a,b,c) which satisfies the condition and the extra condition that gcd(a,b,c) = 1. I will not write the proof of this fact here. Unfortunately, this generator is not perfect, as some non-primitive triples such as (9, 12, 15) cannot be generated by Euclid’s Formula. Instead, we must use the generalized generator for Pythagorean triplets: (d(m^2-n^2), d(2mn), d(m^2+n^2)) for positive integers d, m, and n such that m>n.

Using this generator, the sum of the sides of a general Pythagorean triple is d(2m^2 + 2mn) = 2dm(m+n). Thus, we have reduced this diophantine equation to finding integer solutions to 2dm(m+n) = 1000 such that m>n.

I attempted to solve this Diophantine equation by iterating through the possible values of d and then iterating through the possible values of m. It is worth noting that m < m+n < 2m by the definition of the Pythagorean triplet. Thus, m^2 < 1000/(2d) < 2m^2. Using this fact, we can bound the possible values of m.

Here is an implementation of this method in Python 2.7:

'''
Author: Walker Kroubalkian
General Pythagorean Triple Approach to Project Euler Problem #8
'''

import time
from math import sqrt

def listFactorsGivenFactorization(p,e):
    total = []
    if(len(p) >= 1):
        for x in range(e[0]+1):
            total.append(p[0]**x)
    else:
        return [1]
    previousTotal = listFactorsGivenFactorization(p[1:],e[1:])
    actualTotal = []
    for a in total:
        for b in previousTotal:
            actualTotal.append(a*b)
    return actualTotal

def listFactors(n):
    primes = []
    exponents = []
    c = 2
    while(c<=sqrt(n)):
        if(n%c==0):
            primes.append(c)
            t = 0
            while(n%c==0):
                n/=c
                t+=1
            exponents.append(t)
        c+=1
    if(n>1):
        primes.append(n)
        exponents.append(1)
    return sorted(listFactorsGivenFactorization(primes,exponents))

def projectEulerProblemNine(n):
    if(n%2==1):
        return -1
    maxProduct = -1
    possibleD = listFactors(n/2)
    n/=2
    for d in possibleD:
        newProduct = n/d
        lowerBound = int(sqrt(newProduct/2))
        upperBound = int(sqrt(newProduct))
        for m in range(lowerBound,upperBound+1):
            if(m>0 and newProduct%m==0):
                o = newProduct/m - m
                if(0<o<m):
                    possible = d*d*d*(m*m-o*o)*(2*m*o)*(m*m+o*o)
                    if(possible>maxProduct):
                        maxProduct = possible
    return maxProduct

start = time.time()
print projectEulerProblemNine(1000)
print ("--- %s seconds ---" % (time.time()-start))

'''
Prints

31875000
--- 4.91142272949e-05 seconds ---

for input of n=1000
'''

As shown above, even with a very sloppy implementation of this method, using math can make solving complex problems far more efficient.

Thanks for reading! See you tomorrow.