BZOJ 2957 楼房重建

发表于 2017-01-15 | 更新于 2018-06-17

Link

Solution

可以看出，一个点对于它左边的区间的答案是不会产生影响的，而对右边的区间的影响是遮住了所有与原点连线的斜率比它小的。于是脑中模模糊糊地浮现出一棵线段树，那棵树上的每个节点长着一棵Treap。插入一个新节点，把它右边的区间中比它左边区间中斜率最大值大，而又比新插入的点斜率值小的点从右区间的答案中删去。删除一个节点，把它右边区间中比它的斜率小但比它左边区间中最大值大的点加入答案。 $\uparrow$ 口胡做法不一定对

最后发现，正解的做法非常巧妙。大致思路差不多，但在更新答案的时候只需要在线段树上走即可。首先，把新值赋到线段树的对应节点上。回溯的时候考虑已经知道了左半区间和右半区间的答案，如何得到整个区间的答案。思路就是，左半区间的答案不会改变，可以直接累加，而右半区间需要把比左半区间最大值小的点删掉。对于右半区间 $[mid,r]$ 递归处理，设左半区间的最大值为 $ms$ 。假如 $[mid,r]$ 左半区间的最大值 $< ms$ ，那左半区间就别想了，只需要统计右半区间有多少 $> ms$ 的点。 $==ms$ ，那左半区间就别想了，直接返回右半区间的答案对答案的贡献，即pos->cnt-pos->son[0]->cnt $> ms$ ，那就不会影响右半区间的答案，只需要统计左半区间有多少 $> ms$ 的点，加上右半区间的答案即可。很神奇的思路。每次操作只会走一条到叶节点的路径。单点的修改，把所有的路径长度累加，是 $1+2+...+\log n$ ，也就是 $\dfrac {\log n(1+\log n)}2$ ，即总 $O(\log^2 n)$ 再加上线段树本来的复杂度，总 $O(n\log^2 n)$ ，常数很小，跑得很快。

Tips

依旧是分析问题性质是最重要的，通过左右区间的影响关系的彻底分析，才得出这种很好写跑得快的方法。

Code

//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=100000+10;
typedef long double ld;
using std::max;
namespace iSeg
{
	const int SIZE=MAXN<<2;
	struct Node
	{
		Node *son[2];
		ld ms;int cnt;
		void up();
		bool isleaf();
		Node(){son[0]=son[1]=0;ms=cnt=0;}
	}*POOL=(Node*)malloc(SIZE*sizeof(Node)),*ME=POOL,*null=new(ME++) Node,*root=null;
	Node *newNode(ld slope=0)
	{
		static Node* cur;cur=ME++;
		cur->son[0]=cur->son[1]=null;cur->ms=slope;
		return cur;
	}
	bool Node::isleaf(){return son[0]==son[1] && son[0]==null;}
	int Count(Node *pos,ld s)
	{
		if(pos==null) return 0;
		if(pos->isleaf()) return pos->ms>s;
		if(pos->son[0]->ms<s) return Count(pos->son[1],s); 
		else if(pos->son[0]->ms==s) return pos->cnt-pos->son[0]->cnt;
		else return pos->cnt-pos->son[0]->cnt+Count(pos->son[0],s);
	}
	void Node::up()
	{
		ms=max(son[0]->ms,son[1]->ms);
		cnt=son[0]->cnt+Count(son[1],son[0]->ms);
	}
	void Edit(Node *&pos,int L,int R,int goal,ld s)
	{
		if(pos==null) pos=newNode();
		if(L==R) pos->ms=s,pos->cnt=1;
		else
		{
			int MID=(L+R)>>1;
			if(goal<=MID) Edit(pos->son[0],L,MID,goal,s);
			else Edit(pos->son[1],MID+1,R,goal,s);
			pos->up();
		}
	}
}using namespace iSeg;
int main()
{
	int n,m;get(n),get(m);
	for(int i=1;i<=m;i++)
	{
		int x,y;get(x),get(y);
		Edit(root,1,n,x,(ld)y/x);
		printf("%d\n",root->cnt);
	}
	return 0;
}
/* AC Record(Bugs) *
 * 
 * * * * * * * * * */

BZOJ 1094 粒子运动

发表于 2017-01-15 | 更新于 2018-06-18

Link

Solution

~~模拟即可~~ ~~写了一上午还好意思这么说~~ 在外层循环枚举点对，转化成求点 $a,b$ 在整个过程中的最近距离。可以根据粒子撞壁把整个过程划分成 $O(K)$ 个阶段。假设一个时刻， $a$ 撞了壁，转向， $b$ 正在去撞壁的路上。用它们的速度算出来谁会先撞壁，求出来时间 $t$ 。在 $t$ 时间内，两个点的轨迹都是一个直线，所以可以用二次函数求出在这一端 $t$ 时间内两点距离的最小值。如果用直线的一般式方程，化出的式子比较麻烦。而用直线的点法式方程，可以很简单地求出最值。 = $\sqrt{(v_a-v_b)^2t^2+2(v_a-v_b)(P_a-P_b)t+(P_a-P_b)^2}$ 撞壁的时候需要计算转向，画个图感受感受就会了。。

Tips

没错，OI题也可以用二次函数求最值。细节很多。。难写。。

Code

