[数论][莫比乌斯反演][杜教筛] 51Nod 1220 约数之和

Description

求∑i=1n∑i=1nd(ij)

Solution

∑i=1n∑j=1nd(ij)A(n)B(n)S(n)=======∑i=1n∑j=1n∑a∣i∑b∣j[(a,b)=1]ajb∑k=1nμ(k)kA(⌊nk⌋)B(⌊nk⌋)∑a=1na⌊na⌋∑b=1n∑j=1⌊nb⌋j∑i=1n(μ⋅id)(i)∑i=1n[(μ∗1)⋅id](i)−∑i=2nS(⌊ni⌋)i1−∑i=2nS(⌊ni⌋)i

A,B分块搞，μ⋅id杜教筛。

#include <bits/stdc++.h>
using namespace std;

typedef long long ll;
const int MOD = 1000000007;
const int N = 1010101;

inline char get(void) {
static char buf[100000], *S = buf, *T = buf;
if (S == T) {
T = (S = buf) + fread(buf, 1, 100000, stdin);
if (S == T) return EOF;
}
return *S++;
}
template<typename T>
inline void read(T &x) {
static char c; x = 0; int sgn = 0;
for (c = get(); c < '0' || c > '9'; c = get()) if (c == '-') sgn = 1;
for (; c >= '0' && c <= '9'; c = get()) x = x * 10 + c - '0';
if (sgn) x = -x;
}

int prime
, vis
;
int mu
, pre
;
int n, lim, ans, Pcnt;
map<int, int> mp;

inline int Mod(int x) {
return (x % MOD + MOD) % MOD;
}
inline void Add(int &x, int a) {
x += a; while (x >= MOD) x -= MOD;
}
inline int Sum1(int n) {
return (ll)n * (n + 1) / 2 % MOD;
}
inline int Sum1(int l, int r) {
return Mod(Sum1(r) - Sum1(l - 1));
}
inline void Pre(int n) {
mu[1] = 1; int x;
for (int i = 2; i <= n; i++) {
if (!vis[i]) {
mu[i] = -1; prime[++Pcnt] = i;
}
for (int j = 1; j <= Pcnt && (x = prime[j] * i) <= n; j++) {
vis[x] = 1;
if (i % prime[j]) mu[x] = -mu[i];
else {
mu[x] = 0; break;
}
}
}
for (int i = 1; i <= n; i++)
Add(pre[i], Mod(mu[i] * i + pre[i - 1]));
}
inline int S(int n) {
if (n <= lim) return pre
;
if (mp.count(n)) return mp
;
int res = 1, pos;
for (int i = 2; i <= n; i = pos + 1) {
pos = n / (n / i);
Add(res, MOD - (ll)S(n / i) * Sum1(i, pos) % MOD);
}
return mp
= res;
}
inline int A(int n) {
int res = 0, pos;
for (int a = 1; a <= n; a = pos + 1) {
pos = n / (n / a);
Add(res, (ll)Sum1(a, pos) * (n / a) % MOD);
}
return res;
}
inline int B(int n) {
int res = 0, pos;
for(int b = 1; b <= n; b = pos + 1) {
pos = n / (n / b);
Add(res, (ll)Sum1(n / b) * (pos - b + 1) % MOD);
}
return res;
}

int main(void) {
freopen("1.in", "r", stdin);
freopen("1.out", "w", stdout);
read(n); Pre(lim = pow(n, 0.666));
int pos;
for (int i = 1; i <= n; i = pos + 1) {
pos = n / (n / i);
Add(ans, (ll)Mod(S(pos) - S(i - 1)) * A(n / i) % MOD * B(n / i) % MOD);
}
cout << ans << endl;
return 0;
}