[Codeforces1036E]Covered Points(计算几何求线段交点)

[Codeforces1036E]Covered Points(计算几何求线段交点)

题面

给出(n)条起点和终点都是整点的线段,问这些线段能够覆盖多少个整点.保证线段两两不共线
(n leq 1000)

分析

先不考虑线段相交的情况,只考虑一条线段(AB)上有多少个整点.设(vec{u}=vec{AB}),那么答案就是(gcd(|vec{u}_x|,|vec{u}_y|))
线段的交点会被重复算,因此考虑计算线段的交点被多少条线段覆盖,把重复算的减掉.我们两两计算线段交点,用一个map记录整交点出现的次数,如果点(P)出现了(a)次,即有(a)对线段交点为(P),因此线段个数(x)满足(frac{x(x-1)}{2}=a),解得(x=frac{1+sqrt{1+8a}}{2}). 那么就从答案中减去(x-1)即可

如何判断线段相交和求交点呢,可以用直线两点式爆算。但是有很优美的向量解法,且没有精度误差。

线段相交的判断

快速排斥实验:

以线段(AB)和(CD)为对角线做矩形,如果这两个矩形不相交,那么这两条线段也一定不相交

(点的标号有点问题,但是应该能看出来)

把线段投影到(x)轴上,那么显然有(max(A.x,B.x)geq min(C.x,D.x) 且 max(C.x,D.x) geq min(A.x,B.x))
y轴同理
这个判断被称为快速排斥实验.但是矩形相交线段就一定相交吗?下面给出了一个反例

因此,我们还要用到跨立实验.

---

跨立实验:

如果点(C,D)在直线(AB)的不同侧,则称线段(CD)跨立直线(AB).两线段相交,当且仅当它通过快速排斥实验,且互相跨立。

可以用叉积实现.根据叉积符号的定义知,((vec{AB} imes vec{AD})(vec{AB} imes vec{AC})<0)

bool is_cross(point A,point B,point C,point D){
	if(max(A.x,B.x)>=min(C.x,D.x)&&max(C.x,D.x)>=min(A.x,B.x)&&max(A.y,B.y)>=min(C.y,D.y)&&max(C.y,D.y)>=min(A.y,B.y)){//快速排斥实验
        //跨立实验
		if(sgn(cross(C-A,B-A))*sgn(cross(D-A,B-A))<=0&&sgn(cross(A-C,D-C))*sgn(cross(B-C,D-C))<=0) return 1;//sgn(x)表示x的符号(-1,0,1),不要直接乘,防止溢出 
		else return 0; 
	}else return 0;
}

向量法求线段交点

如果上一步我们已经知道两线段相交,那么就可以利用叉积求交点。

根据叉积的几何意义知道,(S_{Delta ADB}=frac{1}{2}|vec{AB} imes vec{AD}|,S_{Delta ACB}=frac{1}{2}|vec{AB} imes vec{AC}|)

做出三角形的高,显然(frac{CH_2}{H_1H_2}=frac{S_{Delta ACB}}{S_{Delta ADB}})

根据相似三角形得(vec{CP}=frac{CH_2}{H_1H_2+CH_2}vec{CD}=frac{|vec{AB} imes vec{AC}|}{|vec{AB} imes vec{AD}|+|vec{AB} imes vec{AC}|}vec{CD})

(vec{CP})知道了,(P)的坐标也就知道了.

实现上还需要注意一个细节.容易发现除了最后一步的除法之外没有精度误差。那么判断整点只需要在最后一步判余数是否为0即可,不要强行除了之后再取整判断,否则会被卡精度

代码

#include<iostream>
#include<cstdio>
#include<cstring> 
#include<algorithm>
#include<cmath>
#include<map>
#define maxn 1000
#define eps 1e-7
using namespace std;
typedef long long ll;
typedef long long ll; 
ll gcd(ll a,ll b){
	if(b==0) return a;
	else return gcd(b,a%b);
}
struct Vector{
	ll x;
	ll y;
	Vector(){
		
	}
	Vector(ll _x,ll _y){
		x=_x;
		y=_y;
	}
	friend Vector operator + (Vector p,Vector q){
		return Vector(p.x+q.x,p.y+q.y);
	}
	friend Vector operator - (Vector p,Vector q){
		return Vector(p.x-q.x,p.y-q.y);
	}
	friend Vector operator * (Vector p,ll k){
		return Vector(p.x*k,p.y*k);
	} 
	friend Vector operator / (Vector p,ll k){
		return Vector(p.x/k,p.y/k);
	} 
};
typedef Vector point;
inline ll dot(Vector p,Vector q){
	return p.x*q.x+p.y*q.y;
}
inline ll cross(Vector p,Vector q){
	return p.x*q.y-p.y*q.x;
}
inline ll sgn(ll x){
	if(x==0) return 0;
	else if(x>0) return 1;
	else return -1;
}
bool is_cross(point A,point B,point C,point D){
	if(max(A.x,B.x)>=min(C.x,D.x)&&max(C.x,D.x)>=min(A.x,B.x)&&max(A.y,B.y)>=min(C.y,D.y)&&max(C.y,D.y)>=min(A.y,B.y)){
		if(sgn(cross(C-A,B-A))*sgn(cross(D-A,B-A))<=0&&sgn(cross(A-C,D-C))*sgn(cross(B-C,D-C))<=0) return 1;//不要直接乘,防止溢出 
		else return 0; 
	}else return 0;
}

int n;
point A[maxn+5],B[maxn+5];
map< pair<ll,ll>,ll>vis;
void get_cross(point A,point B,point C,point D){
	ll S1=abs(cross(C-A,B-A));
	ll S2=abs(cross(D-A,B-A));
	Vector v=D-C;
	if(v.x*S1%(S1+S2)!=0||v.y*S1%(S1+S2)!=0) return;
	v=v*S1/(S1+S2);
	point P=C+v;
	vis[make_pair(P.x,P.y)]++;
}
ll get_ans(ll a){//x*(x-1)/2=a
	ll x=(1+sqrt(1+8*a))/2;
	return floor(x);
}
int main(){
	ll ans=0;
	scanf("%d",&n);
	for(int i=1;i<=n;i++){
		scanf("%lld %lld %lld %lld",&A[i].x,&A[i].y,&B[i].x,&B[i].y);
		Vector v=B[i]-A[i];
		ans+=gcd(abs(v.x),abs(v.y))+1;
	}
//	printf("debug: %lld
",ans);
	for(int i=1;i<=n;i++){
		for(int j=i+1;j<=n;j++){
			if(is_cross(A[i],B[i],A[j],B[j])){
				get_cross(A[i],B[i],A[j],B[j]);
			}
		}
	}
	for(map< pair<ll,ll>, ll>::iterator it=vis.begin();it!=vis.end();it++){
		ll tim=it->second;
		tim=get_ans(tim);
		ans-=tim-1;
	}
	printf("%lld
",ans);
}