推进 project euler 510 的上限！

282842712474

问题的原始来源是
https://projecteuler.net/problem=510

最终，这个问题由@wayne 找出了通解：
$$\begin{aligned}
r_a =& k p^2 (p+q)^2\\
r_b =& k q^2 (p+q)^2\\
r_c =& k p^2 q^2
\end{aligned}$$
推导过程就不发了，看上去是几何题，本质上是数论题，其实过程也不繁琐，有兴趣的朋友请教@wayne 吧（推卸责任 ）

问题是要求：
$$S(n)=\sum_{0 < r_c \leq r_b \leq r_a \leq n} (r_a + r_b + r_c)$$
原题要求的$n$是$10^9$。可是我用python写了个脚本，测试能够给出正确结果，并且用时不到1s，因此，这个上限有很大提升空间。

现在的擂台是：找出你能计算的最大的$S(n)$，当然，在合理的时间内。至于这个“合理”嘛，见仁见智吧，暂定为1分钟吧。当然，具体情况具体分析。

282842712474

我暂时算到了$10^{15}$

\begin{array}{c|cc}
n & \text{时间(秒)} & S(n) \\
\hline
  9 & 0.0759450905081973 &315306518862563689 \\
10 & 0.257001835605271 &31530956879912196121 \\
11 & 0.8157740317934493 &3153105325203283282791 \\
12 & 2.6701689522824696 &315310838250195945316697 \\
13 & 8.891065307607503 &31531093463422331742627898 \\
14 & 28.55497651464791 &3153109651585438043478518292 \\
15 & 95.10288827428485 &315310974814380507221575440565
\end{array}

我的代码（没什么好注释的，很清晰很简单）：

1. from gmpy2 import *
   
2. import time
   
3. 
   
4. st = time.clock()
   
5. bmax = 10**9
   
6. nmax = int(bmax**0.25)
   
7. 
   
8. s = 0
   
9. for n in range(1, nmax+1):
   
10.   m = 1
    
11.   while m <= n:
    
12.     if gcd(m,n) == 1:
    
13.       k = bmax//(n**2 * (m+n)**2)
    
14.       s = s + k*(k+1)//2 * ((m**2+n**2)*(m+n)**2+m**2 * n**2)
    
15.     m = m+1
    
16. 
    
17. ed = time.clock()
    
18. print('result:',s)
    
19. print('time:',ed-st)

复制代码

282842712474

282842712474 发表于 2015-8-31 17:33
我暂时算到了$10^{15}$

\begin{array}{c|cc}

从计算结果可以看出，貌似这个结果还有个极限
$$\lim_{n\to\infty} \frac{S(n)}{10^{2n}} = 0.3153109...$$

如果没法很好地表示这个极限的话，那么这个擂台也等价于求这个极限的位数。

wayne

写个GMP/C版本的，应该能打破你的记录

282842712474

wayne 发表于 2015-8-31 18:15
写个GMP/C版本的，应该能打破你的记录

这应该只是常数级别优化的，优化空间不大。

最精妙的是从算法上优化

另外，我的计算机配置偏低，换台主流的计算机都能打破我的记录了

282842712474

wayne 发表于 2015-8-31 18:15
写个GMP/C版本的，应该能打破你的记录

自己先来个Python下的GMP的版本，运行时间如下：
\begin{array}{c|cc}  
n & \text{时间(秒)} & S(n) \\  
\hline  
14 & 15.843699989331057 & 3153109651585438043478518292 \\  
15 & 55.72865253989227 & 315310974814380507221575440565 \\
16 & 214.97727446812584 & 31531097786817229071240053411627
\end{array}
速度大约有近一倍的提升。

Python的GMP，也就是一开始就引入的gmpy2这个库，跟C/C++的类似~~C/C++版本的，交给大家来完成啦。

代码很简单，仅仅是改了一下：

1. from gmpy2 import *
   
2. import time
   
3. 
   
4. st = time.clock()
   
5. bmax = mpz(10**14)
   
6. nmax = mpz(bmax**0.25)
   
7. 
   
8. s = mpz(0)
   
9. for n in range(1, nmax+1):
   
10.   m = mpz(1)
    
