Comet OJ - Contest #13 (D,E,F)

D - 「火鼠的皮衣 -不焦躁的内心-」

给定$n,a,b,p$，求下面的式子；
$$\sum _{i=0}^{\lfloor \frac{n}{2}\rfloor }{a}^i{b}^{n-2i}\left(\genfrac{}{}{0ex}{}{n}{2i}\right)\mathrm{m}\mathrm{o}\mathrm{d}p$$

数据组数$T\le 10000$，$1\le n,a,b,p\le 10^{18}$，不保证$p$是质数。

Sol

设$a^{\prime }=\sqrt{a}(\mathrm{m}\mathrm{o}\mathrm{d}p)$，则原式可以化成：
$$\sum _{i=0}^n{a}^{\prime i}{b}^{n-i}\left(\genfrac{}{}{0ex}{}{n}{i}\right)[\text{i is even}]$$

也就是$(a^{\prime }+b{)}^n$的展开式中，$a^{\prime }$的次数为偶数的项的和。

设$f_{i,0/1}$表示$(a^{\prime }+b{)}^i$中$a^{\prime }$次数为偶数/奇数的项的和，那么$f_{i+1,d}=f_{i,d}\cdot b+f_{i,1-d}\cdot a^{\prime }$。

矩阵快速幂就可以求出$f_{n,0}$。

由于$a$在模$p$的意义下可能不存在二次剩余。但是我们在运算中会用到的数一定都可以表示成$x+a^{\prime }\cdot y$的形式，记录下$(x,y)$就可以进行加、减、乘运算。最后得到的答案$y$一定等于$0$，输出$x$就可以了。

还有一个问题是，这里如果用long double的那种$O(1)$快速乘精度会炸，$O(\log p)$的慢速乘会TLE，要用int128来做乘法。

#include <cstdio>
#include <iostream>
#include <algorithm>
#include <cstring>
#define ll long long
using namespace std;
template <class T>
inline void rd(T &x) {
	x=0; char c=getchar(); int f=1;
	while(!isdigit(c)) { if(c=='-') f=-1; c=getchar(); }
	while(isdigit(c)) x=x*10-'0'+c,c=getchar(); x*=f;
}
ll mod,W,n,a,b;
ll add(ll x) { return x>=mod?x-mod:x; }
ll dec(ll x) { return x<0?x+mod:x; }
ll Mul(ll x,ll y) {
	return (__int128_t)x*y%mod;
//	if(!y||(1ll<<62)/y>=x) return x*y%mod;
//	ll tmp=(long double)x*y/mod;
//	return ((x*y-tmp*mod)%mod+mod)%mod;
//	ll res=0;
//	while(y) {
//		if(y&1) res=add(res+x);
//		x=add(x+x),y>>=1;
//	}
//	return res;
}
struct Num {
	ll x,y;
	Num(ll x=0,ll y=0): x(x),y(y) {}
	friend Num operator +(Num A,Num B) { return Num(add(A.x+B.x),add(A.y+B.y)); }
	friend Num operator -(Num A,Num B) { return Num(dec(A.x-B.x),dec(A.y-B.y)); }
	friend Num operator *(Num A,Num B) { return Num(add(Mul(A.x,B.x)+Mul(W,Mul(B.y,A.y))),add(Mul(A.x,B.y)+Mul(A.y,B.x))); }
};
struct Mat {
	Num a[2][2];
	Mat() { memset(a,0,sizeof(a)); }
	void print() { for(int i=0;i<2;++i,puts("")) for(int j=0;j<2;++j) printf("(%lld,%lld) ",a[i][j].x,a[i][j].y); }
};
Mat MulMat(Mat A,Mat B) {
	Mat C;
	for(int i=0;i<2;++i)
		for(int k=0;k<2;++k) if(A.a[i][k].x||A.a[i][k].y)
			for(int j=0;j<2;++j)
				C.a[i][j]=C.a[i][j]+(A.a[i][k]*B.a[k][j]);
	return C;
}
Mat Pow(Mat x,ll y) {
	Mat res; for(int i=0;i<2;++i) res.a[i][i]=Num(1,0);
	while(y) {
		if(y&1) res=MulMat(res,x);
		x=MulMat(x,x),y>>=1;
	}
	return res;
}
int main() {
	int T; rd(T);
	while(T--) {
		rd(n),rd(a),rd(b),rd(mod);
		W=a;
		Mat M;
		for(int i=0;i<2;++i) for(int j=0;j<2;++j) M.a[i][j]=i==j?Num(b,0):Num(0,1);
		M=Pow(M,n);
		printf("%lld\n",M.a[0][0].x);
	}
	return 0;
}

