projectEuler 454 -- Diophantine reciprocals

282842712474

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 ...

$\psi(x,y)$是$O(\sqrt(x))$的，那$\mu(d)$呢，我怎么感觉它的计算量也不低呀

sunwukong

由 13#，

\(n^2=st,s \lt n \lt t\le L-n\)

设 \(s=k*a^2\),其中，\(k\) 无平方素因子，由整数分解定理可知，这个分解是唯一的。

由  \(n^2=k*a^2*t\)

得，\(t=k*b^2,n=k*a*b\)，由 \(s\lt t\)，得 \(a\lt b\)，
由  \(t\le L-n\)，得  \(k*b*(a+b)\le L\)

所以，原问题转化为

给定 \(L\) ，求满足

\(\mu(k)\not=0\)
\(k*b*(a+b)\le L\)
\(a\lt b\)

的解 \(\left(k,a,b \right)\) 的个数。

因为 \(b*(a+b)>=2*(1+2)=6\)，所以 \(k\le L/6\)
所以
\[
\begin{split}
f(L)=\sum_{\begin{subarray}{cc} k=1 \\ \mu(k)\not=0 \end{subarray}}^{L/6} \sum_{b=2}^{\sqrt{L/k}} \min\left( b-1,\lfloor \frac{L}{kb} \rfloor-b \right)
\end{split}
\]

注： \(\mu(x)\) 是莫比乌斯函数 http://zh.wikipedia.org/wiki/默比乌斯函数

282842712474

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 ...

刚刚写完了$O(L^{3/4})$的公式的代码，运行测试$f(10^9)$的时间为24s左右（由上次的190s到了24s）。
但是$f(10^{10})$运行了半个钟都没有结果，而运行$f(2\times 10^9)$需要75s
我觉得由于$d|n$的存在，需要分解n，两个公式其实都不止$O(N)$

$O(L^{3/4})$算法的代码如下（还是python)

1. import time
   
2. import os
   
3. start=time.clock()
   
4. import math
   
5. 
   
6. #主要计算函数（最后一个求和号）
   
7. def psi(l, n, d):
   
8.   sum=0
   