//Code by Lucida
#include<bits/stdc++.h>
#define red(x) scanf("%d",&x)
#define fred(x) scanf("%lf",&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;}
typedef double ld;
const int MAXN=100+10;
typedef double ld;
const ld eps=1e-5,INF=1e100,pi=acos(-1.0);
using std::min;
using std::pair;
using std::make_pair;
int fcmp(ld x)
{
	if(-eps<x && x<eps) return 0;
	return x<0?-1:1;
}
template <class T> T sqr(T x){return x*x;}
template <class T> T abs(T x){return x>0?x:-x;}
struct vec
{
	ld x,y;
	vec(){}
	vec(ld _x,ld _y):x(_x),y(_y){}
	ld norm(){return sqrt(x*x+y*y);}//sqrt(y*y)???
	vec& operator +=(vec a){x+=a.x,y+=a.y;return *this;}
	vec& operator -=(vec a){x-=a.x,y-=a.y;return *this;}
	vec& operator *=(ld t){x*=t,y*=t;return *this;}
	vec& operator /=(ld t){x/=t,y/=t;return *this;}
};
typedef vec point;
vec operator +(vec a,vec b){return vec(a.x+b.x,a.y+b.y);}
vec operator -(vec a,vec b){return vec(a.x-b.x,a.y-b.y);}
vec operator *(vec a,ld t){return vec(a.x*t,a.y*t);}
vec operator /(vec a,ld t){return vec(a.x/t,a.y/t);}
ld inner(vec a,vec b){return a.x*b.x+a.y*b.y;}
ld outer(vec a,vec b){return a.x*b.y-a.y*b.x;}
struct circle
{
	point o;ld r;
	circle(){}
	circle(point _o,ld _r):o(_o),r(_r){}
};
struct line
{
	point a;vec v;
	line();
	line(point _a,vec _v):a(_a),v(_v){}
};
ld dist(point a,line l){return abs(outer(a-l.a,l.v))/l.v.norm();}
ld dist(point a,point b){return sqrt(sqr(a.x-b.x)+sqr(a.y-b.y));}
vec rotate(vec a,ld rad){return vec(a.x*cos(rad)-a.y*sin(rad),a.x*sin(rad)+a.y*cos(rad));}
point cross(line a,line b){return a.a+a.v*outer(a.a-b.a,b.v)/outer(b.v,a.v);}
pair<point,point> cross(line a,circle c)
{
	point p=cross(a,line(c.o,rotate(a.v,pi/2)));
	ld t=sqrt(sqr(c.r)-sqr(dist(c.o,a)))/a.v.norm();//t只是个长度 需要比norm才是比值。。naive
	return make_pair(p+a.v*t,p-a.v*t);
}
point inflex(line a,circle c)//有向直线的交点
{
	pair<point,point> res=cross(a,c);
	return fcmp(inner(a.v,res.first-a.a))<0?res.second:res.first;//a.v......
}
//ld sd(line a,line b){return min(dist(a.a,b.a),dist(a.a+a.v,b.a+b.v));} naive
ld sd(line a,line b,ld t)
{
	ld A=inner(a.v-b.v,a.v-b.v),B=2*inner(a.v-b.v,a.a-b.a),C=inner(a.a-b.a,a.a-b.a);
	ld axis=0;
	if(fcmp(A))
	{
		axis=-B/(2*A);
		if(fcmp(axis)<0) axis=0;
		else if(fcmp(axis-t)>0) axis=t;
	}
	else if(fcmp(B))
		axis=fcmp(B)<0?0:t;
	return sqrt(A*axis*axis+B*axis+C);
}
int K;
bool coincid(point a,point b){return fcmp(a.x-b.x)==0 && fcmp(a.y-b.y)==0;}
void turn(line &l,circle c)
{
	line axis=line(l.a,rotate(l.a-c.o,pi/2));
	line bot=line(l.a+l.v,c.o-l.a);
	point p=cross(axis,bot);
	l.v=l.a+l.v+(p-(l.a+l.v))*2-l.a;
}
ld calc(line li,line lj,circle c)
{
	ld res=INF;
	for(int i=0,j=0;i<K && j<K;)//已经碰撞了i j次
	{
		point pi=inflex(li,c),pj=inflex(lj,c);
		ld ti=dist(li.a,pi)/li.v.norm(),tj=dist(lj.a,pj)/lj.v.norm(),t; 
		chkmn(res,sd(li,lj,t=min(ti,tj)));
		if(coincid(li.a+=li.v*t,pi)) turn(li,c),i++;
		if(coincid(lj.a+=lj.v*t,pj)) turn(lj,c),j++;
	}
	return res;
}
int main()
{
	//freopen("input","r",stdin);
	static vec v[MAXN];
	static circle c;
	static point p[MAXN];
	fred(c.o.x),fred(c.o.y),fred(c.r);
	ld ANS=INF;int n;red(n),red(K);
	for(int i=1;i<=n;i++)
		fred(p[i].x),fred(p[i].y),fred(v[i].x),fred(v[i].y);
	for(int i=1;i<=n-1;i++)
		for(int j=i+1;j<=n;j++)
			chkmn(ANS,calc(line(p[i],v[i]),line(p[j],v[j]),c));
	printf("%.3lf\n",ANS);
	return 0;
}

BZOJ 3990 排序

发表于 2017-01-15 | 更新于 2018-06-18

Link

唉，还是too young。

Solution

操作序列的顺序不影响最终的结果。于是就枚举操作序列，最终乘上步数的阶乘。及时排除掉不可能排好的情况。

从小到大DFS，有点归并的感觉。

Code