---

E - 「燕的子安贝 -永命线-」

平面上有$n$个点。你现在要画一条直线，使得到这条直线的距离不超过$d$的点尽可能多。

输出可能的最多的点数。

$n\le 2000$，坐标的绝对值不超过$10000$，$0\le d\le 10000$

Sol

转化一下问题，以每个点为圆心，画半径为$d$的圆，问画一条直线至多能够交多少个圆。

那么画的直线一定与某个圆相切。枚举直线与哪个圆相切，考虑直线绕这个圆周旋转且始终与圆相切的过程中，其它的圆与这个直线相交的时间是连续的区间，求出被区间覆盖次数最大的点就可以了。

我的代码写丑了，精度有点炸，特判了$d=0$才过的。

时间复杂度$O(n^2\log n)$。

#include <cstdio>
#include <iostream>
#include <algorithm>
#include <cstring>
#include <cmath>
#include <vector>
#define PB push_back
#define ll long long
#define db long double
using namespace std;
template <class T>
inline void rd(T &x) {
	x=0; char c=getchar(); int f=1;
	while(!isdigit(c)) { if(c=='-') f=-1; c=getchar(); }
	while(isdigit(c)) x=x*10-'0'+c,c=getchar(); x*=f;
}
const int N=2010;
db eps=1e-12,Pi=acos(-1.0);
int sgn(db x) { return x<-eps?-1:(x>eps); }
struct Point {
	db x,y;
	Point(db x=0,db y=0): x(x),y(y) {}
	friend Point operator +(Point A,Point B) { return Point(A.x+B.x,A.y+B.y); }
	friend Point operator -(Point A,Point B) { return Point(A.x-B.x,A.y-B.y); }
	friend Point operator *(Point A,db B) { return Point(A.x*B,A.y*B); }
}p[N];
int n,d;
struct item {
	db ag; int t;
	item(db ag=0,int t=0): ag(ag),t(t) {};
	friend bool operator <(item A,item B) { return sgn(A.ag-B.ag)?A.ag<B.ag:A.t>B.t; }
}que[N*8];
int tot;
void add(db l,db r,int t) {
	if(sgn(r-l)>=0) que[++tot]=item(l,t),que[++tot]=item(r,-t);
	else add(l,Pi,t),add(-Pi,r,t);
}
db dis2(Point A) { return A.x*A.x+A.y*A.y; }
vector<db> vec;
void sol() {
	int ans=0;
	for(int i=1;i<=n;++i) {
		vec.clear();
		int cnt0=0;
		for(int j=1;j<=n;++j) if(i!=j) {
			if(sgn(p[i].x-p[j].x)==0) cnt0++;
			else vec.PB((p[j].y-p[i].y)/(p[j].x-p[i].x));
		}
		ans=max(ans,cnt0+1);
		sort(vec.begin(),vec.end());
		for(int l=0,r;l<vec.size();l=r+1) {
			r=l;
			while(r+1<vec.size()&&sgn(vec[r+1]-vec[l])==0) r++;
			ans=max(ans,r-l+2);
		}
	}
	printf("%d",ans);
}
int main() {
	int ans=0;
	rd(n),rd(d);
	for(int i=1;i<=n;++i) cin>>p[i].x>>p[i].y;
	if(d==0) {
		sol();
		return 0;
	}
	for(int i=1;i<=n;++i) {
		tot=0;
		for(int j=1;j<=n;++j) if(j!=i) {
			Point A=p[i],B=p[j];
			db l=atan2(-(B.x-A.x),B.y-A.y),r=l+Pi; if(r>Pi) r-=2*Pi;
			
			if(sgn(dis2(A-B)-4ll*d*d)>0) {
				db d1=sqrt(dis2((B-A)*0.5));
				db ag0=atan2(B.y-A.y,B.x-A.x);
				db delt=acos(d/d1);
				db a1=ag0-delt; if(-Pi>a1) a1+=2*Pi;
				db a2=ag0+delt; if(a2>Pi) a2-=2*Pi;
				add(l,a1,1),add(a2,r,1);
			}
			else add(l,r,1);
		}
		
		sort(que+1,que+tot+1);
		int cnt=0;
		for(int j=1;j<=tot;++j) cnt+=que[j].t,ans=max(ans,cnt+1);
	}
	printf("%d",ans);
	return 0;
}