11.   n = mpz(n)
    
12.   while m <= n:
    
13.     if gcd(m,n) == 1:
    
14.       k = bmax//(n**2 * (m+n)**2)
    
15.       s = s + k*(k+1)//2 * ((m**2+n**2)*(m+n)**2+m**2 * n**2)
    
16.     m = m+1
    
17. 
    
18. ed = time.clock()
    
19. print('result:',s)
    
20. print('time:',ed-st)

复制代码

上述代码的瓶颈是for循环，众所周知for循环效率不高，改成用map函数效率更高，有50%以上的速度提升：

1. from gmpy2 import *
   
2. import time
   
3. 
   
4. st = time.clock()
   
5. bmax = mpz(10**15)
   
6. nmax = mpz(bmax**0.25)
   
7. 
   
8. global s
   
9. s = mpz(0)
   
10. 
    
11. def jisuan(n):
    
12.   global s
    
13.   m = mpz(1)
    
14.   n = mpz(n)
    
15.   while m <= n:
    
16.     if gcd(m,n) == 1:
    
17.       k = bmax//(n**2 * (m+n)**2)
    
18.       s = s + k*(k+1)//2 * ((m**2+n**2)*(m+n)**2+m**2 * n**2)
    
19.     m = m+1
    
20. 
    
21. list(map(jisuan, range(1, nmax+1)))
    
22. 
    
23. ed = time.clock()
    
24. print('result:',s)
    
25. print('time:',ed-st)

复制代码

\begin{array}{c|cc}  
n & \text{时间(秒)} & S(n) \\  
\hline  
14 & 13.480884810981898 & 3153109651585438043478518292 \\  
15 & 37.99569617586795 & 315310974814380507221575440565 \\
16 & 118.63811014725877 & 31531097786817229071240053411627
\end{array}

aimisiyou

问题来源网址怎么打不开?

wayne

17:  3153109788339240470789175907861511
18:  315310979139296541791878037384668861
19:  31531097923586107532823734833133478270
20:  3153109792663973490168308890337939863948
21:  315310979276053840222630577715772302609714
22:  31531097927910748493556944524769918311060947
23:  3153109792800731344476163185764035193616778666
\begin{array}{c|cc}  
n & S(n) \\  
\hline  
17 & 3153109788339240470789175907861511 \\  
18  & 315310979139296541791878037384668861 \\
19  & 31531097923586107532823734833133478270 \\
20 & 3153109792663973490168308890337939863948\\
21 & 315310979276053840222630577715772302609714\\
22 & 31531097927910748493556944524769918311060947\\
23 &  3153109792800731344476163185764035193616778666
\end{array}

该代码只适用于int64的范围，即10^18次幂，对比了python，效率只提高了固定的比例，再往下计算，由于for会比较大，就没有计算的必要了，除非结构有本质的改变...

1. #include<gmp.h>
   
2. #include<stdio.h>
   
3. #include<stdlib.h>
   
4. #include<stdint.h>
   
5. #include<math.h>
   
6. #include<sys/time.h>
   
7. #include<unistd.h>
   
8. int64_t gcd(int64_t u, int64_t v){
   
9. while ( v != 0) {
   
10.         int64_t r = u % v;
    
11.         u = v;
    
12.         v = r;
    
13.     }
    
14.     return u;       
    
15. }
    
16. 
    
