「BZOJ」DZY Loves Math

DZY Loves Math 系列.


「BZOJ 3309」DZY Loves Math

Description

定义 f(n)

f(n)=\max_{p | n}\left\{\nu_{p}(n)\right\}

给定 a,b ,求

\sum_{i=1}^{a}\sum_{j=1}^{b}f\left(\gcd(i, j)\right)

每个测试点 T 组数据。

Constraints

1\le T\le 10^{4},1\le a,b\le 10^{7}

Analysis

下文中记号 ' \ast ' 均表示 Dirichlet Convolution.

先推一波式子。

a<b ,有:

\begin{aligned} &\sum_{i=1}^{a}\sum_{j=1}^{b}f\left(\gcd(i, j)\right)\\ =&\sum_{d=1}^{a}f(d)\sum_{i=1}^{a}\sum_{j=1}^{b}\left[(i, j)=d\right]\\ =&\sum_{d=1}^{a}f(d)\sum_{i=1}^{\left\lfloor\frac{a}{d}\right\rfloor}\sum_{j=1}^{\left\lfloor\frac{b}{d}\right\rfloor}\left[(i, j)=1\right]\\ =&\sum_{d=1}^{a}f(d)\sum_{d'=1}^{\left\lfloor\frac{a}{d}\right\rfloor}\mu\left({d}'\right)\left\lfloor\frac{a}{dd'}\right\rfloor\left\lfloor\frac{b}{dd'}\right\rfloor\\ =&\sum_{s=1}^{c}\left\lfloor\frac{a}{s}\right\rfloor\left\lfloor\frac{b}{s}\right\rfloor\sum_{d | s}f(d)\mu\left(\frac{s}{d}\right) \end{aligned}

g=f\ast\mu ,则需要对 g(n) 快速求值:

g(n)=\sum_{d | n}f(d)\mu\left(\frac{n}{d}\right)

考虑到对于所有有贡献的 d \mu\left(\frac{n}{d}\right)\ne 0 ,则可知

\forall p | \frac{n}{d}:\nu_{p}\left(\frac{n}{d}\right)=1

\forall p | d:0\le\nu_{p}(n)-\nu_{p}(d)\le 1

可以发现 f(d) 只可能有 f(n) f(n)-1 两种取值。

定义 h(n)

h(n)=\sum_{p | n}\left[\nu_{p}(n)=f(n)\right]

亦即 n 所有质因数中次数最大的质因数个数。

h(n)=\omega(n) 时,考虑其中 f(n)-1 的系数,易知其为 (-1)^{\omega(n)}
考虑其中 f(n) 的系数:

\begin{aligned} &\sum_{i=0}^{\omega(n)-1}(-1)^{i}\binom{\omega(n)}{i}\\ =&\left(\left(1-1\right)^{\omega(n)}-(-1)^{\omega(n)}\right)\\ =&(-1)^{\omega(n)+1} \end{aligned}

故当 h(n)=\omega(n) 时,有:

\begin{aligned} g(n)=&(-1)^{\omega(n)}\left(f(n)-1\right)+(-1)^{\omega(n)+1}f(n)\\ =&(-1)^{\omega(n)+1} \end{aligned}

h(n)<\omega(n) 时,考虑其中 f(n)-1 的系数:

\begin{aligned} &(-1)^{h(n)}\sum_{i=0}^{\omega(n)-h(n)}(-1)^{i}\binom{\omega(n)-h(n)}{i}\\ =&(-1)^{h(n)}\left(1-1\right)^{\omega(n)-h(n)}\\ =&0 \end{aligned}

考虑其中 f(n) 的系数:

\begin{aligned} &\sum_{i=0}^{h(n)-1}(-1)^{i}\binom{h(n)}{i}\sum_{j=0}^{\omega(n)-h(n)}(-1)^{j}\binom{\omega(n)-h(n)}{j}\\ =&\sum_{i=0}^{h(n)-1}(-1)^{i}\binom{h(n)}{i}\left(1-1\right)^{\omega(n)-h(n)}\\ =&0 \end{aligned}

故当 h(n)<\omega(n) 时,有:

g(n)=0

则:

g(n)=\begin{cases} 0 & h(n)=\omega(n) \\ (-1)^{\omega(n)+1} & h(n)<\omega(n) \end{cases}

对于每个数,设其最小质因子为 p ,记录 \nu_{p}(n) p^{\nu_{p}(n)} 即可线性筛出 g(n)

时间复杂度 O(n)

Code

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
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
#include<cstdio>

using i64=long long;
constexpr int Maxn=10000000;
constexpr int U=10000000;

template<typename _Tp>
inline _Tp Min(const _Tp&x,const _Tp&y)
{return x<y?x:y;}

struct IO{
#define ISIZE 16384
#define OSIZE 65536

char IBUF[ISIZE],*IS,*IT;
char OBUF[OSIZE+128],*OS,*OT;

IO():IS(IBUF),IT(IBUF),OS(OBUF),OT(OBUF+OSIZE){}
~IO(){Flush();}

inline char gc(){return(IS==IT)&&(IT=(IS=IBUF)+fread(IBUF,1,ISIZE,stdin),IS==IT)?EOF:*IS++;}
template<typename _Tp>
inline IO&operator>>(_Tp&x){
_Tp s=0;char c;
if(IS+10<IT){
do c=*IS++;while(c<48||c>57);
do s=s*10+c-48,c=*IS++;while(c>47&&c<58);
}else{
do c=gc();while(c<48||c>57);
do s=s*10+c-48,c=gc();while(c>47&&c<58);
}return(x=s,*this);
}

inline void Flush(){fwrite(OBUF,1,OS-OBUF,stdout);OS=OBUF;}
inline void pc(const char&c){if(OS>OT)Flush();*OS++=c;}
inline IO&operator<<(const char&c){return(pc(c),*this);}
template<typename _Tp>
inline IO&operator<<(const _Tp&s){
if(OS>OT)Flush();
static char OStack[30],*Top=OStack;
if(!s)*OS++='0';
else{
_Tp x=s;
for(_Tp y;x;x=y){
y=x/10;
*++Top=x-y*10+48;
}for(;Top!=OStack;)
*OS++=*Top--;
}return*this;
}

#undef ISIZE
#undef OSIZE
}io;

