题目大意
求 ∑ i = 1 n ∑ j = 1 n i ∗ j ∗ g c d ( i , j ) \sum_{i=1}^n\sum_{j=1}^ni*j*gcd(i,j) ∑i=1n∑j=1ni∗j∗gcd(i,j)
莫比乌斯反演
在莫比乌斯反演那里学到了这个式子: n = ∑ d ∣ n φ ( d ) n=\sum_{d|n}\varphi(d) n=∑d∣nφ(d)
那么将其代入上式可以得到: ∑ i = 1 n ∑ j = 1 n i ∗ j ∑ d ∣ g c d ( i , j ) φ ( d ) \sum_{i=1}^n\sum_{j=1}^ni*j\sum_{d|gcd(i,j)}\varphi(d) ∑i=1n∑j=1ni∗j∑d∣gcd(i,j)φ(d)
那么我们可以枚举 d ∈ [ 1 , n ] d \in [1,n] d∈[1,n],然后设 i = x ∗ d , j = y ∗ d i=x*d,j=y*d i=x∗d,j=y∗d,枚举 x , y x,y x,y并替换 i i i,我们可以得到:
∑ d = 1 n φ ( d ) ∗ d 2 ∑ x = 1 ⌊ n d ⌋ x ∑ y = 1 ⌊ n d ⌋ y \sum_{d=1}^n \varphi(d) *d^2\sum_{x=1}^{\lfloor \frac{n}{d}\rfloor} x \sum_{y=1}^{\lfloor \frac{n}{d}\rfloor} y ∑d=1nφ(d)∗d2∑x=1⌊dn⌋x∑y=1⌊dn⌋y
令 g ( x ) = x ( x + 1 ) 2 g(x)=\frac{x(x+1)}{2} g(x)=2x(x+1),那么上式等价于:
∑ d = 1 n φ ( d ) ∗ d 2 [ g ( ⌊ n d ⌋ ) ] 2 \sum_{d=1}^n \varphi(d) *d^2 [g(\lfloor \frac{n}{d}\rfloor)]^2 ∑d=1nφ(d)∗d2[g(⌊dn⌋)]2
如果我们能预处理出 ∑ d = 1 n φ ( d ) ∗ d 2 \sum_{d=1}^n \varphi(d) *d^2 ∑d=1nφ(d)∗d2,但是 n ≤ 1 e 10 n\leq 1e^{10} n≤1e10,显然无法线性筛,此时就需要杜教筛了
杜教筛
对于 ∑ d = 1 n φ ( d ) ∗ d 2 \sum_{d=1}^n \varphi(d) *d^2 ∑d=1nφ(d)∗d2,是比较简单的杜教筛,令 f = φ ∗ i d 2 , g = i d 2 f=\varphi*id_2,g=id_2 f=φ∗id2,g=id2,由狄利克雷卷积得 h = f ∗ g = ∑ d ∣ n ( φ ( d ) ∗ d 2 ) ∗ ( n d ) 2 h=f*g=\sum_{d|n}(\varphi(d)*d^2)*(\frac{n}{d})^2 h=f∗g=∑d∣n(φ(d)∗d2)∗(dn)2
展开得到 ∑ d ∣ n φ ( d ) ∗ n 2 \sum_{d|n}\varphi(d)*n^2 ∑d∣nφ(d)∗n2,根据莫比乌斯反演的结论 n = ∑ d ∣ n φ ( d ) n=\sum_{d|n}\varphi(d) n=∑d∣nφ(d),代入得 h = n 3 h=n^3 h=n3
设 f ( n ) f(n) f(n)的前缀和为 S ( n ) S(n) S(n),那么根据杜教筛的常用套路:
S ( n ) = ∑ i = 1 n i 3 − ∑ i = 2 n φ ( i ) ∗ i 2 S(n)=\sum_{i=1}^ni^3-\sum_{i=2}^n\varphi(i)*i^2 S(n)=∑i=1ni3−∑i=2nφ(i)∗i2
预处理 n 2 3 n^{\frac{2}{3}} n32以内的 f f f的前缀和,然后套用这个公式,至于 ∑ i = 1 n i 3 \sum_{i=1}^ni^3 ∑i=1ni3,显然不难得到: ( n ( n + 1 ) 2 ) 2 (\frac{n(n+1)}{2})^2 (2n(n+1))2
代码
一开始线性递推求 ∑ i = 1 n i 3 \sum_{i=1}^ni^3 ∑i=1ni3,但是不知道为什么样例都不对,调了好久也不知道为什么
//
// Created by Happig on 2020/9/27
//
#include <bits/stdc++.h>
#include <unordered_map>
#include <unordered_set>
using namespace std;
#define fi first
#define se second
#define pb push_back
#define ins insert
#define Vector Point
#define ENDL "\n"
#define lowbit(x) (x&(-x))
#define mkp(x, y) make_pair(x,y)
#define mem(a, x) memset(a,x,sizeof a);
typedef long long ll;
typedef long double ld;
typedef unsigned long long ull;
typedef pair<int, int> pii;
typedef pair<ll, ll> pll;
typedef pair<double, double> pdd;
const double eps = 1e-8;
const double pi = acos(-1.0);
const int inf = 0x3f3f3f3f;
const double dinf = 1e300;
const ll INF = 1e18;
ll Mod = 1e9 + 7;
const int maxn = 5e6 + 10;
vector<int> prime;
bitset<maxn> vis;
int N;
ll sum[maxn], phi[maxn];
unordered_map<ll, ll> mp;
ll inv;
ll qkp(ll x, ll n, ll p) {
x %= p;
ll ans = 1;
while (n) {
if (n & 1) ans = ans * x % p;
x = x * x % p;
n >>= 1;
}
return ans;
}
void init(ll n) {
sum[1] = phi[1] = 1;
vis.reset(), prime.clear();
//N = pow(n, 2 / 3);
inv = qkp(6, Mod - 2, Mod);
for (int i = 2; i < maxn; i++) {
if (!vis[i]) {
prime.push_back(i);
phi[i] = i - 1;
}
for (int j = 0; j < prime.size() && 1LL * i * prime[j] < maxn; j++) {
vis[i * prime[j]] = 1;
if (i % prime[j]) {
phi[i * prime[j]] = phi[i] * phi[prime[j]] % Mod;
} else {
phi[i * prime[j]] = phi[i] * prime[j] % Mod;
break;
}
}
ll x = 1LL * i * i % Mod;
sum[i] = (sum[i - 1] + phi[i] * x % Mod) % Mod;
//pre[i] = (pre[i - 1] + x * i % Mod) % Mod;
}
}
ll g(ll x) {
x %= Mod; //注意
return x * (x + 1) / 2 % Mod;
}
ll f(ll x) {
x %= Mod; //注意
return x * (x + 1) % Mod * (x + x + 1) % Mod * inv % Mod;
}
ll getSum(ll n) {
if (n < maxn) return sum[n];
if (mp[n]) return mp[n];
ll x = g(n);
ll ans = x * x % Mod;
for (ll l = 2, r; l <= n; l = r + 1) {
r = n / (n / l);
if (r > n) r = n;
ans = (ans - (f(r) - f(l - 1) + Mod) % Mod * getSum(n / l) % Mod + Mod) % Mod;
}
return mp[n] = ans;
}
int main() {
//freopen("in.txt","r",stdin);
//freopen("out.txt","w",stdout);
ios_base::sync_with_stdio(0), cin.tie(0), cout.tie(0);
ll n;
cin >> Mod >> n;
init(n);
ll ans = 0;
for (ll l = 1, r; l <= n; l = r + 1) {
r = n / (n / l);
if (r > n) r = n;
ll res = g(n / l) * g(n / l) % Mod;
ans += (getSum(r) - getSum(l - 1) + Mod) % Mod * res % Mod;
ans %= Mod;
}
cout << ans << ENDL;
return 0;
}