//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;
using std::swap;
const int N=12,MAXN=20,MAXL=1<<13;
LL fac[MAXN];int a[MAXL],n;
void swap_n(int *x,int n,int *y)
{
	for(int i=1;i<=n;i++)
		swap(*x++,*y++);
}
bool check(int d)
{
	for(int i=1;i+d<=(1<<n);i+=d*2)
		if(a[i]+d!=a[i+d]) return 0;
	return 1;
}
LL DFS(int cur,int cnt)//进行完cnt，要进行cur
{
	if(cur==n) return fac[cnt];
	int stack[5],top=0;
	int d=(1<<cur);
	for(int i=1;i+d<=(1<<n);i+=d*2)
	{
		if(a[i]+d!=a[i+d])
		{
			if(top==4) return 0;
			stack[++top]=i;
			stack[++top]=i+d;
		}
	}
	if(top==0) return DFS(cur+1,cnt);
	LL res=0;
	if(top==2)
	{
		swap_n(a+stack[1],d,a+stack[2]);
		res=DFS(cur+1,cnt+1);
		swap_n(a+stack[1],d,a+stack[2]);
	}
	else
	{
		for(int i=1;i<=top-1;i++)
			for(int j=i+1;j<=top;j++)
			{
				swap_n(a+stack[i],d,a+stack[j]);
				if(check(d))
					res+=DFS(cur+1,cnt+1);
				swap_n(a+stack[i],d,a+stack[j]);
			}
	}
	return res;
}
int main()
{
	fac[0]=1;for(int i=1;i<=N;i++) fac[i]=fac[i-1]*i;
//	freopen("input","r",stdin);
	get(n);
	for(int i=1;i<=(1<<n);i++) get(a[i]);
	printf("%lld\n",DFS(0,0));
	return 0;
}
/* AC Record(Bugs)
 * WA
 * 看题解没研究清楚 实现不同
 */

BZOJ 1018 堵塞的交通

发表于 2017-01-15 | 更新于 2018-06-17

Link

Solution

线段树维护这几种可能性：

Tips

做过道路修建之后觉得这个很水。随便一写，WA。于是找到数据，发现把Y的判成了N。脑子一转，发现这道题和那道还不一样：这道的答案也有可能来自一端出去再回来。于是改，WA。发现把N判成了Y。看看代码，发现合并信息的时候少写了一个条件。于是改，WA。发现把Y判成了N。脑子一转，发现答案可能来自两端同时转出去再转回来。于是改，AC。汗。。细节问题一定要想好。。细节多的题一定要对拍。。对拍的时候一定要KISS，以覆盖到各种情况。。

Code

//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=100000+10;
using std::swap;
struct Data
{
	int a[2][2],tor[2],sf[2];
	int *operator [](int x){return a[x];}
	Data(){memset(a,0,sizeof(a));tor[0]=tor[1]=sf[0]=sf[1]=0;}
};
Data Merge(Data a,Data b)
{
	Data res;
	res[0][0]=(a[0][0]&b[0][0]&a.tor[0]) | (a[0][1]&b[1][0]&a.tor[1]);
	res[1][0]=(a[1][0]&b[0][0]&a.tor[0]) | (a[1][1]&b[1][0]&a.tor[1]);
	res[0][1]=(a[0][0]&b[0][1]&a.tor[0]) | (a[0][1]&b[1][1]&a.tor[1]);
	res[1][1]=(a[1][0]&b[0][1]&a.tor[0]) | (a[1][1]&b[1][1]&a.tor[1]);
	res.tor[0]=b.tor[0],res.tor[1]=b.tor[1];
	res.sf[0]=a.sf[0] | (b.sf[0] & a[0][0] & a.tor[0] & a[1][1] & a.tor[1]);
	res.sf[1]=b.sf[1] | (a.sf[1] & b[0][0] & a.tor[0] & b[1][1] & a.tor[1]);
	return res;
}
namespace iSeg
{
	const int SIZE=MAXN<<2;
	struct Node
	{
		Node *son[2];
		Data v;
		void up(){v=Merge(son[0]->v,son[1]->v);}
		Node(){son[0]=son[1]=0;}
	}*POOL=(Node*)malloc(SIZE*sizeof(Node)),*ME=POOL,*null=new(ME++) Node,*root=null;
	Node *newNode()
	{
		static Node *cur;cur=ME++;
		cur->son[0]=cur->son[1]=null;
		return cur;
	}
	void Build(Node *&pos,int L,int R)
	{
		pos=newNode();
		if(L==R)
			 pos->v[0][0]=pos->v[1][1]=1;
		else
		{
			int MID=(L+R)>>1;
			Build(pos->son[0],L,MID);
			Build(pos->son[1],MID+1,R);
			pos->up();
		}
	}
	Data Access(Node *pos,int L,int R,int l,int r)
	{
		if(L==l && R==r) return pos->v;
		else
		{
			int MID=(L+R)>>1;
			if(r<=MID) return Access(pos->son[0],L,MID,l,r);
			else if(MID+1<=l) return Access(pos->son[1],MID+1,R,l,r);
			else return Merge(Access(pos->son[0],L,MID,l,MID),Access(pos->son[1],MID+1,R,MID+1,r));
		}
	}
	void Edit(Node *pos,int L,int R,int x1,int x2,int y1,int y2,int d)
	{
		if(L==R)
		{
			if(y1==y2) pos->v[0][1]=pos->v[1][0]=pos->v.sf[0]=pos->v.sf[1]=d;
			else pos->v.tor[x1]=d;
		}
		else
		{
			int MID=(L+R)>>1;
			if(y1<=MID) Edit(pos->son[0],L,MID,x1,x2,y1,y2,d);
			else Edit(pos->son[1],MID+1,R,x1,x2,y1,y2,d);
			pos->up();
		}
	}
}using namespace iSeg;
int main()
{
	int n;get(n);
	Build(root,1,n);
	char opt[10];int x1,y1,x2,y2,lc,rc;Data res;
	while(scanf("%s",opt),opt[0]!='E')
	{
		get(x1),get(y1),get(x2),get(y2);
		x1--,x2--;
		if(y1>y2 || (y1==y2 && x1>x2)) swap(x1,x2),swap(y1,y2);
		switch(opt[0])
		{
			case 'A':
				//static int cnt;if(++cnt==91) puts("");
				res=Access(root,1,n,y1,y2);
				lc=Access(root,1,n,1,y1).sf[1];
				rc=Access(root,1,n,y2,n).sf[0];
				puts(res[x1][x2] || (lc && res[x1^1][x2]) || (rc && res[x1][x2^1]) || (lc && rc && res[x1^1][x2^1])?"Y":"N");
				break;
			case 'O':Edit(root,1,n,x1,x2,y1,y2,1);break;
			case 'C':Edit(root,1,n,x1,x2,y1,y2,0);break;
		}
	}
	return 0;
}
/* AC Record(Bugs) *
 * 
 * * * * * * * * * */