bool notp[U+1];
int pri[U/7];
int cnt;

int nu1[Maxn+1];
int pp1[Maxn+1];
int g[Maxn+1];

void Sieve(){
for(int i=2;i<=U;++i){
if(!notp[i])
pri[++cnt]=pp1[i]=i,nu1[i]=g[i]=1;
for(int j=1;j<=cnt;++j){
int pn=pri[j],v=i*pn;
if(v>U)break;
notp[v]=true;
if(i%pn)
nu1[v]=1,pp1[v]=pn,
g[v]=(nu1[i]==1?-g[i]:0);
else{
nu1[v]=nu1[i]+1;
pp1[v]=pp1[i]*pn;
int d=i/pp1[i];
g[v]=(d==1?1:(nu1[v]==nu1[d]?-g[d]:0));
break;
}
}
}
}

i64 calc(const int&x,const int&y){
i64 Ans=0;
for(int d=1,nxd,L=Min(x,y);d<=L;d=nxd+1){
nxd=Min(x/(x/d),y/(y/d));
Ans+=(i64)(g[nxd]-g[d-1])*(x/d)*(y/d);
}return Ans;
}

int main(){
Sieve();
for(int i=1,v=0;i<=U;++i)
v+=g[i],g[i]=v;
int T;
for(io>>T;T;--T){
int a,b;
io>>a>>b;
io<<calc(a,b)<<'\n';
}return 0;
}

「BZOJ 3462」DZY Loves Math Ⅱ

Description

对于正整数 S,n ,称 \left(p_{1},p_{2},...,p_{k}\right) k 为任意正整数)为 n S -质数拆分,当且仅当:

  1. p_{1}+p_{2}+...+p_{k-1}+p_{k}=n
  2. p_{1},p_{2},...,p_{k-1},p_{k} 均为质数
  3. p_{1}\le p_{2}\le ...\le p_{k-1}\le p_{k}
  4. \operatorname{lcm}\left(p_{1},p_{2},...,p_{k}\right)=S

现给定一个正整数 S Q 次询问对于一个正整数 n ,其有多少种不同的质数拆分。

答案对 10^{9}+7 取模。

Constraints

2\le S\le 2\times 10^{6},1\le n\le 10^{18},1\le Q\le 10^{5}

Analysis

首先显然当 \exists p:\nu_{p}(S)>1 时答案必为 0

发现 S 的每个质因数 p 都必须出现,则方案不同即每种质因子的出现次数不同。
p_{i} 出现次数为 c_{i} ,问题转化为求满足 \sum_{i=1}^{k}c_{i}p_{i}=n,1\le c_{i} \left(c_{1},c_{2},...,c_{k}\right) 的方案数。

易知对于任意 p_{i} ,其必为 S 的质因数,故 c_{i} 可以表示成 a_{i}\frac{S}{p_{i}}+b_{i} 的形式,其中 b_{i}\equiv c_{i}\pmod{\frac{S}{p_{i}}}

有:

\begin{aligned} n=&\sum_{i=1}^{k}c_{i}p_{i}\\ =&\sum_{i=1}^{k}\left(a_{i}\frac{S}{p_{i}}+b_{i}\right)p_{i}\\ =&\sum_{i=1}^{k}\left(a_{i}S+b_{i}p_{i}\right)\\ =&S\sum_{i=1}^{k}a_{i}+\sum_{i=1}^{k}b_{i}p_{i} \end{aligned}

故问题被转化成了如下形式:

对于每个满足 zS\le n z ,求出:

  1. \sum_{i=1}^{k}a_{i}=z 的方案数,设其为 ans_{1}
  2. \sum_{i=1}^{k}b_{i}p_{i}=n-zS 的方案数,设其为 ans_{2}

则答案即为 \sum_{z\le\frac{n}{S}}ans_{1}\cdot ans_{2}

ans_{1} 的计算是组合数的经典模型,此处不作赘述。

ans_{1}=\binom{z+k-1}{k-1}

对于 ans_{2} ,直接 DP 即可。
由于

2\times 3\times 5\times 7\times 11\times 13\times 17=510510<2\times 10^{6}

2\times 3\times 5\times 7\times 11\times 13\times 17\times 19=9699690>2\times 10^{6}

S 最多有 7 个质因子,亦即 k\le 7

复杂度 O\left(k^{2}S+Qk\right)

Code

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
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
#include<cstdio>

using i64=long long;

constexpr int p=1000000007;
constexpr int MaxS=2000000;
constexpr int Maxk=7;

