发表于 | 更新于

Link

Solution

考虑NN的一个约数dd会对答案产生多大的贡献，即有多少i≤N,gcd(i,N)==di\le N,gcd(i,N)==d 稍加分析就知道是ϕ(Nd)\phi(\dfrac Nd)个。 ii从11枚举到N\sqrt N，如果是约数就O(N)O(\sqrt N)地计算ϕ(i)\phi(i)和ϕ(Ni)\phi(\dfrac Ni)。

这样不会T。。从11到NN的约数个数和是∑[Ni]\sum [\dfrac Ni]，所以NN的约数个数就是∑[Ni]−∑[N−1i]\sum[\dfrac Ni]-\sum[\dfrac {N-1}i]，取值不会很大的。。 求靠谱一些的靠谱性证明。。。

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
//Code by Lucida
#include<bits/stdc++.h>
#define red(x) scanf("%d",&x)
#define llred(x) scanf("%lld",&x)
//#define debug(x) std::cout<<#x<<'='<<x<<std::endl
template <class T> inline bool chkmx(T &a,const T &b){return a<b?a=b,1:0;}
template <class T> inline bool chkmn(T &a,const T &b){return a>b?a=b,1:0;}
typedef long long LL;
const int N=65537,MAXN=N+10;
int phi[MAXN],prime[MAXN],pcnt;
void Prep()
{
	phi[1]=1;
	static bool notp[MAXN];
	for(int i=2;i<=N;i++)
	{
		if(!notp[i])
		{
			prime[++pcnt]=i;
			phi[i]=i-1;
		}
		for(int j=1;j<=pcnt && i*prime[j]<=N;j++)//保证prime[j]是最小的质因子
		{
			notp[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]*phi[prime[j]];
		}
	}
}
LL Phi(LL x)
{
	if(x<=N) return phi[x];
	LL mid=(LL)sqrt(x+0.5);
	LL res=x;
	for(int i=1;i<=pcnt && prime[i]<=x;i++) 
		if(x%prime[i]==0)
		{
			res=res*(prime[i]-1)/prime[i];
			while(x%prime[i]==0) x/=prime[i];
		}
	if(x>1) res=res*(x-1)/x;
	return res;
}
int main()
{
	//freopen("input","r",stdin);
	Prep();
	LL n;llred(n);
	LL ANS=0,mid=(LL)sqrt(n+0.5);
	for(int i=1;i<=mid;i++)
		if(n%i==0)
			ANS+=Phi(n/i)*i,ANS+=Phi(i)*(n/i);
	printf("%lld\n",ANS);
	return 0;
}

发表于 | 更新于

Link

Solution

结论很神。。但知道结论了也没推出来。。

∑i=1Nd(i)=∑i=1N[Ni]\sum_{i=1}^N d(i)=\sum_{i=1}^N[\frac Ni]

(AHOI2005约数探究) 这个是把每个数对答案的贡献的累计得出的结论。 继续这个思路，有可能猜到

∑i=1N∑j=1Md(ij)=∑i=1N∑j=1M[Ni][Mj]\sum_{i=1}^N\sum_{j=1}^Md(ij)=\sum_{i=1}^N\sum_{j=1}^M[\frac Ni][\frac Mj]

也就是按照统计对答案的贡献这个思路，只不过多用个乘法原理。 然后再完善一下就是最后的结论

∑i=1N∑j=1Md(ij)=∑i=1N∑j=1M[Ni][Mj][gcd(i,j)==1]\sum_{i=1}^N\sum_{j=1}^Md(ij)=\sum_{i=1}^N\sum_{j=1}^M[\frac Ni][\frac Mj][\gcd(i,j)==1]

(说的好像真的能做出来似的) 然后开始反演

F(n)=∑i=1N∑j=1M∑n∣gcd(i,j)[Ni][Mj]=∑i=1[Nn]∑j=1[Mn][Nin][Mjn]=∑n∣df(d)F(n)=\sum_{i=1}^N\sum_{j=1}^M\sum_{n|\gcd(i,j)}[\frac Ni][\frac Mj]=\sum_{i=1}^{[\frac Nn]}\sum_{j=1}^{[\frac Mn]}[\frac N{in}][\frac M{jn}]=\sum_{n|d}f(d)

第三个式子的i,ji,j枚举的是倍数。。

f(n)=∑n∣dμ(dn)F(d)f(n)=\sum_{n|d}\mu(\frac dn)F(d)Ans=f(1)=∑dμ(d)∑i=1[Nd]∑j=1[Md][Nid][Mjd]Ans=f(1)=\sum_{d}\mu(d)\sum_{i=1}^{[\frac Nd]}\sum_{j=1}^{[\frac Md]}[\frac N{id}][\frac M{jd}]

下一步！！ 注意到右边的那个东西是两个关于[Nd][\dfrac Nd]和[Md][\dfrac Md]的函数！！于是把它们提取出来处理一个前缀和就可以加速了！

Tips

所以，要分块加速，不一定要把相应项的分母换成同一个，还可以根据求和指标等等方面判断是否能化成一个函数的两个取值。

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
//Code by Lucida
#include<bits/stdc++.h>
#define get(x) scanf("%d",&x)
//#define debug(x) std::cout<<#x<<'='<<x<<std::endl
template <class T> inline bool chkmx(T &a,const T &b){return a<b?a=b,1:0;}
template <class T> inline bool chkmn(T &a,const T &b){return a>b?a=b,1:0;}
const int N=50000,MAXN=N+10;
typedef long long LL;
using std::swap;
using std::min;
int prime[MAXN],pcnt,mu[MAXN],Smu[MAXN],g[MAXN];
LL Sg[MAXN];
void Sift()
{
	static bool notp[MAXN];
	mu[1]=1;
	for(int i=2;i<=N;i++)
	{
		if(!notp[i])
		{
			prime[++pcnt]=i;
			mu[i]=-1;
		}
		for(int j=1;j<=pcnt && prime[j]*i<=N;j++)
		{
			notp[i*prime[j]]=1;
			if(i%prime[j]==0)
			{
				mu[i*prime[j]]=0;
				break;
			}
			else
				mu[i*prime[j]]=-mu[i];
		}
	}
	for(int i=1;i<=N;i++)
		for(int j=i;j<=N;j+=i)
			g[j]++;
	for(int i=1;i<=N;i++)
	{
		Smu[i]=Smu[i-1]+mu[i];
		Sg[i]=Sg[i-1]+g[i];
	}
}
LL Calc(int n,int m)
{
	if(n>m) swap(n,m);
	LL res=0;
	for(int i=1,last;i<=n;i=last+1)
	{
		last=min(n/(n/i),m/(m/i));
		res+=(LL)(Smu[last]-Smu[i-1])*Sg[n/i]*Sg[m/i];
	}
	return res;
}
int main()
{
	Sift();
	freopen("input","r",stdin);
	int T;get(T);
	while(T--)
	{
		int n,m;get(n),get(m);
		printf("%lld\n",Calc(n,m));
	}
	return 0;
}
/* AC Record(Bugs) *
 * 
 * * * * * * * * * */

