对于各种积性函数,都可以通过两种方法进行计算:
- 分解质因数
- 线性筛
那么,我们就可以得到一些常用的积性函数值。
typedef pair<int, int> P;
int prime[MAXN], cnt, phi[MAXN], mu[MAXN], mnsum[MAXN], a[MAXN];
bool isp[MAXN];
P f[MAXN];//约数和
void Euler(int n) {
mu[1] = 1;
phi[1] = 1;
f[1] = P(1, 1);
for(int i = 2; i <= n; i++) {
f[i].second = i;
if(!isp[i]) {
prime[++cnt] = a[i] = i;
phi[i] = i - 1;
mu[i] = -1;
f[i].first = mnsum[i] = i + 1;
}
for(int j = 1; j <= cnt && i * prime[j] <= n; j++) {
isp[i * prime[j]] = 1;
if(i % prime[j] == 0) {
phi[i * prime[j]] = phi[i] * prime[j];
mu[i * prime[j]] = 0;
mnsum[i * prime[j]] = mnsum[i] + a[i] * prime[j];
f[i * prime[j]].first = f[i].first / mnsum[i] * mnsum[i *
prime[j]]; a[i * prime[j]] = a[i] * prime[j];
break;
}
phi[i * prime[j]] = phi[i] * (prime[j] - 1);
mu[i * prime[j]] = -mu[i];
f[i * prime[j]].first = f[i].first * (prime[j] + 1);
a[i * prime[j]] = prime[j];
mnsum[i * prime[j]] = prime[j] + 1;
}
}
}
先来一道简单题:bzoj2005
题目大意:求
ΣniΣmjgcd(i,j)∗2−1
然后推一推式子
ΣniΣmjgcd(i,j)∗2−1
=
−n∗m+2∗ΣniΣmjgcd(i,j)
=
−n∗m+2∗ΣniΣmjΣd|gcd(i,j)ϕ(d)
=
−n∗m+2∗Σmax(n,m)dϕ(d)∗⌊nd⌋⌊md⌋
然后枚举d这道题就可以A掉了,复杂度为O(n)。
#include <cstdio>
#include <algorithm>
#include <cmath>
using namespace std;
const int MAXN = 100003;
int phi[MAXN], prime[MAXN], cnt;
bool isp[MAXN];
void Euler(int n) {
phi[1] = 1;
for(int i = 2; i <= n; i++) {
if(!isp[i]) {
prime[++cnt] = i;
phi[i] = i - 1;
}
for(int j = 1; j <= cnt && i * prime[j] <= n; j++) {
isp[i * prime[j]] = 1;
if(i % prime[j] == 0) {
phi[i * prime[j]] = phi[i] * prime[j];
break;
}
else phi[i * prime[j]] = phi[i] * (prime[j] - 1);
}
}
}
long long ans;
int main() {
int n, m, k;
scanf("%d%d", &n, &m);
k = max(n, m);
Euler(k);
ans -= (long long)n * m;
for(int d = 1; d <= k; d++) {
ans += 2LL * phi[d] * (n / d) * (m / d);
}
printf("%lld\n", ans);
return 0;
}
来看另一道题:bzoj1101
题目大意:对于给定的整数N,M和d,有多少正整数对x,y,满足x<=N,y<=M,并且gcd(x,y)=d。
等价于x<=N/d,y<=M/d,互质的x,y的对数。
那么我们另n=N/d, m = M/d
于是原题就变成了求
ΣniΣmje(gcd(i,j))
其中
e(x)=(x==1)=Σd|xμ(d)
然后根据套路我们再来推一波式子
ΣniΣmje(gcd(i,j))
=
ΣniΣmjΣd|gcd(i,j)μ(d)
=
Σmin(n,m)dμ(d)∗⌊nd⌋⌊md⌋
根据上一题的做法,每次询问复杂度O(n),显然不能胜任。
但是我们观察式子,显然
⌊nd⌋⌊md⌋
的取值只有
n√+m−−√
个,那么我们可以预处理
μ
的前缀和,并分块完成。
#include <cstdio>
#include <algorithm>
using namespace std;
const int MAXN = 50003;
int mu[MAXN], prime[MAXN], cnt;
bool isp[MAXN];
void Euler(int n) {
mu[1] = 1;
for(int i = 2; i <= n; i++) {
if(!isp[i]) {
prime[++cnt] = i;
mu[i] = -1;
}
for(int j = 1; j <= cnt && i * prime[j] <= n; j++) {
isp[i * prime[j]] = 1;
if(i % prime[j] == 0) {
mu[i * prime[j]] = 0;
break;
}
mu[i * prime[j]] = -mu[i];
}
}
}
int cal(int n, int m) {
if(n > m) swap(n, m);
int res = 0, last;
for(int i = 1; i <= n; i = last + 1) {
last = min(n / (n / i), m / (m / i));
res += (n / i) * (m / i) * (mu[last] - mu[i - 1]);
}
return res;
}
int main() {
Euler(50000);
for(int i = 2; i <= 50000; i++) mu[i] += mu[i - 1];
int T, a, b, d;
scanf("%d", &T);
while(T--) {
scanf("%d%d%d", &a, &b, &d); a /= d; b /= d;
printf("%d\n", cal(a, b));
}
return 0;
}
还有一题:bzoj3529
题目大意:另F(i)表示i的约数和,q次给定n,m,a,求
ΣniΣmjandF(gcd(i,j))<=aF(gcd(i,j))mod231
有个a的限制,好像式子并不能推了啊。所以先把它去掉。
另
g(i)=ΣnxΣmye(gcd(x,y)==i)
那么可以得到
g(i)=Σi|dμ(di)∗⌊nd⌋⌊md⌋
于是
ans=Σmin(n,m)iF(i)∗g(i)
展开化简得到
Σmin(n,m)d⌊nd⌋⌊md⌋Σi|dF(i)μ(di)
然后对
Σi|dF(i)μ(di)
求个前缀和
枚举每个i更新i的倍数即可
现在考虑将a离线处理,询问按a排序
用树状数组维护前缀和即可。
#include <cstdio>
#include <algorithm>
using namespace std;
typedef pair<int, int> P;
const int MAXN = 100003;
int prime[MAXN], cnt, mu[MAXN], mnsum[MAXN], a[MAXN];
bool isp[MAXN];
P f[MAXN];
void Euler(int n) {
mu[1] = 1;
f[1] = P(1, 1);
for(int i = 2; i <= n; i++) {
f[i].second = i;
if(!isp[i]) {
prime[++cnt] = a[i] = i;
mu[i] = -1;
f[i].first = mnsum[i] = i + 1;
}
for(int j = 1; j <= cnt && i * prime[j] <= n; j++) {
isp[i * prime[j]] = 1;
if(i % prime[j] == 0) {
mu[i * prime[j]] = 0;
mnsum[i * prime[j]] = mnsum[i] + a[i] * prime[j];
f[i * prime[j]].first = f[i].first / mnsum[i] * mnsum[i * prime[j]];
a[i * prime[j]] = a[i] * prime[j];
break;
}
mu[i * prime[j]] = -mu[i];
f[i * prime[j]].first = f[i].first * (prime[j] + 1);
a[i * prime[j]] = prime[j];
mnsum[i * prime[j]] = prime[j] + 1;
}
}
}
int bit[MAXN];
void add(int x, int v) {
for(; x <= 100000; x += x & -x) bit[x] += v;
}
int sum(int x) {
int res = 0;
for(; x; x -= x & -x) res += bit[x];
return res;
}
struct Node {
int n, m, a, id;
bool operator < (const Node &x) const {
return a < x.a;
}
inline void read(int i) {
id = i; scanf("%d%d%d", &n, &m, &a);
}
}q[MAXN];
int ans[MAXN];
int cal(int n, int m) {
if(n > m) swap(n, m);
int res = 0, last;
for(int i = 1; i <= n; i = last + 1) {
last = min(n / (n / i), m / (m / i));
res += (n / i) * (m / i) * (sum(last) - sum(i - 1));
}
return res & 0x7fffffff;
}
int main() {
Euler(100000);
sort(f + 1, f + 100001);
int T;
scanf("%d", &T);
for(int i = 1; i <= T; i++) q[i].read(i);
sort(q + 1, q + T + 1);
int cur = 1, N, M, A;
for(int Q = 1; Q <= T; Q++) {
N = q[Q].n, M = q[Q].m, A = q[Q].a;
while(cur <= 100000 && f[cur].first <= A) {
for(int i = f[cur].second; i <= 100000; i += f[cur].second)
add(i, f[cur].first * mu[i / f[cur].second]);
cur++;
}
ans[q[Q].id] = cal(N, M);
}
for(int i = 1; i <= T; i++) printf("%d\n", ans[i]);
return 0;
}
来一道lcm的题:bzoj2154
题目大意:求
ΣniΣmjlcm(i,j)mod20101009
这里跟之前的区别在于,gcd(i,j)在分母上,那么考虑反演
另
f(x)=1x
F(x)=f
X
μ
这里X是Dirichlet积
即
F(x)=Σd|xf(d)∗μ(xd)
那么由莫比乌斯反演我们有
f=F
X
1
那么
=
ΣniΣmji∗j∗f(gcd(i,j))
=
ΣniΣmji∗j∗Σd|gcd(i,j)F(d)
另
sum(x,y)=ΣxiΣyji∗j
于是可以进一步化简
ans=ΣndF(d)∗d∗d∗sum(nd,md)
#include <iostream>
#include <cstdio>
#include <algorithm>
using namespace std;
typedef long long ll;
const int mod = 20101009, MAXN = 10000003;
int mu[MAXN], prime[MAXN], sum[MAXN];
int cnt;
bool isp[MAXN];
void getmu(int n) {
mu[1] = 1;
for(int i = 2; i <= n; i++) {
if(!isp[i]) {
mu[i] = -1;
prime[++cnt] = i;
}
for(int j = 1; j <= cnt && i * prime[j] <= n; j++) {
isp[i * prime[j]] = 1;
if(i % prime[j] == 0) {
mu[i * prime[j]] = 0;
break;
}
mu[i * prime[j]] = -mu[i];
}
}
}
ll n, m, ans;
ll query(ll x, ll y) { return (x * (x + 1) / 2 % mod) * (y * (y + 1) / 2 % mod) % mod; }
ll F(ll x, ll y) {
ll res = 0, last;
for(ll i = 1; i <= min(x, y); i = last + 1) {
last = min(x / (x / i), y / (y / i));
res = (res + (sum[last] - sum[i - 1]) * query(x / i, y / i) % mod) % mod;
}
return res;
}
int main() {
cin>>n>>m;
getmu(min(n, m));
for(ll i = 1; i <= min(n, m); i++) sum[i] = (sum[i - 1] + (i * i * mu[i]) % mod) % mod;
ll last;
for(ll d = 1; d <= min(n, m); d = last + 1) {
last = min(n / (n / d), m / (m / d));
ans = (ans + (last - d + 1) * (d + last) / 2 % mod * F(n / d, m / d) % mod) % mod;
}
ans = (ans + mod) % mod;
cout<<ans<<endl;
return 0;
}
P.S.:推式子治好了我多年的公式恐惧症啊
总结一下:当推式子陷入僵局的时候,尝试将一些东西处理出来,通过分块、反演等方法降低复杂度。对于式子本身的特性进行观察,有时也可以发现优化的地方。
在推式子过程中,将有关变量放在一块处理是化简的常用技巧。