struct IO{
#define ISIZE 16384
#define OSIZE 65536

char IBUF[ISIZE],*IS,*IT;
char OBUF[OSIZE+128],*OS,*OT;

IO():IS(IBUF),IT(IBUF),OS(OBUF),OT(OBUF+OSIZE){}
~IO(){Flush();}

inline char gc(){return(IS==IT)&&(IT=(IS=IBUF)+fread(IBUF,1,ISIZE,stdin),IS==IT)?EOF:*IS++;}
template<typename _Tp>
inline IO&operator>>(_Tp&x){
_Tp s=0;char c;
if(IS+20<IT){
do c=*IS++;while(c<48||c>57);
do s=s*10+c-48,c=*IS++;while(c>47&&c<58);
}else{
do c=gc();while(c<48||c>57);
do s=s*10+c-48,c=gc();while(c>47&&c<58);
}return(x=s,*this);
}

inline void Flush(){fwrite(OBUF,1,OS-OBUF,stdout);OS=OBUF;}
inline void pc(const char&c){if(OS>OT)Flush();*OS++=c;}
inline IO&operator<<(const char&c){return(pc(c),*this);}
template<typename _Tp>
inline IO&operator<<(const _Tp&s){
if(OS>OT)Flush();
static char OStack[30],*Top=OStack;
if(!s)*OS++='0';
else{
_Tp x=s;
for(_Tp y;x;x=y){
y=x/10;
*++Top=x-y*10+48;
}for(;Top!=OStack;)
*OS++=*Top--;
}return*this;
}

#undef ISIZE
#undef OSIZE
}io;

template<typename x_tp,typename y_tp>
inline void inc(x_tp&x,const y_tp&y)
{x+=y;(p<=x)&&(x-=p);}
template<typename x_tp,typename y_tp>
inline void dec(x_tp&x,const y_tp&y)
{x-=y;(x<0)&&(x+=p);}
template<typename x_tp,typename y_tp>
inline x_tp Min(const x_tp&x,const y_tp&y)
{return x<y?x:y;}
template<typename _Tp>
inline void swap(_Tp&x,_Tp&y)
{_Tp z=x;x=y;y=z;}

int inv[Maxk+1]=
{0,1,500000004,333333336,250000002,400000003,166666668,142857144};

inline int C(const i64&n,const int&k){
int Ans=1;
for(int i=1;i<=k;++i)
Ans=(i64)Ans*((n-i+1)%p)%p*inv[i]%p;
return Ans;
}

int cnt[2][MaxS*(Maxk+2)];
int div[Maxk+1];

int main(){
int S,Q;
int divcnt=0;
int divsum=0;

io>>S>>Q;

/*Division*/
int s=S;
for(int i=2;i*i<=s;++i)if(!(s%i)){
div[++divcnt]=i;
divsum+=i;
s/=i;
if(!(s%i)){
for(;Q;--Q)
io<<'0'<<'\n';
return 0;
}
}if(s!=1)
div[++divcnt]=s,divsum+=s;

int *las=cnt[0];
int *cur=cnt[1];
/*DP*/
las[0]=1;
for(int i=1,L=S*divcnt-divsum;i<=divcnt;++i,swap(las,cur))
for(int j=0,d=div[i];j<d;++j)
for(int k=j,tmp=0;k<=L;k+=d){
inc(tmp,las[k]);
if(S+j<=k)
dec(tmp,las[k-S]);
cur[k]=tmp;
}

/*Answer*/
for(;Q;--Q){
i64 n;
io>>n;

if(n<divsum)
io<<'0'<<'\n';
else{
n-=divsum;
int Ans=0;
i64 d=n/S;
i64 r=n-d*S;
for(int i=Min((i64)divcnt,d);~i;--i)
inc(Ans(i64)C(d-i+divcnt-1,divcnt-1)*las[i*S+r]%p);
io<<Ans<<'\n';
}
}

return 0;
}

「BZOJ 3481」DZY Loves Math Ⅲ

Description

给定整数 P,Q ,求满足方程 xy\equiv Q\pmod{P} 的整数解 (x, y) 的数量。
其中 0\le x,y<P P,Q p_{1}p_{2}...p_{n},q_{1}q_{2}...q_{n} 的形式给出。

答案对 10^{9}+7 取模。

Constraints

1\le n\le 10,0\le q_{i}\le 10^{18},1\le p_{i}\le 10^{18},2\le P

Analysis

考虑固定 x 的值,计算合法的 y 的数量,由裴蜀定理可知其有解当且仅当 \left(x,P\right) | Q
此时 xy\pmod{P} 的循环节长度为 \frac{P}{\left(x,P\right)} ,故共有 \left(x,P\right) 个合法的 y
答案即为

\begin{aligned} &\sum_{x=0}^{P-1}\left[\left(x,P\right) | Q\right]\left(x,P\right)\\ =&\sum_{d=1}^{P-1}d\left[d | Q\right]\sum_{x=0}^{P-1}\left[\left(x,P\right)=d\right]\\ =&\sum_{d | Q}d\sum_{x=0}^{P-1}\left[\left(x,P\right)=d\right]\\ =&\sum_{d |(p, q)}d\sum_{x=0}^{\frac{P}{d}-1}\left[\left(x,\frac{P}{d}\right)=1\right]\\ =&\sum_{d |(p, q)}d\varphi\left(\frac{P}{d}\right) \end{aligned}

又发现 f(n)=n\varphi\left(\frac{P}{n}\right) 为积性函数,故可以进一步推导:

\begin{aligned} &\sum_{d |(p, q)}d\varphi\left(\frac{P}{d}\right)\\ =&\prod_{p |(p, q)}\sum_{i=0}^{\nu_{p}\left(\gcd(p, q)\right)}p^{i}\varphi\left(p^{\nu_{p}(p)-i}\right)\\ =&\prod_{p |(p, q)}\sum_{i=0}^{\min\left(\nu_{p}(p),\nu_{p}\left(Q\right)\right)}p^{i}p^{\nu_{p}(p)-i-1}\left(p-1\right)\\ =&\prod_{p |(p, q)}\min\left(\nu_{p}(p),\nu_{p}\left(Q\right)\right)p^{\nu_{p}(p)-1}\left(p-1\right)\\ \end{aligned}

