倒数丢番图方程

Sat, 25th January 2020Edit on Github丢番图倒数

wayne于2014年3月转载projectEuler 454问题:
设对于给定的 LL，等式 1X+1Y=1N,(0<X<Y⩽L)\frac1X + \frac1Y = \frac1N, (0 \lt X \lt Y \leqslant L) 成立总共有 f(L)f(L) 种情况 . 举几个简单的例子, 有 f(6)=1,f(12)=3,f(1000)=1069f(6) =1, f(12)=3, f(1000)=1069, 求f(1012)f(10^{12}).
参考链接

公式推导

wayne指出：
x+y(x,y)∣(x,y)\frac{x+y}{(x,y)}\mid (x,y)
不过,在此基础上,可以导出这种一般性的解表达: x=km(m+n),y=kn(m+n),(m,n)=1,m<nx=km(m+n) ,y=kn(m+n), (m,n)=1, m\lt n, 即L=kn(m+n)L=kn(m+n)

但要计算出f(1012)f(10^{12}),还需要进一步的挖掘,优化.
LwinsG给出进一步公式推导过程:
直接计算$ \sum{n=1}^{\sqrt{L}} \sum_{m=1,(n,m)=1}^{n-1} [\frac{L}{n(n+m)}] 只需要只需要O(L)$的时间，我认为这是可以承受的。不过我们也可试着优化一下：

∑n=1L∑m=1,(n,m)=1n−1[Ln(n+m)]=∑n=1L∑s=n+1,(n,s)=12n−1[Lns]=∑n=1L∑s=n+12n−1[Lns][1(n,s)]=∑n=1L∑s=n+12n−1[Lns]∑d∣(n,s)μ(d)=∑n=1L∑d∣n∑n<s≤2n−1,d∣s[Lns]μ(d)=∑n=1L∑d∣nμ(d)∑nd<s′≤2n−1d[Lnds′]=∑n=1L∑d∣nμ(d)(ψ([Lnd],2n−1d)−ψ([Lnd],nd))\sum_{n=1}^{\sqrt{L}} \sum_{ m=1,(n,m)=1 }^{n-1} [ \frac{L}{n(n+m)} ] =\sum_{n=1}^{\sqrt{L}} \sum_{s=n+1 , (n,s)=1 }^{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_{ n < s \leq 2n-1, d | s } \left[ \frac{L}{ns} \right] \mu(d) =\sum_{n=1}^{\sqrt{L}} \sum_{d | n} \mu(d) \sum_{\frac{n}{d} < s^{\prime} \leq \frac{2n-1}{d}} \left[ \frac{L}{nds^{\prime}} \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)

计算ψ(x,y)=∑n=1y[xn]\psi(x,y) = \sum_{n=1}^{y} \left[ \frac{x}{n} \right]显然是O(x)O(\sqrt{x})的。 于是使用该公式可以做到在O(L3/4)O(L^{3/4})的复杂度内计算。

cn8888给出了这个丢番图方程的几何意义: Diophrecip
倪举鹏表示:
还有种几何关系满足这个等式：平面过一点的三射线中间射线与两边射线夹角都是60度。一条直线截这3射线为三段，这三段长度满足这个等式.

sunwukong认为:
因为
1x+1y=1n\frac1x+\frac1y=\frac1n，0<x<y≤L0\lt x\lt y\leq L
所以 n<x<y≤Ln\lt x\lt y\leq L
设 x=n+sx=n+s，y=n+ty=n+t，则1≤s<t≤L−n1\leq s\lt t\leq L-n
1n+s+1n+t=1n  ⟺  n2=st\frac1{n+s}+\frac1{n+t}=\frac1{n} \iff n^2=st
所以
n2=st≥1×2=2⇒n≥2n^2=st\geq1\times 2=2 \Rightarrow n\geq 2
s=n2t≥n2L−ns=\frac{n^2}t\geq \frac{n^2}{L-n}
n2L−n≤s<n<t≤L−n\frac{n^2}{L-n}\leq s\lt n\lt t\leq L-n
1n=1x+1y>1L+1L=2L⇒n<L2\frac1{n}=\frac1x+\frac1y\gt \frac1L+\frac1L=\frac2L \Rightarrow n\lt \frac{L}{2}
所以
f(L)f(L)
=∑n=2⌊L−12⌋∑s≥n2L−n,s∣n2n−11=\sum_{n=2}^{\lfloor \frac{L-1}{2} \rfloor} \sum_{ s\geq\frac{n^2}{L-n} , s| n^2}^{n-1} 1
=∑n=2⌊L−12⌋∑t=n+1,t∣n2L−n1=\sum_{n=2}^{\lfloor \frac{L-1}{2} \rfloor} \sum_{ t=n+1, t|n^2 }^{L-n} 1
=∑n=2⌊L−12⌋(∑u≥n2L−n,u∣n2L−n1)−12=\sum_{n=2}^{\lfloor \frac{L-1}{2} \rfloor} \frac{(\sum_{ u\geq\frac{n^2}{L-n} , u| n^2}^{L-n} 1)-1}2

关键的是
给定 nn，LL，求满足
t∣n2t| n^2，n<t≤L−nn\lt t\leq L-n
的tt的数目
这一步的复杂度与n2n^2的质因数分解相同

kastin指出:
其实就是涉及大数分解，虽然大数分解本身没有什么高效率的算法，但是这个题目只是统计解的个数，而非找到每一组解。 联想到给出不大于xx的素数个数π(x)\pi(x)的算法目前最高效的算法中，并不需要找出所有小于x的素数然后统计个数（见M. DELEGLISE AND J. RIVAT.COMPUTING  (x): THE MEISSEL, LEHMER, LAGARIAS,MILLER, ODLYZKO METHOD) 因而，n2n^2满足一定限制要求的分解方案数，应该有相应不需要找出每种具体方案的高效算法。

开始计算

wayne试着开始计算

L=10^4;phi=Sum[Floor[#1/n],{n,1,#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]}]

但是他发现这份Mathematica代码计算10710^7就要40秒

Lwins_G测试后认为Mathematica计算ψ(x,y)\psi(x,y)时所消耗的时间是平凡的，为O(y)O(y)，这会导致复杂度退化为

∑n=1L∑d∣nO(nd)=O(L)\sum\limits_{n=1}^{\sqrt{L}}\sum\limits_{d\mid n} O(\frac nd)=O(L)

所以需要自己重新实现。

282842712474直接使用Lwins_G的公式用phyton计算出f(107)=30093331,f(108)=349446716f(10^7)=30093331,\quad f(10^8)=349446716，前者用时1.7秒，后者用时17秒。 f(109)=3979600400f(10^9)=3979600400，用时196.748994145秒... 按这个速度101210^{12}预计需要3天时间才能计算出来。
对应的python代码
然后他继续实现O(L34)O(L^{\frac34})复杂度的代码, 测试运行f(109)f(10^9)提升到24秒。 最终他实现了O(L)O(\sqrt{L})复杂度的算法,这份代码用时101秒计算出f(1012)=5435004633092f(10^{12})=5435004633092.

Github