论坛 › 编程擂台 › 推进 project euler 510 的上限！

282842712474 发表于 2015-8-31 17:29:30

推进 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 吧（推卸责任:lol ）

问题是要求：
$$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 发表于 2015-8-31 17:33:15

我暂时算到了$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}

我的代码（没什么好注释的，很清晰很简单）：
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)

282842712474 发表于 2015-8-31 17:37:21

本帖最后由 282842712474 于 2015-8-31 17:38 编辑

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

\begin{array}{c|cc}


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

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

wayne 发表于 2015-8-31 18:15:08

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

282842712474 发表于 2015-8-31 19:19:20

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

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

最精妙的是从算法上优化

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

282842712474 发表于 2015-8-31 19:51:40

本帖最后由 282842712474 于 2015-8-31 21:33 编辑

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++版本的，交给大家来完成啦。

代码很简单，仅仅是改了一下：
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%以上的速度提升：
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}

aimisiyou 发表于 2015-8-31 20:06:18

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

wayne 发表于 2015-8-31 23:04:42

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会比较大，就没有计算的必要了，除非结构有本质的改变...
#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){
                structtimeval 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,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');
        }
       
}

wayne 发表于 2015-9-1 11:38:16

对比python+gmp 和C+gmp 代码的效率,结论是：同样的算法，C的执行效率是python的三倍，
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



附带任意精度的版本：

#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) {
      structtimeval 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, 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');
    }

}

wayne 发表于 2015-9-2 09:38:33

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

查看完整版本: 推进 project euler 510 的上限！