P,Q 质因数分解后快速幂计算即可。

时间复杂度 O\left(np^{\frac{1}{4}}\log^{3}{p}+\log{P}\right)

Code

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
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
#include<cstdio>
#include<tr1/unordered_map>

using i64=long long;
using u64=unsigned long long;

constexpr int Maxd=100000;
constexpr int test_time=8;

i64 p;
template<typename _Tp>
inline void setp(const _Tp&np)
{p=np;}
template<typename _Tp>
inline _Tp Abs(const _Tp&v)
{return v<0?-v:v;}
template<typename x_tp,typename y_tp>
inline void inc(x_tp&x,const y_tp&y)
{(x>=p-y)?(x+=y-p):(x+=y);}
template<typename x_tp,typename y_tp>
inline void dec(x_tp&x,const y_tp&y)
{x-=y;(x<0)&&(x+=p);}
template<typename x_tp,typename y_tp>
inline x_tp pls(const x_tp&x,const y_tp&y)
{return(x>=p-y)?(x-p+y):(x+y);}
template<typename x_tp,typename y_tp>
inline x_tp mns(const x_tp&x,const y_tp&y)
{return(x<y)?(x-y+p):(x-y);}
template<typename x_tp,typename y_tp>
inline x_tp mtp(const x_tp&_x,const y_tp&_y){
x_tp Ans=0;
for(x_tp x=_x,y=_y;y;y>>=1,x=pls(x,x))
if(y&1)Ans=pls(Ans,x);
return Ans;
}
template<typename x_tp,typename y_tp>
inline x_tp fpow(const x_tp&s,const y_tp&t){
x_tp Ans=1;
x_tp v=s;
for(y_tp n=t;n;n>>=1,v=mtp(v,v))
if(n&1)Ans=mtp(Ans,v);
return Ans;
}
template<typename _Tp>
inline _Tp gcd(const _Tp&x,const _Tp&y)
{return y?gcd(y,x%y):x;}

struct random_generator{
u64 _seed1;
u64 _seed2;
random_generator():_seed1((u64)(new unsigned char*)),_seed2((u64)(new int*)){}

template<typename _Tp>
inline _Tp random(){
_seed1^=(_seed2<<13)^(_seed2>>17);
_seed2=(_seed2^((_seed1>>15)|(_seed1<<21)))-(_seed1&_seed2);
return Abs((_Tp)(_seed1+_seed2))%(p-1)+1;
}
}rg;

#define random(_Tp) rg.random<_Tp>()

template<typename _Tp>
bool Evedence(const _Tp&v,const _Tp&n,const _Tp&r,const _Tp&k){
i64 x=fpow(v,r);
i64 y=x;
for(int i=1;i<=k;++i,y=x){
x=mtp(x,x);
if(x==1&&y!=1&&y!=n-1)
return true;
}return x!=1;
}

template<typename _Tp>
bool Miller_Rabin(const _Tp&n){
if(n==2)
return true;
if((~n&1)||n<=1)
return false;
setp(n);
_Tp v=n-1;
int k=0;
for(;~v&1;v>>=1)++k;
for(int Time=test_time;Time;--Time)
if(Evedence(random(_Tp),n,v,(i64)k))
return false;
return true;
}

template<typename _Tp>
_Tp rho(const _Tp&n){
if(~n&1)
return 2;
if(0<n){
setp(n);
_Tp c=random(_Tp),d=1,x=random(_Tp),y=x;
for(_Tp i=1,k=2;d==1;++i){
x=pls(mtp(x,x),c);
d=gcd(n,Abs(x-y));
if(i==k)
y=x,k<<=1;
}return d;
}return 1;
}

std::tr1::unordered_map<i64,int>idx;
i64 divp[Maxd*10];
int nup[Maxd+1];
int nuq[Maxd+1];
int cnt;

bool Flag=true;

template<typename _Tp>
inline void insp(const _Tp&v){
if(Flag){
if(idx.count(v))
return(void)(++nup[idx[v]]);
divp[++cnt]=v;
idx[v]=cnt;
nup[cnt]=1;
}else if(idx.count(v)){
int idv=idx[v];
if(nuq[idv]<nup[idv])
++nuq[idv];
}
}

template<typename _Tp>
void ins(const _Tp&v){
if(v!=1){
if(Miller_Rabin(v))
return insp(v);
i64 d;
do d=rho(v);while(d==v);
ins(d);ins(v/d);
}
}

int main(){
int n;
i64 v;

scanf(" %d",&n);
for(int i=1;i<=n;++i)
scanf(" %lld",&v),ins(v);
Flag=false;
for(int i=1;i<=n;++i){
scanf(" %lld",&v);
if(!v){
for(int j=1;j<=cnt;++j)
nuq[j]=nup[j];
break;
}ins(v);
}

i64 Ans=1;
setp(1000000007);
for(int i=1;i<=cnt;++i)
Ans=mtp(Ans,mtp(pls(mtp(nuq[i]+1,(divp[i]-1)%p),(nup[i]==nuq[i])),fpow(divp[i]%p,nup[i]-1)));
printf("%lld",Ans);

return 0;
}

「BZOJ 3512」DZY Loves Math Ⅳ