发表于 | 更新于

Link

Solution

多条直径相交的形态一定是 类似树网的核，求出直径之后，从直径的一端开始走。记一个点pp到直径左、右端点的距离为dlp,drpdl_p,dr_p,设当前走到的边的左右两点分别为pl,prp_l,p_r，如果从plp_l不走直径的边能走出一条长度为drpldr_{p_l}的路径，或者prp_r能走出dlprdl_{p_r}那这个边就不是直径的公共边。 走到第一条公共边开始累加答案，到第一条非公共边break，就可以解出来了。

Tips

树的直径的性质分析用到了很多反证法。有些结论正向思考无法下手，不妨试试能不能把自己举出的反例否掉，以尝试得出反证的思路。

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
//Code by Lucida
#include<bits/stdc++.h>
#define get(x) scanf("%d",&x)
//#define debug(x) std::cout<<#x<<'='<<x<<std::endl
template <class T> inline bool chkmx(T &a,const T &b){return a<b?a=b,1:0;}
template <class T> inline bool chkmn(T &a,const T &b){return a>b?a=b,1:0;}
const int MAXN=2e5+10;
using std::reverse;
typedef long long LL;
struct Edge{int to,pre,v;}e[MAXN<<1];int id[MAXN],ec=1,n;
void Adde(int f,int t,int v)
{
	e[++ec].to=t;e[ec].pre=id[f];id[f]=ec;e[ec].v=v;
	e[++ec].to=f;e[ec].pre=id[t];id[t]=ec;e[ec].v=v;
}
int trd[MAXN],pree[MAXN],us[MAXN],ut,dct;
LL d[MAXN],dlen[MAXN],mxd[MAXN],dpf[MAXN],dsf[MAXN];
int BFS(int s)
{
	static int Q[MAXN],he,ta,res;
	Q[he=ta=1]=s;us[s]=ut;d[s]=0,res=s,pree[s]=0;
	while(he<=ta)
	{
		int cur=Q[he++];
		for(int i=id[cur];i;i=e[i].pre)
		{
			int u=e[i].to;
			if(us[u]!=ut)
			{
				us[u]=ut;
				d[u]=d[cur]+e[i].v;
				pree[u]=i;
				if(d[u]>d[res]) res=u;
				Q[++ta]=u;
			}
		}
	}
	return res;
}
void getd()
{
	++ut;int st=BFS(1);
	++ut;int nd=BFS(st);
	++ut;
	for(int i=pree[nd];i;i=pree[e[i^1].to])
	{
		trd[++dct]=e[i].to;
		dlen[dct]=e[i].v;
		us[trd[dct]]=ut;
	}
	reverse(trd+1,trd+1+dct);
	reverse(dlen+1,dlen+1+dct);
	for(int i=1;i<=dct;i++)
	{
		mxd[i]=d[BFS(trd[i])];
		dpf[i]=dpf[i-1]+dlen[i];
	}
	for(int i=dct;i;i--) dsf[i]=dsf[i+1]+dlen[i+1];
}
int main()
{
	get(n);
    for(int i=1;i<=n-1;i++)
    {
		int a,b,c;get(a),get(b),get(c);
		Adde(a,b,c);
    }
	getd();
	int ANS=0,p;LL mind;
	for(p=2;p<=dct;p++)
	{
		if(dpf[p]!=mxd[p] && dsf[p-1]!=mxd[p-1])	
			break;
	}
	for(;p<=dct;p++)
	{
		if(dpf[p]==mxd[p] || dsf[p-1]==mxd[p-1]) break;
		else ANS++;
	}
	printf("%lld\n%d\n",dpf[dct],ANS);
	return 0;
}
/* AC Record(Bugs) *
 * WA 多了1？什么情况 [BFS()]..naive
 * * * * * * * * * */

发表于 | 更新于

Link

Solution

最小割+最小割可行边判断+退流

Tips

被坑。 论Dinic与ISAP的效率。 学网络流算法的时候听说ISAP更高效，就直接没看Dinic。做的几道网络流题也全都没有因为算法被卡过。 但是，今天被卡飞了。 ISAP写完，TLE。本地测试，最大的一个点需要4s。卡常数，最大的一个点还是3s。于是各种办法都没用。 发现网上过掉的程序几乎都是Dinic。学了Dinic。写出来调对之后很快就通过了。 以下纯属一个蒟蒻的自言自语 d[S]<N || num[x]==0。这是ISAP的终止条件。当对网络流进行局部调整的时候，如果两个点之间不联通，可能尝试绕出很多条增广路，才会判定出它们不联通。 加上计数器，计算一次ISAP不联通时的迭代次数，发现单次ISAP，1400个点的图最多居然出现了15W次无效的循环才结束。 而Dinic就每次O(n)O(n)一遍BFS就保证了在不联通时就能及时结束。 So，整体一遍的用ISAP，对图进行局部调整的用Dinic。 欢迎拍砖 欢迎吐槽

Code

