The problem in PE 108 is to find such an n for a equation \( \frac{1}{x} + \frac{1}{y} = \frac{1}{n} \), where x, y, n are all integers, and where number of disctinct solutions exceeds 1000.
It was easy, I solved it using Wolfram Mathematica:
countSolutions[n_] :=
Length[Solve[1/x + 1/y == 1/n && x > 0 && y > 0 && x >= y, {x, y},
Integers]]
Do[cnt = countSolutions[n];
If[cnt > 1000, Print[{n, cnt}]; Abort[]], {n, 1, 10^10}]
Slow (several hours), but doable.
PE 110 problem harder -- the bar is raised up to \( 4 \cdot 10^6 \) solutions. No bruteforce is possible.
In the thread on PE about problem 108, I've found a mention about OEIS A018892. I hooked on this clue and tried not to read other solutions and hints. OK, what is this?
A018892 Number of ways to write 1/n as a sum of exactly 2 unit fractions
a(n) = (tau(n^2)+1)/2. Number of elements in the set {(x,y): x|n, y|n, x<=y, gcd(x,y)=1}
In Wolfram Mathematica, it can be expressed as:
A018892[n_] := (Length[Divisors[n^2]] + 1)/2
The problem is now different. Find such a (minimal) number for which A018892 function will exceed \( 4 \cdot 10^6 \). In other words, find a (minimal) number for which Divisors[] function will exceed \( 2 \cdot 4 \cdot 10^6 - 1 \).
Let's ask ourselves, how to find number of divisors? Oh, it's easy, just factorize that number. For example: 12348. Run msieve with it. It will say:
Wed Jul 17 00:56:20 2024 factoring 12348 (5 digits) Wed Jul 17 00:56:20 2024 p1 factor: 2 Wed Jul 17 00:56:20 2024 p1 factor: 2 Wed Jul 17 00:56:20 2024 p1 factor: 3 Wed Jul 17 00:56:20 2024 p1 factor: 3 Wed Jul 17 00:56:20 2024 p1 factor: 7 Wed Jul 17 00:56:20 2024 p1 factor: 7 Wed Jul 17 00:56:20 2024 p1 factor: 7
Indeed, \( 2\cdot 2\cdot 3\cdot 3\cdot 7\cdot 7\cdot 7==12348 \). Another way of representing it is: \( 2^ 2\cdot 3^2 \cdot 7^3 \).
This expression can be viewed as: \( p1^{c1} \cdot p2^{c2} \cdot ...\), where p1, p2 - prime numbers and c1, c2 -- coefficients or powers or exponents.
Number of divisors is: \( (c1+1)(c2+1)(c3+1)... \) For our number, coefficients are 2, 2 and 3. \( (2+1)(2+1)(3+1)=36 \).
Yes, 36 divisors, let's check in Wolfram Mathematica:
In[]:= Length@Divisors[12348] Out[]= 36
Now another question, how can you find a number that has 4099948 solutions? Or \( 2 \cdot 4099948-1=8199895 \) divisors? Ah, it's easy. Factor this number:
In[]:= FactorInteger[2*4099948-1]
Out[]= {{5, 1}, {11, 1}, {29, 1}, {53, 1}, {97, 1}}
There are 5 primes: 5, 11, 29, 53 and 97. But in fact, these are only coefficients to another expression.
We only have to find p1, p2, p3, p4 and p5 for the expression: \( p1^{5-1} \cdot p2^{11-1} \cdot p3^{29-1} \cdot p4^{53-1} \cdot p5^{97-1} \). Let's plug first 5 primes:
In[]:= p1=2; p2=3; p3=5; p4=7; p5=11; In[]:= x=p1^(5-1) * p2^(11-1) * p3^(29-1) * p4^(53-1) * p5^(97-1) Out[]= 291936574071356139797183542300802080273263700534750232097183894934260994188740440655397\ 62795433862023628730959481433237126302432907576845964140240483880043029785156250000 In[]:= Length@Divisors[x] Out[]= 8199895
In other words, number of solutions of equation \( \frac{1}{x} + \frac{1}{y} = \frac{1}{n} \) is 4099948, where n is:
In[]:= Sqrt[x] Out[]= 5403115527835363060896154958947115524913225008014850942090010084438627775073242187500
So far so good. But this number is too big.
Now the problem is to find such a group of 5 (distinct) primes for which that number will be minimized.
Simplest idea -- just sort that list, let first (smallest) prime have biggest coefficient and last (biggest) prime have smallest coefficient.
#!/usr/bin/python3
import math
# https://math.stackexchange.com/questions/1039519/finding-prime-factors-by-taking-the-square-root
# DIY factoring functions, but for small numbers it is OK
def factor_helper(n):
for i in range(2, int(math.ceil(math.sqrt(n)))+1):
if n % i == 0:
return i
# this is prime
return n
def factor(n):
rt=[]
while True:
p=factor_helper(n)
if p==1:
break
rt.append (p)
n=n//p
return rt
first_primes=[2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97]
smallest_result=-1
smallest_divisors=-1
smallest_solutions=-1
def try_(solutions_we_want, divisors_we_want):
global smallest_result
global smallest_divisors
global smallest_solutions
factors=factor(divisors_we_want)
if max(factors)>100:
# skip it, because we need small results anyway
return
factors=sorted(factors, reverse=True)
rt=1
for p, f in zip(first_primes, factors):
rt*=p**(f-1)
try:
result=int(math.sqrt(rt))
except OverflowError:
# skip it, because we need small results anyway
return
print (divisors_we_want, result)
if (smallest_result==-1) or (result<smallest_result):
smallest_result=result
smallest_divisors=divisors_we_want
smallest_solutions=solutions_we_want
#for solutions_we_want in range(1000, 2000): # for PE 108
for solutions_we_want in range(4*10**6, 4*10**6+100000): # for PE 110
divisors_we_want=(solutions_we_want*2)-1;
try_ (solutions_we_want, divisors_we_want)
print (f"{smallest_result=}")
print (f"{smallest_solutions=}")
print (f"{smallest_divisors=}")
And that worked. For PE 108:
smallest_result=180180 smallest_solutions=1013 smallest_divisors=2025
For PE 110:
smallest_result=9350130049860600 smallest_solutions=4018613 smallest_divisors=8037225
And both 180180 and 9350130049860600 numbers are happens to be there:
A126098 Where records occur in A018892.
As if Project Euler author took idea right from here. Maybe it was so. But the code at OEIS is simpler than mine.
Yes, I know about these lousy Disqus ads. Please use adblocker. I would consider to subscribe to 'pro' version of Disqus if the signal/noise ratio in comments would be good enough.