BZOJ 1295 最长距离

发表于 2017-01-15 | 更新于 2018-06-17

神思路啊！！

Link

Solution

最大值最大？？根本没有思路。 $x+y=c$ , $x^2+y^2=c^2-2xy$ ？走扁矩形的对角线求曼哈顿距离的最长路再转换？？换个思路，变成枚举。 $O(n^2)$ 地枚举所有点对，看它们的距离能否成为答案。如何判断能否成为答案？？最短路求出最少经过多少障碍点与 $T$ 比较即可。因为偷懒用了Shortest Path Faster Algorithm，所以复杂度为 $O(nke)$ ( $n$ 为点数)，其中 $k\in(2,n)$ 。

Tips

暴力出奇迹。用给出的题目的性质找不到答案，就寻找答案的性质(单峰->三分/单调->二分/总量不大->枚举所有可能性)，然后变成判定问题。

Code

//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;}
using std::pair;
using std::queue;
const int MAXN=30+5,d[][2]={{0,1},{0,-1},{1,0},{-1,0}};
template<class T> inline T sqr(T x){return x*x;}
int map[MAXN][MAXN],sp[MAXN][MAXN][MAXN][MAXN],n,m;
void BFS(int x,int y,int sp[MAXN][MAXN])
{
	typedef pair<int,int> Node;
	static int inq[MAXN][MAXN],ut;
	static queue<Node> Q;
	sp[x][y]=map[x][y];Q.push(Node(x,y));inq[x][y]=++ut;
	while(!Q.empty())
	{
		x=Q.front().first,y=Q.front().second;inq[x][y]--;Q.pop();
		for(int p=0;p<4;++p)
		{
			int dx=x+d[p][0],dy=y+d[p][1];
			if(1<=dx && dx<=n && 1<=dy && dy<=m 
			   && chkmn(sp[dx][dy],sp[x][y]+map[dx][dy]) && inq[dx][dy]!=ut)
				Q.push(Node(dx,dy)),inq[dx][dy]=ut;
		}
	}
}
int main()
{
	freopen("input","r",stdin);
	memset(sp,31,sizeof(sp));
	int T,c;get(n),get(m),get(T);
	for(int i=1;i<=n;i++)
		for(int j=1;j<=m;j++)
		{
			while(!isdigit(c=getchar()));
			map[i][j]=c-'0';
		}
	for(int i=1;i<=n;i++)
		for(int j=1;j<=m;j++)
			BFS(i,j,sp[i][j]);
	double ANS=0;
	for(int x1=1;x1<=n;x1++)
		for(int y1=1;y1<=m;y1++)
			for(int x2=1;x2<=n;x2++)
				for(int y2=1;y2<=m;y2++)
					if(sp[x1][y1][x2][y2]<=T)
						chkmx(ANS,sqrt(sqr(x1-x2)+sqr(y1-y2)));
	printf("%.6lf\n",ANS);
	return 0;
}
/* AC Record(Bugs) *
 * 
 * * * * * * * * * */

BZOJ 1951 古代猪文

发表于 2017-01-15 | 更新于 2018-06-18

Link

Solution

求首先，根据费马小定理，当 $p$ 为质数， $gcd(a,p)==1$ 时，这个模数是质数，很好！当时， $G==999911659$ ，答案为 $0$ ，需要特判掉。否则，只要求出 $\sum_{k|N}C_N^K 999911658$ 问题就解决了。注意到 $999911658$ 不是个质数，但它的分解是四个一次的质数相乘，可以直接用Lucas定理，再用CRT合并。

Tips

不管是做数学题还是做OI题，在分析的时候需要有严谨的思维，各个定理的适用条件一定要熟知。数论题能打表的范围内如果可以的话一定要打表，优化的效果不是盖的。

Code

//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 Pow(int base,int v,int P)
{
	int res=1;
	while(v)
	{
		if(v&1)
			res=(LL)res*base%P;
		base=(LL)base*base%P;
		v>>=1;
	}
	return res;
}
void exgcd(int a,int b,LL &x,LL &y)
{
	if(!b){x=1,y=0;}
	else
	{
		exgcd(b,a%b,x,y);
		LL nx=y,ny=x-(a/b)*y;
		x=nx,y=ny;
	}
}
int inv(int a,int n)
{
	LL x,y;
	exgcd(a,n,x,y);
	((x%=n)+=n)%=n;
	return x;
}