ISAP

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
//Code by Lucida
#include<bits/stdc++.h>
#define get(x) scanf("%d",&x)
//#define debug(x) std::cout<<#x<<'='<<x<<std::endl
template <class T> inline bool chkmx(T &a,const T &b){return a<b?a=b,1:0;}
template <class T> inline bool chkmn(T &a,const T &b){return a>b?a=b,1:0;}
const int MAXN=(700+10)<<1,MAXM=MAXN*MAXN,INF=0x1f1f1f1f;
using std::sort;
struct Edge{int to,cap,v,pre;}e[MAXM<<1];int id[MAXN],S,T,N,ec;
int Adde(int f,int t,int cap)
{
	e[++ec].to=t;e[ec].cap=cap;e[ec].v=0;e[ec].pre=id[f];id[f]=ec;
	e[++ec].to=f;e[ec].cap=0;e[ec].v=0;e[ec].pre=id[t];id[t]=ec;
	return ec-1;
}
int pree[MAXN],d[MAXN],num[MAXN],now[MAXN];
int Aug(int T)
{
	int minf=INF;
	for(int i=pree[T];i;i=pree[e[i^1].to])
		chkmn(minf,e[i].cap-e[i].v);
	for(int i=pree[T];i;i=pree[e[i^1].to])
	{
		e[i].v+=minf;
		e[i^1].v-=minf;
	}
	return minf;
}
bool BFS(int S,int T)
{
	static int Q[MAXN],he,ta;
	memset(d,-1,sizeof(d));d[T]=0;
	Q[he=ta=1]=T;
	while(he<=ta)
	{
		int cur=Q[he++];
		for(int i=id[cur];i;i=e[i].pre)
		{
			int u=e[i].to;
			if(e[i^1].cap>e[i^1].v && d[u]==-1)
			{
				d[u]=d[cur]+1;
				Q[++ta]=u;
			}
		}
	}
	return d[S]!=-1;
}
int ISAP(int S=S,int T=T)
{
	if(!BFS(S,T)) return 0;
	memset(num,0,sizeof(num));
	memset(pree,0,sizeof(pree));
	memcpy(now,id,sizeof(id));
	for(int i=0;i<=N;i++) num[d[i]]++;
	int flow=0,cur=S;
	while(d[S]<N)
	{
		if(cur==T)
		{
			flow+=Aug(T);
			cur=S;
		}
		bool adv=0;
		for(int i=now[cur];i;i=e[i].pre)
		{
			int u=e[i].to;
			if(e[i].cap>e[i].v && d[u]+1==d[cur])
			{
				adv=1;
				now[cur]=i;
				pree[cur=u]=i;
				break;
			}
		}
		if(!adv)
		{
			if(--num[d[cur]]==0) break;
			int mind=N-1;
			for(int i=(now[cur]=id[cur]);i;i=e[i].pre)
				if(e[i].cap>e[i].v) chkmn(mind,d[e[i].to]);
			++num[d[cur]=mind+1];
			if(cur!=S) cur=e[pree[cur]^1].to;
		}
	}
	return flow;
}
void CLEAR(){ec=1;memset(id,0,sizeof(id));}
int c[MAXN];
bool cmp(int i,int j){return c[i]<c[j];}
void WORK()
{
	static int a[MAXN],b[MAXN],f[MAXN],d[MAXN],s[MAXN],slen,link[MAXN];
	#define node(x,y) (x+y*n)
	int n;get(n);S=0,T=node(n,1)+1,N=T+1;
	for(int i=1;i<=n;i++) get(a[i]);
	for(int i=1;i<=n;i++) get(b[i]);
	for(int i=1;i<=n;i++) get(c[i]),d[i]=i;
	sort(d+1,d+1+n,cmp);
	int maxlen=0;
	for(int i=1;i<=n;i++)
	{
		f[i]=1;
		for(int j=1;j<i;j++)
			if(a[j]<a[i]) chkmx(f[i],f[j]+1);
		chkmx(maxlen,f[i]);
	}
	for(int i=1;i<=n;i++)
		if(f[i]==1) Adde(S,node(i,0),INF);
		else if(f[i]==maxlen) Adde(node(i,1),T,INF);
	for(int i=1;i<=n;i++)
		for(int j=1;j<i;j++)
			if(a[j]<a[i] && f[j]+1==f[i])
				Adde(node(j,1),node(i,0),INF);
	for(int i=1;i<=n;i++) link[i]=Adde(node(i,0),node(i,1),b[i]);
	printf("%d ",ISAP());slen=0;
	for(int i=1,j;i<=n;i++)
	{
		if(BFS(node(d[i],0),node(d[i],1))) continue;
		s[++slen]=d[i];
		ISAP(node(d[i],0),S);
		ISAP(T,node(d[i],1));
		j=link[d[i]];
		e[j].cap=e[j].v=e[j^1].cap=e[j^1].v=0;
	}
	printf("%d\n",slen);
	sort(s+1,s+1+slen);
	for(int i=1;i<=slen;i++)
		printf("%d%s",s[i],i==slen?"\n":" ");
}
int main()
{
	int T;get(T);
	while(T--) CLEAR(),WORK();
	return 0;
}

Dinic

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
//Code by Lucida
#include<bits/stdc++.h>
#define get(x) scanf("%d",&x)
//#define debug(x) std::cout<<#x<<'='<<x<<std::endl
template <class T> inline bool chkmx(T &a,const T &b){return a<b?a=b,1:0;}
template <class T> inline bool chkmn(T &a,const T &b){return a>b?a=b,1:0;}
const int MAXN=(700+10)<<1,MAXM=MAXN*MAXN,INF=0x1f1f1f1f;
using std::sort;
using std::min;

struct Edge{int to,cap,v,pre;}e[MAXM<<1];int id[MAXN],S,T,N,ec;
int Adde(int f,int t,int cap)
{
	e[++ec].to=t;e[ec].cap=cap;e[ec].v=0;e[ec].pre=id[f];id[f]=ec;
	e[++ec].to=f;e[ec].cap=0;e[ec].v=0;e[ec].pre=id[t];id[t]=ec;
	return ec-1;
}
int pree[MAXN],d[MAXN],now[MAXN];
bool BFS(int S,int T)
{
	static int Q[MAXN],he,ta;
	memset(d,0,sizeof(d));
	Q[he=ta=1]=S;d[S]=1;
	while(he<=ta)
	{
		int cur=Q[he++];
		for(int i=id[cur];i;i=e[i].pre)
		{
			int u=e[i].to;
			if(e[i].cap>e[i].v && !d[u])
			{
				d[u]=d[cur]+1;Q[++ta]=u;
				if(u==T) return 1;
			}
		}
	}
	return 0;
}
int DFS(int pos,int t,int cmf)
{
	if(pos==t || cmf==0) return cmf;
	int flow=0,minf;
	for(int &i=now[pos];i;i=e[i].pre)
	{
		int u=e[i].to;
		if(d[u]==d[pos]+1 && (minf=DFS(u,t,min(cmf,e[i].cap-e[i].v))))
		{
			e[i].v+=minf;e[i^1].v-=minf;
			flow+=minf;cmf-=minf;
			if(!cmf) break;
		}
	}
	return flow;
}
int Dinic(int S,int T)
{
	int flow=0;
	while(BFS(S,T))
	{
		memcpy(now,id,sizeof(id));
		flow+=DFS(S,T,INF);
	}
	return flow;
}
void CLEAR(){ec=1;memset(id,0,sizeof(id));}
int c[MAXN];
bool cmp(int i,int j){return c[i]<c[j];}
void WORK()
{
	static int a[MAXN],b[MAXN],f[MAXN],d[MAXN],s[MAXN],slen,link[MAXN];
	#define node(x,y) (x+y*n)
	int n;get(n);S=0,T=node(n,1)+1,N=T+1;
	for(int i=1;i<=n;i++) get(a[i]);
	for(int i=1;i<=n;i++) get(b[i]);
	for(int i=1;i<=n;i++) get(c[i]),d[i]=i;
	sort(d+1,d+1+n,cmp);
	int maxlen=0;
	for(int i=1;i<=n;i++)
	{
		f[i]=1;
		for(int j=1;j<i;j++)
			if(a[j]<a[i]) chkmx(f[i],f[j]+1);
		chkmx(maxlen,f[i]);
	}
	for(int i=1;i<=n;i++)
		if(f[i]==1) Adde(S,node(i,0),INF);
		else if(f[i]==maxlen) Adde(node(i,1),T,INF);
	for(int i=1;i<=n;i++)
		for(int j=1;j<i;j++)
			if(a[j]<a[i] && f[j]+1==f[i])
				Adde(node(j,1),node(i,0),INF);
	for(int i=1;i<=n;i++) link[i]=Adde(node(i,0),node(i,1),b[i]);
	printf("%d ",Dinic(S,T));slen=0;
	for(int i=1,j;i<=n;i++)
	{
		if(BFS(node(d[i],0),node(d[i],1))) continue;
		s[++slen]=d[i];
		Dinic(node(d[i],0),S);
		Dinic(T,node(d[i],1));
		j=link[d[i]];
		e[j].cap=e[j].v=e[j^1].cap=e[j^1].v=0;
	}
	printf("%d\n",slen);
	sort(s+1,s+1+slen);
	for(int i=1;i<=slen;i++)
		printf("%d%s",s[i],i==slen?"\n":" ");
}
int main()
{
//	freopen("input","r",stdin);
	int T;get(T);
	while(T--) CLEAR(),WORK();
	return 0;
}

