倒数丢番图方程

wayne于2014年3月转载projectEuler 454问题:
设对于给定的 LL,等式 1X+1Y=1N,(0<X<YL)\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)\left. \frac{x+y}{(x,y)}\middle|(x,y) \right.
不过,在此基础上,可以导出这种一般性的解表达: 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\limits
{n=1}^{\sqrt{L}} \sum\limits_{m=1,(n,m)=1}^{n-1} \left[ \frac{L}{n(n+m)} \right] 只需要只需要O(L)$的时间,我认为这是可以承受的。不过我们也可试着优化一下:

n=1Lm=1,(n,m)=1n1[Ln(n+m)]=n=1Ls=n+1,(n,s)=12n1[Lns]=n=1Ls=n+12n1[Lns][1(n,s)]=n=1Ls=n+12n1[Lns]d(n,s)μ(d)=n=1Ldnn<s2n1,ds[Lns]μ(d)=n=1Ldnμ(d)nd<s2n1d[Lnds]=n=1Ldnμ(d)(ψ([Lnd],2n1d)ψ([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' \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)

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

cn8888给出了这个丢番图方程的几何意义:

倪举鹏表示:
还有种几何关系满足这个等式:平面过一点的三射线中间射线与两边射线夹角都是60度。一条直线截这3射线为三段,这三段长度满足这个等式.

sunwukong认为:
因为
1x+1y=1n\frac1x+\frac1y=\frac1n0<x<yL0\lt x\lt y\leq L
所以 n<x<yLn\lt x\lt y\leq L
x=n+sx=n+sy=n+ty=n+t,则1s<tLn1\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=st1×2=2n2n^2=st\geq1\times 2=2 \Rightarrow n\geq 2
s=n2/tn2Lns=n^2/t\geq \frac{n^2}{L-n}
n2Lns<n<tLn\frac{n^2}{L-n}\leq s\lt n\lt t\leq L-n
1n=1x+1y>1L+1L=2Ln<L2\frac1{n}=\frac1x+\frac1y\gt \frac1L+\frac1L=\frac2L \Rightarrow n\lt \frac{L}{2}
所以
f(L)f(L)
=n=2L12sn2Ln,sn2n11=\sum_{n=2}^{\lfloor \frac{L-1}{2} \rfloor} \sum_{ s\geq\frac{n^2}{L-n} , s| n^2 }^{n-1} 1
=n=2L12t=n+1,tn2Ln1=\sum_{n=2}^{\lfloor \frac{L-1}{2} \rfloor} \sum_{ t=n+1 , t| n^2 }^{L-n} 1
=n=2L12(un2Ln,un2Ln1)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

关键的是
给定 nnLL,求满足
tn2t| n^2n<tLnn\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}]&amp;;
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秒

LwinsG测试后认为Mathematica计算ψ(x,y)\psi(x,y)时所消耗的时间是平凡的,为O(y)\mathcal{O}(y),这会导致复杂度退化为 $ \sum{n=1}^{\sqrt{L}} \sum{d | n} \mathcal{O} \left( \frac{n}{d} \right) = \mathcal{O}(L) $ 所以需要自己重新实现。 282842712474直接[使用LwinsG的公式](https://bbs.emath.ac.cn/forum.php?mod=redirect&goto=findpost&ptid=5423&pid=54248&fromuid=20)用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.

© 2021 Blog Boost Starter