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 (distincs) 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.