namespace iLucas
{
	const int dv[]={0,2,3,4679,35617},dc=4,tM[]={0,499955829,333303886,289138806,877424796};
	const int MODU=999911658,SIZE=1e5;
	int table[5][SIZE];
	void Init()
	{
		for(int i=1,*fac;i<=dc;i++)
		{
			fac=table[i];fac[0]=1;
			for(int j=1;j<=dv[i];j++)
				fac[j]=(LL)fac[j-1]*j%dv[i];
		}
	}
	int __C(int n,int m,int idx)
	{
		if(n<m) return 0;
		static int P,*fac;
		P=dv[idx],fac=table[idx];
		return (LL)fac[n]*inv(fac[m],P)%P*inv(fac[n-m],P)%P;
	}
	int Lucas(int n,int m,int idx)
	{
		int res=1,P=dv[idx];
		while(n && m)
		{
			res=(LL)res*__C(n%P,m%P,idx)%P;
			n/=P,m/=P;
		}
		return res;
	}
	int A[5];
	int CRT()
	{
		int res=0;
		for(int i=1;i<=dc;i++) 
			res=((LL)A[i]*tM[i]%MODU+res)%MODU;
		return res;
	}
	int C(int n,int m)
	{
		
		for(int i=1;i<=dc;i++) A[i]=Lucas(n,m,i);
		return CRT();
	}
}
int C(int n,int m)
{
	if(n<m) return 0;
	return iLucas::C(n,m);
	
}
int Calc(int n)
{
	const int P=999911658;
	int res=0;
	for(int i=1;i*i<=n;i++)
	{
		if(n%i==0)
		{
			res=((LL)res+C(n,i))%P;
		if(i*i!=n)
			res=((LL)res+C(n,n/i))%P;
		}
	}
	return res;
}
int main()
{
	const int MODU=999911659;
	iLucas::Init();
	int n,G;get(n),get(G);
	if(G==MODU) puts("0");
	else printf("%d\n",Pow(G,Calc(n),MODU));
	return 0;
}
/* AC Record(Bugs) *
 * TLE fac table
 * WA if i*i!=n!!!
 * WA if n==G!!!
 * * * * * * * * * */

BZOJ 3199 Escape

发表于 2017-01-15 | 更新于 2018-06-17

Link

Solution

可以发现每个亲戚管辖的区域是由与其它亲戚连线的中垂线分界的半平面，分界线上是两人管辖，分界线交点处是3个及以上的人管辖。如果不考虑交点的话，就是把相邻的半平面连一条权为1的边，跑最短路。然后就被交点弄晕了。去翻题解，发现根本没有人单独考虑了交点的情况。换个思路一想，因为贪心地不走交点至少不会更差，而且也不会出现只能走交点的情况。剩下就好说了。

Tips

建立半平面的时候方向随便选的，查了一节课；onleft写错，查了一节课。还是要练习静态查错被边界情况困扰的时候，想想这种边界是不是必须要被考虑（同SDOI2015道路修建）。

Code

//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=600+10,MAXM=MAXN*MAXN;
using std::sort;
using std::queue;