---

F - 「蓬莱的弹枝 -七色的弹幕-」

你要维护一个序列$a$，支持以下操作：

1）将区间$[l,r]$内的数向左偏移，具体来说就是$a_i=a_{i+1}(l\le i<r)$，$a_r=a_l$。
2）将区间$[l,r]$内的数全部$+1$。
3）给出一个下标$x$，查询与$a_x$值相同的元素的下标与$x$的差的绝对值的最小值。

$n\le 10^5,q\le 10^5$，初始$a_i\le 10^5$。

Sol

对序列分块，对每个块维护：

• 加标记$tag$。
• 每个元素的真实值-加标记的值$a[]$。
• 真实值-加标记的每个值的出现次数$f[]$。
• 下标偏移值：对这个块记一个下标$s$，表示序列中这一块的元素从左到右依次是$a[s],a[s+1]\cdots a[m],a[1],\cdots a[s-1]$，其中$m$是块大小。

区间加的时候这些东西都很好维护。

区间偏移的时候，对$l,r$所属的块进行暴力，而对于中间的每一个块，相当于是删掉了最左边的元素，将其元素的下标都$-1$，然后在最右边插入一个元素，可以通过维护$s$在$O(1)$的时间内更新。

查询的时候，先找出左边/右边第一个出现了$a_x$的块，然后在块内暴力。

时间复杂度$O(n\sqrt{n})$。