发表于 | 更新于

Link And This

求

∑i<=N∑j<=Mlcm(i,j)\sum_{i<=N}\sum_{j<=M}lcm(i,j)

Solution

Ans=∑i<=N∑j<=Mlcm(i,j)=∑i<=N∑j<=Mijgcd(i,j)Ans=\sum_{i<=N}\sum_{j<=M}lcm(i,j)=\sum_{i<=N}\sum_{j<=M}\dfrac {ij}{\gcd(i,j)}

改成枚举d=gcd(i,j)d=\gcd(i,j)

Ans=∑d∑i<=N∑j<=Mijd[gcd(i,j)==d]=∑d∑i<=N∑j<=Mijd[gcd([id],[jd])==1]Ans=\sum_d\sum_{i<=N}\sum_{j<=M}\dfrac {ij}d[\gcd(i,j)==d]=\sum_d\sum_{i<=N}\sum_{j<=M}\dfrac {ij}d[\gcd([\dfrac id],[\dfrac jd])==1]

令x=[id],y=[jd]x=[\dfrac id],y=[\dfrac jd]，则

Ans=∑d∑x<=[Nd]∑j<=[Md]xyd[gcd(x,y)==1]=∑dd{∑x<=[Nd]∑j<=[Md]xy[gcd(x,y)==1]}Ans=\sum_d\sum_{x<=[\frac Nd]}\sum_{j<=[\frac Md]}xyd[\gcd(x,y)==1]=\sum_dd\{\sum_{x<=[\frac Nd]}\sum_{j<=[\frac Md]}xy[\gcd(x,y)==1]\}

好了，现在的问题是：求

∑i<=a∑j<=bij[gcd(i,j)==1]\sum_{i<=a}\sum_{j<=b}ij[\gcd(i,j)==1]

定义F(n)=∑i<=a∑j<=bij[n∣gcd(i,j)],f(n)=∑i<=a∑j<=bij[gcd(i,j)==n]F(n)=\sum_{i<=a}\sum_{j<=b}ij[n|\gcd(i,j)],f(n)=\sum_{i<=a}\sum_{j<=b}ij[\gcd(i,j)==n]，则

F(n)=∑n∣df(d)F(n)=\sum_{n|d}f(d)f(n)=∑n∣dμ(dn)F(d)f(n)=\sum_{n|d}\mu(\dfrac dn)F(d)

要求的是

f(1)=∑dμ(d)F(d)f(1)=\sum_{d}\mu(d)F(d)

考虑求F(n)F(n)，