typedef double ld;
const ld eps=1e-6;
const ld pi=acos(-1.0);
inline int fcmp(ld x)
{
	if(-eps<x && x<eps) return 0;
	return x<0?-1:1;
}
struct vec
{
	ld x,y;
	vec(){}
	vec(ld x,ld y):x(x),y(y){}
};
typedef vec poi;
vec operator +(vec a,vec b){return vec(a.x+b.x,a.y+b.y);}
vec operator -(vec a,vec b){return vec(a.x-b.x,a.y-b.y);}
vec operator *(vec a,ld b){return vec(a.x*b,a.y*b);}
vec operator /(vec a,ld b){return vec(a.x/b,a.y/b);}
ld inner(vec a,vec b){return a.x*b.x+a.y*b.y;}
ld outer(vec a,vec b){return a.x*b.y-a.y*b.x;}
vec Rotate(vec a,ld rad){return vec(a.x*cos(rad)-a.y*sin(rad),a.x*sin(rad)+a.y*cos(rad));}
struct line
{
	poi a;vec v;int id;ld ang;
	line(){}
	line(poi a,vec v):a(a),v(v){ang=atan2(v.y,v.x);}
	bool operator <(const line& a) const 
	{
		return ang<a.ang;//atan2(v.y,v.x)<atan2(a.v.y,a.v.x);
	}
};
poi p[MAXN];line dvd[MAXN];int dcnt;
line NormalBisector(poi a,poi b)
{
	return line((a+b)/2,Rotate(b-a,pi/2));
}
int width,height;
void BuildEdge(int cur,int n)
{
	dcnt=0;
	for(int i=1;i<=n;i++)
		if(i!=cur)
		{
			dvd[++dcnt]=NormalBisector(p[cur],p[i]);
			dvd[dcnt].id=i;
		}
	dvd[++dcnt]=line(poi(0,0),vec(1,0));dvd[dcnt].id=0;
	dvd[++dcnt]=line(poi(0,height),vec(0,-1));dvd[dcnt].id=0;
	dvd[++dcnt]=line(poi(width,0),vec(0,1));dvd[dcnt].id=0;
	dvd[++dcnt]=line(poi(width,height),vec(-1,0));dvd[dcnt].id=0;
}
line Ql[MAXN];poi Qp[MAXN];int he,ta;
bool onleft(poi a,line b){return fcmp(outer(a-b.a,b.v))<0;}
bool paral(vec a,vec b){return fcmp(outer(a,b))==0;}
bool paral(line a,line b){return paral(a.v,b.v);}
poi Cross(line a,line b)
{
	ld t=outer(a.a-b.a,b.v)/outer(b.v,a.v);
	return a.a+a.v*t;
}
ld Dist(poi a,poi b){return sqrt(inner(a-b,a-b));}
bool online(poi a,line l){return paral(a-l.a,l.v);}
void Intersect()
{
	sort(dvd+1,dvd+1+dcnt);
	#define dist(a,b) (b-a+1)
	Ql[he=ta=1]=dvd[1];
	for(int i=2;i<=dcnt;i++)
	{
		while(dist(he,ta)>=2 && !onleft(Qp[ta-1],dvd[i])) 
			ta--;
		while(dist(he,ta)>=2 && !onleft(Qp[he],dvd[i])) 
			he++;//没写。。:w
		Ql[++ta]=dvd[i];
		if(paral(Ql[ta-1],Ql[ta]))
		{
			ta--;
			//if((!onleft(Qp[ta-1],dvd[i]))!=(onleft(dvd[i].a,Ql[ta]))) assert((!onleft(Qp[ta-1],dvd[i]))==(onleft(dvd[i].a,Ql[ta])));
			//if(!onleft(Qp[ta-1],dvd[i]))交点不一定存在。 
			if(onleft(dvd[i].a,Ql[ta]))
				Ql[ta]=dvd[i];
		}
		if(dist(he,ta)>=2)//需要判断
			Qp[ta-1]=Cross(Ql[ta],Ql[ta-1]);
	}
	while(dist(he,ta)>=2 && !onleft(Qp[ta-1],Ql[he])) ta--; 
	#undef dist
}
struct Edge{int to,pre,v;}e[MAXM<<1];int id[MAXN],ec;
void Adde(int f,int t,int v)
{
	e[++ec].to=t;e[ec].v=v;e[ec].pre=id[f];id[f]=ec;
}
int SP(int S)
{
	static queue<int> Q;
	static int dis[MAXN];static bool inq[MAXN];
	memset(dis,31,sizeof(dis));memset(inq,0,sizeof(inq));
	Q.push(S);inq[S]=1;dis[S]=0;
	while(!Q.empty())
	{
		int cur=Q.front();Q.pop();inq[cur]=0;
		for(int i=id[cur];i;i=e[i].pre)
		{
			int u=e[i].to;
			if(chkmn(dis[u],dis[cur]+e[i].v) && !inq[u])
				Q.push(u),inq[u]=1;
		}
	}
	return dis[0];
}
void CLEAR()
{
	memset(id,0,sizeof(id));ec=0;
}
void WORK()
{
	CLEAR();
	int n,x0,y0;get(n);
	get(width),get(height),get(x0),get(y0);
	if(n==0) {puts("0");return;}
	poi S=poi(x0,y0);int block=1;
	for(int i=1;i<=n;i++)
	{
		int x,y;get(x),get(y);
		p[i].x=x,p[i].y=y;
		if(fcmp(Dist(S,p[block])-Dist(S,p[i]))>0) 
			block=i;
	}
	for(int i=1;i<=n;i++)
	{
		BuildEdge(i,n);
		Intersect();
		for(line *ptr=Ql+he;ptr<=Ql+ta;++ptr)
			Adde(i,ptr->id,1);
	}
	printf("%d\n",SP(block));
}
int main()
{
	int T;get(T);
	while(T--)WORK();
	return 0;
}
/* AC Record(Bugs) *
 * WA atan2 
 * * * * * * * * * */

BZOJ 3533 向量集

发表于 2017-01-15 | 更新于 2018-06-18

Link

Solution

好一道数据结构套计算几何题啊。。（虽然是基础数据结构套基础计算几何。。但对我这种蒟蒻就已经有超强杀伤力了。。）维护一个支持插入的向量序列，输出区间中的向量与给定向量最大的内积。首先很容易看出答案向量的终点在凸包上。还有一个看起来对但不会证的结论：给定向量与终点在凸包上的向量内积值 $v$ 是单峰的。 $y>0$ ， $v$ 在上凸壳上有最大值 $y<0$ ， $v$ 在下凸壳上有最大值于是用线段树维护序列，每个节点维护区间的凸包。但是凸包的合并类似Graham(Andrew)，是 $O(n)$ 的，每次都合并肯定合不起。但是：一个区间只有塞满了之后才会被查询。于是在区间塞满的时候再合并。这样合并的总复杂度是 $O(n\log n)$ 的。查询的时候在凸包上三分。一开始想得是要把查询的区间在线段树上对应的 $\log n$ 个区间的凸包合并出来得到整个区间的凸包再三分。单次查询的复杂度就成了 $O(n\log n)$ ? 然后发现网上给出的做法都是分别在每个节点上三分，然后取最大值， $O(\log^2 n)$ 解决问题。

Tips

如果信息合并的复杂度高，就根据查询的情况来减少合并的次数。这个在每个节点上三分的做法可以类比区间第K大的线段树套平衡树做法，只是这个可能会产生一些无用的解，但不影响正确性和复杂度。

Code

//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;
const int MAXN=4e5+10;

