题目描述##
求
[sumlimits_{i=1}^{n} sumlimits_{j=1}^{n} i*j*gcd(i,j) pmod{p}
]
(n<=10^{10}),(p)是质数
题解##
推导很长就省略啦，，
有空补回来
最后推得这个式子：
[sumlimits_{T = 1}^{n} (frac{lfloor frac{n}{T}
floor * (lfloor frac{n}{T}
floor + 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[n];
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[n] = 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;
}