TeX parse error: Undefined control sequence \(

好了，

TeX parse error: Undefined control sequence \(

好了,

Ans=∑pp{∑x<=[Np]∑j<=[Mp]xy[gcd(x,y)==1]}Ans=\sum_pp\{\sum_{x<=[\frac Np]}\sum_{j<=[\frac Mp]}xy[\gcd(x,y)==1]\}

TeX parse error: Undefined control sequence \(

令T=pdT=pd，则

TeX parse error: Undefined control sequence \(

为什么我想推O(n)O(n)的却直接推出了O(n)O(\sqrt n)的 不过好歹算是自己推出来了，，虽然用了半个小时。。主要是变量重名把自己绕晕了，就重新推。。推反演式子的时候多用几个字母。。重名晕了就GG了。。 前面分块，后面前缀和。。 而后面的前缀和是一个积性函数。

积性函数

一个数论函数ff，如果在gcd(x,y)==1\gcd(x,y)==1的时候有f(xy)==f(x)f(y)f(xy)==f(x)f(y)，则称为积性函数。

  1. 积性函数的积是积性函数，这个很显然。
  2. 积性函数的约数和是积性函数 对于积性函数ff，定义g(n)=∑d∣nf(d)g(n)=\sum_{d|n} f(d) 那么g(x)g(y)=∑d1∣xf(d1)∑d2∣yf(d2)g(x)g(y)=\sum_{d_1|x}f(d_1)\sum_{d_2|y}f(d_2) 因为gcd(x,y)==1gcd(x,y)==1，所以∑d1∣xd1∑d2∣xd2\sum_{d_1|x}d_1\sum_{d_2|x}d_2依次相乘展开就是xyxy的所有约数，所以∑d1∣xf(d1)∑d2∣xf(d2)\sum_{d_1|x}f(d_1)\sum_{d_2|x}f(d_2)依次相乘展开的和就是∑d∣xyf(d)=g(xy)\sum_{d|xy}f(d)=g(xy)

有了第二个结论之后，h(T)=T∑d∣Tμ(d)dh(T)=T\sum_{d|T}\mu(d)d就是个积性函数了。因为μ\mu是个积性函数，而idid也是个积性函数(id(n)=nid(n)=n)。 积性函数就可以用线性筛筛了 

1
2
3
4
5
6
7
8
...
if(!notp[i])  prime[++pcnt]=i;
for(int j=1;j<=pcnt && prime[j]*i<=N;j++)
{
    notp[prime[j]*i]=1;
    if(i%prime[j]==0) break;
}
...

如果 ，说明h(i×primej)==h(i)h(primej)h(i\times prime_j)==h(i)h(prime_j) 否则，对于这个函数，i×primeji\times prime_j比ii多出来的因数一定都是primej2×xprime_j^2\times x的形式，它们的μ\mu都是00，所以不会对答案产生贡献。所以对答案的影响只是前面的系数ii变成了i×primeji\times prime_j，即h(i×primej)=i×primej∑d∣iμ(d)dh(i\times prime_j)=i\times prime_j\sum_{d|i}\mu(d)d

Code

BZOJ 2154

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
//Code by Lucida
#include<bits/stdc++.h>
#define red(x) scanf("%d",&x)
//#define debug(x) std::cout<<#x<<'='<<x<<std::endl
template <class T> inline bool chkmx(T &a,const T &b){return a<b?a=b,1:0;}
template <class T> inline bool chkmn(T &a,const T &b){return a>b?a=b,1:0;}
const int P=20101009,N=1e7,MAXN=N+10;
using std::swap;
using std::min;

typedef long long LL;
int prime[MAXN],pcnt,mu[MAXN];
LL Smud2[MAXN];
void Shake(int N)
{
	static bool notp[MAXN];
	mu[1]=1;
	for(int i=2;i<=N;i++)
	{
		if(!notp[i])
		{
			prime[++pcnt]=i;
			mu[i]=-1;
		}
		for(int j=1;j<=pcnt && (LL)prime[j]*i<=N;j++)
		{
			notp[prime[j]*i]=1;
			if(i%prime[j]==0)
			{
				mu[i*prime[j]]=0;
				break;
			}
			else
				mu[i*prime[j]]=-mu[i];
		}
	}
	//Smud2[1]=mu[1]*1*1;
	Smud2[0]=0;
	for(int i=1;i<=N;i++)
		Smud2[i]=(Smud2[i-1]+(LL)mu[i]*i*i%P)%P;
}
inline LL S(LL x,LL y){return ( (x*(x+1)/2%P) * (y*(y+1)/2%P) )%P;}
inline LL SMUD2(int l,int r)
{
	LL res=Smud2[r]-Smud2[l-1];
	((res%=P)+=P)%=P;
	return res;
}
LL F(int n,int m)
{
	if(n>m) swap(n,m);
	LL res=0;
	for(int i=1,last;i<=n;i=last+1)
	{
		last=min(n/(n/i),m/(m/i));
		(res+=(SMUD2(i,last)*S(n/i,m/i))%P)%=P;
	}
	return res;
}
inline LL Sum(LL l,LL r){return ((l+r)*(r-l+1)/2)%P;}
int main()
{
	int n,m;red(n),red(m);
	if(n>m) swap(n,m);
	Shake(n);
	LL ANS=0;
	for(int i=1,last;i<=n;i=last+1)//last=i+1..
	{
		last=min(n/(n/i),m/(m/i));
		(ANS+=(Sum(i,last)*F(n/i,m/i))%P)%=P;
	}
	((ANS%=P)+=P)%=P;
	printf("%lld\n",ANS);
	return 0;
}

BZOJ 2693

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
//Code by Lucida
#include<bits/stdc++.h>
#define red(x) scanf("%d",&x)
//#define debug(x) std::cout<<#x<<'='<<x<<std::endl
template <class T> inline bool chkmx(T &a,const T &b){return a<b?a=b,1:0;}
template <class T> inline bool chkmn(T &a,const T &b){return a>b?a=b,1:0;}
const int N=10000000,MAXN=N+10,P=1e8+9,MAXT=10000+10;
using std::min;
using std::swap;
using std::max;
typedef long long LL;
int prime[MAXN],pcnt,mu[MAXN];
LL H[MAXN],SH[MAXN];
void Shake(int N)
{
	static bool notp[MAXN];
	mu[1]=H[1]=1;
	for(int i=2;i<=N;i++)
	{
		if(!notp[i])
		{
			prime[++pcnt]=i;
			mu[i]=-1;
			H[i]=(((LL)i*(1*mu[1]+i*mu[i])%P)+P)%P;
		}
		for(int j=1;j<=pcnt && (LL)prime[j]*i<=N;j++)
		{
			notp[prime[j]*i]=1;
			if(i%prime[j]==0)
			{
				mu[i*prime[j]]=0;
				H[i*prime[j]]=H[i]*prime[j]%P;
				break;
			}
			else
			{
				mu[i*prime[j]]=-mu[i];
				H[i*prime[j]]=H[i]*H[prime[j]]%P;
			}
		}
	}
	/*for(int i=1;i<=N;i++)
		for(int j=i;j<=N;j+=i)
			H[j]+=((LL)j*mu[i]*i)%P,(H[j]+=P)%=P;
	*/
	SH[0]=0;
	for(int i=1;i<=N;i++)
		SH[i]=(SH[i-1]+H[i])%P;
}
LL S(LL x,LL y)
{
	return ( ((x*(x+1)/2)%P)*((y*(y+1)/2)%P) )%P;
}
LL HSum(LL l,LL r)
{
	LL res=SH[r]-SH[l-1];
	((res%=P)+=P)%=P;
	return res;
}
LL Calc(LL n,LL m)
{
	LL res=0;
	if(n>m) swap(n,m);
	for(int i=1,last;i<=n;i=last+1)
	{
		last=min(n/(n/i),m/(m/i));
		(res+=HSum(i,last)*S(n/i,m/i))%=P;
	}
	return res;
}
int main()
{
	//freopen("input","r",stdin);
	static int Q[MAXT][2];
	int T;red(T);
	int n=0;
	for(int i=1;i<=T;i++)
	{
		red(Q[i][0]),red(Q[i][1]);
		chkmx(n,max(Q[i][0],Q[i][1]));
	}
	Shake(n);
	for(int i=1;i<=T;i++)
		printf("%lld\n",Calc(Q[i][0],Q[i][1]));
	return 0;
}

发表于 | 更新于

Link

想到了b=0是BSGS，然后化式子，结果化错了，以为不会，就弃疗了。。

Solution

化出式子， 同乘a−1a-1，BSGS，真好哈。 然后写，写完WA。本地测90。 想不透。研究数据无果。事实证明我太naive。 这个xix_i每次都要 ，怎么能把分母带着呢？？ 于是把1a−1\dfrac 1{a-1}换成inv(a−1)inv(a-1)，重新划划式子，改上，WA。 无奈，又查了查，发现BSGS写错了。mid=Pmid=\sqrt P，应该取上整，这样imid−jimid-j才能取遍[0,P][0,P]。之前写BSGS从没出过这种问题，说明还算是比较隐蔽的。。要好好记着。

Tips

思考题还是要慢点，不要急躁。 用高中数学方法化出来的式子不一定能直接套到模型里，尤其是很可能不适用于在同余系。

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
//Code by Lucida
#include<bits/stdc++.h>
#define get(x) scanf("%d",&x)
//#define debug(x) std::cout<<#x<<'='<<x<<std::endl
template <class T> inline bool chkmx(T &a,const T &b){return a<b?a=b,1:0;}
template <class T> inline bool chkmn(T &a,const T &b){return a>b?a=b,1:0;}
typedef long long LL;
int P;
int add(int a,int b){return ((LL)a+b)%P;}
int mul(int a,int b){return (LL)a*b%P;}
int Pow(int x,int v)
{
	int res=1;
	while(v)
	{
		if(v&1)
			res=mul(res,x);
		x=mul(x,x);
		v>>=1;
	}
	return res;
}
int inv(int x)
{
	return Pow(x,P-2);
}
int exgcd(int a,int b,LL &x,LL &y)
{
	if(!b)
	{
		x=1,y=0;
		return a;
	}
	else
	{
		int d=exgcd(b,a%b,x,y);
		LL nx=y,ny=x-(a/b)*y;
		x=nx,y=ny;
		return d;
	}
}
LL ExGCD(int a,int b,int c)//ax+by=c
{
	LL x,y;
	int d=exgcd(a,b,x,y);
	if(c%d) return -2;
	x*=c/d;int k=b/d;
	((x%=k)+=k)%=k;
	return x;
}

namespace iHash
{
	const int SIZE=1e6,P=1e5+7;
	struct Node
	{
		int val,key;Node *pre;
		Node(int val,int key,Node* pre):val(val),key(key),pre(pre){}
	}*POOL=(Node*)malloc(SIZE*sizeof(Node)),*ME=POOL,*head[P];
	Node *newNode(int val,int key,Node *pre){return new(ME++) Node(val,key,pre);}
	void Insert(int val,int key)
	{
		int hs=key%P;
		head[hs]=newNode(val,key,head[hs]);
	}
	int Find(int key)
	{
		int hs=key%P;
		for(Node *ptr=head[hs];ptr;ptr=ptr->pre)
			if(ptr->key==key) return ptr->val;
		return -1;
	}
	void Clear()
	{
		ME=POOL;memset(head,0,sizeof(head));
	}
}
LL BSGS(int a,int b,int P)//return -2!!
{
	if(a%P==0) return -2;
	iHash::Clear();
	int mid=(int)sqrt(P)+1;
	int aib=b;
	for(int i=0;i<=mid;i++)//val key
	{
		if(iHash::Find(aib)==-1)
			iHash::Insert(i,aib);
		aib=mul(aib,a);
	}
	int asp=Pow(a,mid),akj=asp;
	for(int i,j=1;j<=mid;j++)
	{
		if((i=iHash::Find(akj))!=-1)
			return mid*j-i;
		akj=mul(akj,asp);
	}
	return -2;
}

int main()
{
	//freopen("input","r",stdin);
	int T;get(T);int a,b,X1,t;
	while(T--)
	{
		get(P),get(a),get(b),get(X1),get(t);
		if(X1==t) puts("1");
		else if(a==0) printf("%d\n",b==t?2:-1);
		else if(a==1) printf("%lld\n",ExGCD(b,P,t-X1)+1);
		else
		{
			/*int x=X1;
			for(int i=1;i<=100;i++)
				printf("f(%d)=%d\n",i,x=add(mul(x,a),b));*/
			int binvam1=mul(b,inv(a-1)),y=add(t,binvam1),k=add(X1,binvam1);
			printf("%lld\n",BSGS(a,mul(y,inv(k)),P)+1);
		}
	}
	return 0;
}
/* AC Record(Bugs) *
 * WA Naive.
 * * * * * * * * * */

发表于 | 更新于

Link

Solution

良心出题人。。写一个LCT再写一个裸主席树再写一个暴力就给85pts。。 正解启发式。启发式真是好，好写又不慢。。 如果没有Link操作，直接树上主席树（见CLJ论文）。 对于Link操作，采用把小树合到大树里面去的方法。小树所有节点的主席树都暴力重建。 查询的时候要求LCA。因为是动态地连接，只能用倍增（可以不用吗？？）。因此在之前的DFS中搜到每个点都预处理一下fa[][]fa[][]数组。 对于一个点，重建一次的代价是O(logn)O(\log n)(主席树+倍增)。一个点至多被重构O(logn)O(\log n)次：因为每次重构，这个点所在的子树的大小都至少扩大了一倍，所以最多logn\log n次操作就能这个点与其它所有的点连成一棵树了。类似并查集的按秩合并。 一开始WA了，然后RE了。因为fa[][]fa[][]数组没有彻底清零，写成了 

1
if(!(fa[pos][i]=fa[fa[pos][i-1]][i-1])) break;

Tips

倍增LCA数组重建之后一定要彻底清零。。

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
//Code by Lucida
#include<bits/stdc++.h>
#define get(x) scanf("%d",&x)
//#define debug(x) std::cout<<#x<<'='<<x<<std::endl
template <class T> inline bool chkmx(T &a,const T &b){return a<b?a=b,1:0;}
template <class T> inline bool chkmn(T &a,const T &b){return a>b?a=b,1:0;}
//重构一个点是log n 一个点最多被重构log n次
using std::sort;
using std::lower_bound;
using std::swap;
const int MAXN=8e4+10,LOG=20;
namespace iSeg
{
	const int TL=1;int TR;
	const int SIZE=MAXN*LOG*LOG;
	struct Node
	{
		Node *son[2];int cnt;
		Node(Node *lc=0,Node *rc=0){son[0]=lc,son[1]=rc;cnt=0;}
		void up(){cnt=son[0]->cnt+son[1]->cnt;}
	}*POOL=(Node*)malloc(SIZE*sizeof(Node)),*ME=POOL,*null=new(ME++) Node,*root[MAXN];
	Node *newNode()	{return new(ME++) Node;}
	Node *newNode(Node *old)
	{
		static Node *cur;
		cur=newNode();*cur=*old;return cur;
	}
	void Edit(Node *&pos,int L,int R,int goal,int v)
	{
		pos=newNode(pos);
		if(L==R) pos->cnt+=v;
		else
		{
			int MID=(L+R)>>1;
			if(goal<=MID) Edit(pos->son[0],L,MID,goal,v);
			else Edit(pos->son[1],MID+1,R,goal,v);
			pos->up();
		}
	}
	int GetK(Node *p1,Node *p2,Node *lca,Node *lcaf,int L,int R,int K)//可以改成非递归
	{
		if(L==R) return L;
		else
		{
			int sum=p1->son[0]->cnt+p2->son[0]->cnt-lca->son[0]->cnt-lcaf->son[0]->cnt,MID=(L+R)>>1;
			if(K<=sum) return GetK(p1->son[0],p2->son[0],lca->son[0],lcaf->son[0],L,MID,K);//写成了sum<=K
			else return GetK(p1->son[1],p2->son[1],lca->son[1],lcaf->son[1],MID+1,R,K-sum);//写成了K
		}
	}
}using namespace iSeg;
struct Edge{int to,pre;}e[MAXN<<1];int ec,id[MAXN],w[MAXN];
void Adde(int f,int t)
{
	e[++ec].to=t;e[ec].pre=id[f];id[f]=ec;
	e[++ec].to=f;e[ec].pre=id[t];id[t]=ec;
}
int fa[MAXN][LOG],dep[MAXN];
int LCA(int p1,int p2)
{
	if(dep[p1]<dep[p2]) swap(p1,p2);
	for(int j=LOG-1;~j;j--)
		if(dep[fa[p1][j]]>=dep[p2]) p1=fa[p1][j];//写成了<=
	//0,1的dep不能相等。。否则会跳过头。。
	if(p1==p2) return p1;
	for(int j=LOG-1;~j;j--)
		if(fa[p1][j]!=fa[p2][j]) p1=fa[p1][j],p2=fa[p2][j];
	return fa[p1][0];
}
int ud[MAXN],ut,sz[MAXN],rt[MAXN];
void DFS(int pos,int r)
{
	for(int j=1;j<LOG;j++)
		fa[pos][j]=fa[fa[pos][j-1]][j-1];
	ud[pos]=ut;
	rt[pos]=r;sz[r]++;
	dep[pos]=dep[fa[pos][0]]+1;
	root[pos]=root[fa[pos][0]];
	Edit(root[pos],TL,TR,w[pos],1);
	for(int i=id[pos];i;i=e[i].pre)
	{
		int u=e[i].to;
		if(ud[u]==ut) continue;
		fa[u][0]=pos;
		DFS(u,r);
	}
}
int main()
{
	null->son[0]=null->son[1]=null;root[0]=null;

	static int t[MAXN];
	scanf("%*d");
	int n,m,q;get(n),get(m),get(q);TR=n;
	for(int i=1;i<=n;i++) get(w[i]),t[i]=w[i];
	sort(t+1,t+1+n);
	for(int i=1;i<=n;i++) w[i]=lower_bound(t+1,t+1+n,w[i])-t;//写成了-w.....上天了...
	for(int i=1;i<=m;i++)
	{
		int x,y;get(x),get(y);
		Adde(x,y);
	}
	ut++;
	for(int i=1;i<=n;i++)
		if(ud[i]!=ut)
			DFS(i,i);
	int last=0;
	#define decode(x) (x^=last)
	for(int i=1;i<=q;i++)
	{
		char opt;int x,y,K;
		while(!isalpha(opt=getchar()));
		get(x),decode(x),get(y),decode(y);
		if(opt=='L')//Link写的有问题??
		{
			if(sz[rt[x]]>sz[rt[y]]) swap(x,y);
			sz[rt[x]]=0;fa[x][0]=y;
			++ut;DFS(x,rt[y]);
			Adde(x,y);
		}
		else
		{
			get(K),decode(K);
			int lca=LCA(x,y);
			last=t[GetK(root[x],root[y],root[lca],root[fa[lca][0]],TL,TR,K)];
			printf("%d\n",last);
		}
	}
	return 0;
}
/* AC Record(Bugs) *
 * RE  WA once. 
 * * * * * * * * * */

发表于 | 更新于

Link

Solution

找出一条从11到nn，max(ai)+max(bi)max(a_i)+max(b_i)最小的道路 两个量不好一起算， 可以通过枚举aia_i，求出相应最小的bib_i来做。 用LCT维护，边转点，在点上记录边的bib_i。 按照aia_i排个序，然后依次枚举边ee，设 , 为起点，终点 如果 和 不联通，那就连上 否则如果 到 的道路上bib_i最大值大于当前的beb_e，就切掉 否则就continue 然后如果没有continue的话，把这个边连上，更新答案。

Tips

很早就做过，复习拿出来看看打算写写新的下传标记的方法。 于是就惨了 在 里是只有 需要下传，但在LCT里，按照那个原则，在 和 都要下传。也就是说只要 的父子关系要变动，之前就必须要 。于是就很愉快地debug了大半上午。

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
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
//Code by Lucida
#include<bits/stdc++.h>
#define get(x) scanf("%d",&x)
#define put(x) printf("%d",x)
#define getc() getchar()
#define putc(x) putchar(x)
#define puts(x) printf("%s",x)

template <class T> inline bool chkmx(T &a,const T &b){return a<b?a=b,1:0;}
template <class T> inline bool chkmn(T &a,const T &b){return a>b?a=b,1:0;}

using std::max;
using std::min;
using std::sort;
using std::swap;
using std::pair;
using std::make_pair;

const int INF=0x1f1f1f1f,MAXN=50000+11,MAXM=100000+11;
//修改父子关系 必须down().
namespace LCT
{
	const int SIZE=MAXN+MAXM;
	typedef pair<int,int> Data;
	struct Node
	{
		Node *son[2],*fa;
		bool rev;
		Data self,imax;
		int kind(){return fa->son[1]==this;}
		bool isrt(){return fa->son[0]!=this && fa->son[1]!=this;}
		Node()
		{
			rev=0;
			son[0]=son[1]=fa=0;
			self=imax=make_pair(-INF,0);
		}
		void Rev()
		{
			rev^=1;
			swap(son[0],son[1]);
		}
		void down()
		{
			if(rev)
			{
				son[0]->Rev();
				son[1]->Rev();
				rev=0;
			}
		}
		void up(){imax=max(self,max(son[0]->imax,son[1]->imax));}
	}*POOL=(Node*)malloc(SIZE*sizeof(Node)),*ME=POOL,*null=new(ME++) Node,*node[SIZE];
	
	Node *newNode(Data v=make_pair(-INF,0))
	{
		static Node *cur;cur=ME++;
		cur->self=cur->imax=v;
		cur->rev=0;
		cur->son[0]=cur->son[1]=cur->fa=null;
		return cur;
	}
	void Trans(Node *pos)
	{
		static Node *f,*grf;
		static int d;
		f=pos->fa,grf=f->fa;
		f->down(),pos->down();
		d=pos->kind();//有可能变化！！
		if(!f->isrt())
			grf->son[f->kind()]=pos;
		pos->fa=grf;
		f->son[d]=pos->son[!d],pos->son[!d]->fa=f;
		pos->son[!d]=f,f->fa=pos;
		f->up();
		pos->up();
	}
	void Splay(Node *pos)
	{
		while(!pos->isrt())
		{
			Node *f=pos->fa;
			if(!f->isrt())
				Trans(pos->kind()==f->kind()?f:pos);//isrt..
			Trans(pos);
		}
	}
	void Access(Node *pos)
	{
		Node *pred=null;
		while(pos!=null)
		{
			Splay(pos);
			pos->down(),pos->son[1]=pred,pos->up();
			pred=pos;
			pos=pos->fa;
		}
	}
	void BeRoot(Node *pos)
	{
		Access(pos);
		Splay(pos);
		pos->Rev();
	}
	Node *FindRoot(Node *pos)
	{
		Access(pos);Splay(pos);
		while(pos->son[0]!=null) 
			pos=pos->son[0];
		Splay(pos);
		return pos;
	}
	void Link(Node *p1,Node *p2)
	{
		BeRoot(p1);
		p1->fa=p2;
		Splay(p2);
	}
	void Cut(Node *p1,Node *p2)
	{
		BeRoot(p1);
		Access(p2);
		Splay(p2);
		p2->down();p2->son[0]=null;
		p1->fa->down();p1->fa=null;
		p2->up();
	}

	Data Query(Node *p1,Node *p2)
	{
		BeRoot(p1);
		Access(p2);
		Splay(p2);
		return p2->imax;
	}
}using namespace LCT;
struct Edges{int from,to,va,vb;}edges[MAXM];
bool operator <(const Edges &e1,const Edges &e2)
{
	if(e1.va!=e2.va) return e1.va<e2.va;
	return e1.vb<e2.vb;
}
int main()
{
	freopen("input","r",stdin);
	int n,m;get(n),get(m);
	#define edge(i) (i+n)
	#define point(i) (i)
	for(int i=1;i<=m;++i)
		get(edges[i].from),get(edges[i].to),get(edges[i].va),get(edges[i].vb);
	sort(edges+1,edges+1+m);
	for(int i=1;i<=n;++i)
		node[point(i)]=newNode(make_pair(-INF,-i));
	int Ans=INF;
	for(int i=1;i<=m;++i)
	{
		int from=edges[i].from,to=edges[i].to,va=edges[i].va,vb=edges[i].vb;
		if(FindRoot(node[point(from)])==FindRoot(node[point(to)]))
		{
			Data res;
			if((res=Query(node[point(from)],node[point(to)])).first>vb)
			{
				int idx=res.second;
				Cut(node[point(edges[idx].from)],node[edge(idx)]);
				Cut(node[point(edges[idx].to)],node[edge(idx)]);
			}
			else continue;
		}
	
		Node *curEdge=node[edge(i)]=newNode(make_pair(vb,i));
		Link(node[point(from)],curEdge);
		Link(node[point(to)],curEdge);

		if(FindRoot(node[point(1)])==FindRoot(node[point(n)]))
			chkmn(Ans,va+Query(node[point(1)],node[point(n)]).first);
	}
	put(Ans==INF?-1:Ans),putc('\n');
	return 0;
}

发表于 | 更新于

Link

Solution

这种算法都是玄学 每次先random_shuffle来随机改变序列（好像并不能增加随机性？不管了反正本来就是玄学）

每次random地改动一个元素所在的集合，如果更优直接改掉； 否则，有一定的概率改掉，直观地想下，这个概率应该需要逐渐减小的（否则不就等于每次都无视之前所做的努力吗。。），由此就引入了温度的概念。 从初始温度开始模拟退火，每次循环的温度降低一些，接受不优解的概率是TT0\dfrac TT_0，通过随机数来实现。

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
//Code by Lucida
#include<bits/stdc++.h>
#define get(x) scanf("%d",&x)
//#define debug(x) std::cout<<#x<<'='<<x<<std::endl
template <class T> inline bool chkmx(T &a,const T &b){return a<b?a=b,1:0;}
template <class T> inline bool chkmn(T &a,const T &b){return a>b?a=b,1:0;}
const int MAXN=30;
using std::min_element;
using std::swap;
using std::random_shuffle;
template <class T> inline T sqr(T x){return x*x;}
int n,m,a[MAXN];double ANS=1e100,avg;
void Try()
{
	random_shuffle(a+1,a+1+n);
	static double sum[MAXN];static int bn[MAXN];
	memset(sum,0,sizeof(sum));
	double res=0;
	for(int i=1;i<=n;i++) sum[bn[i]=rand()%m+1]+=a[i];
	for(int i=1;i<=m;i++) res+=sqr(sum[i]-avg);
	int P=2e5;
	double Temper=P;
	while(Temper>1e-1)
	{
		Temper*=0.9;
		int cs,now,to;
		do	
		{
			cs=rand()%n+1,now=bn[cs],to=bn[cs];
			if(Temper>500) to=min_element(sum+1,sum+1+m)-sum;
			else to=rand()%m+1;
		}while(now==to);
		double temp=res;
		temp-=sqr(sum[now]-avg),temp-=sqr(sum[to]-avg);
		sum[now]-=a[cs],sum[to]+=a[cs];
		temp+=sqr(sum[now]-avg),temp+=sqr(sum[to]-avg);
		if(temp<res || rand()%P<Temper) bn[cs]=to,res=temp;
		else sum[now]+=a[cs],sum[to]-=a[cs];
	}
	chkmn(ANS,res);
}
int main()
{
	srand(522133279);
	get(n),get(m);
	for(int i=1;i<=n;i++) 
	{
		get(a[i]);
		avg+=a[i];
	}
	avg/=m;
	for(int loop=1;loop<=10000;loop++) Try();
	printf("%.2lf\n",sqrt(ANS/m));
	return 0;
}
/* AC Record(Bugs) *
 * 
 * * * * * * * * * */

发表于 | 更新于

学习了

  • 基础计算几何
  • 基础莫比乌斯反演
  • 基础概率与期望
  • 虚树构造
  • 扩展Floyd
  • Lucas和扩展Lucas
  • FHQ Treap
  • 可持久化化Trie的一些基础用法
  • 大体了解了爬山，模拟退火
  • 大体了解了一点XOR方程组

复习了

  • 基础AC自动机
  • BSGS
  • 线性筛法的一些用处

做的题

  • 乱做了SDOI的题目

总结 这20天做的题目还是有一半以上不能独立做出来，相比之下码力的提高却是很大的。 以后做题还是要多总结，多思考，培养分析问题的能力，补齐自己的短板。 始终要牢记我的目标，不要被盲目地刷Ranklist，让我做的每一道题都有意义、有价值。

从现在到过年的这几天 精做NOI的题目