推进 project euler 510 的上限！

公子相

问题的原始来源是
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分钟吧。当然，具体情况具体分析。

azanuvahejih

我暂时算到了$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 * import time st = time.clock() bmax = 10**9 nmax = int(bmax**0.25) s = 0 for n in range(1, nmax+1):   m = 1   while m <= n:     if gcd(m,n) == 1:       k = bmax//(n**2 * (m+n)**2)       s = s + k*(k+1)//2 * ((m**2+n**2)*(m+n)**2+m**2 * n**2)     m = m+1 ed = time.clock() print('result:',s) print('time:',ed-st)

复制代码

ipeuric

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

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

Marleslots

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

油条大哥大

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

最精妙的是从算法上优化

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

ojaupobcisi

自己先来个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 * import time st = time.clock() bmax = mpz(10**14) nmax = mpz(bmax**0.25) s = mpz(0) for n in range(1, nmax+1):   m = mpz(1)   n = mpz(n)   while m <= n:     if gcd(m,n) == 1:       k = bmax//(n**2 * (m+n)**2)       s = s + k*(k+1)//2 * ((m**2+n**2)*(m+n)**2+m**2 * n**2)     m = m+1 ed = time.clock() print('result:',s) print('time:',ed-st)

复制代码

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

1. from gmpy2 import * import time st = time.clock() bmax = mpz(10**15) nmax = mpz(bmax**0.25) global s s = mpz(0) def jisuan(n):   global s   m = mpz(1)   n = mpz(n)   while m <= n:     if gcd(m,n) == 1:       k = bmax//(n**2 * (m+n)**2)       s = s + k*(k+1)//2 * ((m**2+n**2)*(m+n)**2+m**2 * n**2)     m = m+1 list(map(jisuan, range(1, nmax+1))) ed = time.clock() print('result:',s) 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}

isotaztiqotew

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

urinowiq

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> #include<stdio.h> #include<stdlib.h> #include<stdint.h> #include<math.h> #include<sys/time.h> #include<unistd.h> int64_t gcd(int64_t u, int64_t v){ while ( v != 0) {         int64_t r = u % v;         u = v;         v = r;     }     return u;         } int main(int c,char** v) {         if(c>1){                 struct  timeval start,end;                 gettimeofday(&start,NULL);                 mpz_t max_tmp,maxT,ans,tmp;                 mpz_inits(max_tmp,maxT,ans,tmp,'\0');                 mpz_set_str(max_tmp,v[1],10);                 int64_t max = mpz_get_ui(max_tmp);                 mpz_root(maxT,max_tmp,4);                 int64_t qmax = mpz_get_ui(maxT);                 mpz_fdiv_q_ui(max_tmp,max_tmp,4);                 mpz_root(maxT,max_tmp,4);                 int64_t pmax = mpz_get_ui(maxT);                 for(int64_t p=1;p<=pmax;++p)                 {                         int64_t q=p;                         while(q<=qmax){                                 if(gcd(p,q)==1){                                         uint64_t k=max/q/q/(p+q)/(p+q);                                         uint64_t a=p*p*(p+q)*(p+q);                                         uint64_t b=q*q*(p+q)*(p+q);                                         uint64_t c=p*p*q*q;                                         mpz_mul_ui(tmp,tmp,0);                                         mpz_add_ui(tmp,tmp,a);                                         mpz_add_ui(tmp,tmp,b);                                         mpz_add_ui(tmp,tmp,c);                                         mpz_mul_ui(tmp,tmp,k);                                         mpz_mul_ui(tmp,tmp,k+1);                                         mpz_fdiv_q_ui(tmp,tmp,2);                                         mpz_add(ans,ans,tmp);                                         //gmp_printf("%lld,%lld: %Zd, %Zd\n",p,q,tmp,ans);                                 }                                 q++;                         }                 }                 gettimeofday(&end,NULL);                 double diff = (end.tv_sec-start.tv_sec)+ (end.tv_usec-start.tv_usec)*1.0/1000000;                 gmp_printf("time:%fs,   %lld: %Zd\n",diff,max,ans);                 mpz_clears(max_tmp,maxT,ans,tmp,'\0');         }          }

复制代码

Zpasarbola

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

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

复制代码

附带任意精度的版本：

1. #include<gmp.h> #include<stdio.h> #include<stdlib.h> #include<stdint.h> #include<math.h> #include<sys/time.h> #include<unistd.h> int main(int cnt, char** v) {     if (cnt > 1) {         struct  timeval start, end;         gettimeofday(&start, NULL);         mpz_t max, pmax, qmax, p, q, max_tmp, ans, gcd, tmp, a, b, c, k;         mpz_inits(max, pmax, qmax, '\0');         mpz_inits(p, q, max_tmp, '\0');         mpz_inits(ans, gcd, tmp, '\0');         mpz_inits(a, b, c, k, '\0');         mpz_set_str(max, v[1], 10);         //int64_t max = mpz_get_ui(max_tmp);         mpz_root(qmax, max, 4);         //int64_t qmax = mpz_get_ui(maxT);         mpz_fdiv_q_ui(max_tmp, max, 4);         mpz_root(pmax, max_tmp, 4);         //int64_t pmax = mpz_get_ui(maxT);         mpz_set_str(p, "1", 10);         while (mpz_cmp(p, pmax) <= 0)         {             mpz_set(q, p);             while (mpz_cmp(q, qmax) <= 0) {                 mpz_gcd(gcd, p, q);                 if (mpz_cmp_ui(gcd, 1) == 0) {                     // int g = mpz_get_ui(gcd);                     // if(g==1){                     // uint64_t k=max/q/q/(p+q)/(p+q);                     // uint64_t a=p*p*(p+q)*(p+q);                     // uint64_t b=q*q*(p+q)*(p+q);                     // uint64_t c=p*p*q*q;                     mpz_set_str(tmp, "0", 10);                     mpz_set_str(max_tmp, "0", 10);                     //mpz_set_str(k,"1",10);                     mpz_set_str(a, "1", 10);                     mpz_set_str(b, "1", 10);                     mpz_set_str(c, "1", 10);                     mpz_add(max_tmp, p, q);                     mpz_fdiv_q(k, max, q);                     mpz_fdiv_q(k, k, q);                     mpz_fdiv_q(k, k, max_tmp);                     mpz_fdiv_q(k, k, max_tmp);                     mpz_mul(a, p, p);                     mpz_mul(a, a, max_tmp);                     mpz_mul(a, a, max_tmp);                     mpz_mul(b, q, q);                     mpz_mul(b, b, max_tmp);                     mpz_mul(b, b, max_tmp);                     mpz_mul(c, p, p);                     mpz_mul(c, c, q);                     mpz_mul(c, c, q);                     mpz_add(tmp, tmp, a);                     mpz_add(tmp, tmp, b);                     mpz_add(tmp, tmp, c);                     mpz_mul(tmp, tmp, k);                     mpz_add_ui(k, k, 1);                     mpz_mul(tmp, tmp, k);                     mpz_fdiv_q_ui(tmp, tmp, 2);                     mpz_add(ans, ans, tmp);                     // }                 }                 mpz_add_ui(q, q, 1);             }             mpz_add_ui(p, p, 1);         }         gettimeofday(&end, NULL);         double diff = (end.tv_sec - start.tv_sec) + (end.tv_usec - start.tv_usec) / 1000000.;         gmp_printf("time:%fs, %Zd\n", diff, ans);         mpz_clears(max, pmax, qmax, p, q, max_tmp, ans, gcd, tmp, a, b, c, k, '\0');     } }

复制代码

asomiajht

挂机算了足足一晚上，推到$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$