17. int main(int c,char** v) {
    
18.         if(c>1){
    
19.                 struct  timeval start,end;
    
20.                 gettimeofday(&start,NULL);
    
21.                 mpz_t max_tmp,maxT,ans,tmp;
    
22.                 mpz_inits(max_tmp,maxT,ans,tmp,'\0');
    
23.                 mpz_set_str(max_tmp,v[1],10);
    
24. 
    
25.                 int64_t max = mpz_get_ui(max_tmp);
    
26.                 mpz_root(maxT,max_tmp,4);
    
27.                 int64_t qmax = mpz_get_ui(maxT);
    
28.                 mpz_fdiv_q_ui(max_tmp,max_tmp,4);
    
29.                 mpz_root(maxT,max_tmp,4);
    
30.                 int64_t pmax = mpz_get_ui(maxT);
    
31.                 for(int64_t p=1;p<=pmax;++p)
    
32.                 {
    
33.                         int64_t q=p;
    
34.                         while(q<=qmax){
    
35.                                 if(gcd(p,q)==1){
    
36.                                         uint64_t k=max/q/q/(p+q)/(p+q);
    
37.                                         uint64_t a=p*p*(p+q)*(p+q);
    
38.                                         uint64_t b=q*q*(p+q)*(p+q);
    
39.                                         uint64_t c=p*p*q*q;
    
40.                                         mpz_mul_ui(tmp,tmp,0);
    
41.                                         mpz_add_ui(tmp,tmp,a);
    
42.                                         mpz_add_ui(tmp,tmp,b);
    
43.                                         mpz_add_ui(tmp,tmp,c);
    
44.                                         mpz_mul_ui(tmp,tmp,k);
    
45.                                         mpz_mul_ui(tmp,tmp,k+1);
    
46.                                         mpz_fdiv_q_ui(tmp,tmp,2);
    
47.                                         mpz_add(ans,ans,tmp);
    
48.                                         //gmp_printf("%lld,%lld: %Zd, %Zd\n",p,q,tmp,ans);
    
49.                                 }
    
50.                                 q++;
    
51.                         }
    
52.                 }
    
53.                 gettimeofday(&end,NULL);
    
54.                 double diff = (end.tv_sec-start.tv_sec)+ (end.tv_usec-start.tv_usec)*1.0/1000000;
    
55.                 gmp_printf("time:%fs,   %lld: %Zd\n",diff,max,ans);
    
56.                 mpz_clears(max_tmp,maxT,ans,tmp,'\0');
    
57.         }
    
58.        
    
59. }

复制代码

wayne

对比python+gmp 和C+gmp 代码的效率,结论是：同样的算法，C的执行效率是python的三倍，

1. 9 0.011497s, 0.01612                315306518862563689
   
2. 10 0.034581s, 0.052466                        31530956879912196121
   
3. 11 0.095133s, 0.164893                         3153105325203283282791
   
4. 12 0.234725s, 0.49664                        315310838250195945316697
   
5. 13 0.599286s, 1.577465                        31531093463422331742627898
   
6. 14 1.796593s, 4.937906                        3153109651585438043478518292
   
7. 15 5.764869s, 15.912529                315310974814380507221575440565
   
8. 16 17.976543s, 50.962873                31531097786817229071240053411627
   
9. 17 57.295435s, 159.026102                3153109788339240470789175907861511
   
10. 18 182.901621s, *****                        315310979139296541791878037384668861
    
11. 19 571.448696s, ***                                31531097923586107532823734833133478270
    
12. 20 1829.964781s, ****                        3153109792663973490168308890337939863948
    
13. 21 5759.034041s, *****                        315310979276053840222630577715772302609714
    
14. 22                                                                 31531097927910748493556944524769918311060947
    
15. 23                                                                3153109792800731344476163185764035193616778666
    

复制代码

附带任意精度的版本：

1. #include<gmp.h>
   
2. #include<stdio.h>
   
3. #include<stdlib.h>
   
4. #include<stdint.h>
   
5. #include<math.h>
   
6. #include<sys/time.h>
   
7. #include<unistd.h>
   
8. 
   
