projectEuler 454 -- Diophantine reciprocals

cn8888

几何意义: 1/A+1/B=1/C
要是有人能把这个图片用latex实现,那就把这个图片删除了吧,
不过删之前先发短信息告诉我一下

cn8888

几何意义就是三个纵坐标都是整数,然后还满足上面的图形要求

sunwukong

因为
\(\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

可否使用多线程？

wayne

282842712474 发表于 2014-6-23 22:41
可否使用多线程？

这个是求和性质的，所以可以使用多线程提升速度。但也只是常数优化。

倪举鹏

算这个的精确值困难的话，算算估计值呢。就像n!   当然有公式算精确值最好了

kastin

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

有没有已知的数据进行对比？
我用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 ... 2667&fromuid=70
的第一个公式算的。

补充：
算了$f(10^9)=3979600400$，用时196.748994145秒...

补充内容 (2014-6-26 08:38):
我看错了...200000秒是差不多3天了~

282842712474

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的整数型来说已经溢出了）。

1. import time
   
2. start=time.clock()
   
3. import math
   
4. l=10000000 #上限
   
5. s=0 #统计个数
   
6. rl=int(math.sqrt(l))
   
7. rrl=int(math.sqrt(int(math.sqrt(l))))
   
8. 
   
9. #弄出个素数表来
   
10. prime=[i for i in range(1,rl+1)]#定义整数表
    
11. j=2
    
12. while j<=rrl:
    
13.   if prime[j-1] != 0:
    
14.     m=int(rl/j)
    
15.     for i in range(j,m+1):
    
16.       prime[i*j-1]=0
    
17.     j=j+1
    
18.   else:
    
19.     j=j+1
    
20. prime.sort() #从小到大重新排列（让0位于前面）
    
21. z=prime.count(0) #统计0的个数
    
22. prime=prime[z+1:rl]
    
23. #素数表生成完毕
    
24. 
    
25. for m in range(1,int(rl/1.4)):
    
26.   #找出与m互素的数，存在p里
    
27.   p=[i for i in range(1,rl+1)]
    
28.   m2=m
    
29.   for i in prime:
    
30.     if i>m:
    
31.       break
    
32.     if m%i==0:
    
33.       m=m/i
    
34.       for j in range(m2//i+1,int(rl/i)+1): #!!!
    
35.         p[i*j-1]=0
    
36.   p.sort() #从小到大重新排列（让0位于前面）
    
37.   z=p.count(0) #统计0的个数
    
38.   p=p[z+m2:rl]
    
39.   #找出与m互数的完毕。
    
40.   for n in p:
    
41.     s=s+l//((n+m2)*n)
    
42. print(s)
    
43. end=time.clock()
    
44. print(end-start)

复制代码