[BZOJ4916]神犇和蒟蒻
试题描述
很久很久以前,有一只神犇叫yzy;
很久很久之后,有一只蒟蒻叫lty;
输入
请你读入一个整数 (N); (1 le N le 10^9), (A)、(B) 模 (10^9+7);
输出
请你输出一个整数 (A=sum_{i=1}^N{mu (i^2)});
请你输出一个整数 (B=sum_{i=1}^N{varphi (i^2)});
输入示例
1
输出示例
1
1
数据规模及约定
见“输入”
题解
知道莫比乌斯函数定义的人都知道 (A) 恒为 (1)。
然后我们把欧拉函数展开一下就会发现 (varphi(i^2) = i cdot varphi(i))。
那么就是杜教筛 (i cdot varphi(i)) 的前缀和。这个东西很神奇,你会发现将 (i cdot varphi(i)) 和 (i) 狄利克雷卷积一下就出来了。
[sum_{i=1}^N sum_{j|i} j cdot varphi(j) cdot frac{i}{j} \
= sum_{i=1}^N i sum_{j|i} varphi(j) \
= sum_{i=1}^N i^2
]
目前还不知道这能干什么,令 (F(N) = sum_{i=1}^N i cdot varphi(i)),接着往下推就知道了
[sum_{i=1}^N sum_{j|i} j cdot varphi(j) cdot frac{i}{j} \
= sum_{j=1}^N { j cdot varphi(j) sum_{i=1}^{lfloor frac{N}{j}
floor} i } \
= sum_{i=1}^N { i sum_{j=1}^{lfloor frac{N}{i}
floor} j cdot varphi(j) } \
= sum_{i=1}^N i cdot F left( lfloor frac{N}{i}
floor
ight)
]
于是有
[sum_{i=1}^N i^2 = sum_{i=1}^N i cdot F left( lfloor frac{N}{i}
floor
ight) \
F(N) = sum_{i=1}^N i^2 - sum_{i=2}^N i cdot F left( lfloor frac{N}{i}
floor
ight)
]
解决!
#include <iostream>
#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <cctype>
#include <algorithm>
#include <map>
using namespace std;
#define rep(i, s, t) for(int i = (s), mi = (t); i <= mi; i++)
#define dwn(i, s, t) for(int i = (s), mi = (t); i >= mi; i--)
int read() {
int x = 0, f = 1; char c = getchar();
while(!isdigit(c)){ if(c == '-') f = -1; c = getchar(); }
while(isdigit(c)){ x = x * 10 + c - '0'; c = getchar(); }
return x * f;
}
#define maxn 1000010
#define MOD 1000000007
#define inv6 166666668
#define LL long long
bool vis[maxn];
int prime[maxn], cp, phi[maxn], sum[maxn];
void init() {
int n = maxn - 1;
phi[1] = sum[1] = 1;
rep(i, 2, n) {
if(!vis[i]) prime[++cp] = i, phi[i] = i - 1;
for(int j = 1; j <= cp && i * prime[j] <= n; j++) {
vis[i*prime[j]] = 1;
if(i % prime[j] == 0) {
phi[i*prime[j]] = phi[i] * prime[j];
break;
}
phi[i*prime[j]] = phi[i] * (prime[j] - 1);
}
sum[i] = sum[i-1] + (LL)phi[i] * i % MOD;
if(sum[i] >= MOD) sum[i] -= MOD;
}
return ;
}
map <int, int> hSum;
LL pre(int n) { return ((LL)n * (n + 1) >> 1) % MOD; }
int calc(int n) {
if(n < maxn) return sum[n];
if(hSum.count(n)) return hSum[n];
int ans = (LL)n * (n + 1) % MOD * (n << 1 | 1) % MOD * inv6 % MOD;
for(int i = 2; i <= n; ) {
int r = min(n / (n / i), n);
ans -= (LL)calc(n / i) * (pre(r) - pre(i - 1) + MOD) % MOD;
if(ans < 0) ans += MOD;
i = r + 1;
}
return hSum[n] = ans;
}
int main() {
init();
int n = read();
printf("1
%d
", calc(n));
return 0;
}