题目700：欧拉币
本帖最后由 永恒的蓝色梦想 于 2020-4-24 08:17 编辑https://s1.ax1x.com/2020/04/24/J0NebV.png
https://s1.ax1x.com/2020/04/24/J0NnET.png 没人吗？ 冰河星云 发表于 2020-5-6 13:49
没人吗？
这么难的题…… 本帖最后由 guosl 于 2021-3-4 21:21 编辑
答案：1517926517777556（用时：4.18425秒）
要用高精度数值计算和一次同余方程解的理论来求。如果要加快计算速度还可以用平行编程来实现。
#include <iostream>
#include <stack>
#include <vector>
#include <algorithm>
#include <gmp.h>
#include <gmpxx.h>
#include <omp.h>
using namespace std;
const int nN = 50000000;
const mpz_class MOD = mpz_class("4503599627370517");
const mpz_class uEu = mpz_class("1504170715041707");
struct stNode
{
unsigned long long v;
unsigned long long n;
bool operator<(const stNode &nd) const
{
return (n < nd.n);
}
};
stack<stNode> sk;
vector<stNode> nd;
mpz_class qpow(mpz_class a, mpz_class p)
{
if (p == 1)
return a;
mpz_class r;
if (p % 2 == 0)
{
r = qpow(a, p / 2);
r *= r;
}
else
{
r = qpow(a, p - 1);
r *= a;
}
r %= MOD;
return r;
}
int main(void)
{
double t = omp_get_wtime();
nd.resize(nN);
nd.v = uEu.get_ui();
nd.n = 1;
sk.push(nd);
unsigned long long pp;
int kk;
mpz_class u = qpow(uEu, MOD - 2);
#pragma omp parallel shared(nd,pp,kk,u)
{
#pragma omp for
for (int i = 0; i < nN; ++i)
{
mpz_t integ;
mpz_init(integ);
mpz_mul_si(integ, uEu.get_mpz_t(), i + 1);
mpz_mod(integ, integ, MOD.get_mpz_t());
nd.v = mpz_get_ui(integ);
mpz_clear(integ);
nd.n = (i + 1);
}
#pragma omp single
{
for (int i = 0; i < nN; ++i)
{
stNode &d = sk.top();
if (nd.v < d.v)
sk.push(nd);
}
pp = nd.n;
kk = sk.top().v;
nd.clear();
nd.resize(kk);
}
#pragma omp for
for (int i = 1; i < kk; ++i)
{
mpz_t integ;
mpz_init(integ);
mpz_mul_si(integ, u.get_mpz_t(), i);
mpz_mod(integ, integ, MOD.get_mpz_t());
nd.n = mpz_get_ui(integ);
mpz_clear(integ);
nd.v = i;
}
}
sort(nd.begin(), nd.end());
for (int i = 0; i < (int)nd.size(); ++i)
{
if (nd.n > pp && nd.v < sk.top().v)
sk.push(nd);
}
unsigned long long Sum = 0;
while (!sk.empty())
{
Sum += sk.top().v;
sk.pop();
}
t = omp_get_wtime() - t;
cout << Sum << endl << t << endl;
return 0;
}
页:
[1]