题面好长,看了好久的。其实要求的东西很简单,看懂题之后秒出思路了。

由于我以前封装过《区间加,维护区间gcd》的板子,所以直接贴一波板子,然后写+调只花了40min,然而喜提Tle on test 22。此后就是漫长的卡常,无果,期间一度怀疑是出题人特意卡掉了这种复杂度。然后突然想起cf上好像scanf巨慢,于是换成关同步流cin就过了。。。总共花了1h

现在看之前封装的《区间加,维护区间gcd》依托答辩!(不过正确性还是有保证的)打算近期重新封装下。

先贴个代码在这:(好像确实有不用线段树复杂度更优的做法,之后再看)

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
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
#include<bits/stdc++.h>
using namespace std;
using ll=long long;

struct ADD_GCD//要求支持:a数组区间加,查询区间gcd
{
const int n;
const vector<ll>a;//原数组a是差分数组d的前缀和

vector<ll>d;//差分数组时刻保持正确

ADD_GCD(const vector<ll> &init) : n(init.size()-1) , a(init)
{
assert(a[0]==0);
d.resize(n+1) , d[0]=0;
for(int i=1;i<=n;++i)d[i]=a[i]-a[i-1];

function<void(int,int,int)>build = [&](int rt,int l,int r)
{
if(l==r)return info[rt]={d[l]} , void();

int mid=l+r>>1;
build(rt<<1,l,mid) , build(rt<<1|1,mid+1,r);
gengxin(rt);
};
info.resize(4<<__lg(n)) , build(1,1,n);

BIT.resize(n+1) , BIT[0]=0;
for(int i=1;i<=n;++i)BIT[i]=a[i]-a[i-lowbit(i)];
}


struct INFO
{
ll g;
INFO operator+(const INFO &a)const
{
return {__gcd(g,a.g)};
}
};
static constexpr INFO E={0};//单位元
vector<INFO>info;//线段树维护差分数组d(单点赋值,查询区间gcd)
void gengxin(int rt)
{
info[rt]=info[rt<<1]+info[rt<<1|1];
}

void ddxg(int rt,int l,int r,int x,const INFO &C)
{
if(l==r)return info[rt]=C , void();

int mid=l+r>>1;
if(x<=mid)ddxg(rt<<1,l,mid,x,C);
else ddxg(rt<<1|1,mid+1,r,x,C);
gengxin(rt);
}
void ddxg(int x,ll C)//d[x]=C
{
assert(d[x]==C) , ddxg(1,1,n,x,{d[x]});
}

INFO chaxun(int rt,int l,int r,int L,int R)
{
if(L<=l && r<=R)return info[rt];
INFO res=E; int mid=l+r>>1;
if(L<=mid)res=res+chaxun(rt<<1,l,mid,L,R);
if(R>mid)res=res+chaxun(rt<<1|1,mid+1,r,L,R);
return res;
}
ll chaxun(int L,int R)//gcd(d[L],...,d[R])
{
return chaxun(1,1,n,L,R).g;
}


inline int lowbit(int x){return x&-x;}
vector<ll>BIT;//BIT维护差分数组d
void xg(int l,int r,ll del)//原数组a区间加,a[l~r]+=del
{
int x=l;
d[x]+=del , ddxg(x,d[x]);
while(x<=n)BIT[x]+=del , x+=lowbit(x);

x=r+1;if(x>n)return;
d[x]-=del , ddxg(x,d[x]);
while(x<=n)BIT[x]-=del , x+=lowbit(x);
}
ll cx(int x)//查询原数组a的单点值a[x]
{
ll res=0;
while(x)res+=BIT[x] , x-=lowbit(x);
return res;
}


void op1(int l,int r,ll del)
{
xg(l,r,del);
}
ll op2(int l,int r)
{
//return assert(l==1 && r==n) , abs(chaxun(1,n));
//这题只会查全部的gcd,卡常一下

ll al=cx(l) , res;

if(l==r)res=al;
else res=__gcd(al,chaxun(l+1,r));

return abs(res);//abs很关键,保证gcd非负!
}
};


constexpr int N=5e5+5;
struct E{int to,nex;ll w;}e[N<<1];
int head[N],tot;
void ae(int u,int v,ll w)
{
++tot , e[tot]={v,head[u],w} , head[u]=tot;
}