#include <cstdio>
#include <iostream>
#include <algorithm>
#include <cstring>
#define ll long long
using namespace std;
template <class T>
inline void rd(T &x) {
	x=0; char c=getchar(); int f=1;
	while(!isdigit(c)) { if(c=='-') f=-1; c=getchar(); }
	while(isdigit(c)) x=x*10-'0'+c,c=getchar(); x*=f;
}
const int N=100010,M=500,C=100000;
int n;
struct BLOCK {
	int a[M+10],f[N*3];
	int tg,s;
	int m;
	void init(int *b,int _m) {
		m=_m;
		for(int i=0;i<m;++i) f[C+(a[i]=b[i])]++;
		tg=0,s=0;
	}
	void add(int l,int r) { // [l, r)
		if(l==0&&r==m) return (void)(tg++);
		for(int i=(s+l)%m;i!=(s+r)%m;i=(i+1)%m)
			f[C+a[i]]--,f[C+(++a[i])]++;
	}
	int flip(int y) {
		int x=a[s]+tg;
		f[C+a[s]]--;
		f[C+(a[s]=y-tg)]++;
		s=(s+1)%m;
		return x;
	}
	int flipl(int l,int r,int x) {
		int c=(s+l)%m,R=(s+r)%m,tmp=a[c]+tg;
		f[C+a[c]]--;
		while(c!=R) a[c]=a[(c+1)%m],c=(c+1)%m;
		f[C+(a[c]=x-tg)]++;
		return tmp;
	}
	void flipi(int l,int r) {
		int tmp=a[(s+l)%m],i,R=(s+r)%m;
		for(i=(s+l)%m;i!=R;i=(i+1)%m) a[i]=a[(i+1)%m];
		a[i]=tmp;
	}
	int Qval(int p) { return a[(s+p)%m]+tg; }
	int Qex(int x) { return f[C+(x-tg)]; }
	int Qex_posr(int x) {
		for(int i=s;;i=(i+1)%m) if(a[i]+tg==x) return (i-s+m)%m;
	}		
	int Qex_posl(int x) {
		for(int i=(s-1+m)%m;;i=(i-1+m)%m) if(a[i]+tg==x) return (i-s+m)%m;
	}
	int Qex_pos(int x,int l,int r,int d) {
		if(l==r) return -1;
		for(int i=(s+l+m)%m;i!=(s+r+m)%m;i=(i+d+m)%m) if(a[i]+tg==x) return (i-s+m)%m;
		return -1;
	}
}A[210];
int a[N],b[N],s[N];
void flip(int l,int r) {
	if(b[l]==b[r]) return A[b[l]].flipi(l-s[b[l]],r-s[b[r]]);
	int tmp=A[b[l]].Qval(l-s[b[l]]);
	tmp=A[b[r]].flipl(0,r-s[b[r]],tmp);
	for(int i=b[r]-1;i>b[l];--i) tmp=A[i].flip(tmp);
	A[b[l]].flipl(l-s[b[l]],M-1,tmp);
}
void add(int l,int r) {
	if(b[l]==b[r]) return A[b[l]].add(l-s[b[l]],r-s[b[r]]+1);
	A[b[l]].add(l-s[b[l]],M);
	A[b[r]].add(0,r-s[b[r]]+1);
	for(int i=b[l]+1;i<b[r];++i) A[i].tg++;
}
int query(int x) {
	int v=A[b[x]].Qval(x-s[b[x]]);
	int R=A[b[x]].Qex_pos(v,x-s[b[x]]+1,0,1);
	int L=A[b[x]].Qex_pos(v,x-s[b[x]]-1,-1,-1);
	if(L==-1) {
		int p=b[x]-1;
		while(p>=0&&!A[p].Qex(v)) p--;
		if(p>=0) L=s[p]+A[p].Qex_posl(v);
		else L=-1e8;
	}
	else L+=s[b[x]];
	if(R==-1) {
		int p=b[x]+1;
		while(p<=b[n-1]&&!A[p].Qex(v)) p++;
		if(p<=b[n-1]) R=s[p]+A[p].Qex_posr(v);
		else R=1e8;
	}
	else R+=s[b[x]];
	int ans=min(x-L,R-x);
	return ans<=n?ans:-1;
}
int main() {
	int q; rd(n),rd(q);
	for(int i=0;i<n;++i) rd(a[i]),b[i]=i/M;
	for(int p=0;p<n;p+=M) A[p/M].init(a+p,min(M,n-p)),s[p/M]=(p/M)*M;
	while(q--) {
		int op,l,r; rd(op);
		if(op==1) rd(l),rd(r),flip(--l,--r);
		else if(op==2) rd(l),rd(r),add(--l,--r);
		else rd(l),printf("%d\n",query(--l));
	}
	return 0;
}
 ans<=n?ans:-1;
}
int main() {
	int q; rd(n),rd(q);
	for(int i=0;i<n;++i) rd(a[i]),b[i]=i/M;
	for(int p=0;p<n;p+=M) A[p/M].init(a+p,min(M,n-p)),s[p/M]=(p/M)*M;
	while(q--) {
		int op,l,r; rd(op);
		if(op==1) rd(l),rd(r),flip(--l,--r);
		else if(op==2) rd(l),rd(r),add(--l,--r);
		else rd(l),printf("%d\n",query(--l));
	}
	return 0;
}