Problem #108 concerns the Diophantine equation 1/x + 1/y = 1/n for arbitrary integers n. The question reads:
Project Euler Problem 108: Diophantine reciprocals I In the following equation x, y, and n are positive integers. 1/x + 1/y = 1/n For n = 4 there are exactly three distinct solutions: 1/5 + 1/20 = 1/4 1/6 + 1/12 = 1/4 1/8 + 1/8 = 1/4 What is the least value of n for which the number of distinct solutions exceeds one-thousand? NOTE: This problem is an easier version of Problem 110; it is strongly advised that you solve this one first.
I believe the most natural approach to solving this problem with knowledge of math is Simon’s Favorite Factoring Trick. Here’s my solution:
Solution #1: SFFT Approach
We can rearrange the given equation to xy = nx+ny. By Simon’s Favorite Factoring Trick, this equation can be factored as (x-n)(y-n) = n^2. It follows that the number of ordered solutions (x,y) is the number of factors of n^2. Thus, it suffices to find the smallest value of n such that n^2 has at least 1999 factors. I did this by using the Sieve of Eratosthenes to find all primes less than an arbitrary upper bound and then using this list of primes to check each potential value of n until a sufficient value was found. Here’s an implementation of this approach in Python 2.7:
'''
Author: Walker Kroubalkian
Factoring Approach to Project Euler Problem #108
'''
import time
from math import sqrt
def countFactorSquare(n):
a = n*n
t = 0
for x in range(1,n):
if(a%x==0):
t+=1
return t+1
def sieveEratosthenes(n):
myPrimes = []
primePossible = [True]*(n+1)
primePossible[0] = False
primePossible[1] = False
for (i,possible) in enumerate(primePossible):
if possible:
for x in range(i*i, (n+1), i):
primePossible[x] = False
myPrimes.append(i)
return myPrimes
def projectEulerProblemOneHundredEight(n,m):
myPrimes = sieveEratosthenes(n)
primes = []
exponents = []
for a in range(n+1):
primes.append([])
exponents.append([])
for x in myPrimes:
for y in range(x,n+1,x):
primes[y].append(x)
exponents[y].append(1)
power = x*x
while(power<=n):
for y in range(power,n+1,power):
l = len(exponents[y])
exponents[y][l-1]+=1
power*=x
totals = []
for x in range(n+1):
myProd = 1
for a in exponents[x]:
myProd*=(2*a+1)
myProd = (myProd+1)/2
totals.append(myProd)
for x in range(n+1):
if(totals[x]>m):
return x
return -1
start = time.time()
print projectEulerProblemOneHundredEight(1000000,1000)
print ("--- %s seconds ---" % (time.time()-start))
'''
Prints
180180
--- 3.04166793823 seconds ---
for input of n = 1000000, m = 1000.
'''
And with that, we’re done. Even with an arbitrary upper bound for primes, this solution took over 3 seconds to run. I may come back to this problem eventually to look for a more efficient approach.
Thanks for reading! See you tomorrow.
The Inner Machinations of Walker 