using std::vector;
using std::max;
struct point
{
	int x,y;
	point(){}
	point(int x,int y):x(x),y(y){}
	bool operator <(const point &a)const
	{
		if(x!=a.x) return x<a.x;
		return y<a.y;
	}
};
point operator -(const point &a,const point &b){return point(a.x-b.x,a.y-b.y);}
LL inner(const point &a,const point &b){return (LL)a.x*b.x+(LL)a.y*b.y;}
LL outer(const point &a,const point &b){return (LL)a.x*b.y-(LL)a.y*b.x;}
typedef vector<point> hull;
const int UP_HULL=0,DOWN_HULL=1;
void Merge(hull& a,hull& b,hull &dest,int kind)
{
	static point stack[MAXN];int top=0;
	int na=a.size()-1,nb=b.size()-1;
	dest.clear();
	for(int pa=0,pb=0;pa<=na || pb<=nb;)
	{
		point cur;
		if(pa>na) cur=b[pb++];
		else if(pb>nb) cur=a[pa++];
		else cur=a[pa]<b[pb]?a[pa++]:b[pb++];
		while(top>=2 && outer(stack[top]-stack[top-1],cur-stack[top])*(kind==UP_HULL?1:-1)>=0) top--;
		stack[++top]=cur;
	}
	for(int i=1;i<=top;i++) dest.push_back(stack[i]);
}
LL Trisect(hull& a,const point& p)
{
	int L=0,R=a.size()-1,M1,M2;
	while(L!=R)
	{
		int M1=L+(R-L)/3,M2=R-(R-L)/3;
		if(inner(a[M1],p)>inner(a[M2],p))
			R=M2-1;
		else
			L=M1+1;
	}
	return inner(a[L],p);
}
namespace iSeg
{
	const int SIZE=MAXN<<2;
	struct Node
	{
		hull convex[2];
		Node *son[2];
		Node(Node *lc=0,Node *rc=0)
		{
			son[0]=lc,son[1]=rc;
			convex[0].clear(),convex[1].clear();
		}
		void up()
		{
			Merge(son[0]->convex[0],son[1]->convex[0],convex[0],0);
			Merge(son[0]->convex[1],son[1]->convex[1],convex[1],1);
		}
	}*POOL=(Node*)malloc(SIZE*sizeof(Node)),*ME=POOL,*null=new(ME++) Node;
	Node *newNode(){return new(ME++) Node(null,null);}
	struct Seg
	{
		Node *root;
		Seg(){root=null;}
		int TL,TR;
		void Init(int l,int r){TL=l;TR=r;}
		void Edit(Node *&pos,int L,int R,int goal,const point &p)
		{
			if(pos==null) pos=newNode();
			if(L==R)
			{
				pos->convex[UP_HULL].push_back(p);
				pos->convex[DOWN_HULL].push_back(p);
			}
			else
			{
				int MID=(L+R)>>1;
				if(goal<=MID) Edit(pos->son[0],L,MID,goal,p);
				else Edit(pos->son[1],MID+1,R,goal,p);
				if(goal==R) pos->up();
			}
		}

		LL Access(Node *pos,int L,int R,int l,int r,const point &p)//一个区间只会合并一次 n log n
		{
			assert(pos!=null);
			if(L==l && R==r) return Trisect(pos->convex[p.y>0?UP_HULL:DOWN_HULL],p);
			else
			{
				int MID=(L+R)>>1;
				if(r<=MID) return Access(pos->son[0],L,MID,l,r,p);
				else if(MID+1<=l) return Access(pos->son[1],MID+1,R,l,r,p);
				else return max(Access(pos->son[0],L,MID,l,MID,p),Access(pos->son[1],MID+1,R,MID+1,r,p));
			}
		}
		void Add(point p)
		{
			static int cnt=0;
			++cnt;
			Edit(root,TL,TR,cnt,p);
		}
		LL Answer(int l,int r,point p){return Access(root,TL,TR,l,r,p);}
	}seg;
}using namespace iSeg;


int main()
{
	char dkind;int n;
	get(n);while(!isalpha(dkind=getchar()));
	seg.Init(1,n);
	#define decode(x) ( x=(dkind=='E')?x:(x^(last & 0x7fffffff)) )
	LL last=0;
	for(int i=1;i<=n;i++)
	{
		char opt;int x,y,l,r;
		while(!isalpha(opt=getchar()));
		get(x),decode(x),get(y),decode(y);
		if(opt=='A') seg.Add(point(x,y));
		else
		{
			get(l),decode(l),get(r),decode(r);
			printf("%lld\n",last=seg.Answer(l,r,point(x,y)));
		}
	}
	return 0;
}
/* AC Record(Bugs) *
 * 
 * * * * * * * * * */

胡扯扩展Floyd算法

发表于 2017-01-15 | 更新于 2018-06-17

Problem

给定一个有向图，求一个至少经过2个点的最小环。

Solution

Floyd算法的代码如下

for(int k=1;k<=n;k++)
	for(int i=1;i<=n;i++)
		for(int j=1;j<=n;j++)
        	chkmn(sp[i][j],sp[i][k]+sp[k][j]);

外层循环执行完 $k==X$ 的循环，表示两两节点之间只走 $1...X$ 的最短路已经被求出。考虑一个环，上面一定有一个编号最大的点以上图为例，在第66次循环之前，所有点对经过 $1...65$ 的最短路已经求出，由于最短路中一定没有点 $66$ ，所以可以放心地尝试把 $66$ 号点加进去而不必担心出现重边。于是算法就出来了：额外记录原图的边 $edge[][]$ ，在 $sp[][]$ 上求最短路。第 $k$ 次循环之前，枚举每对点 $i,j< k$ ，把 $sp[i][j]$ 连到 $edge[j][k]$ 和 $edge[k][i]$ 上，用这个环的大小来更新答案。