bool isu[N];
struct subtree_dfn
{
int mn,mx;
void merge(const subtree_dfn &a)
{
if(a.mn==-1)return assert(a.mx==-1) , void();
else
{
if(mn==-1)
{
assert(mx==-1);
mn=a.mn , mx=a.mx;
}
else
{
assert(mx!=-1);
mn=min(mn,a.mn) , mx=max(mx,a.mx);
}

assert(mn!=-1 && mx!=-1);
}
}
}dfn[N];
int dfn_table[N],js_dfn;
ll shen[N];
void get_dfn(int u,int fa)//只有关键点才++dfn
{
if(isu[u])
{
++js_dfn;
dfn[u]={js_dfn,js_dfn};
dfn_table[js_dfn]=u;
}
else dfn[u]={-1,-1};

for(int i=head[u];i;i=e[i].nex)
{
int v=e[i].to;
if(v==fa)continue;

shen[v]=shen[u]+e[i].w , get_dfn(v,u);
dfn[u].merge(dfn[v]);
}
}

ll smdis,ans=1e18;
struct OP
{
int l,r;
ll del;

bool ok(){return l<=r;}
void cao(ADD_GCD &QJ,int op=1)
{
if(op>0)QJ.op1(l,r,del);
else QJ.op1(l,r,-del);
}
};
void dfs(int u,int fa,ADD_GCD &QJ)
{
ll g=QJ.op2(1,QJ.n);
assert(smdis%g==0);
ans=min(ans,smdis/g);

for(int i=head[u];i;i=e[i].nex)
{
int v=e[i].to;
if(v==fa)continue;

int qi=dfn[v].mn , mo=dfn[v].mx;
int jsv = (qi==-1) ? 0 : mo-qi+1;
int outv = QJ.n-jsv;

ll del=e[i].w*(outv-jsv);
smdis+=del;

stack<OP>stk;
if(jsv==0)
{
OP op={1,QJ.n,+e[i].w};
assert(op.ok());
op.cao(QJ,1) , stk.push(op);
}
else
{
OP op={qi,mo,-e[i].w};
assert(op.ok());
op.cao(QJ,1) , stk.push(op);

op={1,qi-1,+e[i].w};
if(op.ok())op.cao(QJ,1) , stk.push(op);

op={mo+1,QJ.n,+e[i].w};
if(op.ok())op.cao(QJ,1) , stk.push(op);
}

dfs(v,u,QJ);

while(!stk.empty())
{
OP op=stk.top();stk.pop();
op.cao(QJ,-1);
}
smdis-=del;
}
}