9.   for s in xrange(n//d+1, (2*n-1)//d + 1):
   
10.     sum = sum + l//(n*d*s)
    
11.   return sum
    
12. 
    
13. #生成一个素数表的莫比乌斯表
    
14. def mobi(p):
    
15.   length = len(p)
    
16.   if length == 1:
    
17.     mu = [1,-p[0]]
    
18.   else:
    
19.     mu2=mobi(p[0:length-1])
    
20.     mu=mu2[:]
    
21.     for i in mu2:
    
22.       mu.append(-1*p[length-1]*i)
    
23.   return mu
    
24. 
    
25. #开始
    
26. l=2000000000 #上限
    
27. s=0 #统计个数
    
28. rl=int(math.sqrt(l))
    
29. rrl=int(math.sqrt(int(math.sqrt(l))))
    
30. rrrl=int(math.sqrt(int(math.sqrt(int(math.sqrt(l))))))
    
31. 
    
32. #生成rrl内的素数表
    
33. prime=[i for i in xrange(1,rrl+1)]#定义整数表
    
34. j=2
    
35. while j<=rrrl:
    
36.   if prime[j-1] != 0:
    
37.     m=int(rrl/j)
    
38.     for i in xrange(j,m+1):
    
39.       prime[i*j-1]=0
    
40.     j=j+1
    
41.   else:
    
42.     j=j+1
    
43. prime.sort() #从小到大重新排列（让0位于前面）
    
44. z=prime.count(0) #统计0的个数
    
45. prime=prime[z+1:rl]
    
46. #素数表生成完毕
    
47. 
    
48. prime.append(rl) #往素数表里边添加一个大数，避免溢出错误
    
49. 
    
50. for n in xrange(2,rl+1):
    
51.   p=[] #n的素因子表
    
52.   mu=[] #莫比乌斯表
    
53.   nn=n
    
54.   i=0
    
55.   rn=int(math.sqrt(n))
    
56.   #开始生成n的素因子表
    
57.   while prime[i] <= rn:
    
58.     if n%prime[i]==0:
    
59.       n=n//prime[i]
    
60.       if n%prime[i] != 0:
    
61.         p.append(prime[i])
    
62.         i=i+1
    
63.     else:i=i+1
    
64.   #下面这句的意思是，如果n的一个因子不能被前面的素数整除，那么这个因子就是素数。
    
65.   if n>rn:
    
66.     p.append(n)
    
67.   #n的素因子表p生成完毕
    
68.   mu=mobi(p) #生成莫比乌斯表
    
69.   for d in mu:
    
70.       s=s+d//abs(d)*psi(l,nn,abs(d))
    
71. 
    
72. print(s)
    
73. end=time.clock()
    
74. print(end-start)
    

复制代码

282842712474

我发觉主要的问题是，我想不出$psi(x,y)$的$O(\sqrt{x})$算法...求指教。

282842712474

终于做出来了。$f(10^{12})=5435004633092$，用时101秒！
稍后会整理代码出来。

=========================

http://spaces.ac.cn/index.php/archives/2665/

1. import time
   
2. start=time.clock()
   
3. import math
   
4. 
   
5. #主要计算函数（最后一个求和号）
   
6. def psi(x, y1, y2):
   
7.   s=0
   
8.   rx=int(math.sqrt(x))
   
9.   if y2<= rx:
   
10.     for d in xrange(y1, y2+1):
    
11.       s=s+x//d
    
12.   
    
13.   elif y1<=rx:
    
14.     for d in xrange(y1, rx+1):
    
15.       s=s+x//d
    
16.     p=x//rx-1
    
17.     q=x//y2
    
18.     d=0
    
19.     q1=x//(p-d+1)
    
20.     while d < p-q:
    
21.       q2=q1
    
22.       q1=x//(p-d)
    
23.       s=s+(q1-q2)*(p-d)
    
24.       d=d+1
    
25.     s=s+(y2-x//(q+1))*q
    
26.   
    
27.   else:
    
28.     p=x//y1-1
    
29.     q=x//y2
    
30.     if p<q:
    
31.       s=(1+y2-y1)*q
    
32.     else:
    
33.       d=0
    
34.       q1=x//(p-d+1)
    
35.       while d < p-q:
    
36.         q2=q1
    
37.         q1=x//(p-d)
    
38.         s=s+(q1-q2)*(p-d)
    
39.         d=d+1
    
40.       s=s+(y2-x//(q+1))*q
    
41.       s=s+(x//(p+1)-y1+1)*(p+1)
    
42.   return s
    
43. 
    
44. #生成一个素数表的莫比乌斯表
    
45. def mobi(p):
    
46.   length = len(p)
    
47.   mu=[1,-p[0]]
    
48.   for i in xrange(1,length):
    
49.     mu2=mu[:]
    
50.     for j in mu2:
    
51.       mu.append(-1*j*p[i])
    
52.   return mu
    
53. 
    
54. #开始
    
55. l=1000000000000 #上限
    
56. s=0 #统计个数
    
57. rl=int(math.sqrt(l))
    
58. rrl=int(math.sqrt(int(math.sqrt(l))))
    
59. 
    
60. #生成rl内的素数表
    
61. prime=[i for i in range(1,rl+1)]#定义整数表
    
62. j=2
    
63. while j<=rrl:
    
64.   if prime[j-1] != 0:
    
65.     m=int(rl/j)
    
66.     for i in range(j,m+1):
    
67.       prime[i*j-1]=0
    
68.     j=j+1
    
69.   else:
    
70.     j=j+1
    
71. prime.sort() #从小到大重新排列（让0位于前面）
    
72. z=prime.count(0) #统计0的个数
    
73. prime=prime[z+1:rl]
    
74. #素数表生成完毕
    
75. 
    
76. prime.append(rl) #往素数表里边添加一个大数，避免溢出错误
    
77. 
    
78. for n in xrange(2,rl+1):
    
79.   p=[] #n的素因子表
    
80.   mu=[] #莫比乌斯表
    
81.   nn=n
    
82.   i=0
    
83.   rn=int(math.sqrt(n))
    
84.   #开始生成n的素因子表
    
85.   while prime[i] <= rn:
    
86.     if n%prime[i]==0:
    
87.       n=n//prime[i]
    
88.       if n%prime[i] != 0:
    
89.         p.append(prime[i])
    
90.         i=i+1
    
91.     else:i=i+1
    
92.   #下面这句的意思是，如果n的一个因子不能被前面的素数整除，那么这个因子就是素数。
    
93.   if n>rn:
    
94.     p.append(n)
    
95.   #n的素因子表p生成完毕
    
96.   mu=mobi(p) #生成莫比乌斯表
    
97.   for d in mu:
    
98.     ll=l//(nn*abs(d))
    
99.     s=s+d//abs(d)*(psi(ll,nn//abs(d)+1,(2*nn-1)//abs(d)))
    
100. 
     
101. print(s)
     
102. end=time.clock()
     
103. print(end-start)
     
104. 
     

复制代码