9. int main(int cnt, char** v) {
   
10.     if (cnt > 1) {
    
11.         struct  timeval start, end;
    
12.         gettimeofday(&start, NULL);
    
13. 
    
14.         mpz_t max, pmax, qmax, p, q, max_tmp, ans, gcd, tmp, a, b, c, k;
    
15.         mpz_inits(max, pmax, qmax, '\0');
    
16.         mpz_inits(p, q, max_tmp, '\0');
    
17.         mpz_inits(ans, gcd, tmp, '\0');
    
18.         mpz_inits(a, b, c, k, '\0');
    
19. 
    
20.         mpz_set_str(max, v[1], 10);
    
21.         //int64_t max = mpz_get_ui(max_tmp);
    
22.         mpz_root(qmax, max, 4);
    
23.         //int64_t qmax = mpz_get_ui(maxT);
    
24.         mpz_fdiv_q_ui(max_tmp, max, 4);
    
25.         mpz_root(pmax, max_tmp, 4);
    
26.         //int64_t pmax = mpz_get_ui(maxT);
    
27. 
    
28.         mpz_set_str(p, "1", 10);
    
29.         while (mpz_cmp(p, pmax) <= 0)
    
30.         {
    
31.             mpz_set(q, p);
    
32.             while (mpz_cmp(q, qmax) <= 0) {
    
33.                 mpz_gcd(gcd, p, q);
    
34.                 if (mpz_cmp_ui(gcd, 1) == 0) {
    
35.                     // int g = mpz_get_ui(gcd);
    
36.                     // if(g==1){
    
37.                     // uint64_t k=max/q/q/(p+q)/(p+q);
    
38.                     // uint64_t a=p*p*(p+q)*(p+q);
    
39.                     // uint64_t b=q*q*(p+q)*(p+q);
    
40.                     // uint64_t c=p*p*q*q;
    
41.                     mpz_set_str(tmp, "0", 10);
    
42.                     mpz_set_str(max_tmp, "0", 10);
    
43.                     //mpz_set_str(k,"1",10);
    
44.                     mpz_set_str(a, "1", 10);
    
45.                     mpz_set_str(b, "1", 10);
    
46.                     mpz_set_str(c, "1", 10);
    
47.                     mpz_add(max_tmp, p, q);
    
48. 
    
49.                     mpz_fdiv_q(k, max, q);
    
50.                     mpz_fdiv_q(k, k, q);
    
51.                     mpz_fdiv_q(k, k, max_tmp);
    
52.                     mpz_fdiv_q(k, k, max_tmp);
    
53. 
    
54.                     mpz_mul(a, p, p);
    
55.                     mpz_mul(a, a, max_tmp);
    
56.                     mpz_mul(a, a, max_tmp);
    
57. 
    
58.                     mpz_mul(b, q, q);
    
59.                     mpz_mul(b, b, max_tmp);
    
60.                     mpz_mul(b, b, max_tmp);
    
61. 
    
62.                     mpz_mul(c, p, p);
    
63.                     mpz_mul(c, c, q);
    
64.                     mpz_mul(c, c, q);
    
65. 
    
66.                     mpz_add(tmp, tmp, a);
    
67.                     mpz_add(tmp, tmp, b);
    
68.                     mpz_add(tmp, tmp, c);
    
69.                     mpz_mul(tmp, tmp, k);
    
70.                     mpz_add_ui(k, k, 1);
    
71.                     mpz_mul(tmp, tmp, k);
    
72.                     mpz_fdiv_q_ui(tmp, tmp, 2);
    
73.                     mpz_add(ans, ans, tmp);
    
74.                     // }
    
75.                 }
    
76.                 mpz_add_ui(q, q, 1);
    
77.             }
    
78.             mpz_add_ui(p, p, 1);
    
79.         }
    
80.         gettimeofday(&end, NULL);
    
81.         double diff = (end.tv_sec - start.tv_sec) + (end.tv_usec - start.tv_usec) / 1000000.;
    
82.         gmp_printf("time:%fs, %Zd\n", diff, ans);
    
83.         mpz_clears(max, pmax, qmax, p, q, max_tmp, ans, gcd, tmp, a, b, c, k, '\0');
    
84.     }
    
85. 
    
86. }

复制代码

wayne

挂机算了足足一晚上，推到$S(10^23)=3153109792800731344476163185764035193616778666$了，耗时58777.435285s，在楼上已经更新

另外，看了下PE的论坛，有人给出了${S(n)}/{10^{2n}}$的极限值，\[-{{1800\,{\rm Li}_4\left(\frac{1}{2}\right)}\over{\pi^4}}-{{3465\,\log 2\,\zeta\left(3
\right)}\over{2\,\pi^4}}+{{135\,\zeta\left(3\right)}\over{2\,\pi^4}}
-{{75\,\log ^42}\over{\pi^4}}+{{75\,\log ^22}\over{\pi^2}}+{{1305
}\over{64}}\]

$0.3153109792805197234817566743674355895396730753785832114718406388865685408449715027332238093190782468$