Description

给定 n,m ,求

\sum_{i=1}^{n}\sum_{j=1}^{m}\varphi\left(ij\right)

答案对 10^{9}+7 取模。

Constraints

1\le n\le 10^{5},1\le m\le 10^{9}

Analysis

对于 \varphi\left(kn\right) 的取值,考虑 k n 的质因子组成可知:

\mu(k)\ne 0 时,有:

\begin{aligned} \varphi\left(kn\right)=&\varphi(n)\varphi\left(\frac{k}{\left(n,k\right)}\right)\left(n,k\right)\\ =&\varphi(n)\varphi\left(\frac{k}{\left(n,k\right)}\right)\sum_{d |\left(n,k\right)}\varphi(d)\\ =&\varphi(n)\varphi\left(\frac{k}{\left(n,k\right)}\right)\sum_{d |\left(n,k\right)}\varphi\left(\frac{\left(n,k\right)}{d}\right)\\ =&\varphi(n)\sum_{d |\left(n,k\right)}\varphi\left(\frac{k}{d}\right) \end{aligned}

\mu(k)=0 时,设

s=\max_{d | k,\mu(d)\ne 0}d

有:

\varphi\left(kn\right)=\frac{\varphi\left(sn\right)k}{s}

S\left(n,k\right)=\sum_{i=1}^{n}\varphi\left(ki\right)

\mu(k)\ne 0 时:

\begin{aligned} S\left(n,k\right)=&\sum_{i=1}^{n}\varphi\left(ki\right)\\ =&\sum_{i=1}^{n}\varphi(i)\sum_{d |\left(i,k\right)}\varphi\left(\frac{k}{d}\right)\\ =&\sum_{d | k}\varphi\left(\frac{k}{d}\right)\sum_{i=1}^{\left\lfloor\frac{n}{d}\right\rfloor}\varphi\left(id\right)\\ =&\sum_{d | k}\varphi\left(\frac{k}{d}\right)S\left(\left\lfloor\frac{n}{d}\right\rfloor,d\right) \end{aligned}

\mu(k)=0 时,有:

\begin{aligned} S\left(n,k\right)=&\sum_{i=1}^{n}\varphi\left(ki\right)\\ =&\sum_{i=1}^{n}\frac{\varphi\left(si\right)k}{s}\\ =&\frac{k}{s}\sum_{i=1}^{n}\varphi\left(si\right)\\ =&\frac{k}{s}S(n, s) \end{aligned}

于是可以得到一个递推式:

S\left(n,k\right)=\frac{k}{s}\sum_{d | s}\varphi\left(\frac{s}{d}\right)S\left(\left\lfloor\frac{n}{d}\right\rfloor,d\right)

递归计算即可,在 k=1 时使用杜教筛。

复杂度分析:
对于最底层的所有 S\left(n,1\right) ,可以通过杜教筛在 O\left(n^{\frac{2}{3}}\right) 的复杂度内求出。
对于 S\left(n,k\right) ,其在第一层就只有 O\left(k\sqrt{n}\right) 种取值,上界为 O\left(\sqrt{k}\right)

故总复杂度为 O\left(n\left(\sqrt{n}+\sqrt{m}\right)+m^{\frac{2}{3}}\right)

Code

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
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
#include<cstdio>
#include<cstring>
#include<tr1/unordered_map>

using i64=long long;

constexpr int p=1000000007;
constexpr int inv_2=(p+1)>>1;
constexpr int U=1000000;

template<typename x_tp,typename y_tp>
inline void inc(x_tp&x,const y_tp&y)
{(x>=p-y)?(x+=y-p):(x+=y);}
template<typename x_tp,typename y_tp>
inline x_tp pls(const x_tp&x,const y_tp&y)
{return x+(x>=p-y)?(y-p):y;}

template<typename x_tp,typename y_tp>
inline void dec(x_tp&x,const y_tp&y)
{x-=y;(x<0)&&(x+=p);}
template<typename x_tp,typename y_tp>
inline x_tp mns(const x_tp&x,const y_tp&y)
{return(x<y)?(x-y+p):(x-y);}

template<typename x_tp,typename y_tp>
inline x_tp fpow(const x_tp&x,const y_tp&t){
x_tp Ans=1,v=x;
for(y_tp n=t;n;n>>=1,v*=v)
if(n&1)Ans*=v;
return Ans;
}

template<typename idx_T,typename val_T,typename sum_T,idx_T _size>
struct Multiplicative_Function{

val_T val[_size+1];
sum_T sums[_size+1];
std::tr1::unordered_map<idx_T,sum_T>sumb;

sum_T(*S)(const idx_T&idx);

Multiplicative_Function():S(nullptr){
memset(val,0,sizeof val);
memset(sums,0,sizeof sums);
}
Multiplicative_Function(sum_T(*_S)(const idx_T&idx)):S(_S){
memset(val,0,sizeof val);
memset(sums,0,sizeof sums);
}

inline val_T&operator[](const idx_T&idx){return val[idx];}
inline void setS(sum_T(*_S)){S=_S;}

inline void init(){
sum_T v=(sum_T)val[1];
sums[1]=v;
for(int i=2;i<=_size;++i)
inc(v,val[i]),sums[i]=v;
}

sum_T calc(const idx_T&idx){
if(idx<=_size)
return sums[idx];
if(sumb.count(idx))
return sumb[idx];
sum_T v=S(idx);
for(idx_T d=2,nxd;nxd<idx;d=nxd+1){
nxd=(idx_T)idx/(idx/d);
dec(v,(sum_T)(nxd-d+1)*calc(idx/d)%p);
}return sumb[idx]=v%p;
}

inline sum_T operator()(const idx_T&idx)
{return calc(idx);}
};

