## [Math] Project Euler 108 and 110

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.

###### (the post first published at 20240717.)

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.