Decription
给定nn,求
Solution
∑i=1n∑j=1nijgcd(i,j)∑i=1n∑j=1nijgcd(i,j)
=∑d3∑i=1⌊nd⌋∑j=1⌊nd⌋ij[gcd(i,j=1)]=∑d3∑i=1⌊nd⌋∑j=1⌊nd⌋ij[gcd(i,j=1)]
设S(n)=∑ni=1∑nj=1ij[gcd(i,j=1)]S(n)=∑i=1n∑j=1nij[gcd(i,j=1)]
则原式
=∑d3G(⌊nd⌋)=∑d3G(⌊nd⌋)
S(n)=∑i=1n∑j=1nij[gcd(i,j=1)]S(n)=∑i=1n∑j=1nij[gcd(i,j=1)]
=∑i=1ni∑j=1nj[gcd(i,j=1)]=∑i=1ni∑j=1nj[gcd(i,j=1)]
=2∑i=1ni∑j=1ij[gcd(i,j=1)]−1=2∑i=1ni∑j=1ij[gcd(i,j=1)]−1
=2∑i=1niiϕ(i)+[i=1]2−1=2∑i=1niiϕ(i)+[i=1]2−1
=∑i=1ni2iϕ(i)=∑i=1ni2iϕ(i)
设f(i)=i2ϕ(i),g(i)=i2f(i)=i2ϕ(i),g(i)=i2,则S(n)=∑ni=1f(i)S(n)=∑i=1nf(i)
最后用杜教筛求出S(n)S(n)
S(n)=∑i=1n(f∗g)(i)−∑i=2ng(i)S(i)=∑i=1ni3−∑i=2ng(i)S(i)=(n∗(n+1)2)2−∑i=2ng(i)S(i)S(n)=∑i=1n(f∗g)(i)−∑i=2ng(i)S(i)=∑i=1ni3−∑i=2ng(i)S(i)=(n∗(n+1)2)2−∑i=2ng(i)S(i)
#include <bits/stdc++.h>
using namespace std;
const int maxn = 5000005, Max = 5000000;
typedef long long lint;
int p, inv1, inv2;
lint n;
inline int pow(int x, int k)
{
int ret = 1;
while (k) {
if (k & 1) ret = (lint)ret * x % p;
x = (lint)x * x % p; k >>= 1;
}
return ret;
}
int s[maxn], phi[maxn], pr[maxn], pcnt;
bool isp[maxn];
void prepare(int n)
{
memset(isp, -1, sizeof(isp));
isp[1] = 0; phi[1] = 1;
for (int i = 2; i <= n; ++i) {
if (isp[i]) {phi[i] = i - 1; pr[++pcnt] = i;}
for (int j = 1; i * pr[j] <= n && j <= pcnt; ++j) {
isp[i * pr[j]] = 0;
if (i % pr[j]) phi[i * pr[j]] = (lint)phi[i] * phi[pr[j]] % p;
else phi[i * pr[j]] = (lint)phi[i] * pr[j] % p;
}
}
for (int i = 1; i <= n; ++i) {
s[i] = (lint)i * i % p * phi[i] % p;
s[i] += s[i - 1];
if (s[i] >= p) s[i] -= p;
}
}
bool vis[maxn];
lint S[maxn];
inline lint calc1(lint x)
{
x %= p;
return (lint)x * (x + 1) % p * (2 * x + 1) % p * inv1 % p;
}
inline lint calc2(lint x)
{
x %= p;
x = x * (x + 1) % p * inv2 % p;
return x * x % p;
}
lint dfs(lint x)
{
if (x < Max) return s[x];
int m = n / x;
if (vis[m]) return S[m];
vis[m] = true;
lint sum = calc2(x);
for (lint i = 2, j; i <= x; i = j + 1) {
j = x / (x / i);
(sum -= (calc1(j) - calc1(i - 1)) * dfs(x / j) % p) %= p;
}
if (sum < 0) sum += p;
return S[m] = sum;
}
int main()
{
scanf("%d%lld", &p, &n); inv1 = pow(6, p - 2); inv2 = pow(2, p - 2);
prepare(Max);
dfs(n);
lint ans = 0;
for (lint i = 1, j; i <= n; i = j + 1) {
j = n / (n / i);
if (n / i <= Max) ans += (lint)(calc2(j) - calc2(i - 1) + p) % p * s[n / i] % p;
else ans += (lint)(calc2(j) - calc2(i - 1) + p) % p * S[i] % p;
if (ans >= p) ans -= p;
}
printf("%lld\n", ans);
return 0;
}