Code

for(int k=1;k<=n;k++)
{
    for(int i=1;i<k;i++)
        for(int j=i+1;j<k;j++)
            chkmn(res,sp[i][j]+dist[j][k]+dist[k][i]);
    for(int i=1;i<=n;i++)
        for(int j=1;j<=n;j++)
            chkmn(sp[i][j],sp[i][k]+sp[k][j]);
}

BZOJ 1390 Fence

发表于 2017-01-15 | 更新于 2018-06-18

Link

Solution

好题！首先，围栏的凸包之外的树一定是圈不到的，就不管它们，把答案加上。而凸包里面的是一定要圈的，因为要圈一个最多得多建立3个栏杆，建栏杆的花费肯定比舍弃树的代价小。然而，圈凸包不一定是最优的，最优的一定是能圈起所有凸包内树的最小的环。下一步转化就非常厉害了：对于环上的任意一条边，一定满足凸包内的树都在它的同侧。那就枚举每对围栏（有向，即 $p_i,p_j$ 需要枚举和），用叉积判断是否所有的树都在它的左侧。如果是，就把 $p_i$ 到 $p_j$ 连一条长度为 $1$ 的边。最后用扩展Floyd求最小环即可。

Tips

这个解法依然是通过观察答案性质得到的。

Code

//Code by Lucida
#include<bits/stdc++.h>
//#define red(x) scanf("%d",&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::sort;
using std::swap;
const int MAXN=100+10,INF=0x1f1f1f1f;
struct vec
{
	int x,y;
	vec(){}
	vec(int _x,int _y):x(_x),y(_y){}
	void operator +=(vec a){x+=a.x,y+=a.y;}
	void operator -=(vec a){x-=a.x,y-=a.y;}
};
typedef vec point;
vec operator -(vec a,vec b){return vec(a.x-b.x,a.y-b.y);}
int inner(vec a,vec b){return a.x*b.x+a.y*b.y;}
int outer(vec a,vec b){return a.x*b.y-a.y*b.x;}
bool cmp(point a,point b)
{
	if(outer(a,b)!=0)
		return outer(a,b)>0;
	return inner(a,a)<inner(b,b);
}
int Graham(point *p,int n,point *hull)
{
	int mi=1;
	for(int i=1;i<=n;i++)
		if(p[i].x<p[mi].x || (p[i].x==p[mi].x && p[i].y<p[mi].y)) mi=i;
	swap(p[1],p[mi]),p[0]=p[1];
	for(int i=1;i<=n;i++) p[i]-=p[0];
	sort(p+2,p+1+n,cmp);
	for(int i=1;i<=n;i++) p[i]+=p[0];
	p[n+1]=p[1];
	static point S[MAXN];int top=0;
	for(int i=1;i<=n+1;i++)
	{
		while(top>=2 && outer(p[i]-S[top],S[top]-S[top-1])>0) top--;//>=0
		S[++top]=p[i];
	}
	top--;for(int i=1;i<=top;i++) hull[i]=S[i];
	return top;
}
bool cover(point a,point b,point c){return outer(b-a,c-a)==0 && (inner(c-a,c-a)<inner(b-a,b-a));}
bool cover(point *p,int n,point a)
{
	int wn=0;
	#define suc(i) (i+1>n?1:i+1)
	for(int i=1;i<=n;i++)
	{
		if(cover(p[suc(i)],p[i],a)) 
			return 1;
		if(outer(p[suc(i)]-p[i],a-p[i])>0 && a.y>p[i].y && a.y<=p[suc(i)].y) 
			wn++;
		if(outer(p[suc(i)]-p[i],a-p[i])<0 && a.y>p[suc(i)].y && a.y<=p[i].y) 
			wn--;
	}
	#undef suc
	return wn!=0;
}
bool check(point a,point b,point *t,int n)
{
	for(int i=1;i<=n;i++)
		if(outer(b-a,t[i]-a)<0) return 0;
	return 1;	
}
int main()
{
	static point p[MAXN],tr[MAXN],t[MAXN],hull[MAXN];
	//freopen("input","r",stdin);
	int n,m;red(n),red(m);
	for(int i=1;i<=n;i++)
		red(p[i].x),red(p[i].y);
	for(int i=1;i<=m;i++)
		red(tr[i].x),red(tr[i].y);
	int hc=Graham(p,n,hull);
	int tc=0,ANS=0;
	for(int i=1;i<=m;i++)
		if(cover(hull,hc,tr[i]))
			t[++tc]=tr[i];
		else
			ANS+=111;
	if(tc)
	{
		static int sp[MAXN][MAXN],dist[MAXN][MAXN];
		memset(sp,31,sizeof(sp));
		for(int i=1;i<=n;i++)//有向。。。
			for(int j=1;j<=n;j++)
				if(i!=j && check(p[i],p[j],t,tc))
					sp[i][j]=1;
		for(int i=1;i<=n;i++) sp[i][i]=0;
		memcpy(dist,sp,sizeof(sp));
		int res=INF;
		for(int k=1;k<=n;k++)
		{
			for(int i=1;i<k;i++)
				for(int j=i+1;j<k;j++)
					chkmn(res,sp[i][j]+dist[j][k]+dist[k][i]);
			for(int i=1;i<=n;i++)
				for(int j=1;j<=n;j++)
					chkmn(sp[i][j],sp[i][k]+sp[k][j]);
		}
		ANS+=res*20;
	}
	printf("%d\n",ANS);
	return 0;
}

Lucida

It's been a while.

222 日志

127 标签