inline int Sphi(const int&idx)
{return(i64)idx*(idx+1)%p*inv_2%p;}

Multiplicative_Function<int,int,int,U>phi(Sphi);

bool notp[U+1];
int pri[U/7];
int ms[U+1];
int cnt;

void Sieve(const int&_U){
phi[1]=1;ms[1]=1;
for(int i=2;i<=_U;++i){
if(!notp[i])
pri[++cnt]=i,phi[i]=i-1,ms[i]=i;
for(int j=1,cp=phi[i],np;j<=cnt;++j){
np=pri[j];
int v=np*i;
if(_U<v)break;
notp[v]=true;
if(i%np)
phi[v]=cp*(np-1),
ms[v]=ms[i]*np;
else{
phi[v]=cp*np;
ms[v]=ms[i];
break;
}
}
}
}

inline int _phi(const int&n)
{return mns(phi.sums[n],phi.sums[n-1]);}

std::tr1::unordered_map<int,int>_S[U+1];

int S(const int&n,const int&k){
if(!n)
return 0;
if(k==1)
return phi(n)%p;
if(_S[k].count(n))
return _S[k][n];
int Ans=0;
for(int d=1,v;d*d<=k;++d)
if(!(k%d)){
v=k/d;
inc(Ans,(i64)_phi(v)*S(n/d,d)%p);
if(d!=v)
inc(Ans,(i64)_phi(d)*S(n/v,v)%p);
}
return _S[k][n]=Ans%p;
}

int main(){
Sieve(U);
phi.init();

int n,m;
scanf(" %d %d",&n,&m);

int Ans=0;
for(int k=1;k<=n;++k)
inc(Ans,(i64)k/ms[k]*S(m,ms[k])%p);
printf("%d",Ans);

return 0;
}

「BZOJ 3560」DZY Loves Math Ⅴ

Description

给定 n 个正整数 a_{1},a_{2},...,a_{n} ,求

\sum_{i_{1} | a_{1}}\sum_{i_{2} | a_{2}}...\sum_{i_{n} | a_{n}}\varphi\left(i_{1}i_{2}...i_{n}\right)

答案对 10^{9}+7 取模。

Constraints

1\le n\le 10^{5},1\le a_{i}\le 10^{7}

Analysis

\varphi 函数的积性可知

\begin{aligned} &\sum_{i_{1} | a_{1}}\sum_{i_{2} | a_{2}}...\sum_{i_{n} | a_{n}}\varphi\left(i_{1}i_{2}...i_{n}\right)\\ =&\prod_{p | \operatorname{lcm}\{a_{n}\}}\left(\sum_{k_{1}=0}^{\nu_{p}\left(a_{1}\right)}\sum_{k_{2}=0}^{\nu_{p}\left(a_{2}\right)}...\sum_{k_{n}=0}^{\nu_{p}\left(a_{n}\right)}\varphi\left(p^{\sum_{i=1}^{n}k_{i}}\right)\right)\\ =&\prod_{p | \operatorname{lcm}\{a_{n}\}}\left(1+\frac{p-1}{p}\sum_{k_{1}=0}^{\nu_{p}\left(a_{1}\right)}\sum_{k_{2}=0}^{\nu_{p}\left(a_{2}\right)}...\sum_{k_{n}=0}^{\nu_{p}\left(a_{n}\right)}\left[\sum_{i=1}^{n}k_{i}\ne 0\right]p^{\sum_{i=1}^{n}k_{i}}\right)\\ =&\prod_{p | \operatorname{lcm}\{a_{n}\}}\left(1+\frac{p-1}{p}\left(\left(\sum_{k_{1}=0}^{\nu_{p}\left(a_{1}\right)}\sum_{k_{2}=0}^{\nu_{p}\left(a_{2}\right)}...\sum_{k_{n}=0}^{\nu_{p}\left(a_{n}\right)}p^{\sum_{i=1}^{n}k_{i}}\right)-1\right)\right)\\ =&\prod_{p | \operatorname{lcm}\{a_{n}\}}\left(1+\frac{p-1}{p}\left(\left(\prod_{i=1}^{n}\sum_{k=0}^{\nu_{p}\left(a_{i}\right)}p^{k}\right)-1\right)\right) \end{aligned}

对于每个 p 维护 p^{k} 的前缀和即可。

复杂度 O(n \log{n})

Code

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
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
#include<cstdio>
#include<algorithm>

using i64=long long;
constexpr int p=1000000007;
constexpr int Maxn=100000;
/*
2*3*5*7*11*13*17*19 = 9699690
2*3*5*7*11*13*17*19*23 = 223092870
9699690 < 10^{7} < 223092870
ω(a_{i})<=8
*/
constexpr int Maxd=800000;
/* sqrt 10^{7} = 3162 */
constexpr int U=3161;

template<typename x_tp,typename y_tp>
inline x_tp pls(const x_tp&x,const y_tp&y)
{return(x>=p-y)?(x-p+y):(x+y);}

template<typename _Tp,class __Tp>
inline _Tp fpow(const _Tp&x,const __Tp&t){
_Tp v=x,Ans=1;
for(__Tp n=t;n;n>>=1,v=(i64)v*v%p)
if(n&1)Ans=(i64)Ans*v%p;
return Ans;
}
template<typename _Tp>
inline _Tp inv(const _Tp&v)
{return fpow(v,p-2);}
template<typename _Tp>
inline _Tp sqr(const _Tp&v)
{return v*v;}

