几何意义: 1/A+1/B=1/C
要是有人能把这个图片用latex实现,那就把这个图片删除了吧,
不过删之前先发短信息告诉我一下
几何意义就是三个纵坐标都是整数,然后还满足上面的图形要求
因为
\(\frac1x+\frac1y=\frac1n\)，\(0\lt x\lt y\leq L\)
所以 \(n\lt x\lt y\leq L\)
设 \(x=n+s\)，\(y=n+t\)，则\(1\leq s\lt t\leq L-n\)
\(\frac1{n+s}+\frac1{n+t}=\frac1{n} \iff n^2=st\)
所以
\(n^2=st\geq1\times 2=2 \Rightarrow n\geq 2\)
\(s=n^2/t\geq \frac{n^2}{L-n}\)
\(\frac{n^2}{L-n}\leq s\lt n\lt t\leq L-n\)
\(\frac1{n}=\frac1x+\frac1y\gt \frac1L+\frac1L=\frac2L \Rightarrow n\lt \frac{L}{2}\)
所以
\(f(L)\)
\(=\sum_{n=2}^{\lfloor \frac{L-1}{2} \rfloor} \sum_{\begin{subarray}{cc} s\geq\frac{n^2}{L-n} \\ s\,|\, n^2 \end{subarray}}^{n-1} 1\)
\(=\sum_{n=2}^{\lfloor \frac{L-1}{2} \rfloor} \sum_{\begin{subarray}{cc} t=n+1 \\ t\,|\, n^2 \end{subarray}}^{L-n} 1\)
\(=\sum_{n=2}^{\lfloor \frac{L-1}{2} \rfloor} ((\sum_{\begin{subarray}{cc} u\geq\frac{n^2}{L-n} \\ u\,|\, n^2 \end{subarray}}^{L-n} 1)-1)/2\)
关键的是
给定 \(n\)，\(L\)，求满足
\(t\,|\, n^2\)，\(n\lt t\leq L-n\)
的\(t\)的数目
这一步的复杂度与\(n^2\)的质因数分解相同
可否使用多线程？
282842712474 发表于 2014-6-23 22:41
可否使用多线程？
这个是求和性质的，所以可以使用多线程提升速度。但也只是常数优化。
算这个的精确值困难的话，算算估计值呢。就像n! 当然有公式算精确值最好了
sunwukong 发表于 2014-5-23 13:17
因为
\(\frac1x+\frac1y=\frac1n\)，\(0\lt x\lt y\leq L\)
所以 \(n\lt x\lt y\leq L\)
其实就是涉及大数分解，虽然大数分解本身没有什么高效率的算法，但是这个题目只是统计解的个数，而非找到每一组解。
联想到给出不大于`x`的素数个数`\pi(x)`的算法目前最高效的算法中，并不需要找出所有小于x的素数然后统计个数（见M. DELEGLISE AND J. RIVAT.COMPUTING (x): THE MEISSEL, LEHMER, LAGARIAS,MILLER, ODLYZKO METHOD)
因而，`n^2`满足一定限制要求的分解方案数，应该有相应不需要找出每种具体方案的高效算法。
还有种几何关系满足这个等式：平面过一点的三射线中间射线与两边射线夹角都是60度。一条直线截这3射线为三段，这三段长度满足这个等式
本帖最后由 282842712474 于 2014-6-25 20:02 编辑
有没有已知的数据进行对比？
我用python算得$f(10^7)=30093331,\quad f(10^8)=349446716$，前者用时1.7秒，后者用时17秒。
照$O(N)$复杂度算，$f(10^12)$需要200000秒，也就是六七个小时左右（实际可能更多）。使用多个线程，可以把时间减半吧~
我是直接用
http://bbs.emath.ac.cn/forum.php?mod=redirect&goto=findpost&ptid=5423&pid=52667&fromuid=70
的第一个公式算的。
补充：
算了$f(10^9)=3979600400$，用时196.748994145秒...
补充内容 (2014-6-26 08:38):
我看错了...200000秒是差不多3天了~
本帖最后由 282842712474 于 2014-6-25 20:06 编辑
282842712474 发表于 2014-6-25 19:45
有没有已知的数据进行对比？
我用python算得$f(10^7)=30093331,\quad f(10^8)=349446716$，前者用时1.7秒 ...
下面的代码在python2.7和python3均可以运行，17秒的结果是在pypy下运行得到
由于python是号称运行速度最慢的编程语言，如果转为C++，应该有优化空间。但是我不知道C++上面如何处理大整数（$10^12$对于C的整数型来说已经溢出了）。
import time
start=time.clock()
import math
l=10000000 #上限
s=0 #统计个数
rl=int(math.sqrt(l))
rrl=int(math.sqrt(int(math.sqrt(l))))
#弄出个素数表来
prime=#定义整数表
j=2
while j<=rrl:
if prime != 0:
m=int(rl/j)
for i in range(j,m+1):
prime=0
j=j+1
else:
j=j+1
prime.sort() #从小到大重新排列（让0位于前面）
z=prime.count(0) #统计0的个数
prime=prime
#素数表生成完毕
for m in range(1,int(rl/1.4)):
#找出与m互素的数，存在p里
p=
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=0
p.sort() #从小到大重新排列（让0位于前面）
z=p.count(0) #统计0的个数
p=p
#找出与m互数的完毕。
for n in p:
s=s+l//((n+m2)*n)
print(s)
end=time.clock()
print(end-start)