projectEuler 454 -- Diophantine reciprocals

wayne

设对于给定的 \(L\)，等式 \(\D\frac1X + \frac1Y = \frac1N, (0 \lt X \lt Y \leqslant L)\) 成立总共有 \(f(L)\) 种情况 .
举几个简单的例子, 有 \(f(6) =1, f(12)=3, f(1000)=1069\), 求\(f(10^{12})\)

参考链接:  http://projecteuler.net/problem=454

Lwins_G

\[\frac1x+\frac1y=\frac1n \iff x+y\,|\, xy \iff \left. (x,y)\,\middle|\,\frac{x+y}{(x,y)} \right.\]

wayne

不错的条件,   能否继续前进,导出计数公式
ps,好像写反了吧:\[\left. \frac{x+y}{(x,y)}\,\middle|\,(x,y) \right.\]
不过,在此基础上,可以导出这种一般性的解表达:   \(x=km(m+n) ,y=kn(m+n), (m,n)=1, m<n\), 即\(L=kn(m+n)\)

但要计算出$f(10^12)$,还需要进一步的挖掘,优化

Lwins_G

直接计算\[ \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] \]只需要`\mathcal{O}(L)`的时间，我认为这是可以承受的。不过我们也可试着优化一下：
\[\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}
\]
计算`\displaystyle \psi(x,y) = \sum_{n=1}^{y} \left[ \frac{x}{n} \right]`显然是`\mathcal{O}(\sqrt{x})`的。
于是使用该公式可以做到在`\mathcal{O}(L^{3/4})`的复杂度内计算。

wayne

写了个Mathematica函数, 计算10^7花了40秒钟,再算下去不太现实呢,或许C/C++语言实现要快一些

1. L=10^4;phi=Sum[Floor[#1/n],{n,1,#2}]&;
   
2. Sum[MoebiusMu[j](phi[Floor[L/(j i)],(2i-1)/j]-phi[Floor[L/(j i)],i/j]),{i,1,Sqrt[L]},{j,Divisors[i]}]

复制代码

Lwins_G

wayne 发表于 2014-3-21 13:34
写了个Mathematica函数, 计算10^7花了40秒钟,再算下去不太现实呢,或许C/C++语言实现要快一些

经测试，Mathematica计算`\psi(x,y)`时所消耗的时间是平凡的，为`\mathcal{O}(y)`，这会导致复杂度退化为
\[ \sum_{n=1}^{\sqrt{L}} \sum_{d | n} \mathcal{O} \left( \frac{n}{d} \right) = \mathcal{O}(L) \]

wayne

Lwins_G 发表于 2014-3-21 19:18
经测试，Mathematica计算`\psi(x,y)`时所消耗的时间是平凡的，为`\mathcal{O}(y)`，这会导致复杂度退化 ...

可以将两个相减的合成一个, 但这也只是常数级的优化,时间砍了一半:

1. L=10^6;Sum[MoebiusMu[j]Sum[Floor[L/(j i n)],{n,i/j,(2i-1)/j}],{i,1,Sqrt[L]},{j,Divisors[i]}]-L

复制代码

Ruthles

Lwins_G 发表于 2014-3-21 12:01
直接计算\[ \sum_{n=1}^{\sqrt{L}} \sum_{\begin{subarray}{cc} m=1 \\ (n,m)=1 \end{subarray}}^{n-1} \le ...

(x,y)中x是10^12的级别的。

cn8888

1/a+1/b=1/c
等价于
(c-a)*(c-b)=c^2
这样问题是不是简单了???????????

cn8888

(c-a)*(c-b)=c^2写成(a-c)*(b-c)=c^2这种形式更好