题目700：欧拉币

永恒的蓝色梦想

                               
登录/注册后可看大图

                               
登录/注册后可看大图

冰河星云

没人吗？

永恒的蓝色梦想

冰河星云 发表于 2020-5-6 13:49
没人吗？

这么难的题……

guosl

答案：1517926517777556（用时：4.18425秒）
要用高精度数值计算和一次同余方程解的理论来求。如果要加快计算速度还可以用平行编程来实现。

1. #include <iostream>
   
2. #include <stack>
   
3. #include <vector>
   
4. #include <algorithm>
   
5. #include <gmp.h>
   
6. #include <gmpxx.h>
   
7. #include <omp.h>
   
8. using namespace std;
   
9. 
   
10. const int nN = 50000000;
    
11. const mpz_class MOD = mpz_class("4503599627370517");
    
12. const mpz_class uEu = mpz_class("1504170715041707");
    
13. 
    
14. struct stNode
    
15. {
    
16.   unsigned long long v;
    
17.   unsigned long long n;
    
18.   bool operator<(const stNode &nd) const
    
19.   {
    
20.     return (n < nd.n);
    
21.   }
    
22. };
    
23. 
    
24. stack<stNode> sk;
    
25. vector<stNode> nd;
    
26. 
    
27. mpz_class qpow(mpz_class a, mpz_class p)
    
28. {
    
29.   if (p == 1)
    
30.     return a;
    
31.   mpz_class r;
    
32.   if (p % 2 == 0)
    
33.   {
    
34.     r = qpow(a, p / 2);
    
35.     r *= r;
    
36.   }
    
37.   else
    
38.   {
    
39.     r = qpow(a, p - 1);
    
40.     r *= a;
    
41.   }
    
42.   r %= MOD;
    
43.   return r;
    
44. }
    
45. 
    
46. int main(void)
    
47. {
    
48.   double t = omp_get_wtime();
    
49.   nd.resize(nN);
    
50.   nd[0].v = uEu.get_ui();
    
51.   nd[0].n = 1;
    
52.   sk.push(nd[0]);
    
53.   unsigned long long pp;
    
54.   int kk;
    
55.   mpz_class u = qpow(uEu, MOD - 2);
    
56. #pragma omp parallel shared(nd,pp,kk,u)
    
57.   {
    
58. #pragma omp for
    
59.     for (int i = 0; i < nN; ++i)
    
60.     {
    
61.       mpz_t integ;
    
62.       mpz_init(integ);
    
63.       mpz_mul_si(integ, uEu.get_mpz_t(), i + 1);
    
64.       mpz_mod(integ, integ, MOD.get_mpz_t());
    
65.       nd[i].v = mpz_get_ui(integ);
    
66.       mpz_clear(integ);
    
67.       nd[i].n = (i + 1);
    
68.     }
    
69. #pragma omp single
    
70.     {
    
71.       for (int i = 0; i < nN; ++i)
    
72.       {
    
73.         stNode &d = sk.top();
    
74.         if (nd[i].v < d.v)
    
75.           sk.push(nd[i]);
    
76.       }
    
77.       pp = nd[nN - 1].n;
    
78.       kk = sk.top().v;
    
79.       nd.clear();
    
80.       nd.resize(kk);
    
81.     }
    
82. #pragma omp for
    
83.     for (int i = 1; i < kk; ++i)
    
84.     {
    
85.       mpz_t integ;
    
86.       mpz_init(integ);
    
87.       mpz_mul_si(integ, u.get_mpz_t(), i);
    
88.       mpz_mod(integ, integ, MOD.get_mpz_t());
    
89.       nd[i - 1].n = mpz_get_ui(integ);
    
90.       mpz_clear(integ);
    
91.       nd[i - 1].v = i;
    
92.     }
    
93.   }
    
94.   sort(nd.begin(), nd.end());
    
95.   for (int i = 0; i < (int)nd.size(); ++i)
    
96.   {
    
97.     if (nd[i].n > pp && nd[i].v < sk.top().v)
    
98.       sk.push(nd[i]);
    
99.   }
    
100.   unsigned long long Sum = 0;
     
101.   while (!sk.empty())
     
102.   {
     
103.     Sum += sk.top().v;
     
104.     sk.pop();
     
105.   }
     
106.   t = omp_get_wtime() - t;
     
107.   cout << Sum << endl << t << endl;
     
108.   return 0;
     
109. }
     

复制代码