洛谷P3768 简单的数学题 【莫比乌斯反演 + 杜教筛】

题目描述

求
\[\sum\limits_{i=1}^{n} \sum\limits_{j=1}^{n} i*j*gcd(i,j) \pmod{p}\]

\(n<=10^{10}\),\(p\)是质数

题解

推导很长就省略啦，，
有空补回来

最后推得这个式子：
\[\sum\limits_{T = 1}^{n} (\frac{\lfloor \frac{n}{T} \rfloor * (\lfloor \frac{n}{T} \rfloor + 1)}{2})^2 * T^2 * \varphi(T)\]
前边分块，后边杜教筛
杜教筛的\(g(n)\)取\(g(n) = n^2\)

#include<iostream>
#include<cstdio>
#include<cmath>
#include<map>
#include<cstring>
#include<algorithm>
#define LL long long int
#define Redge(u) for (int k = h[u],to; k; k = ed[k].nxt)
#define REP(i,n) for (int i = 1; i <= (n); i++)
#define BUG(s,n) for (int i = 1; i <= (n); i++) cout<<s[i]<<' '; puts("");
using namespace std;
const int maxn = 5000005,maxm = 100005,INF = 1000000000;
typedef map<LL,LL> Map;
Map _f;
LL P,N,v6,v2;
LL p[maxn],pi,phi[maxn],f[maxn];
int isn[maxn];
LL qpow(LL a,LL b){
LL ans = 1;
for (; b; b >>= 1,a = a * a % P)
if (b & 1) ans = ans * a % P;
return ans;
}
void init(LL n){
v6 = qpow(6,P - 2);
v2 = qpow(2,P - 2);
N = (LL)pow(n,2.0 / 3.0);
phi[1] = 1;
for (LL i = 2; i < N; i++){
if (!isn[i]) p[++pi] = i,phi[i] = (i - 1) % P;
for (LL j = 1; j <= pi && i * p[j] < N; j++){
isn[i * p[j]] = true;
if (i % p[j] == 0){
phi[i * p[j]] = phi[i] * p[j] % P;
break;
}
phi[i * p[j]] = phi[i] * (p[j] - 1) % P;
}
}
for (LL i = 1; i < N; i++) f[i] = (f[i - 1] + i * i % P * phi[i] % P) % P;
}
LL sum(LL n){
n %= P;
LL tmp = n * (n + 1) % P * v2 % P;
return tmp * tmp % P;
}
LL sum2(LL n){
n %= P;
return n * (n + 1) % P * (2 * n % P + 1) % P * v6 % P;
}
LL S(LL n){
if (n < N) return f
;
Map::iterator it;
if ((it = _f.find(n)) != _f.end())
return it->second;
LL ans = n % P * ((n + 1) % P) % P * v2 % P;
ans = ans * ans % P;
for (LL i = 2,nxt; i <= n; i = nxt + 1){
nxt = n / (n / i);
ans = (ans - (sum2(nxt) - sum2(i - 1)) % P * S(n / i) % P) % P;
}
ans = (ans + P) % P;
return _f
 = ans;
}
int main(){
LL n,ans = 0;
cin >> P >> n;
init(n);
for (LL i = 1,nxt; i <= n; i = nxt + 1){
nxt = n / (n / i);
ans = (ans + sum(n / i) * ((S(nxt) - S(i - 1)) % P) % P) % P;
}
ans = (ans + P) % P;
cout << ans << endl;
return 0;
}