[SPOJ] DIVCNT2 - Counting Divisors (square) 积性函数前缀和

75 阅读 0 评论 50 点赞

我是靠谱客的博主温暖豆芽，这篇文章主要介绍[SPOJ] DIVCNT2 - Counting Divisors (square) 积性函数前缀和，现在分享给大家，希望可以做个参考。

题面
题目要我们求

\sum i = 1 n σ 0 (i 2)

把

i 分解质因数

pα11∗pα22...pαkk
就是求

\sum i = 1 n \prod j = 1 k (2 α j + 1)

考虑每个乘式选

2αj 还是

1 ，枚举集合就是

\sum i = 1 n \sum S \subseteq {1, 2, . . ., k} 2 | S | * \prod j α j

把

S 集合看成选

αj>0的，把

∏jαj 看成枚举每个

>0 的

αj ，事实上就是枚举

i 的约数，每个约数的贡献是

2ω(d)，其中

ω(i) 表示

i 的质因数个数，所以就是

\sum i = 1 n \sum d | n 2 ω (d)

其实把

ω(d) 看成每个质因数选或不选，就是

d 的无平方因子约数的个数，所以式子变成

\sum i = 1 n \sum d | n \sum p | d μ 2 (p)

改成枚举

μ2(i) ，就是

\sum i = 1 n μ 2 (i) \sum j = 1 ⌊ n i ⌋ σ 0 (j)

于是利用整除分块，我们只要知道

μ2(i) 和

σ0(i) 的前缀和，按照杜教筛的复杂度分析，就能在预处理前

n23 项的情况下，

O(n23) 求出。
前者就是

1 n 的无平方因子数个数，考虑容斥

O(n√)

\sum i = 1 n μ 2 (i) = \sum i = 1 n \sqrt μ (i) * ⌊ n i 2 ⌋

后者就是sb套路，整除分块

O(n√)

\sum i = 1 n σ 0 (i) = \sum i = 1 n ⌊ n i ⌋

代码：

复制代码

#include<iostream>
#include<cstdio>
#include<cstring>
#include<cmath>
#define ll long long
using namespace std;
const int maxn=50000000;
ll n;
int pri[10000000],num,mu[maxn+5],d[maxn+5],sm[maxn+5],sd[maxn+5],ci[maxn+5];
bool flag[maxn+5];
void getpri(int n)
{
memset(flag,1,sizeof(flag));
flag[1]=0;mu[1]=1;d[1]=1;
for(int i=2;i<=n;i++)
{
if(flag[i]) pri[++num]=i,mu[i]=-1,ci[i]=1,d[i]=2;
for(int j=1;j<=num&&i*pri[j]<=n;j++)
{
int v=i*pri[j];
flag[v]=0;ci[v]=1;
if(i%pri[j]==0) {ci[v]+=ci[i];mu[v]=0;d[v]=d[i]*(ci[v]+1)/ci[v];break;}
mu[v]=-mu[i];d[v]=(d[i]<<1);
}
}
}
ll qsd(ll m)
{
if(m<=maxn) return sd[m];
ll re=0;
for(ll l=1,r=1,k;l<=m;l=r+1)
{
k=m/l;r=m/k;
re+=(ll)(r-l+1)*k;
}
return re;
}
ll qsm(ll m)
{
if(m<=maxn) return sm[m];
ll gen=round(sqrt(m)+1),re=0;
for(ll i=1;i<=gen;i++)
re+=mu[i]*(m/(i*i));
return re;
}
ll solve()
{
ll re=0;
for(ll l=1,r=1,k;l<=n;l=r+1)
{
k=n/l;r=n/k;
re+=qsd(k)*(qsm(r)-qsm(l-1));
}
return re;
}
int main()
{
int ca,lim=maxn;
scanf("%d",&ca);
if(ca>800) lim=10000;
getpri(lim);
for(int i=1;i<=lim;i++)
sm[i]=sm[i-1]+mu[i]*mu[i],sd[i]=sd[i-1]+d[i];
while(ca--)
{
scanf("%lld",&n);
printf("%lldn",solve());
}
return 0;
}

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
#include<iostream>
#include<cstdio>
#include<cstring>
#include<cmath>
#define ll long long
using namespace std;
const int maxn=50000000;
ll n;
int pri[10000000],num,mu[maxn+5],d[maxn+5],sm[maxn+5],sd[maxn+5],ci[maxn+5];
bool flag[maxn+5];
void getpri(int n)
{
memset(flag,1,sizeof(flag));
flag[1]=0;mu[1]=1;d[1]=1;
for(int i=2;i<=n;i++)
{
if(flag[i]) pri[++num]=i,mu[i]=-1,ci[i]=1,d[i]=2;
for(int j=1;j<=num&&i*pri[j]<=n;j++)
{
int v=i*pri[j];
flag[v]=0;ci[v]=1;
if(i%pri[j]==0) {ci[v]+=ci[i];mu[v]=0;d[v]=d[i]*(ci[v]+1)/ci[v];break;}
mu[v]=-mu[i];d[v]=(d[i]<<1);
}
}
}
ll qsd(ll m)
{
if(m<=maxn) return sd[m];
ll re=0;
for(ll l=1,r=1,k;l<=m;l=r+1)
{
k=m/l;r=m/k;
re+=(ll)(r-l+1)*k;
}
return re;
}
ll qsm(ll m)
{
if(m<=maxn) return sm[m];
ll gen=round(sqrt(m)+1),re=0;
for(ll i=1;i<=gen;i++)
re+=mu[i]*(m/(i*i));
return re;
}
ll solve()
{
ll re=0;
for(ll l=1,r=1,k;l<=n;l=r+1)
{
k=n/l;r=n/k;
re+=qsd(k)*(qsm(r)-qsm(l-1));
}
return re;
}
int main()
{
int ca,lim=maxn;
scanf("%d",&ca);
if(ca>800) lim=10000;
getpri(lim);
for(int i=1;i<=lim;i++)
sm[i]=sm[i-1]+mu[i]*mu[i],sd[i]=sd[i-1]+d[i];
while(ca--)
{
scanf("%lld",&n);
printf("%lldn",solve());
}
return 0;
}