int main()
{
ios::sync_with_stdio(0);
int n,k;cin>>n>>k;
if(k==1)
{
cout<<"0\n";return 0;
}

for(int i=1;i<=k;++i)
{
int u;cin>>u;
isu[u]=true;
}
for(int i=1;i<n;++i)
{
int u,v;ll w;
cin>>u>>v>>w;
ae(u,v,w) , ae(v,u,w);
}
js_dfn=0 , shen[1]=0 , get_dfn(1,0);
smdis=0;
for(int u=1;u<=n;++u)
{
if(isu[u])smdis+=shen[u];
}

vector<ll>a(k+1,0);//注意,是维护到k个点距离的gcd
assert(js_dfn==k);
for(int i=1;i<=k;++i)//枚举dfn,dfn是连续的,所以能维护
{
int u=dfn_table[i];
a[i]=shen[u];
}
ADD_GCD QJ(a);

dfs(1,0,QJ);
cout<<2*ans<<"\n";
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
#include<bits/stdc++.h>

template<class S,class OP=S(*)(S,S)>
class SegmentTree
{
private:
OP op;
public:
SegmentTree(OP _op=OP()) : op(_op) {}
S f(S a,S b)
{
return op(a,b);
}
};

int mul(int a,int b)
{
return a*b;
}
struct Max
{
int operator()(int a, int b)
{
return std::max(a,b);
}
};

void mytest()
{
std::function<int(int,int)> add = [](int a,int b){return a+b;};
SegmentTree<int,decltype(add)>t1(add);
std::cout<<t1.f(3,5)<<"\n";//8

SegmentTree<int>t2(mul);
std::cout<<t2.f(3,5)<<"\n";//15

SegmentTree<int,Max>t3;
std::cout<<t3.f(3,5)<<"\n";//5
}
int main()
{
mytest();
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
#include<bits/stdc++.h>
using namespace std;

int main()
{
int const a[]={1,2,3};

cout<<a[0]<<" "<<a[1]<<" "<<a[2]<<"\n";//output:1 2 3
//++a[0];//[Error] increment of read-only location 'a[0]'

/*正确的,但不好复现
{
const int * pp=a;
cout<<pp<<" -> "<<*pp<<"\n";//output:0x73fde0 -> 1

++pp , cout<<pp<<" -> "<<*pp<<"\n";//output:0x73fde4 -> 2
//int * mypp=pp-1;//[Error] invalid conversion from 'const int*' to 'int*' [-fpermissive]

int * ppp = (int *)0x73fde0;
cout<<ppp<<" -> "<<*ppp<<"\n";//output:0x73fde0 -> 1

++(*ppp);
cout<<ppp<<" -> "<<*ppp<<"\n";//output:0x73fde0 -> 2
cout<<a<<" -> "<<a[0]<<"\n";//output:0x73fde0 -> 2

//可以看到,a[0]确实被改变了!!!
}
*/


{
//int * const p=a;//[Error] invalid conversion from 'const int*' to 'int*' [-fpermissive]
//无法声明!

int b[]={100,200,300};
int * const p=b;
cout<<*p<<"\n";//output:100
*p=1000 , cout<<*p<<" = "<<b[0]<<"\n";//output:1000 = 1000
//++p;//[Error] increment of read-only variable 'p'

//结论:p不可以变,但可以改变p指向的东西
}


{
int const *p=a;
cout<<p<<" -> "<<*p<<"\n";//output:0x73fde0 -> 1
//*p=100;//[Error] assignment of read-only location '* p'
++p , cout<<p<<" -> "<<*p<<"\n";//output:0x73fde4 -> 2
//*p=100;//[Error] assignment of read-only location '* p'

//结论:p可以随意变,但p指向的区域永远不能被修改(即:并不只是声明时指向的区域)
}


{
int const * const p = &a[0];
cout<<p<<" -> "<<*p<<"\n";//0x73fde0 -> 1
//++p;//[Error] increment of read-only variable 'p'
//++(*p);//[Error] increment of read-only location '*(const int*)p'

//结论:p被const修饰,不可变;*p被const修饰,所以永远无法通过p改变p指向的区域
}
return 0;
}

①概述:

作用:快速求某些积性函数f(x)的前缀和

对f(x)的要求:

由于要求1,故绝大多数能用min25筛的情况下,f(p)都是关于p的简单多项式(简单多项式的前缀和可以快速求出)

比如:

复杂度:



②前置知识:

1.埃氏筛

2.线性筛

复杂度线性 if(i%p[j]==0)break;

i是p[j]的倍数时,i*p[j+1]这个合数肯定也能被p[j]筛掉,所以不用继续筛下去了

(因为i中含有p[j],而p[j]肯定小于p[j+1])



③正式开始:

part1:质数部分 (dp)

需要引入 完全积性函数 F(x),且满足质数处有F(p)=f(p) g(n,j),如何状态转移?

其实就是考虑第j轮埃氏筛,会筛掉哪些数!

g(30,1) : 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

g(30,2) : 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

g(30,3) : 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

g(n,j-1) -> g(n,j),减小的值为最小质因子为p[j]的合数的值!

设被筛掉的合数为x,有:


part2:加上合数部分的贡献,求出答案

质数处的和已经在part1求出。

容易想到,可以通过枚举最小质因子的幂次的方式,来统计合数部分的贡献(且只用到积性函数的性质就够了) 再简单解释一下:



④复杂度分析

part1:(dp模拟埃氏筛)

part2:

线



⑤代码实现

洛谷模板题:

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
#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
const ll mod=1e9+7;
ll ksm(ll a,ll b)
{
ll res=1;
while(b)
{
if(b&1)res=res*a%mod;
a=a*a%mod , b>>=1;
}
return res;
}
ll ny(ll x){return ksm(x,mod-2);}


namespace MIN25//f1(x)=x , f2(x)=x*x
{
const ll N=1e10+7 , __N=1e5+5;
ll n,__n;//筛到__n(根号n)即可
//因为n以内,没有 最小质因子>根号n 的合数

void read_n(){scanf("%lld",&n) , __n=sqrt(n);}

bool vis[__N];
ll p[__N];int jsp;//根号n以内的p的个数
ll sp1[__N/*jsp*/],sp2[__N/*jsp*/];//F在质数处的前缀和
void xxs()//线性筛预处理
{
for(int i=2;i<=__n;++i)
{
if(!vis[i])
{
p[++jsp]=i;
sp1[jsp] = (sp1[jsp-1]+i)%mod;
sp2[jsp] = (sp2[jsp-1]+1LL*i*i)%mod;
}

for(int j=1;j<=jsp && p[j]*i<=__n;++j)
{
vis[p[j]*i]=true;
if(i%p[j]==0)break;
}
}
}

ll w[__N*2/*tot*/];int tot;//所有整除分块的值(返回n/l),tot为块号
int id1[__N*2/*n/l的值*/],id2[__N*2/*相当于(n/l)的值*/];//w的逆映射,返回块号tot
inline int get_id(ll val)//val是n的整除分块的值(n/l)
{
if(val<=__n)return id1[val];
else return id2[n/val];
}
ll g1[__N*2/*tot*/],g2[__N*2/*tot*/];//初始为g(n,0) (对n的每个整除分块的值都要算g的)

void work_g()
{
for(int j=1;j<=jsp;++j)//g(n/... , j),只保存n/...这一维
{
for(int i=1;i<=tot && p[j]*p[j]<=w[i];++i)//每次要修改所有n的整除分块的g
//注意顺序,w数组的值是递减的!
{
ll tmp=w[i]/p[j];
int tmptot = get_id(tmp);

g1[i] -= p[j]/*f1*/*(g1[tmptot]-sp1[j-1]);
g1[i]%=mod;
g2[i] -= p[j]*p[j]/*f2*/%mod*(g2[tmptot]-sp2[j-1]);
g2[i]%=mod;
}
}
}
ll S(ll i,int j)//i是n/... (的值)
{
if(p[j]>=i)return 0;
int itot = get_id(i);
ll res = (g2[itot]-g1[itot]) - (sp2[j]-sp1[j]);//质数部分
res%=mod;

for(int k=j+1;p[k]*p[k]<=i && k<=jsp;++k)
{
ll pe=p[k];//p[k]的e次方
for(int e=1;pe<=i;++e,pe=pe*p[k])//pe这里先别取模!要作为循环条件!
{
ll x=pe%mod;
res += (x*(x-1)%mod)/*f(pe)*/ * (S(i/pe,k)+(e>1));
res%=mod;
}
}

return res;
}

ll work__min25()
{
read_n() , xxs();

ll inv6=ny(6);
auto calc = [inv6](ll i,int op){
i%=mod;
if(op==1)return i*(i+1)/2%mod;
return i*(i+1)%mod*(i+i+1)%mod*inv6%mod;
};

for(ll l=1,r;l<=n;l=r+1)
{
r=min(n,n/(n/l));

w[++tot]=n/l;
if(w[tot]<=__n)id1[w[tot]] = tot;
else id2[n/w[tot]] = tot;

g1[tot]=calc(w[tot],1)-1 , g2[tot]=calc(w[tot],2)-1;
//g1,g2记得减掉F(1)的值 !
}

work_g();
ll res=S(n,0)+1;
res=(res%mod+mod)%mod;
return res;
}
}

int main()
{
printf("%lld\n",MIN25::work__min25());
return 0;
}

Welcome to Hexo! This is your very first post. Check documentation for more info. If you get any problems when using Hexo, you can find the answer in troubleshooting or you can ask me on GitHub.

Quick Start

Create a new post

1
$ hexo new "My New Post"

More info: Writing

Run server

1
$ hexo server

More info: Server

Generate static files

1
$ hexo generate

More info: Generating

Deploy to remote sites

1
$ hexo deploy

More info: Deployment

0%