template<typename _Tp>
struct Divisor{
_Tp v;int k;
Divisor():v(0),k(0){}
Divisor(const _Tp&_v,const int&_k):v(_v),k(_k){}
inline bool operator<(const Divisor<_Tp>&rhs)
const{return v<rhs.v||(v==rhs.v&&k<rhs.k);}
};
Divisor<int>divsor[Maxd+2];
int dcnt;
#define v(i) (divsor[i].v)
#define k(i) (divsor[i].k)
inline void insd(const int&np,const int&nk)
{divsor[++dcnt]=Divisor<int>(np,nk);}

bool notp[U+1];
int pri[U/7];
int pcnt;

void Sieve(const int&_U){
for(int i=2;i<=_U;++i){
if(!notp[i])
pri[++pcnt]=i;
for(int j=1,v;j<=pcnt;++j){
if((v=i*pri[j])>_U)break;
notp[v]=true;
if(!(i%pri[j]))break;
}
}
}

struct IO{
#define ISIZE 16384

char IBUF[ISIZE],*IS,*IT;

IO():IS(IBUF),IT(IBUF){}
~IO(){}

inline char gc(){return(IS==IT)&&(IT=(IS=IBUF)+fread(IBUF,1,ISIZE,stdin),IS==IT)?EOF:*IS++;}
template<typename _Tp>
inline IO&operator>>(_Tp&x){
_Tp s=0;char c;
if(IS+10<IT){
do c=*IS++;while(c<48||c>57);
do s=s*10+c-48,c=*IS++;while(c>47&&c<58);
}else{
do c=gc();while(c<48||c>57);
do s=s*10+c-48,c=gc();while(c>47&&c<58);
}return(x=s,*this);
}

#undef ISIZE
}io;

void Decompose(const int&x){
int v=x;
for(int i=1;i<=pcnt;++i){
int np=pri[i];
if(sqr(np)>v)break;
if(!(v%np)){
int nk=1;
for(v/=np;!(v%np);v/=np)
++nk;
insd(np,nk);
}
}if(v!=1)insd(v,1);
}

int calc(const int&L,const int&R){
/*2^{23} < 10^{7} < 2^{24}*/
static int spw[24]={1};
int np=v(L);
for(int i=1,pw=1,_=k(R);i<=_;++i)
pw=(i64)pw*np%p,spw[i]=pls(spw[i-1],pw);

int v=1;
for(int i=L;i<=R;++i)
v=(i64)v*spw[k(i)]%p;
return pls((i64)(v-1LL)*inv(np)%p*(np-1LL)%p,1);
}

int main(){
Sieve(U);
int n;
io>>n;
for(int i=1,a;i<=n;++i)
io>>a,Decompose(a);

int Ans=1;
std::sort(divsor+1,divsor+dcnt+1);
for(int L=1,R=1;R<=dcnt;++R)
if(v(R)!=v(R+1))
Ans=(i64)Ans*calc(L,R)%p,L=R+1;
return printf("%d",Ans),0;
}

「BZOJ 3561」DZY Loves Math Ⅵ

Description

给定 n,m ,求

\sum_{i=1}^{n}\sum_{j=1}^{m}\operatorname{lcm}(i, j)^{\gcd(i, j)}

答案对 10^{9}+7 取模。

Constraints

1\le n,m\le 5\times 10^{5}

Analysis

这恐怕是 DZY Loves Math 系列中最简单的一道了...

n\le m ,有:

