分类算法下的文章

Codeforces 809E. Surprise me!

May 25, 2017

###Description
给一棵树，每个点点权 $$a_i$$，保证 $$a_i$$各不相同，现在随机选两个点 $$u,v$$，求 $$f(u,v)=\phi(a_u\cdot a_v)\cdot dis(u,v)$$ 的期望，$$mod\ 10^9+7$$。

链接：[http://codeforces.com/contest/809/problem/E](http://codeforces.com/contest/809/problem/E "http://codeforces.com/contest/809/problem/E")

###Solution
首先，有一个结论：$$\phi(a\cdot b)=\frac{\phi(a)\cdot\phi(b)\cdot g}{\phi(g)},g=gcd(a,b)$$
证明：考虑找一个 $$x$$，使得 $$x|b,gcd(\frac{b}{x},a)=1,gcd(x,\frac{b}{x})=1$$ 且 $$a$$ 包含 $$x$$ 中所有质因子。
那么 $$\phi(a\cdot b)=\phi(a)\cdot x\cdot\phi(\frac{b}{x})=\frac{\phi(a)\cdot\phi(b)\cdot x}{\phi(x)}$$
显然，$$x$$ 是 $$g$$ 的倍数，所以 $$\frac{x}{g}=\frac{\phi(x)}{\phi(g)}$$，所以$$\phi(a\cdot b)=\frac{\phi(a)\cdot\phi(b)\cdot g}{\phi(g)}$$

有了这个结论之后可以考虑点分，那么设 $$dis(x)$$ 表示 $$x$$ 到当前根的距离。
设 $$h(x)=\frac{x}{\phi(x)}$$，那么$$f(u,v)=\phi(a_u)\cdot\phi(a_v)\cdot h(g)\cdot(dis(u)+dis(v)),g=gcd(a_u,a_v)$$
考虑当 $$u$$ 确定时，需要知道 $$\sum\phi(a_v)\cdot h(gcd(a_u,a_v))$$ 和 $$\sum\phi(a_v)\cdot h(gcd(a_u,a_v))\cdot dis(v)$$。
设 $$g(x)=h(x)-\sum[y|x,y\neq x]g(y)$$，那么 $$h(gcd(a_u,a_v))=\sum[x|a_u,x|a_v]g(x)$$
这样就可以直接枚举约数算贡献了。

###Code
```c++
#include<bits/stdc++.h>
typedef long long ll;
inline void repr(int&a,int b){if(a<b)a=b;}
#define mset(a,b) memset(a,b,sizeof(a))
#define fo0(i,n) for(int i=0,i##end=n;i<i##end;i++)
#define fo1(i,n) for(int i=1,i##end=n;i<=i##end;i++)
#define fo(i,a,b) for(int i=a,i##end=b;i<=i##end;i++)
#define foe(i,x)for(__typeof(x.end())i=x.begin();i!=x.end();++i)

const int N=200007,P=1000000007;

struct edge
{
	int to;edge*ne;
}_e[N*2],*e=_e,*p[N];

inline void add(int a,int b)
{
	*e=(edge){b,p[a]};p[a]=e++;
}

std::vector<int>v[N];
int n,pr[N],pm,phi[N],inv[N],f[N],w[N],sz[N],rsz,nro,nsz,g[N],h[N],q[N],qe,ans;
bool vis[N];

inline void dfs1(int x,int fa)
{
	sz[x]=1;
	for(edge*i=p[x];i;i=i->ne)
		if(i->to!=fa&&!vis[i->to])
		{
			dfs1(i->to,x);
			sz[x]+=sz[i->to];
		}
}

inline void dfs2(int x,int fa)
{
	int t=rsz-sz[x];
	for(edge*i=p[x];i;i=i->ne)
		if(i->to!=fa&&!vis[i->to])
		{
			repr(t,sz[i->to]);
			dfs2(i->to,x);
		}
	if(t<nsz)nsz=t,nro=x;
}

inline void set(int x,int y,int z)
{
	if(!~g[x])q[qe++]=x,g[x]=h[x]=0;
	g[x]=(g[x]+(ll)y*f[x])%P;
	h[x]=(h[x]+(ll)z*f[x])%P;
}

inline void get(int x,int&y,int&z)
{
	if(~g[x])y=g[x],z=h[x];
	else y=z=0;
}

inline void dfs3(int x,int fa,int dis)
{
	int a=phi[w[x]],b=(ll)a*dis%P;
	foe(i,v[w[x]])set(*i,a,b);
	for(edge*i=p[x];i;i=i->ne)
		if(i->to!=fa&&!vis[i->to])
			dfs3(i->to,x,dis+1);
}

inline void dfs4(int x,int fa,int dis)
{
	int a=phi[w[x]],b=(ll)a*dis%P,A,B;
	foe(i,v[w[x]])
	{
		get(*i,A,B);
		ans=(ans+(ll)a*B+(ll)b*A)%P;
	}
	for(edge*i=p[x];i;i=i->ne)
		if(i->to!=fa&&!vis[i->to])
			dfs4(i->to,x,dis+1);
}

inline void work(int x)
{
	dfs1(x,0);
	rsz=nsz=sz[x],nro=x;
	dfs2(x,0);
	vis[x=nro]=1;
	foe(i,v[w[x]])set(*i,phi[w[x]],0);
	for(edge*i=p[x];i;i=i->ne)
		if(!vis[i->to])
		{
			dfs4(i->to,x,1);
			dfs3(i->to,x,1);
		}
	fo0(i,qe)g[q[i]]=h[q[i]]=-1;
	qe=0;
	for(edge*i=p[x];i;i=i->ne)
		if(!vis[i->to])work(i->to);
}

int main()
{
	fo1(i,N-1)for(int j=i;j<N;j+=i)v[j].push_back(i);
	fo(i,2,N-1)
	{
		if(!vis[i])pr[pm++]=i,phi[i]=i-1;
		for(int j=0;i*pr[j]<N;j++)
		{
			vis[i*pr[j]]=1;
			if(i%pr[j]==0)
			{
				phi[i*pr[j]]=phi[i]*pr[j];
				break;
			}
			phi[i*pr[j]]=phi[i]*(pr[j]-1);
		}
	}
	phi[1]=1;
	inv[1]=1;
	fo(i,2,N-1)inv[i]=(ll)(P-P/i)*inv[P%i]%P;
	fo1(i,N-1)
	{
		f[i]=(ll)i*inv[phi[i]]%P;
		foe(j,v[i])if(i^*j)(f[i]+=P-f[*j])%=P;
	}
	scanf("%d",&n);
	fo1(i,n)scanf("%d",w+i);
	fo0(i,n-1)
	{
		int x,y;
		scanf("%d%d",&x,&y);
		add(x,y);
		add(y,x);
	}
	mset(vis,0);
	mset(g,0xff);
	mset(h,0xff);
	work(1);
	ans=(ll)ans*inv[n]%P*inv[n-1]%P*2%P;
	printf("%d",ans);
}
```

BZOJ 2694 & 4659 Lcm

February 28, 2017

###Description
![fa(1).jpg][1]

###Solution
$$\sum_i^n\sum_j^m[\mu(gcd(i,j))\neq 0]lcm(i,j)$$
$$=\sum_k^n|\mu(k)|\sum_i^{\frac{n}{k}}\sum_j^{\frac{n}{k}}[gcd(i,j)=1]i\cdot j\cdot k$$
$$=\sum_k^n|\mu(k)|\cdot k\sum_d^{\frac{n}{k}}\mu(d)\cdot d^2\sum_i^{\frac{n}{k\cdot d}}\sum_j^{\frac{m}{k\cdot d}}i\cdot j$$
$$=\sum_k^n\sum_d[d|k]\mu(d)\cdot d^2\cdot|\mu(\frac{k}{d})|\cdot\frac{p}{d}\cdot\frac{\frac{n}{k}\cdot(\frac{n}{k}+1)}{2}\cdot\frac{\frac{m}{k}\cdot(\frac{m}{k}+1)}{2}$$
最后的式子，前面是积性函数，后面根号分块

###Code
```c++
#define _GLIBCXX_IOSTREAM
#include<bits/stdc++.h>

#define N 4000001

int pri[N],pm,f[N],ps[N],pc[N];
bool np[N];

int main()
{
	for(int i=2;i<N;i++)
	{
		if(!np[i])pri[pm++]=i,f[i]=1,ps[i]=i,pc[i]=1;
		int rf;
		if(pc[i]==1)rf=f[i]*(ps[i]-ps[i]*ps[i]);
		else if(pc[i]==2)rf=-f[i]*ps[i]*ps[i]*ps[i];
		else rf=0;
		for(int j=0;i*pri[j]<N;j++)
		{
			np[i*pri[j]]=1;
			if(i%pri[j])f[i*pri[j]]=rf,ps[i*pri[j]]=pri[j],pc[i*pri[j]]=1;
			else
			{
				f[i*pri[j]]=f[i],ps[i*pri[j]]=pri[j],pc[i*pri[j]]=pc[i]+1;
				break;
			}
		}
		f[i]=rf;
	}
	f[1]=1;
	for(int i=1;i<N;i++)
		f[i]+=f[i-1];
	int T;
	for(scanf("%d",&T);T--;)
	{
		int n,m,ans=0;
		scanf("%d%d",&n,&m);
		for(int i=1,t,p,q;i<=n&&i<=m;i=t+1)
		{
			p=n/i,q=m/i,t=std::min(n/p,m/q);
			ans+=(f[t]-f[i-1])*p*(p+1)*q*(q+1)>>2;
		}
		printf("%d\n",ans&0x3fffffff);
	}
}
```

[1]: /usr/uploads/2017/02/3245901942.jpg
  [2]: /usr/uploads/2017/02/1279320617.png

BZOJ 3926: [Zjoi2015]诸神眷顾的幻想乡

February 20, 2017

##Description
幽香是全幻想乡里最受人欢迎的萌妹子，这天，是幽香的2600岁生日，无数幽香的粉丝到了幽香家门前的太阳花田上来为幽香庆祝生日。 
粉丝们非常热情，自发组织表演了一系列节目给幽香看。幽香当然也非常高兴啦。 
这时幽香发现了一件非常有趣的事情，太阳花田有n块空地。在过去，幽香为了方便，在这n块空地之间修建了n-1条边将它们连通起来。也就是说，这n块空地形成了一个树的结构。 
有n个粉丝们来到了太阳花田上。为了表达对幽香生日的祝贺，他们选择了c中颜色的衣服，每种颜色恰好可以用一个0到c-1之间的整数来表示。并且每个人都站在一个空地上，每个空地上也只有一个人。这样整个太阳花田就花花绿绿了。幽香看到了，感觉也非常开心。 
粉丝们策划的一个节目是这样的，选中两个粉丝A和B（A和B可以相同），然后A所在的空地到B所在的空地的路径上的粉丝依次跳起来（包括端点），幽香就能看到一个长度为A到B之间路径上的所有粉丝的数目（包括A和B）的颜色序列。一开始大家打算让人一两个粉丝（注意：A,B和B,A是不同的，他们形成的序列刚好相反，比如红绿蓝和蓝绿红）都来一次，但是有人指出这样可能会出现一些一模一样的颜色序列，会导致审美疲劳。 
于是他们想要问题，在这个树上，一共有多少可能的不同的颜色序列（子串）幽香可以看到呢？ 
太阳花田的结构比较特殊，只与一个空地相邻的空地数量不超过20个。

###Input
第一行两个正整数n,c。表示空地数量和颜色数量。 
第二行有n个0到c-1之间，由空格隔开的整数，依次表示第i块空地上的粉丝的衣服颜色。（这里我们按照节点标号从小到大的顺序依次给出每块空地上粉丝的衣服颜色）。 
接下来n-1行，每行两个正整数u,v，表示有一条连接空地u和空地v的边。

###Output
一行，输出一个整数，表示答案。

###Sample Input
7 3
0 2 1 2 1 0 0
1 2
3 4
3 5
4 6
5 7
2 5

###Sample Output
30

###HINT
对于所有数据，1<=n<=100000, 1<=c<=10。 
对于15%的数据，n<=2000。 
另有5%的数据，所有空地都至多与两个空地相邻。 
另有5%的数据，除一块空地与三个空地相邻外，其他空地都分别至多与两个空地相邻。 
另有5%的数据，除某两块空地与三个空地相邻外，其他空地都分别至多与两个空地相邻。

###Solution
考虑一条从u到v的路径上的字符串，如果在两端不断加节点，一定可以扩充成一条两个度数为1的节点间的路径上的子串。
那么显然有一种方法是把20个节点两两间的字符串建成广义SAM，然后直接统计，但是这样时间、内存都会超。
建广义SAM时，每次会把last设为初始节点，那么对于这道题，也可以类似优化。
把每个度数为1的节点作为根dfs，每次dfs时记下当前last，然后在处理完一棵子树后还原last。

###Code
```c++
#define _GLIBCXX_IOSTREAM
#include<bits/stdc++.h>

typedef long long ll;

#define N 2000010

struct node
{
	int ma,inc;ll cnt;
	node*link,*nxt[10];
}pool[N*2],*pm=pool+1,*start=pool;

namespace _A{node*lst,*cur,*p,*q,*sq;}
void add(int t)
{
	using namespace _A;
	cur=pm++;
	cur->ma=lst->ma+1;
	for(p=lst;p&&!p->nxt[t];p=p->link)p->nxt[t]=cur;
	if(!p){cur->link=start;goto naive;}
	q=p->nxt[t];
	if(p->ma+1==q->ma){cur->link=q;goto naive;}
	sq=pm++;
	sq->ma=p->ma+1;
	memcpy(sq->nxt,q->nxt,sizeof(q->nxt));
	for(;p&&p->nxt[t]==q;p=p->link)p->nxt[t]=sq;
	sq->link=q->link;
	q->link=sq;
	cur->link=sq;
	naive:lst=cur;
}

struct edge
{
	int to;edge*ne;
}_e[200000],*e=_e,*p[100001];

inline void add(int a,int b)
{
	*e=(edge){b,p[a]};p[a]=e++;
}

int c[100001],deg[100001];

void dfs(int x,int fa)
{
	add(c[x]);
	node*t=_A::lst;
	for(edge*i=p[x];i;i=i->ne)
		if(i->to^fa)
		{
			_A::lst=t;
			dfs(i->to,x);
		}
}

node*q[N*2];

int main()
{
	int n,x,y;
	scanf("%d%*d",&n);
	for(int i=1;i<=n;i++)
		scanf("%d",c+i);
	for(int i=1;i<n;i++)
	{
		scanf("%d%d",&x,&y);
		deg[x]++;
		deg[y]++;
		add(x,y);
		add(y,x);
	}
	for(int i=1;i<=n;i++)if(deg[i]==1)
		_A::lst=start,dfs(i,0);
	for(node*i=pool;i<pm;i++)
		for(int j=0;j<10;j++)if(i->nxt[j])i->nxt[j]->inc++;
	int qe=0;
	for(node*i=pool;i<pm;i++)
		if(!i->inc)i->inc=1e9,q[qe++]=i;
	for(int i=0;i<qe;i++)
		for(int j=0;j<10;j++)if(q[i]->nxt[j])if(!--q[i]->nxt[j]->inc)
			q[i]->nxt[j]->inc=1e9,q[qe++]=q[i]->nxt[j];
	for(int i=qe-1;~i;i--)
	{
		for(int j=0;j<10;j++)if(q[i]->nxt[j])q[i]->cnt+=q[i]->nxt[j]->cnt;
		q[i]->cnt++;
	}
	printf("%lld",start->cnt-1);
}
```

对于各种并查集写法速度的研究

January 24, 2017

2018.3.3：更新了测试代码，加入了更多写法，修复了一些致命错误
测试机CPU为Intel® Celeron® Processor G3900，内存为8GB DDR4，系统为Debian 9。
编译命令为g++ test.cpp -o test -O3。
先放测试代码：
```c++
#include<cstdio>
#include<ctime>
#include<random>
#include<cassert>
#define N 1024
#define T 100000
int fa[N+233],st[N+233],sz[N+233];
int find1(int x)
{
	return x==fa[x]?x:fa[x]=find1(fa[x]);
}
inline int find2(int x)
{
	return x==fa[x]?x:fa[x]=find2(fa[x]);
}
inline int find3(int x)
{
	while(x!=fa[x])x=fa[x]=fa[fa[x]];return x;
}
inline int find4(int x)
{
	int se=0;st[0]=x;
	while(fa[x]!=x)st[++se]=x=fa[x];
	while(se)fa[st[--se]]=x;
	return x;
}
inline int find5(int x)
{
	int t=x,p;
	while(x!=fa[x])x=fa[x];
	while(t!=x)p=fa[t],fa[t]=x,t=p;
	return x;
}
int find6(int x)
{
	return 0>fa[x]?x:fa[x]=find6(fa[x]);
}
inline int find7(int x)
{
	return 0>fa[x]?x:fa[x]=find7(fa[x]);
}
inline int find8(int x)
{
	while(fa[x]>=0&&fa[fa[x]]>=0)x=fa[x]=fa[fa[x]];return fa[x]<0?x:fa[x];
}
long long global_checksum;
unsigned int count,seed;
std::mt19937 ran_(0);
void init_rand(int x)
{
	ran_=std::mt19937(x);
	ran_(),ran_(),ran_();
	count=0;
}
int rand()
{
	if(count&1023)seed=ran_();
	count++;
	seed^=seed<<13;
	seed^=seed>>17;
	seed^=seed<<5;
	return seed;
}
inline long long test(int(*f)(int))
{
	init_rand(19260817);
	long long last,timeusage=0,checksum=0;
	for(int t=0;t<T;t++)
	{
		for(int i=0;i<N;i++)fa[i]=i;
		last=clock();
		for(int i=0,x,y;i<N;i++)
		{
			x=rand(),y=rand();
			x&=N-1,y&=N-1;
			x=f(x),y=f(y);
			if(x!=y)fa[x]=y;
		}
		timeusage+=clock()-last;
		for(int i=0;i<N;i++)checksum=checksum*998244353+(f(i)==f(0));
	}
	assert(global_checksum==checksum||global_checksum==0);
	global_checksum=checksum;
	return timeusage;
}
inline long long test2(int(*f)(int))
{
	init_rand(19260817);
	long long last,timeusage=0,checksum=0;
	for(int t=0;t<T;t++)
	{
		for(int i=0;i<N;i++)fa[i]=i;
		for(int i=0;i<N;i++)sz[i]=1;
		last=clock();
		for(int i=0,x,y;i<N;i++)
		{
			x=rand(),y=rand();
			x&=N-1,y&=N-1;
			x=f(x),y=f(y);
			if(x!=y){if(sz[x]>sz[y])sz[x]+=sz[y],fa[y]=x;else sz[y]+=sz[x],fa[x]=y;}
		}
		timeusage+=clock()-last;
		for(int i=0;i<N;i++)checksum=checksum*998244353+(f(i)==f(0));
	}
	assert(global_checksum==checksum||global_checksum==0);
	global_checksum=checksum;
	return timeusage;
}
inline long long test3(int(*f)(int))
{
	init_rand(19260817);
	long long last,timeusage=0,checksum=0;
	for(int t=0;t<T;t++)
	{
		for(int i=0;i<N;i++)fa[i]=-1;
		last=clock();
		for(int i=0,x,y;i<N;i++)
		{
			x=rand(),y=rand();
			x&=N-1,y&=N-1;
			x=f(x),y=f(y);
			if(x!=y){if(fa[x]<fa[y])fa[x]+=fa[y],fa[y]=x;else fa[y]+=fa[x],fa[x]=y;}
		}
		timeusage+=clock()-last;
		for(int i=0;i<N;i++)checksum=checksum*998244353+(f(i)==f(0));
	}
	assert(global_checksum==checksum||global_checksum==0);
	global_checksum=checksum;
	return timeusage;
}
inline long long test2_1(int(*f)(int))
{
	init_rand(19260817);
	long long last,timeusage=0,checksum=0;
	for(int t=0;t<T;t++)
	{
		for(int i=0;i<N;i++)fa[i]=i;
		for(int i=0;i<N;i++)sz[i]=1;
		last=clock();
		for(int i=0,x,y;i<N;i++)
		{
			x=rand(),y=rand();
			x&=N-1,y&=N-1;
			x=f(x),y=f(y);
			if(x!=y){if(sz[x]>sz[y])fa[y]=x;else{fa[x]=y;if(sz[x]==sz[y])sz[y]++;}}
		}
		timeusage+=clock()-last;
		for(int i=0;i<N;i++)checksum=checksum*998244353+(f(i)==f(0));
	}
	assert(global_checksum==checksum||global_checksum==0);
	global_checksum=checksum;
	return timeusage;
}
inline long long test3_1(int(*f)(int))
{
	init_rand(19260817);
	long long last,timeusage=0,checksum=0;
	for(int t=0;t<T;t++)
	{
		for(int i=0;i<N;i++)fa[i]=-1;
		last=clock();
		for(int i=0,x,y;i<N;i++)
		{
			x=rand(),y=rand();
			x&=N-1,y&=N-1;
			x=f(x),y=f(y);
			if(x!=y){if(fa[x]<fa[y])fa[y]=x;else{if(fa[x]==fa[y])fa[y]--;fa[x]=y;}}
		}
		timeusage+=clock()-last;
		for(int i=0;i<N;i++)checksum=checksum*998244353+(f(i)==f(0));
	}
	assert(global_checksum==checksum||global_checksum==0);
	global_checksum=checksum;
	return timeusage;
}
int main()
{
	printf("CLOCKS_PER_SEC:%d\n",CLOCKS_PER_SEC);
	printf("find1:%lld %lld %lld\n",test(find1),test2(find1),test2_1(find1));
	printf("find2:%lld %lld %lld\n",test(find2),test2(find2),test2_1(find2));
	printf("find3:%lld %lld %lld\n",test(find3),test2(find3),test2_1(find3));
	printf("find4:%lld %lld %lld\n",test(find4),test2(find4),test2_1(find4));
	printf("find5:%lld %lld %lld\n",test(find5),test2(find5),test2_1(find5));
	printf("find6:%lld %lld\n",test3(find6),test3_1(find6));
	printf("find7:%lld %lld\n",test3(find7),test3_1(find7));
	printf("find8:%lld %lld\n",test3(find8),test3_1(find8));
}
```
当N=1024,T=500000时，输出：
```
find1:20089351 15887467 16531550
find2:20097390 15910071 16540963
find3:17256541 13390076 13982165
find4:19644034 15645215 16138683
find5:19427605 15368866 15878426
find6:15846219 16573551
find7:15818575 16527184
find8:14988614 15839679
```
当N=131072,T=5000时，输出：
```
CLOCKS_PER_SEC:1000000
find1:40043295 28478678 29630797
find2:40079319 28487803 29605407
find3:34651392 23403124 24607184
find4:40438077 27544080 28666981
find5:39498635 26400323 27404284
find6:25862629 26859295
find7:25814962 26853594
find8:24541376 25698358
```
当N=2097152,T=500时，输出：
```
CLOCKS_PER_SEC:1000000
find1:240728970 173155213 166865669
find2:240076225 173137194 166896002
find3:257240885 186958166 159086582
find4:240671337 168750996 164743633
find5:210863641 130934764 128655852
find6:107561256 107685991
find7:107630995 107676199
find8:101390000 100197421
```
当N=16777216,T=400时，输出：
```
CLOCKS_PER_SEC:1000000
find1:2537867344 1538207603 1438028041
find2:2542492294 1537816024 1438141065
find3:2878818361 1750536054 1388536804
find4:2497444413 1495104276 1413699601
find5:2212506256 1083761103 1015381467
find6:983906431 926941775
find7:984064763 926552621
find8:892107319 808279658
```
（原结果有误，已删）
最后得出的结果是find8（while非递归+按秩合并+秩和fa记在一个数组上）最快，当N较小时秩选用size，当N较大时秩选用depth。
当不方便非递归时可以写递归的（find6或find7）。
inline似乎没有明显优化。

BZOJ 3160: 万径人踪灭

January 20, 2017

###Description
略。见[BZOJ3160][1]

###Solution
对于“不能是连续的一段”，可以先计算所有答案，然后减去连续一段的，即回文串，这个用回文自动机或者manacher可以得出
考虑以每个位置为中心的对称元素对个数f[i]，那么i位置对答案的贡献是$$2^{f[i]}-1$$
如果把a视为1，b视为-1，然后把原串长度\*4，卷积，那么以i为中心的元素对如果相同，则贡献1，否则贡献-1

###Code
```c++
#include<bits/stdc++.h>
typedef double lf;

#define fo0(i,n) for(int i=0,i##end=n;i<i##end;i++)
#define fo1(i,n) for(int i=1,i##end=n;i<=i##end;i++)

inline int mod(int x,int y){if(x>=y)return x-y;return x;}

#define N 524288

struct num
{
	lf r,i;
	num(){r=i=0;}
	num(lf x){r=x,i=0;}
	num(lf x,lf y){r=x,i=y;}
	num(const num &x){r=x.r,i=x.i;}
	inline num& operator +=(const num &x){r+=x.r,i+=x.i;}
	inline num& operator -=(const num &x){r-=x.r,i-=x.i;}
	inline num& operator *=(const num &x){lf a=r*x.r-i*x.i,b=r*x.i+i*x.r;r=a,i=b;}
	inline num operator +(const num &x){num t=*this;t+=x;return t;}
	inline num operator -(const num &x){num t=*this;t-=x;return t;}
	inline num operator *(const num &x){num t=*this;t*=x;return t;}
};

int id[N];
num tx[N],mf[N];

inline void init_id(int y)
{
	fo0(i,1<<y)id[i]=id[i>>1]>>1|(i&1)<<y-1;
}

#define idft(a,b) dft(a,b,1)

inline void dft(num *s,int n,bool II=0)
{
	fo0(i,1<<n)tx[i]=s[id[i]];
	fo1(p,n)
	{
		lf tmp=M_PI*2/(1<<p);
		if(II)tmp=-tmp;
		num t0(1),t1(cos(tmp),sin(tmp));
		fo0(i,1<<p-1)mf[i]=t0,t0*=t1;
		num X,Y;
		for(int i=0,ie=1<<p,ir=ie/2,pp=ir-1;i<(1<<n);i+=ie)fo0(j,ir)
		{
			X=tx[i|j],Y=tx[i|j|ir]*mf[j];
			tx[i|j]=X+Y,tx[i|j|ir]=X-Y;
		}
	}
	fo0(i,1<<n)s[i]=tx[i];
}

#define N2 100010

struct node
{
	int len,cnt;
	node *nxt[2],*fail;
}t[N2+2],*tm=t+2;

inline void build(char *s,int n)
{
	t[0].len=0,t[1].len=-1;
	node *cur=t+1,*tmp,*t2;
	t[0].fail=t+1;
	for(int i=0;i<n;i++)
	{
		char u=s[i]-'a';
		for(;s[i]!=s[i-1-cur->len];cur=cur->fail);
		if(cur->nxt[u])
			cur=cur->nxt[u],cur->cnt++;
		else
		{
			tmp=cur->nxt[u]=tm++;
			tmp->len=cur->len+2;
			tmp->cnt=1;
			if(tmp->len==1)
				tmp->fail=t;
			else
			{
				for(t2=cur->fail;s[i]!=s[i-1-t2->len];t2=t2->fail);
				tmp->fail=t2->nxt[u];
			}
			cur=tmp;
		}
	}
}

inline void count()
{
	for(node *x=tm-1;x>=t;x--)
		if(x->fail)x->fail->cnt+=x->cnt;
}

#define P 1000000007

inline int cnt()
{
	int ans=0;
	for(node *x=tm-1;x>t+1;x--)
		ans=mod(ans+x->cnt,P);
	return ans;
}

char s[100000];
num p[N];
int po[100001];

int main()
{
	scanf("%s",s);
	int n=strlen(s);
	fo1(i,n/2+5)po[i]=mod(mod(po[i-1]*2+1,P),P);
	int y=0;
	while((1<<(y+1))<n*4)y++;y++;
	fo0(i,n)p[i]=s[i]=='a'?1:-1;
	init_id(y);
	dft(p,y);
	fo0(i,1<<y)p[i]*=p[i];
	idft(p,y);
	lf t=1.0/(1<<y);
	int ans=0;
	fo0(i,n*2-1)
	{
		int cnt=round(p[i].r*t);
		cnt=(std::min(i/2+1,(n*2-i)/2)+(cnt&1?(cnt+1)/2:cnt/2))/2;
		ans=mod(ans+po[cnt],P);
	}
	build(s,n);
	count();
	printf("%d",mod(ans+P-cnt(),P));
}
```

[1]: http://www.lydsy.com/JudgeOnline/problem.php?id=3160

BZOJ 2989: 数列 & 4170: 极光

January 8, 2017

###Description
给定一个长度为n的正整数数列a[i]。
定义2个位置的graze值为两者位置差与数值差的和，即graze(x,y)=|x-y|+|a[x]-a[y]|。
2种操作（k都是正整数）：
1.Modify x k：将第x个数的值修改为k。
2.Query x k：询问有几个i满足graze(x,i)<=k。因为可持久化数据结构的流行，询问不仅要考虑当前数列，还要考虑任意历史版本，即统计任意位置上出现过的任意数值与当前的a[x]的graze值<=k的对数。（某位置多次修改为同样的数值，按多次统计）

###Input
第1行两个整数n,q。分别表示数列长度和操作数。
第2行n个正整数，代表初始数列。
第3--q+2行每行一个操作。

###Output
对于每次询问操作，输出一个非负整数表示答案

###Sample Input
3 5
2 4 3
Query 2 2
Modify 1 3
Query 2 2
Modify 1 2
Query 1 1

###Sample Output
2
3
3

###HINT
N<=60000 修改操作数<=40000 询问<=60000 Max{a[i]}含修改<=100000

###Solution
这道题其实是维护两种操作，向平面上加点，查询菱形内点数
把(x,y)变为(x+y,x-y)，菱形就变成了正方形
假设某个询问是以(x0,y0)为中心，边长为2k的正方形，那么可以拆成两个，一个是x>=x0-k且x<=x0+k且y<=y0+k的点数，另一个是x>=x0-k且x<=x0+k且y<y0-k的点数，把两个询问做差就和原询问相同
那么可以CDQ分治，按y排序，x用树状数组维护

###Code
```c++
#include<bits/stdc++.h>

using namespace std;

int n,m,s[60001],ans[60000],ac,qc,mx;

struct query
{
	int x,y,k,id;char op;
}q[230000],tmp[230000];

inline void add(int x,int y){q[qc].x=x+y,q[qc].y=x-y,qc++,mx=max(mx,x+y);}

int p[160001];
inline void pl(int x){for(;x<=mx;x+=x&-x)p[x]++;}
inline int qr(int x){if(x<=0)return 0;int r=0;for(x=min(x,mx);x;x^=x&-x)r+=p[x];return r;}
inline void cl(int x){for(;x<=mx;x+=x&-x)p[x]=0;}

void work(int l,int r)
{
	if(l==r)return;
	int mid=(l+r)>>1;
	work(l,mid),work(mid+1,r);
	int i1=l,i2=mid+1,ti=l;
	for(;i2<=r;tmp[ti++]=q[i2++])
	{
		for(;i1<=mid&&q[i1].y<=q[i2].y;tmp[ti++]=q[i1++])
			if(!q[i1].op)pl(q[i1].x);
		if(q[i2].op)
		{
			if(q[i2].op==1)
				ans[q[i2].id]+=qr(q[i2].x+q[i2].k)-qr(q[i2].x-q[i2].k-1);
			else
				ans[q[i2].id]+=qr(q[i2].x-q[i2].k-1)-qr(q[i2].x+q[i2].k);
		}
	}
	for(;i1<=mid;tmp[ti++]=q[i1++]);
	for(i1=l;i1<=mid;i1++)if(!q[i1].op)cl(q[i1].x);
	for(i1=l;i1<=r;i1++)q[i1]=tmp[i1];
}

int main()
{
	scanf("%d%d",&n,&m);
	for(int i=1;i<=n;i++)
		scanf("%d",s+i),add(i,s[i]);
	char tmp[10];int a,b;
	while(m--)
	{
		scanf("%s%d%d",tmp,&a,&b);
		if(tmp[0]=='M')
			add(a,s[a]=b);
		else
		{
			int x=a+s[a],y=a-s[a];
			q[qc].x=x,q[qc].y=y+b,q[qc].k=b,q[qc].id=ac,q[qc].op=1,qc++;
			q[qc].x=x,q[qc].y=y-b-1,q[qc].k=b,q[qc].id=ac,q[qc].op=2,qc++;
			ac++;
		}
	}
	work(0,qc-1);
	for(int i=0;i<ac;i++)
		printf("%d\n",ans[i]);
}
```

BZOJ 2958&3269: 序列染色

December 26, 2016

###Description
给出一个长度为N由B、W、X三种字符组成的字符串S，你需要把每一个X染成B或W中的一个。
对于给出的K，问有多少种染色方式使得存在整数a,b,c,d使得:
1<=a<=b<c<=d<=N
Sa,Sa+1,...,Sb均为B
Sc,Sc+1,...,Sd均为W
其中b=a+K-1,d=c+K-1
由于方法可能很多，因此只需要输出最后的答案对109+7取模的结果。

###Input
第一行两个正整数N,K
第二行一个长度为N的字符串S

###Output
一行一个整数表示答案%(109+7)。

###Sample Input
5 2
XXXXX

###Sample Output
4

###Solution
首先搞出B和W的前缀个数，那么可以O(1)判某段区间能否为B或W
然而这样对于超过K个的会算重
对于连续的B和W，考虑用最后K个计算，则不会算重（可以在序列最后加一个X，更好写）
状态可以设计成f[i][j][k]表示考虑前i位，状态为j，这一位选了k，j为0,1,2分别代表BW都没有，只有B，BW都有
转移看代码

###Code
```c++
#include<bits/stdc++.h>

#define mod 1000000007

int n,k,bp[1000010],wp[1000010],f[1000010][3][2];
char s[1000010];

int main()
{
	scanf("%d%d%s",&n,&k,s+1);
	s[++n]='X';
	for(int i=1;i<=n;i++)
	{
		bp[i]=bp[i-1]+(s[i]=='B');
		wp[i]=wp[i-1]+(s[i]=='W');
	}
	f[0][0][0]=1;
	for(int i=1;i<=n;i++)
	{
		if(s[i]=='B'||s[i]=='X')
			for(int j=0;j<3;j++)f[i][j][0]=(f[i-1][j][0]+f[i-1][j][1])%mod;
		if(s[i]=='W'||s[i]=='X')
			for(int j=0;j<3;j++)f[i][j][1]=(f[i-1][j][0]+f[i-1][j][1])%mod;
		if(i<=k)continue;
		if((s[i]=='B'||s[i]=='X')&&!(bp[i-1]-bp[i-k-1]))
		{
			f[i][2][0]=(f[i][2][0]+f[i-k][1][1])%mod;
			f[i][1][0]=(f[i][1][0]-f[i-k][1][1])%mod;
		}
		if((s[i]=='W'||s[i]=='X')&&!(wp[i-1]-wp[i-k-1]))
		{
			f[i][1][1]=(f[i][1][1]+f[i-k][0][0])%mod;
			f[i][0][1]=(f[i][0][1]-f[i-k][0][0])%mod;
		}
	}
	int ans=f[n][2][0];
	if(ans<0)ans+=mod;
	printf("%d\n",ans);
}
```

BZOJ 4398: 福慧双修

December 26, 2016

###Description
菩萨为行，福慧双修，智人得果，不忘其本。
——唐朠立《大慈恩寺三藏法师传》
有才而知进退，福慧双修，这才难得。
——乌雅氏
如何福慧双修？被太后教导的甄嬛徘徊在御花园当中。突然，她发现御花园中的花朵全都是红色和蓝色的。她冥冥之中得到了响应：这就是指导她如何福慧双修的！ 现在御花园可以看作是有N块区域，M条小路，两块区域之间可通过小路连接起来。现在甄嬛站在1号区域，而她需要在御花园中绕一绕，且至少经过1个非1号区 域的区域。但是恰好1号区域离碎玉轩最近，因此她最后还是要回到1号区域。由于太后教导她要福慧双修，因此，甄嬛不能走过任何一条她曾经走过的路。但是， 御花园中来往的奴才们太多了，而且奴才们前行的方向也不一样，因此甄嬛在走某条小路的时候，方向不同所花的时间不一定一样。天色快暗了，甄嬛需要尽快知道 至少需要花多少时间才能学会如何福慧双修。如果甄嬛无法达到目的，输出“-1”。

###Input
第一行仅2个正整数n,m，意义如题。
接下来m行每行4个正整数s,t,v,w，其中s,t为小路所连接的两个区域的编号，v为甄嬛从s到t所需的时间，w为甄嬛从t到s所需的时间。数据保证无重边。

###Output
仅一行，为甄嬛回到1号区域所需的最短时间，若方案不存在，则输出-1

###Sample Input
3 3
1 2 2 3
2 3 1 4
3 1 5 2

###Sample Output
8

###Solution
一条合法路径一定是1->a->...->b->1这样的
考虑a和b一定有至少一个二进制位不同
所以按二进制分组，每次新建一个源点和汇点，分别连向某位为1和0的点，然后跑最短路

###Code
```c++
#include<bits/stdc++.h>

#define mp(a,b) std::make_pair(a,b)

int n,p[40001],dis[40001];
bool vis[40001];
std::priority_queue<std::pair<int,int>,std::vector<std::pair<int,int> >,
    std::greater<std::pair<int,int> > >q;

struct edge
{
	int to,ne,w;
}e[200001];

inline void add(int i,int a,int b,int w)
{
	e[i].to=b,e[i].ne=p[a],e[i].w=w,p[a]=i;
}

inline void yjq(int andv,int anda,int &ans)
{
	memset(dis,0x7f,sizeof(dis));
	memset(vis,0,sizeof(vis));
	for(int j=p[1];j;j=e[j].ne)
		if((e[j].to&andv)==anda)
			q.push(mp(dis[e[j].to]=e[j].w,e[j].to));
	while(!q.empty())
	{
		int k=q.top().second;q.pop();
		if(vis[k])continue;
		vis[k]=1;
		for(int j=p[k];j;j=e[j].ne)
			if(e[j].to==1)
			{
				if((k&andv)!=anda&&ans>dis[k]+e[j].w)ans=std::min(ans,dis[k]+e[j].w);
			}
			else if(dis[e[j].to]>dis[k]+e[j].w)
			{
				q.push(mp(dis[e[j].to]=dis[k]+e[j].w,e[j].to));
			}
	}
}

int main()
{
	int m;
	scanf("%d%d",&n,&m);
	for(int i=0;i<m;i++)
	{
		int a,b,c,d;
		scanf("%d%d%d%d",&a,&b,&c,&d);
		add(i*2+1,a,b,c);
		add(i*2+2,b,a,d);
	}
	int ans=0x7fffffff;
	for(int i=1;i<=n;i<<=1)
	{
		yjq(i,i,ans);
		yjq(i,0,ans);
	}
	printf("%d",ans);
}
```