Finals are approaching, and many students have started reviewing. I am not yet very familiar with several courses I took this semester and intended to use this time to study properly. However, five days ago, while browsing the programming arena of the Mathematics Research Forum, I came across this specific problem:
For a given L, let f(L) be the number of positive integer solutions to the equation \frac{1}{x}+\frac{1}{y}=\frac{1}{n} satisfying 0 < x < y \leq L. For example, f(6)=1, f(12)=3, f(1000)=1069. Find f(10^{12}).
The source of this problem is Project Euler Problem 454: Diophantine reciprocals III. The problem is concise and easy to understand, yet possesses significant depth, perfectly fitting my definition of an ideal problem. Moreover, I have been enjoying learning Python recently and felt eager to try this challenge. Consequently, my last five days vanished, and in the process, I exhausted almost every programming trick I know. Due to continuous testing and execution, my computer generated several times more heat than usual—it has truly worked hard. The final code is, in my opinion, the most brilliant I have written to date. I share it here with everyone for mutual encouragement.
The expression above involves fractions, which is not ideal for programming. Since n=\frac{xy}{x+y}, the problem is equivalent to finding integer solutions where (x+y)|xy (meaning x+y divides xy).
The Rambling Phase
At this stage, the most direct programming technique is a double loop. In Python, it looks like this:
l = 1000 # Set upper limit
s = 0 # Count
for y in xrange(2, l + 1):
for x in xrange(1, y):
if (x * y) % (x + y) == 0:
s = s + 1
print(s)This works fine for small L, but it becomes extremely problematic as L grows (calculating f(10^5) already takes quite a while). Note that we perform \frac{L(L-1)}{2} divisions, meaning the algorithm is at least \mathcal{O}(L^2). This implies that if L increases tenfold, the execution time increases by 100 times! Since we need to calculate for L=10^{12}, the efficiency of this algorithm is unacceptable.
Optimization
The content in this section is largely based on insights from the user Lwins_G on the Mathematics Research Forum.
First, the algorithm can be reduced to \mathcal{O}(L) in one step. It is not difficult to prove that (x+y)|xy if and only if there exist integers m, n, k such that (m, n)=1 and: x=km(m+n), \quad y=kn(m+n) Thus, we only need to find solutions for: kn(m+n) \leq L, \quad (m, n)=1, \quad m < n The formula can be written as: f(L)=\sum_{n=1}^{\sqrt{L}} \sum_{\begin{subarray}{cc} m=1 \\ (n,m)=1 \end{subarray}}^{n-1} \left[ \frac{L}{n(n+m)} \right] This formula is \mathcal{O}(L). Code 1 in the appendix is the algorithm I wrote based on this formula in Python. The trick to using this formula is, given n, to construct a table of integers coprime to n, which is similar to constructing a prime table. If one were to iterate through all m and check if (m, n)=1, efficiency would drop significantly. Running in PyPy, calculating f(10^7)=30093331 took 1.7 seconds, f(10^8)=349446716 took 17 seconds, and f(10^9)=3979600400 took 190 seconds. While this is much faster than the double loop, f(10^{12}) is still a long way off.
Further Optimization
Lwins_G also proposed that the algorithm could be reduced to \mathcal{O}(L^{3/4}) as follows: \begin{split} &\sum_{n=1}^{\sqrt{L}} \sum_{\begin{subarray}{cc} m=1 \\ (n,m)=1 \end{subarray}}^{n-1} \left[ \frac{L}{n(n+m)} \right] \\ =&\sum_{n=1}^{\sqrt{L}} \sum_{\begin{subarray}{cc} s=n+1 \\ (n,s)=1 \end{subarray}}^{2n-1} \left[ \frac{L}{ns} \right] \\ =&\sum_{n=1}^{\sqrt{L}} \sum_{s=n+1}^{2n-1} \left[ \frac{L}{ns} \right] \left[ \frac{1}{(n,s)} \right] \\ =&\sum_{n=1}^{\sqrt{L}} \sum_{s=n+1}^{2n-1} \left[ \frac{L}{ns} \right] \sum_{d | (n,s)} \mu(d) \\ =&\sum_{n=1}^{\sqrt{L}} \sum_{d | n} \sum_{\begin{subarray}{cc} n < s \leq 2n-1 \\ d | s \end{subarray}} \left[ \frac{L}{ns} \right] \mu(d) \\ =&\sum_{n=1}^{\sqrt{L}} \sum_{d | n} \mu(d) \sum_{\frac{n}{d} < s' \leq \frac{2n-1}{d}} \left[ \frac{L}{nds'} \right] \\ =&\sum_{n=1}^{\sqrt{L}} \sum_{d | n} \mu(d) \left( \psi\left( \left[ \frac{L}{nd} \right],\frac{2n-1}{d}\right)- \psi\left(\left[ \frac{L}{nd} \right],\frac{n}{d}\right) \right) \end{split} Where \mu(d) is the Möbius function, and \psi(x,y) = \sum_{n=1}^{y} \left[ \frac{x}{n} \right]. This process is brilliant, cleverly incorporating the Möbius function to (partially) eliminate the (m, n)=1 constraint. Since calculating \psi(x,y) is \mathcal{O}(\sqrt{x}), the final formula is naturally \mathcal{O}(L^{3/4}).
Initially, I didn’t understand why \psi(x,y) is \mathcal{O}(\sqrt{x}). So, I programmed using the formula just before the last equals sign (Code 2). The problem involves the Möbius function, which is best handled by construction. It also involves the d|n constraint, requiring factorization. However, since n only goes up to \sqrt{L}=10^6, it isn’t too large, so trial division is sufficient. This requires preparing a prime table up to 10^3 (which is easy and takes less than half a second). Many values of the Möbius function are zero; we only need to find non-zero terms, which are generated by multiplying distinct prime factors of n. This is done using a recursive construction. The advantage of construction is that it uses multiplication, whereas checking requires division, which is much slower. Furthermore, construction significantly reduces the number of iterations.
The efficiency of Code 2 improved significantly, calculating f(10^9) in just 24 seconds. Strangely, however, when calculating f(10^{10}), it ran for half an hour without a result. I had to optimize further. I realized the main issue was that I didn’t know the \mathcal{O}(\sqrt{x}) algorithm for \psi(x,y) mentioned by Lwins_G; I was just calculating it directly.
Heading Straight for 10^{12}
Yesterday, I tried to optimize some details of Code 2, but the speed didn’t improve noticeably. Today, I continued to think about the \mathcal{O}(\sqrt{x}) algorithm for \psi(x,y). Finally, I figured it out. I had noticed a pattern a few days ago: for \left[\frac{x}{n}\right], if n \leq \sqrt{x}, you have to calculate it directly. But if n > \sqrt{x}, the results start appearing in “batches.” That is, \left[\frac{x}{n}\right] slowly decreases from \left[\sqrt{x}\right] down to 1. By looking at the value of the quotient, we can determine the range of the divisor. For example, \left[\frac{101}{1}\right]=101, \left[\frac{101}{2}\right]=50. We then know that for n from 51 to 101, \left[\frac{101}{n}\right] is always 1. This allows us to calculate 51 divisions in one step.
This optimization is a qualitative leap. Our calculation went from \mathcal{O}(x) to \mathcal{O}(\sqrt{x}), improving the total efficiency from \mathcal{O}(L) to \mathcal{O}(L^{3/4}). One might think the difference between \mathcal{O}(L) and \mathcal{O}(L^{3/4}) is small, but when L=10^{12}, the difference is 1000-fold—three orders of magnitude!
I spent most of the afternoon writing the \psi(x,y) function, even skipping dinner. Using the improved \psi(x,y) algorithm, calculating f(10^9) became trivial, taking less than 3 seconds. However, calculating f(10^{10})=44647347052 took 180 seconds. Regardless, this was the first time I reached 10^{10}, which was exciting. I was puzzled by the "explosion" in time and feared f(10^{11}) was hopeless. But when I tried f(2\times 10^{10}) and f(3\times 10^{10}), the time increase was small; f(3\times 10^{10}) took only 313 seconds, just over double f(10^{10}). This surprised me, so I went straight for f(10^{11})=494986959815. The result: 782 seconds. There was hope!
It seemed that the time for f(10^{10}) was the real baseline, and from there, it followed the \mathcal{O}(L^{3/4}) growth. If so, f(10^{12}) might take about 5000 seconds, or an hour and a half—acceptable.
However, I considered improving the algorithm further because while \psi(x,y) was efficient, it performed redundant work. We only needed the difference between two \psi values. Just as S_n = a_1 + \dots + a_n, to calculate S_{n+k} - S_n, you don’t necessarily need to compute both S_{n+k} and S_n from scratch. I spent most of the evening refining \psi(x,y), resulting in the final version of the code. Testing f(10^{10}) took 6 seconds, and f(10^{11}) took 23 seconds!! This far exceeded my expectations!! So, I headed straight for f(10^{12}).
It was past 9 PM. After starting the code, I turned off the screen and went to shower. I was nervous about unexpected issues—I didn’t know if the time would grow linearly or explode like the jump from f(10^9) to f(10^{10}).
After my shower, I took my time hanging the laundry, sat down, and turned on the screen. The number before me stunned me—
101.629625957 seconds!
f(10^{12})=5435004633092! It was an incredible surprise and a huge thrill! Success! I thought it would take at least dozens of minutes, but it took less than 2 minutes! This is the result of five days of hard work—exhausting almost all my skills! The arena has been conquered! (It seems this algorithm truly achieved \mathcal{O}(L^{3/4}); from 10^{10} to 10^{11} to 10^{12}, the time increased by a factor of 5, which matches \mathcal{O}(L^{3/4}) since 10^{3/4} \approx 5.62. Update: I later tested f(10^{13})=59201396855810, which took only 476.7 seconds, confirming this pattern.)
Below are some values for reference and comparison: \begin{aligned} f(10^2)&=60, & f(10^3)&=1069\\ f(10^4)&=15547, & f(10^5)&=203931\\ f(10^6)&=2524207, & f(10^7)&=30093331\\ f(10^8)&=349446716, & f(10^9)&=3979600400\\ f(10^{10})&=44647347052, & f(10^{11})&=494986959815\\ f(10^{12})&=5435004633092, & f(10^{13})&=59201396855810 \end{aligned}
Code 1: \mathcal{O}(L) Implementation
import time
import math
start = time.clock()
l = 10000000 # Upper limit
s = 0 # Count
rl = int(math.sqrt(l))
rrl = int(math.sqrt(int(math.sqrt(l))))
# Generate prime table
prime = [i for i in range(1, rl + 1)]
j = 2
while j <= rrl:
if prime[j-1] != 0:
m = int(rl / j)
for i in range(j, m + 1):
prime[i * j - 1] = 0
j = j + 1
else:
j = j + 1
prime.sort()
z = prime.count(0)
prime = prime[z + 1:rl]
for m in range(1, int(rl / 1.4)):
# Find numbers coprime to m, store in p
p = [i for i in range(1, rl + 1)]
m2 = m
for i in prime:
if i > m:
break
if m % i == 0:
m = m / i
for j in range(m2 // i + 1, int(rl / i) + 1):
p[i * j - 1] = 0
p.sort()
z = p.count(0)
p = p[z + m2:rl]
for n in p:
s = s + l // ((n + m2) * n)
print(s)
end = time.clock()
print(end - start)Code 2: \mathcal{O}(L^{3/4}) Initial Implementation
import time
import math
start = time.clock()
def psi(l, n, d):
sum = 0
for s in xrange(n // d + 1, (2 * n - 1) // d + 1):
sum = sum + l // (n * d * s)
return sum
def mobi(p):
length = len(p)
if length == 1:
mu = [1, -p[0]]
else:
mu2 = mobi(p[0:length - 1])
mu = mu2[:]
for i in mu2:
mu.append(-1 * p[length - 1] * i)
return mu
l = 2000000000
s = 0
rl = int(math.sqrt(l))
rrl = int(math.sqrt(int(math.sqrt(l))))
rrrl = int(math.sqrt(int(math.sqrt(int(math.sqrt(l))))))
prime = [i for i in xrange(1, rrl + 1)]
j = 2
while j <= rrrl:
if prime[j - 1] != 0:
m = int(rrl / j)
for i in xrange(j, m + 1):
prime[i * j - 1] = 0
j = j + 1
else:
j = j + 1
prime.sort()
z = prime.count(0)
prime = prime[z + 1:rl]
prime.append(rl)
for n in xrange(2, rl + 1):
p = []
nn = n
i = 0
rn = int(math.sqrt(n))
while prime[i] <= rn:
if n % prime[i] == 0:
n = n // prime[i]
if n % prime[i] != 0:
p.append(prime[i])
i = i + 1
else: i = i + 1
if n > rn:
p.append(n)
mu = mobi(p)
for d in mu:
s = s + d // abs(d) * psi(l, nn, abs(d))
print(s)
end = time.clock()
print(end - start)Final Optimized Code
import time
import math
start = time.clock()
def psi(x, y1, y2):
s = 0
rx = int(math.sqrt(x))
if y2 <= rx:
for d in xrange(y1, y2 + 1):
s = s + x // d
elif y1 <= rx:
for d in xrange(y1, rx + 1):
s = s + x // d
p = x // rx - 1
q = x // y2
d = 0
q1 = x // (p - d + 1)
while d < p - q:
q2 = q1
q1 = x // (p - d)
s = s + (q1 - q2) * (p - d)
d = d + 1
s = s + (y2 - x // (q + 1)) * q
else:
p = x // y1 - 1
q = x // y2
if p < q:
s = (1 + y2 - y1) * q
else:
d = 0
q1 = x // (p - d + 1)
while d < p - q:
q2 = q1
q1 = x // (p - d)
s = s + (q1 - q2) * (p - d)
d = d + 1
s = s + (y2 - x // (q + 1)) * q
s = s + (x // (p + 1) - y1 + 1) * (p + 1)
return s
def mobi(p):
length = len(p)
mu = [1, -p[0]]
for i in xrange(1, length):
mu2 = mu[:]
for j in mu2:
mu.append(-1 * j * p[i])
return mu
l = 1000000000000
s = 0
rl = int(math.sqrt(l))
rrl = int(math.sqrt(int(math.sqrt(l))))
rrrl = int(math.sqrt(int(math.sqrt(int(math.sqrt(l))))))
prime = [i for i in xrange(1, rrl + 1)]
j = 2
while j <= rrrl:
if prime[j - 1] != 0:
m = int(rrl / j)
for i in xrange(j, m + 1):
prime[i * j - 1] = 0
j = j + 1
else:
j = j + 1
prime.sort()
z = prime.count(0)
prime = prime[z + 1:rl]
prime.append(rl)
for n in xrange(2, rl + 1):
p = []
nn = n
i = 0
rn = int(math.sqrt(n))
while prime[i] <= rn:
if n % prime[i] == 0:
n = n // prime[i]
if n % prime[i] != 0:
p.append(prime[i])
i = i + 1
else: i = i + 1
if n > rn:
p.append(n)
mu = mobi(p)
for d in mu:
ll = l // (nn * abs(d))
s = s + d // abs(d) * (psi(ll, nn // abs(d) + 1, (2 * nn - 1) // abs(d)))
print(s)
end = time.clock()
print(end - start)Original Address: https://kexue.fm/archives/2665
For more details on reprinting, please refer to: Scientific Space FAQ