\begin{aligned} &\sum_{i=1}^{n}\sum_{j=1}^{m}\operatorname{lcm}(i, j)^{\gcd(i, j)}\\ =&\sum_{d=1}^{n}\sum_{i=1}^{n}\sum_{j=1}^{m}\left[(i, j)=d\right]\left(\frac{ij}{d}\right)^{d}\\ =&\sum_{d=1}^{n}\sum_{i=1}^{\left\lfloor\frac{n}{d}\right\rfloor}\sum_{j=1}^{\left\lfloor\frac{m}{d}\right\rfloor}\left(ijd\right)^{d}\sum_{d'=1}^{n}\left[d' | i\right]\left[d' | j\right]\mu(d')\\ =&\sum_{d=1}^{n}d^{d}\sum_{d'=1}^{n}\mu(d')\sum_{i=1}^{\left\lfloor\frac{n}{dd'}\right\rfloor}\sum_{j=1}^{\left\lfloor\frac{m}{dd'}\right\rfloor}\left(id'jd'\right)^{d}\\ =&\sum_{d=1}^{n}d^{d}\sum_{d'=1}^{n}\mu(d')d'^{2d}\sum_{i=1}^{\left\lfloor\frac{n}{dd'}\right\rfloor}i^{d}\sum_{j=1}^{\left\lfloor\frac{m}{dd'}\right\rfloor}j^{d} \end{aligned}

对于每个 d ,处理出 i^{d} 的前缀和并统计答案即可。

时间复杂度 O(n \log{n})

Code

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
71
72
73
74
75
76
77
78
79
80
81
82
#include<cstdio>

constexpr int p=1000000007;
constexpr int U=500000;
constexpr int Maxn=500000;

template<typename _Tp>
inline void swap(_Tp&x,_Tp&y)
{_Tp z=x;x=y;y=z;}
template<typename _Tp>
inline _Tp sqr(const _Tp&v)
{return(i64)v*v%p;}

template<typename _Tp1,typename _Tp2>
inline void inc(_Tp1&x,const _Tp2&y)
{(x>=p-y)?(x+=y-p):(x+=y);}
template<typename _Tp1,typename _Tp2>
inline _Tp1 pls(const _Tp1&x,const _Tp2&y)
{return(x>=p-y)?(x-p+y):(x+y);}
template<typename _Tp1,typename _Tp2>
inline void dec(_Tp1&x,const _Tp2&y)
{x-=y;(x<0)&&(x+=p);}

template<typename _Tp,class __Tp>
inline _Tp fpow(const _Tp&x,const __Tp&t){
_Tp Ans=1,v=x;
for(_T n=t;n;n>>=1,v=(i64)v*v%p)
if(n&1)Ans=(i64)Ans*v%p;
return Ans;
}

bool notp[U+1];
int pri[U/7];
int cnt;
short mu[U+1]={0,1};

void Sieve(const int&_U){
for(int i=2;i<=_U;++i){
if(!notp[i])
pri[++cnt]=i,mu[i]=-1;
for(int j=1,v;j<=cnt;++j){
if((v=i*pri[j])>_U)break;
notp[v]=true;
if(i%pri[j])
mu[v]=-mu[i];
else break;
}
}
}

int pw[Maxn+1];
int spw[Maxn+1];

int calc(const int&n,const int&m){
for(int i=1;i<=m;++i)
pw[i]=1;

int Ans=0;
for(int d=1;d<=n;++d){
int nx=n/d,mx=m/d;

for(int i=1,v=0,cv=0;i<=mx;++i)
pw[i]=cv=(i64)pw[i]*i%p,
spw[i]=v=pls(v,cv);

int v=0;
for(int _d=1,dv;_d<=nx;++_d)
if(mu[_d])
dv=(i64)sqr(pw[_d])*spw[nx/_d]%p*spw[mx/_d]%p,
~mu[_d]?inc(v,dv):dec(v,dv);

inc(Ans,(i64)fpow(d,d)*v%p);
}return Ans;
}

int main(){
int n,m;
scanf(" %d %d",&n,&m);
if(n>m)swap(n,m);
Sieve(m);
return printf("%d",calc(n,m)),0;
}

「BZOJ 3739」DZY Loves Math Ⅷ

Description

给定 n ,求

\sum_{i=1}^{n}\sum_{j=1}^{i}\mu\left(\operatorname{lcm}(i, j)^{\gcd(i, j)}\right)

T 组询问。

Constraints

1\le n\le 10^{7},1\le T\le 10^{3}

Analysis

容易发现仅有满足 i \perp j 的数对 (i, j) 对答案有贡献:

\begin{aligned} &\sum_{i=1}^{n}\sum_{j=1}^{i}\mu\left(\operatorname{lcm}(i, j)^{\gcd(i, j)}\right)\\ =&\sum_{i=1}^{n}\sum_{j=1}^{i}\left[i \perp j\right]\mu(i)\mu\left(j\right)\\ =&\sum_{i=1}^{n}\sum_{j=1}^{i}\mu(i)\mu\left(j\right)\sum_{d |\gcd(i, j)}\mu(d)\\ =&\sum_{i=1}^{n}\mu(i)\sum_{d | i}\mu(d)\sum_{j=1}^{\left\lfloor\frac{i}{d}\right\rfloor}\mu\left(jd\right) \end{aligned}

S_{d}(n)=\sum_{i=1}^{\left\lfloor\frac{n}{d}\right\rfloor}\mu\left(id\right)

\sum_{i=1}^{n}\sum_{j=1}^{i}\mu\left(\operatorname{lcm}(i, j)^{\gcd(i, j)}\right)=\sum_{i=1}^{n}\mu(i)\sum_{d | i}\mu(d)S_{d}(i)

i 对答案有贡献当且仅当其无平方因子,将其分解质因数并搜索出所有因数更新 S_{d}(i) 即可。

时间复杂度 O\left(n+\frac{6n}{\pi^{2}}n^{\frac{1.066}{\ln\ln{n}}}+T\right)

Code

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<cstdio>
#include<algorithm>

constexpr int U=10000000;

inline int rd(){
static int _;
return scanf("%d",&_),_;
}

unsigned char notp[U+1];
int pri[U/10],
mnd[U+1];
short mu[U+1];
int cnt;

int q[1001],
ans[U+1],
S[U+1];

int dv[10];
int dcnt,del,Sdel;

inline void calc(const int&pos,const int&v){
if(pos==dcnt){
S[v]+=Sdel;
if(mu[v]<0)
del-=S[v];
else del+=S[v];
return;
}calc(pos+1,v);
calc(pos+1,v*dv[pos]);
}

int main(){
const int T=rd();
for(int i=0;i<T;++i)
q[i]=rd();

mu[1]=mnd[1]=1;
const int mxq=*std::max_element(q,q+T);
for(int i=2;i<=mxq;++i){
if(!notp[i])
pri[++cnt]=mnd[i]=i,mu[i]=-1;
for(int j=1,v;j<=cnt&&(v=i*pri[j])<=mxq;++j){
notp[v]=true;
mnd[v]=pri[j];
if(i%pri[j])
mu[v]=-mu[i];
else break;
}
}

for(int i=1;i<=mxq;++i)
if(ans[i]=ans[i-1],mu[i]){
int t=i;
for(dcnt=del=0,Sdel=mu[i];t!=1;t/=mnd[t])
dv[dcnt++]=mnd[t];
calc(0,1);

if(mu[i]<0)
ans[i]-=del;
else ans[i]+=del;
}

for(int i=0;i<T;++i)
printf("%d\n",ans[q[i]]);

return 0;
}