浙大月赛
2014-03-04 00:02
232 查看
A题：暂时没做，有时间补。
B题：素数判定，可以暴力，也可以米勒罗宾法。
米勒罗宾法代码：
/* ***********************************************
Author :xianxingwuguan
Created Time :2014/3/2 12:17:25
File Name :4.cpp
************************************************ */
#pragma comment(linker, "/STACK:102400000,102400000")
#include <stdio.h>
#include <iostream>
#include <algorithm>
#include <sstream>
#include <stdlib.h>
#include <string.h>
#include <limits.h>
#include <string>
#include <time.h>
#include <math.h>
#include <queue>
#include <stack>
#include <set>
#include <map>
using namespace std;
#define INF 0x3f3f3f3f
#define eps 1e-8
#define pi acos(-1.0)
typedef unsigned long long ll;
ll change(ll a,ll b){
ll ans=1;
for(ll i=1;i<a;i++)
ans=ans*b+1;
return ans;
}
const int S=2;
ll mult_mod(ll a,ll b,ll c)
{
a%=c;
b%=c;
ll ret=0;
while(b)
{
if(b&1){ret+=a;ret%=c;}
a<<=1;
if(a>=c)a%=c;
b>>=1;
}
return ret;
}
ll pow_mod(ll x,ll n,ll mod)
{
if(n==1)return x%mod;
x%=mod;
ll tmp=x;
ll ret=1;
while(n)
{
if(n&1)ret=mult_mod(ret,tmp,mod);
tmp=mult_mod(tmp,tmp,mod);
n>>=1;
}
return ret;
}
ll check(ll a,ll n,ll x,ll t)
{
ll ret=pow_mod(a,x,n);
ll last=ret;
for(int i=1;i<=t;i++)
{
ret=mult_mod(ret,ret,n);
if(ret==1 && last!=1 &&last!=n-1)return true;
last=ret;
}
if(ret!=1)return true;
return false;
}
bool judge(ll n)
{
if(n<2)return false;
if(n==2)return true;
if((n&1)==0)return false;
ll x=n-1;
ll t=0;
while((x&1)==0){x>>=1;t++;}
for(int i=0;i<S;i++)
{
ll a=rand()%(n-1)+1;
if(check(a,n,x,t))
return false;
}
return true;
}
int main()
{
//freopen("data.in","r",stdin);
//freopen("data.out","w",stdout);
ll a,b;
while(cin>>a>>b){
//cout<<change(a,b)<<endl;
//cout<<change(a,b)<<endl;
if(judge(change(b,a)))puts("YES");
else puts("NO");
}
return 0;
}
C题:暂时没做，有时间补。
D题：奇偶建图最大流，最大独立点集，
代码：
/* ***********************************************
Author :xianxingwuguan
Created Time :2014/3/3 21:35:47
File Name :F.cpp
************************************************ */
#pragma comment(linker, "/STACK:102400000,102400000")
#include <stdio.h>
#include <iostream>
#include <algorithm>
#include <sstream>
#include <stdlib.h>
#include <string.h>
#include <limits.h>
#include <string>
#include <time.h>
#include <math.h>
#include <queue>
#include <stack>
#include <set>
#include <map>
using namespace std;
#define INF 0x3f3f3f3f
#define eps 1e-8
#define pi acos(-1.0)
typedef long long ll;
const ll maxn=1010;
const ll maxm=400010;
struct Edge{
ll to,next,cap,flow;
Edge(){};
Edge(ll _next,ll _to,ll _cap,ll _flow){
next=_next;to=_to;cap=_cap;flow=_flow;
}
}edge[maxm];
ll head[maxn],tol,gap[maxn],dep[maxn],cur[maxn];
void addedge(ll u,ll v,ll flow){
edge[tol]=Edge(head[u],v,flow,0);head[u]=tol++;
edge[tol]=Edge(head[v],u,0,0);head[v]=tol++;
}
ll Q[maxn];
void bfs(ll start,ll end){
memset(dep,-1,sizeof(dep));
memset(gap,0,sizeof(gap));
gap[0]++;ll front=0,rear=0;
dep[end]=0;Q[rear++]=end;
while(front!=rear){
ll u=Q[front++];
for(ll i=head[u];i!=-1;i=edge[i].next){
ll v=edge[i].to;if(dep[v]==-1&&edge[i].cap)
Q[rear++]=v,dep[v]=dep[u]+1,gap[dep[v]]++;
}
}
}
ll S[maxn];
ll sap(ll start,ll end,ll N){
bfs(start,end);
memcpy(cur,head,sizeof(head));
ll top=0,u=start,ans=0;
while(dep[start]<N){
if(u==end){
ll MIN=INF,id;
for(ll i=0;i<top;i++)
if(MIN>edge[S[i]].cap-edge[S[i]].flow)
MIN=edge[S[i]].cap-edge[S[i]].flow,id=i;
for(ll i=0;i<top;i++)
edge[S[i]].flow+=MIN,edge[S[i]^1].flow-=MIN;
ans+=MIN,top=id,u=edge[S[top]^1].to;
continue;
}
bool flag=0;ll v;
for(ll i=cur[u];i!=-1;i=edge[i].next){
v=edge[i].to;
if(edge[i].cap-edge[i].flow&&dep[v]+1==dep[u]){
flag=1;cur[u]=i;break;
}
}
if(flag){
S[top++]=cur[u];u=v;continue;
}
ll MIN=N;
for(ll i=head[u];i!=-1;i=edge[i].next)
if(edge[i].cap-edge[i].flow&&dep[edge[i].to]<MIN)
MIN=dep[edge[i].to],cur[u]=i;
if(--gap[dep[u]]==0)break;gap[dep[u]=MIN+1]++;
if(u!=start)u=edge[S[--top]^1].to;
}
return ans;
}
ll gcd(ll a,ll b){
if(a==0)return b;
return gcd(b%a,a);
}
ll a[1000],b[1000];
int main()
{
//freopen("data.in","r",stdin);
//freopen("data.out","w",stdout);
ll n,p;
while(~scanf("%lld%lld",&n,&p)){
memset(head,-1,sizeof(head));tol=0;
ll sum=0;
for(ll i=1;i<=n;i++){
ll j,k;
scanf("%lld%lld",&j,&k);
a[i]=j^k;
b[i]=j&k;
sum+=b[i];
}
for(ll i=1;i<=n;i++)
if(a[i]&1){
addedge(0,i,b[i]);
for(ll j=1;j<=n;j++)
if(gcd(a[i]^a[j],p)==1)addedge(i,j,INF);
}
else addedge(i,n+1,b[i]);
printf("%lld\n",sum-sap(0,n+1,n+3));
}
return 0;
}
E题：暂时没做，有时间补。
F题:两种办法，一种是枚举每一个点，找最长边，然后枚举，n^2效率。
另一种:要以每个点为中心先极角排序，然后枚举边
这个时候假设是有向地枚举边
然后另一个点始终在左侧范围内
也就是那个第三个点
这样的话如果当前点到底边的高度比前一个点大，就一直移动,
动不了的时候，再将底边往下移动一条
第一种效率较高，
方法一代码：
/* ***********************************************
Author :xianxingwuguan
Created Time :2014/2/27 10:15:39
File Name :\代码\acm程序\计算几何.cpp
************************************************ */
#pragma comment(linker, "/STACK:102400000,102400000")
#include <stdio.h>
#include <iostream>
#include <algorithm>
#include <sstream>
#include <stdlib.h>
#include <string.h>
#include <limits.h>
#include <string>
#include <time.h>
#include <math.h>
#include <queue>
#include <stack>
#include <set>
#include <map>
using namespace std;
#define INF 0x3f3f3f3f
#define eps 1e-8
#define pi acos(-1.0)
typedef long long ll;
int dcmp(double x){
if(fabs(x)<eps)return 0;
return x>0?1:-1;
}
struct Point{
double x,y;
Point(double _x=0,double _y=0){
x=_x;y=_y;
}
};
Point operator + (Point a,Point b){
return Point(a.x+b.x,a.y+b.y);
}
Point operator - (Point a,Point b){
return Point(a.x-b.x,a.y-b.y);
}
Point operator * (Point a,double p){
return Point(a.x*p,a.y*p);
}
Point operator / (Point a,double p){
return Point(a.x/p,a.y/p);
}
bool operator < (const Point &a,const Point &b){//重载大小排序，
return a.x<b.x||(a.x==b.x&&a.y<b.y);
}
bool operator == (const Point &a,const Point &b){//重载相等
return dcmp(a.x-b.x)==0&&dcmp(a.y-b.y)==0;
}
double Dot(Point a,Point b){//点积
return a.x*b.x+a.y*b.y;
}
double Cross(Point a,Point b){//叉积
return a.x*b.y-a.y*b.x;
}
double Length(Point a){//向量长度
return sqrt(Dot(a,a));
}
double Angle(Point a,Point b){//a->b顺时针为负，逆时针为正。
return acos(Dot(a,b)/Length(a)/Length(b));
}
Point Rotate(Point a,double rad){ //向量逆时针旋转，rad为弧度。
return Point(a.x*cos(rad)-a.y*sin(rad),a.x*sin(rad)+a.y*cos(rad));
}
double Area2(Point a,Point b,Point c){//三角形面积二倍。
return fabs(Cross(b-a,c-a));
}
Point Normal(Point a){//计算向量的单位法线，即左转90度，长度归一化。
double l=Length(a);//用前确保A不是零向量
return Point(-a.y/l,a.x/l);
}
Point GetLineIntersection(Point P,Point v,Point Q,Point w){
Point u=P-Q;//求解直线P+t*v Q+t*w 的交点
double t=Cross(w,u)/Cross(v,w);//用前确保两条直线有唯一交点
return P+v*t;//当且仅当cross(v,w)!=0
}
double DistanceToLine(Point p,Point a,Point b){
Point v1=b-a;//计算点P到直线AB距离
Point v2=p-a;//如果不取绝对值，得的是有向距离
return fabs(Cross(v1,v2))/(Length(v1));
}
double DistanceToSegment(Point p,Point a,Point b){
if(a==b)return Length(p-a);//点到线段AB的距离。
Point v1=b-a,v2=p-a,v3=p-b;
if(dcmp(Dot(v1,v2))<0)return Length(v2);
else if(dcmp(Dot(v1,v3))>0)return Length(v3);
else return fabs(Cross(v1,v2))/Length(v1);
}
Point GetLineProjection(Point p,Point a,Point b){
//求点在直线上投影
//AB直线:A+tv , 投影Q的参数为 t0
//由 Dot(v,P-(A+t0*v))==0 得出
Point v=b-a;
return a+v*(Dot(v,p-a)/Dot(v,v));
}
void getLineGenralEquation(Point p1,Point p2,double &a,double &b,double &c){
a=p2.y-p1.y;//直线两点式转一般式。
b=p1.x-p2.x;
c=-a*p1.x-b*p1.y;
}
bool SegmentProperIntersection(Point a1,Point a2,Point b1,Point b2){
//线段相交判定
double c1=Cross(a2-a1,b1-a1);
double c2=Cross(a2-a1,b2-a1);
double c3=Cross(b2-b1,a1-b1);
double c4=Cross(b2-b1,a2-b1);
return dcmp(c1)*dcmp(c2)<0&&dcmp(c3)*dcmp(c4)<0;
}
bool OnSegment(Point p,Point a1,Point a2){
//判断一个点在线段上
return dcmp(Cross(a1-p,a2-p))==0&&dcmp(Dot(a1-p,a2-p))<0;
}
double ConvexPolygonArea(Point* p,int n){
//从第一个顶点出发把凸多边形分成 n-2 个三角形
double area=0;
for(int i=1;i<n-1;i++)area+=Cross(p[i]-p[0],p[i+1]-p[0]);
return area/2;
}
struct Circle{
Point c;
double r;
Circle(){}
Circle(Point cc,double rr){
c=cc; r=rr;
}
Point point(double a){
//通过圆心角来求坐标
return Point(c.x+cos(a)*r,c.y+sin(a)*r);
}
};
struct Line{
//有向直线,左边就是对应的半平面
Point p,v;
double ang;
Line(){}
Line(Point pp,Point vv){
p=pp;v=vv;
ang=atan2(v.y,v.x);
}
Point point(double a){
return p+(v*a);
}
bool operator<(const Line& L) const{
return ang<L.ang;
}
};
bool OnLeft(Line l,Point pp){
//点Ｐ在有向直线的左边，线上不算
return Cross(l.v,pp-l.p)>0;
}
Point GetLineIntersection(Line a,Line b){
//二直线相交，假设交点唯一存在
Point u=a.p-b.p;
double t=Cross(b.v,u)/Cross(a.v,b.v);
return a.p+a.v*t;
}
int getLineCircleIntersection(Line l,Circle cc,double& t1,double& t2,vector<Point>& sol){
//直线和圆的交点，直线为AB,圆的圆心为C，半径为R，第一种方法是解方程组
//设交点为P＝A+t*(B-A),代入圆方程整理后先得到 (a*t+b)^2+(c*t+d)^2=r^2
//进一步整理后得到一元二次方程 e*t^2+f*t+g=0
//根据判别式的值可以区分3种情况，
//
//函数返回的是交点个数，参数sol存放的是交点本身
//没有清空，可以反复调用函数，把所有的交点放在一个sol里
//---
//另一种方法是几何法
//先求C在AB上的投影P,AB的单位向量V
//交点P－LV，P+LV，由勾股定理算出
double a=l.v.x;
double b=l.p.x-cc.c.x;
double c=l.v.y;
double d=l.p.y-cc.c.y;
double e=a*a+c*c;
double f=2*(a*b+c*d);
double g=b*b+d*d-cc.r*cc.r;
double delta=f*f-4*e*g;//判别式
if(dcmp(delta)<0)return 0;//相离
if(dcmp(delta)==0){
t1=t2=-f/(2*e);
sol.push_back(l.point(t1));
return 1;//相切
}
t1=(-f-sqrt(delta))/(2*e);sol.push_back(l.point(t1));
t2=(-f+sqrt(delta))/(2*e);sol.push_back(l.point(t2));
return 2;//相交
}
double angle(Point v){
//计算向量的极角
return atan2(v.y,v.x);
}
int getCircleCircleIntersection(Circle c1,Circle c2,vector<Point>& sol){ //计算两圆交点
double d=Length(c1.c-c2.c);
if(dcmp(d)==0){
if(dcmp(c1.r-c2.r)==0)return -1;//两圆重合
return 0;
}
if(dcmp(c1.r+c2.r-d)<0)return 0;
if(dcmp(fabs(c1.r-c2.r)-d)>0)return 0;
double a=angle(c2.c-c1.c);//向量C1C2的极角
double da=acos((c1.r*c1.r+d*d-c2.r*c2.r)/(2*c1.r*d));//计算C1C2到C1P1的角
Point p1=c1.point(a-da);Point p2=c1.point(a+da);
sol.push_back(p1);if(p1==p2)return 1;
sol.push_back(p2);
return 2;
}
int getTangents(Point p,Circle c,Point* v){
//过点p到圆C的切线，v[i]是第i条切线的向量，返回切线条数
Point u=c.c-p;
double dist=Length(u);
if(dist<c.r)return 0;
else if(dcmp(dist-c.r)==0){
v[0]=Rotate(u,pi/2);
return 1;
}
else{
double ang=asin(c.r/dist);
v[0]=Rotate(u,-ang);v[1]=Rotate(u,ang);
return 2;
}
}
int getTangents(Circle a,Circle b,Point* aa,Point* bb){
//返回切线条数，-1表示无穷条切线,a[i]和b[i]分别表示第i条切线在圆A和圆B上的切点
int cnt=0;
if(a.r<b.r)swap(a,b), swap(aa,bb);
int d2=(int)((a.c.x-b.c.x)*(a.c.x-b.c.x)+(a.c.y-b.c.y)*(a.c.y-b.c.y));
int rdiff=(int)(a.r-b.r),rsum=(int)(a.r+b.r);
if(d2<rdiff*rdiff)return 0;//内含
double base=atan2(b.c.y-a.c.y,b.c.x-a.c.x);
if(d2==0&&a.r==b.r) return -1;//无穷多条切线
if(d2==rdiff*rdiff){
aa[cnt]=a.point(base);
bb[cnt]=b.point(base);
cnt++;
return 1;//内切，1条
}
//有外公切线
double ang=acos((a.r-b.r)/sqrt(d2));
aa[cnt]=a.point(base+ang);
bb[cnt]=b.point(base+ang);
cnt++;
aa[cnt]=a.point(base-ang);
bb[cnt]=b.point(base-ang);
cnt++;
if(d2==rsum*rsum){
aa[cnt]=a.point(base);
bb[cnt]=b.point(base+pi);
cnt++;//一条内公切线
}
else
if(d2>rsum*rsum){
double ang=acos((a.r+b.r)/sqrt(d2));
aa[cnt]=a.point(base+ang);
bb[cnt]=b.point(pi+base+ang);
cnt++;
aa[cnt]=a.point(base-ang);
bb[cnt]=b.point(pi+base-ang);
cnt++;//两条内公切线
}
return cnt;
}
Circle CircumscribedCircle(Point p1,Point p2,Point p3){
//三角形外接圆
double bx=p2.x-p1.x,by=p2.y-p1.y;
double cx=p3.x-p1.x,cy=p3.y-p1.y;
double d=2*(bx*cy-cx*by);
double ccx=(cy*(bx*bx+by*by)-by*(cx*cx+cy*cy))/d+p1.x;
double ccy=(bx*(cx*cx+cy*cy)-cx*(bx*bx+by*by))/d+p1.y;
Point p=Point(ccx,ccy);
return Circle(p,Length(p1-p));
}
Circle InscribedCircle(Point p1,Point p2,Point p3){
//三角形内切圆
double a=Length(p2-p3),b=Length(p3-p1);
double c=Length(p1-p2);
Point p=(p1*a+p1*b+p3*c)/(a+b+c);
return Circle(p,DistanceToLine(p,p1,p2));
}
double torad(double deg){
//角度转弧度
return deg/180*pi;
}
void get_coord(double r,double lat,double lng,double& x,double& y,double& z){
//经纬度(角度)转化为空间坐标
lat=torad(lat);lng=torad(lng);
x=r*cos(lat)*cos(lng);y=r*cos(lat)*sin(lng);
z=r*sin(lat);
}
//求球面上两点最短距离，走大圆弧
//半径r，弦长d，圆心角 2arcsin(d/2r)
struct Polygon{
int n;
Point p[1000];
Polygon(){ n=0;}
Polygon(int nn,Point* pp){
n=nn;
for(int q=0;q<n;q++) p[q]=pp[q];
}
Point operator[](int a) const{
return p[a];
}
};
int isPointInPolygon(Point p,Polygon po){
//假想有一条向右的射线，统计多边形穿过这条射线正反多少次，把这个数记为绕数
//逆时穿过加1，顺时针穿过减1
int wn=0, n=po.n;
for(int i=0;i<n;i++){
if(OnSegment(p,po[i],po[(i+1)%n]))return -1;//在边界上
int k=dcmp(Cross(po[(i+1)%n]-po[i],p-po[i]));
int d1=dcmp(po[i].y-p.y);
int d2=dcmp(po[(i+1)%n].y-p.y);
if(k>0&&d1<=0&&d2>0)wn++;
if(k<0&&d2<=0&&d1>0)wn--;
}
if(wn!=0) return 1;//内部
return 0;//外部
}
//点在凸多边形内判定更简单，只要判断是否在所有边的左边
//各顶点按逆时针顺序
int ConvexHull(Point* p,int n,Point* ch)
{
//计算凸包，输入点数组p，个数为n，输出点数组ch，函数返回凸包顶点数
//输入不能有重复点，函数执行完后输入点的顺序被破坏
//如果不希望在凸包的边上有输入点，把两个<=改成<
//在精度要求高时建议用dcmp比较
sort(p,p+n);
//n=unique(p,p+n)-p;
int m=0;
for(int i=0;i<n;i++){
while(m>1&&Cross(ch[m-1]-ch[m-2],p[i]-ch[m-2])<=0)m--;
ch[m++]=p[i];
}
int k=m;
for(int i=n-2;i>=0;i--){
while(m>k&&Cross(ch[m-1]-ch[m-2],p[i]-ch[m-2])<=0)m--;
ch[m++]=p[i];
}
if(n>1)m--;
return m;
}
double diameter(Point* o,int n,Point *p){
//求点集直径(最大距离),旋转卡壳
int num=ConvexHull(o,n,p);
int mi=0,ma=0;
for(int q=0;q<num;q++){
if(p[q]<p[mi])mi=q;
if(p[ma]<p[q])ma=q;
}
int v=ma;
double ans=(p[mi].x-p[ma].x)*(p[mi].x-p[ma].x)+(p[mi].y-p[ma].y)*(p[mi].y-p[ma].y);
ans=max(ans,(p[num-1].x-p[ma].x)*(p[num-1].x-p[ma].x)+(p[num-1].y-p[ma].y)*(p[num-1].y-p[ma].y));
ans=max(ans,(p[mi].x-p[ma-1].x)*(p[mi].x-p[ma-1].x)+(p[mi].y-p[ma-1].y)*(p[mi].y-p[ma-1].y));
ans=max(ans,(p[num-1].x-p[ma-1].x)*(p[num-1].x-p[ma-1].x)+(p[num-1].y-p[ma-1].y)*(p[num-1].y-p[ma-1].y));
for(;mi<v;mi++){
for(;;){
int a=dcmp(Cross(p[mi+1]-p[mi],p[ma+1]-p[ma]));
if(a==0){
ans=max(ans,(p[mi].x-p[ma].x)*(p[mi].x-p[ma].x)+(p[mi].y-p[ma].y)*(p[mi].y-p[ma].y));
ans=max(ans,(p[mi+1].x-p[ma].x)*(p[mi+1].x-p[ma].x)+(p[mi+1].y-p[ma].y)*(p[mi+1].y-p[ma].y));
ans=max(ans,(p[mi].x-p[ma+1].x)*(p[mi].x-p[ma+1].x)+(p[mi].y-p[ma+1].y)*(p[mi].y-p[ma+1].y));
ans=max(ans,(p[mi+1].x-p[ma+1].x)*(p[mi+1].x-p[ma+1].x)+(p[mi+1].y-p[ma+1].y)*(p[mi+1].y-p[ma+1].y));
ma++;
break;
}
if(a>0){
ans=max(ans,(p[mi].x-p[ma].x)*(p[mi].x-p[ma].x)+(p[mi].y-p[ma].y)*(p[mi].y-p[ma].y));
ans=max(ans,(p[mi].x-p[ma+1].x)*(p[mi].x-p[ma+1].x)+(p[mi].y-p[ma+1].y)*(p[mi].y-p[ma+1].y));
ma++;
}
if(a<0){
ans=max(ans,(p[mi+1].x-p[ma].x)*(p[mi+1].x-p[ma].x)+(p[mi+1].y-p[ma].y)*(p[mi+1].y-p[ma].y));
ans=max(ans,(p[mi+1].x-p[ma+1].x)*(p[mi+1].x-p[ma+1].x)+(p[mi+1].y-p[ma+1].y)*(p[mi+1].y-p[ma+1].y));
break;
}
}
}
return ans;
}
int HalfplaneIntersection(Line* l,int n,Point* poly){
sort(l,l+n);
int first,last;
Point *p=new Point
;Line *q=new Line
;
q[first=last=0]=l[0];
for(int i=1;i<n;i++){
while(first<last&&!OnLeft(l[i],p[last-1]))last--;
while(first<last&&!OnLeft(l[i],p[first]))first++;
q[++last]=l[i];
if(fabs(Cross(q[last].v,q[last-1].v))<eps){
last--;
if(OnLeft(q[last],l[i].p))q[last]=l[i];
}
if(first<last)p[last-1]=GetLineIntersection(q[last-1],q[last]);
}
while(first<last&&!OnLeft(q[first],p[last-1]))last--;
if(last-first<=1)return 0;
p[last]=GetLineIntersection(q[last],q[first]);
int m=0;
for(int i=first;i<=last;i++)poly[m++]=p[i];
return m;
}
void FindBorder(Point* p,int n){
//寻找多边形边界
int num;
Point ans[10005],z[205];
double x[205];
int a,min1=0;
for(int q=0;q<n;q++)if(p[q]<p[min1])min1=q;
p
=p[0];ans[0]=p[min1];
double s=-pi;num=1;
for(int idx=0;;idx++){
a=0;
for(int q=0;q<n;q++)
if(p[q]==ans[num-1]){
z[a]=p[(q+1)%n];
x[a]=atan2(-p[q].y+z[a].y,-p[q].x+z[a].x)-s-2*pi;
for(;dcmp(x[a])<=0;)x[a]+=2*pi;
a++;z[a]=p[(q+n-1)%n];
x[a]=atan2(-p[q].y+z[a].y,-p[q].x+z[a].x)-s-2*pi;
for(;dcmp(x[a])<=0;)x[a]+=2*pi;
a++;
}
for(int q=0;q<n;q++)
if(OnSegment(ans[num-1],p[q],p[q+1])){
z[a]=p[(q+1)];
x[a]=atan2(-ans[num-1].y+z[a].y,-ans[num-1].x+z[a].x)-s-2*pi;
for(;dcmp(x[a])<=0;)x[a]+=2*pi;
a++;
z[a]=p[q];
x[a]=atan2(-ans[num-1].y+z[a].y,-ans[num-1].x+z[a].x)-s-2*pi;
for(;dcmp(x[a])<=0;)x[a]+=2*pi;
a++;
}
int mi=0;
for(int q=0;q<a;q++)if(x[q]<x[mi])mi=q;
double ma=101000000;
Point f=z[mi],g;
for(int q=0;q<n;q++)
if(SegmentProperIntersection(z[mi],ans[num-1],p[q],p[q+1])){
g=GetLineIntersection(z[mi],z[mi]-ans[num-1],p[q],p[q]-p[q+1]);
if(Length(g-ans[num-1])<ma)ma=Length(g-ans[num-1]),f=g;
}
ans[num]=f;
s=atan2(ans[num-1].y-ans[num].y,ans[num-1].x-ans[num].x);num++;
if(ans[num-1]==ans[0])break;
}
printf("%d\n",num-1);
for(int q=0;q<num-1;q++)
printf("%.4lf %.4lf\n",ans[q].x,ans[q].y);
}
Polygon cutPolygon(Polygon& po,Point a,Point s){
//切割多边形，返回直线左边的部分,要正反调用两次
Point x,c;Polygon ans;
int sum=0;
for(int q=0;q<po.n;q++){
x=po[q];
c=po[(q+1)%po.n];
if(OnLeft(Line(a,s-a),x)||OnSegment(x,a,s))ans.p[sum]=x,sum++;
if(SegmentProperIntersection(x,c,a,s))ans.p[sum]=GetLineIntersection(x,c-x,a,s-a),sum++;
}
ans.n=sum;
return ans;
}
struct Point3{
double x,y,z;
Point3(){
x=0;y=0;z=0;
}
Point3(double xx,double yy,double zz){
x=xx;y=xx;z=zz;
}
bool operator==(Point3 p) const{
return dcmp(x-p.x)==0&&dcmp(y-p.y)==0&&dcmp(z-p.z)==0;
}
};
Point3 operator+(Point3 a,Point3 b){
return Point3(a.x+b.x,a.y+b.y,a.z+b.z);
}
Point3 operator-(Point3 a,Point3 b){
return Point3(a.x-b.x,a.y-b.y,a.z-b.z);
}
Point3 operator*(Point3 a,double p){
return Point3(a.x*p,a.y*p,a.z*p);
}
Point3 operator/(Point3 a,double p){
return Point3(a.x/p,a.y/p,a.z/p);
}
double Dot(Point3 a,Point3 b){
return a.x*b.x+a.y*b.y+a.z*b.y;
}
double Length(Point3 a){
return sqrt(Dot(a,a));
}
double Angle(Point3 a,Point3 b){
return acos(Dot(a,b)/Length(a)/Length(b));
}
//平面用点法式表示(p0,n),p0是平面上一个点，n 是平面法向量
double DistanceToPlane(Point3 &p,Point3 &p0,Point3 &n){
//点p到平面 p0-n 的的距离，n 是单位向量
return fabs(Dot(p-p0,n));
//如果不取绝对值，得到的是有向距离
}
Point3 GetPlaneProjection(Point3 &p,Point3 &p0,Point3 &n){
//点 p 在平面 p0-n 上的投影，n 是单位向量
return p-n*Dot(p-p0,n);
}
Point3 LinePlaneIntersection(Point3 p1,Point3 p2,Point3 p0,Point3 n){
//直线 p1-p2 到平面 p0-n 的交点,分母为 0 对应直线和平面平行，直线在平面上
Point3 v=p2-p1;
double t=(Dot(n,p0-p1)/Dot(n,p2-p1));//判断分母是否为 0
return p1+v*t;//如果是线段，判断 t 是不是在 0 和 1 之间
}
Point3 Cross(Point3 a,Point3 b){
return Point3(a.y*b.z-a.z*b.y,a.z*b.x-a.x*b.z,a.x*b.y-a.y*b.x);
}
double Area2(Point3 a,Point3 b,Point3 c){
return Length(Cross(b-a,c-a));
}
bool PointInTri(Point3 p,Point3 p0,Point3 p1,Point3 p2){
//判断点在三角形内，先判断是否在三角形所在的平面，再用面积关系
double area1=Area2(p,p0,p1);
double area2=Area2(p,p2,p1);
double area3=Area2(p,p0,p2);
return dcmp(area1+area2+area3-Area2(p0,p1,p2))==0;
}
bool TriSegIntersection(Point3 p0,Point3 p1,Point3 p2,Point3 a,Point3 b,Point3 &p){
//判断线段和三角形是否相交,没有考虑线段在平面上的情况
Point3 n=Cross(p1-p0,p2-p0);
if(dcmp(Dot(n,b-a))==0)return false;
else{
double t=Dot(n,p0-a)/Dot(n,b-a);
if(dcmp(t)<0||dcmp(t-1)>0)return false;
p=a+(b-a)*t;//计算交点
return PointInTri(p,p0,p1,p2);//判断交点在三角形内
}
}
double DistanceToLine(Point3 p,Point3 a,Point3 b){
//点 p 到直线 ab 距离
Point3 v1=b-a,v2=p-a;
return Length(Cross(v1,v2))/Length(v1);
}
double DistanceToSegment(Point3 p,Point3 a,Point3 b){
//点 p 到线段 ab 距离
if(a==b)return Length(p-a);
Point3 v1=b-a,v2=p-a,v3=p-b;
if(dcmp(Dot(v1,v2))<0) return Length(v2);
else if(dcmp(Dot(v1,v3))>0) return Length(v3);
else return Length(Cross(v1,v2))/Length(v1);
}
double Volume6(Point3 a,Point3 b,Point3 c,Point3 d){
return Dot(d-a,Cross(b-a,c-a));
}
double TriCir(Point a,Point s,Point d,Circle c,vector<Point> v){
//三角形与圆面积交
//其中圆心是 a ,三角形 abc
double x1=Length(s-a),x2=Length(d-a);
double vv,bb;
if(dcmp(x1-c.r)<=0&&dcmp(x2-c.r)<=0)return Area2(a,s,d)/2;
if(dcmp(x1-c.r)>0&&dcmp(x2-c.r)>0){
if(DistanceToLine(c.c,s,d)>=c.r)return Angle(d-a,s-a)*c.r*c.r/2;
else{
v.clear();
getLineCircleIntersection(Line(s,d-s),c,vv,bb,v);vv=Angle(v[0]-a,v[1]-a);
return min(Angle(v[0]-a,d-a)+Angle(v[1]-a,s-a),Angle(v[1]-a,d-a)+Angle(v[0]-a,s-a))/2*c.r*c.r+sin(Angle(v[0]-a,v[1]-a))/2*c.r*c.r;
}
}
if(dcmp(x1-c.r)>0&&dcmp(x2-c.r)<=0){
v.clear();
getLineCircleIntersection(Line(s,d-s),c,vv,bb,v);
if(OnSegment(v[1],s,d)){
v[0]=v[1];printf("$$");
}
return Area2(a,d,v[0])/2+Angle(v[0]-a,s-a)*c.r*c.r/2;
}
if(dcmp(x1-c.r)<=0&&dcmp(x2-c.r)>0){
v.clear();
getLineCircleIntersection(Line(s,d-s),c,vv,bb,v);
if(OnSegment(v[1],s,d))v[0]=v[1];
return Area2(a,s,v[0])/2+Angle(v[0]-a,d-a)*c.r*c.r/2;
}
}
Point p[600];
double Max(double a,double b){
if(a<=b)return b;
return a;
}
int main(){
//freopen("data.in","r",stdin);
//freopen("data.out","w",stdout);
int n;
while(~scanf("%d",&n)){
for(int i=0;i<n;i++)scanf("%lf%lf",&p[i].x,&p[i].y);
double ans=0.0;
for(int i=0;i<n;i++){
int id=0;
double han=Length(p[i]-p[0]);
for(int j=1;j<n;j++){
if(i==j)continue;
double hh=Length(p[i]-p[j]);
if(hh>han){
han=hh;
id=j;
}
}
for(int j=0;j<n;j++){
if(j==i||j==id)continue;
ans=max(ans,DistanceToLine(p[i],p[id],p[j]));
ans=max(ans,DistanceToLine(p[id],p[i],p[j]));
ans=max(ans,DistanceToLine(p[j],p[i],p[id]));
}
}
printf("%.10lf\n",ans);
}
return 0;
}
方法二代码：
/* ***********************************************
Author :xianxingwuguan
Created Time :2014/2/27 10:15:39
File Name :\代码\acm程序\计算几何.cpp
************************************************ */
#pragma comment(linker, "/STACK:102400000,102400000")
#include <stdio.h>
#include <iostream>
#include <algorithm>
#include <sstream>
#include <stdlib.h>
#include <string.h>
#include <limits.h>
#include <string>
#include <time.h>
#include <math.h>
#include <queue>
#include <stack>
#include <set>
#include <map>
using namespace std;
#define INF 0x3f3f3f3f
#define eps 1e-8
#define pi acos(-1.0)
typedef long long ll;
int dcmp(double x){
if(fabs(x)<eps)return 0;
return x>0?1:-1;
}
struct Point{
double x,y;
Point(double _x=0,double _y=0){
x=_x;y=_y;
}
};
Point operator + (Point a,Point b){
return Point(a.x+b.x,a.y+b.y);
}
Point operator - (Point a,Point b){
return Point(a.x-b.x,a.y-b.y);
}
Point operator * (Point a,double p){
return Point(a.x*p,a.y*p);
}
Point operator / (Point a,double p){
return Point(a.x/p,a.y/p);
}
bool operator < (const Point &a,const Point &b){//重载大小排序，
return a.x<b.x||(a.x==b.x&&a.y<b.y);
}
bool operator == (const Point &a,const Point &b){//重载相等
return dcmp(a.x-b.x)==0&&dcmp(a.y-b.y)==0;
}
double Dot(Point a,Point b){//点积
return a.x*b.x+a.y*b.y;
}
double Cross(Point a,Point b){//叉积
return a.x*b.y-a.y*b.x;
}
double Length(Point a){//向量长度
return sqrt(Dot(a,a));
}
double Angle(Point a,Point b){//a->b顺时针为负，逆时针为正。
return acos(Dot(a,b)/Length(a)/Length(b));
}
Point Rotate(Point a,double rad){ //向量逆时针旋转，rad为弧度。
return Point(a.x*cos(rad)-a.y*sin(rad),a.x*sin(rad)+a.y*cos(rad));
}
double Area2(Point a,Point b,Point c){//三角形面积二倍。
return fabs(Cross(b-a,c-a));
}
Point Normal(Point a){//计算向量的单位法线，即左转90度，长度归一化。
double l=Length(a);//用前确保A不是零向量
return Point(-a.y/l,a.x/l);
}
Point GetLineIntersection(Point P,Point v,Point Q,Point w){
Point u=P-Q;//求解直线P+t*v Q+t*w 的交点
double t=Cross(w,u)/Cross(v,w);//用前确保两条直线有唯一交点
return P+v*t;//当且仅当cross(v,w)!=0
}
double DistanceToLine(Point p,Point a,Point b){
Point v1=b-a;//计算点P到直线AB距离
Point v2=p-a;//如果不取绝对值，得的是有向距离
return fabs(Cross(v1,v2))/(Length(v1));
}
double DistanceToSegment(Point p,Point a,Point b){
if(a==b)return Length(p-a);//点到线段AB的距离。
Point v1=b-a,v2=p-a,v3=p-b;
if(dcmp(Dot(v1,v2))<0)return Length(v2);
else if(dcmp(Dot(v1,v3))>0)return Length(v3);
else return fabs(Cross(v1,v2))/Length(v1);
}
Point GetLineProjection(Point p,Point a,Point b){
//求点在直线上投影
//AB直线:A+tv , 投影Q的参数为 t0
//由 Dot(v,P-(A+t0*v))==0 得出
Point v=b-a;
return a+v*(Dot(v,p-a)/Dot(v,v));
}
void getLineGenralEquation(Point p1,Point p2,double &a,double &b,double &c){
a=p2.y-p1.y;//直线两点式转一般式。
b=p1.x-p2.x;
c=-a*p1.x-b*p1.y;
}
bool SegmentProperIntersection(Point a1,Point a2,Point b1,Point b2){
//线段相交判定
double c1=Cross(a2-a1,b1-a1);
double c2=Cross(a2-a1,b2-a1);
double c3=Cross(b2-b1,a1-b1);
double c4=Cross(b2-b1,a2-b1);
return dcmp(c1)*dcmp(c2)<0&&dcmp(c3)*dcmp(c4)<0;
}
bool OnSegment(Point p,Point a1,Point a2){
//判断一个点在线段上
return dcmp(Cross(a1-p,a2-p))==0&&dcmp(Dot(a1-p,a2-p))<0;
}
double ConvexPolygonArea(Point* p,int n){
//从第一个顶点出发把凸多边形分成 n-2 个三角形
double area=0;
for(int i=1;i<n-1;i++)area+=Cross(p[i]-p[0],p[i+1]-p[0]);
return area/2;
}
struct Circle{
Point c;
double r;
Circle(){}
Circle(Point cc,double rr){
c=cc; r=rr;
}
Point point(double a){
//通过圆心角来求坐标
return Point(c.x+cos(a)*r,c.y+sin(a)*r);
}
};
struct Line{
//有向直线,左边就是对应的半平面
Point p,v;
double ang;
Line(){}
Line(Point pp,Point vv){
p=pp;v=vv;
ang=atan2(v.y,v.x);
}
Point point(double a){
return p+(v*a);
}
bool operator<(const Line& L) const{
return ang<L.ang;
}
};
bool OnLeft(Line l,Point pp){
//点Ｐ在有向直线的左边，线上不算
return Cross(l.v,pp-l.p)>0;
}
Point GetLineIntersection(Line a,Line b){
//二直线相交，假设交点唯一存在
Point u=a.p-b.p;
double t=Cross(b.v,u)/Cross(a.v,b.v);
return a.p+a.v*t;
}
int getLineCircleIntersection(Line l,Circle cc,double& t1,double& t2,vector<Point>& sol){
//直线和圆的交点，直线为AB,圆的圆心为C，半径为R，第一种方法是解方程组
//设交点为P＝A+t*(B-A),代入圆方程整理后先得到 (a*t+b)^2+(c*t+d)^2=r^2
//进一步整理后得到一元二次方程 e*t^2+f*t+g=0
//根据判别式的值可以区分3种情况，
//
//函数返回的是交点个数，参数sol存放的是交点本身
//没有清空，可以反复调用函数，把所有的交点放在一个sol里
//---
//另一种方法是几何法
//先求C在AB上的投影P,AB的单位向量V
//交点P－LV，P+LV，由勾股定理算出
double a=l.v.x;
double b=l.p.x-cc.c.x;
double c=l.v.y;
double d=l.p.y-cc.c.y;
double e=a*a+c*c;
double f=2*(a*b+c*d);
double g=b*b+d*d-cc.r*cc.r;
double delta=f*f-4*e*g;//判别式
if(dcmp(delta)<0)return 0;//相离
if(dcmp(delta)==0){
t1=t2=-f/(2*e);
sol.push_back(l.point(t1));
return 1;//相切
}
t1=(-f-sqrt(delta))/(2*e);sol.push_back(l.point(t1));
t2=(-f+sqrt(delta))/(2*e);sol.push_back(l.point(t2));
return 2;//相交
}
double angle(Point v){
//计算向量的极角
return atan2(v.y,v.x);
}
int getCircleCircleIntersection(Circle c1,Circle c2,vector<Point>& sol){ //计算两圆交点
double d=Length(c1.c-c2.c);
if(dcmp(d)==0){
if(dcmp(c1.r-c2.r)==0)return -1;//两圆重合
return 0;
}
if(dcmp(c1.r+c2.r-d)<0)return 0;
if(dcmp(fabs(c1.r-c2.r)-d)>0)return 0;
double a=angle(c2.c-c1.c);//向量C1C2的极角
double da=acos((c1.r*c1.r+d*d-c2.r*c2.r)/(2*c1.r*d));//计算C1C2到C1P1的角
Point p1=c1.point(a-da);Point p2=c1.point(a+da);
sol.push_back(p1);if(p1==p2)return 1;
sol.push_back(p2);
return 2;
}
int getTangents(Point p,Circle c,Point* v){
//过点p到圆C的切线，v[i]是第i条切线的向量，返回切线条数
Point u=c.c-p;
double dist=Length(u);
if(dist<c.r)return 0;
else if(dcmp(dist-c.r)==0){
v[0]=Rotate(u,pi/2);
return 1;
}
else{
double ang=asin(c.r/dist);
v[0]=Rotate(u,-ang);v[1]=Rotate(u,ang);
return 2;
}
}
int getTangents(Circle a,Circle b,Point* aa,Point* bb){
//返回切线条数，-1表示无穷条切线,a[i]和b[i]分别表示第i条切线在圆A和圆B上的切点
int cnt=0;
if(a.r<b.r)swap(a,b), swap(aa,bb);
int d2=(int)((a.c.x-b.c.x)*(a.c.x-b.c.x)+(a.c.y-b.c.y)*(a.c.y-b.c.y));
int rdiff=(int)(a.r-b.r),rsum=(int)(a.r+b.r);
if(d2<rdiff*rdiff)return 0;//内含
double base=atan2(b.c.y-a.c.y,b.c.x-a.c.x);
if(d2==0&&a.r==b.r) return -1;//无穷多条切线
if(d2==rdiff*rdiff){
aa[cnt]=a.point(base);
bb[cnt]=b.point(base);
cnt++;
return 1;//内切，1条
}
//有外公切线
double ang=acos((a.r-b.r)/sqrt(d2));
aa[cnt]=a.point(base+ang);
bb[cnt]=b.point(base+ang);
cnt++;
aa[cnt]=a.point(base-ang);
bb[cnt]=b.point(base-ang);
cnt++;
if(d2==rsum*rsum){
aa[cnt]=a.point(base);
bb[cnt]=b.point(base+pi);
cnt++;//一条内公切线
}
else
if(d2>rsum*rsum){
double ang=acos((a.r+b.r)/sqrt(d2));
aa[cnt]=a.point(base+ang);
bb[cnt]=b.point(pi+base+ang);
cnt++;
aa[cnt]=a.point(base-ang);
bb[cnt]=b.point(pi+base-ang);
cnt++;//两条内公切线
}
return cnt;
}
Circle CircumscribedCircle(Point p1,Point p2,Point p3){
//三角形外接圆
double bx=p2.x-p1.x,by=p2.y-p1.y;
double cx=p3.x-p1.x,cy=p3.y-p1.y;
double d=2*(bx*cy-cx*by);
double ccx=(cy*(bx*bx+by*by)-by*(cx*cx+cy*cy))/d+p1.x;
double ccy=(bx*(cx*cx+cy*cy)-cx*(bx*bx+by*by))/d+p1.y;
Point p=Point(ccx,ccy);
return Circle(p,Length(p1-p));
}
Circle InscribedCircle(Point p1,Point p2,Point p3){
//三角形内切圆
double a=Length(p2-p3),b=Length(p3-p1);
double c=Length(p1-p2);
Point p=(p1*a+p1*b+p3*c)/(a+b+c);
return Circle(p,DistanceToLine(p,p1,p2));
}
double torad(double deg){
//角度转弧度
return deg/180*pi;
}
void get_coord(double r,double lat,double lng,double& x,double& y,double& z){
//经纬度(角度)转化为空间坐标
lat=torad(lat);lng=torad(lng);
x=r*cos(lat)*cos(lng);y=r*cos(lat)*sin(lng);
z=r*sin(lat);
}
//求球面上两点最短距离，走大圆弧
//半径r，弦长d，圆心角 2arcsin(d/2r)
struct Polygon{
int n;
Point p[1000];
Polygon(){ n=0;}
Polygon(int nn,Point* pp){
n=nn;
for(int q=0;q<n;q++) p[q]=pp[q];
}
Point operator[](int a) const{
return p[a];
}
};
int isPointInPolygon(Point p,Polygon po){
//假想有一条向右的射线，统计多边形穿过这条射线正反多少次，把这个数记为绕数
//逆时穿过加1，顺时针穿过减1
int wn=0, n=po.n;
for(int i=0;i<n;i++){
if(OnSegment(p,po[i],po[(i+1)%n]))return -1;//在边界上
int k=dcmp(Cross(po[(i+1)%n]-po[i],p-po[i]));
int d1=dcmp(po[i].y-p.y);
int d2=dcmp(po[(i+1)%n].y-p.y);
if(k>0&&d1<=0&&d2>0)wn++;
if(k<0&&d2<=0&&d1>0)wn--;
}
if(wn!=0) return 1;//内部
return 0;//外部
}
//点在凸多边形内判定更简单，只要判断是否在所有边的左边
//各顶点按逆时针顺序
int ConvexHull(Point* p,int n,Point* ch)
{
//计算凸包，输入点数组p，个数为n，输出点数组ch，函数返回凸包顶点数
//输入不能有重复点，函数执行完后输入点的顺序被破坏
//如果不希望在凸包的边上有输入点，把两个<=改成<
//在精度要求高时建议用dcmp比较
sort(p,p+n);
//n=unique(p,p+n)-p;
int m=0;
for(int i=0;i<n;i++){
while(m>1&&Cross(ch[m-1]-ch[m-2],p[i]-ch[m-2])<=0)m--;
ch[m++]=p[i];
}
int k=m;
for(int i=n-2;i>=0;i--){
while(m>k&&Cross(ch[m-1]-ch[m-2],p[i]-ch[m-2])<=0)m--;
ch[m++]=p[i];
}
if(n>1)m--;
return m;
}
double diameter(Point* o,int n,Point *p){
//求点集直径(最大距离),旋转卡壳
int num=ConvexHull(o,n,p);
int mi=0,ma=0;
for(int q=0;q<num;q++){
if(p[q]<p[mi])mi=q;
if(p[ma]<p[q])ma=q;
}
int v=ma;
double ans=(p[mi].x-p[ma].x)*(p[mi].x-p[ma].x)+(p[mi].y-p[ma].y)*(p[mi].y-p[ma].y);
ans=max(ans,(p[num-1].x-p[ma].x)*(p[num-1].x-p[ma].x)+(p[num-1].y-p[ma].y)*(p[num-1].y-p[ma].y));
ans=max(ans,(p[mi].x-p[ma-1].x)*(p[mi].x-p[ma-1].x)+(p[mi].y-p[ma-1].y)*(p[mi].y-p[ma-1].y));
ans=max(ans,(p[num-1].x-p[ma-1].x)*(p[num-1].x-p[ma-1].x)+(p[num-1].y-p[ma-1].y)*(p[num-1].y-p[ma-1].y));
for(;mi<v;mi++){
for(;;){
int a=dcmp(Cross(p[mi+1]-p[mi],p[ma+1]-p[ma]));
if(a==0){
ans=max(ans,(p[mi].x-p[ma].x)*(p[mi].x-p[ma].x)+(p[mi].y-p[ma].y)*(p[mi].y-p[ma].y));
ans=max(ans,(p[mi+1].x-p[ma].x)*(p[mi+1].x-p[ma].x)+(p[mi+1].y-p[ma].y)*(p[mi+1].y-p[ma].y));
ans=max(ans,(p[mi].x-p[ma+1].x)*(p[mi].x-p[ma+1].x)+(p[mi].y-p[ma+1].y)*(p[mi].y-p[ma+1].y));
ans=max(ans,(p[mi+1].x-p[ma+1].x)*(p[mi+1].x-p[ma+1].x)+(p[mi+1].y-p[ma+1].y)*(p[mi+1].y-p[ma+1].y));
ma++;
break;
}
if(a>0){
ans=max(ans,(p[mi].x-p[ma].x)*(p[mi].x-p[ma].x)+(p[mi].y-p[ma].y)*(p[mi].y-p[ma].y));
ans=max(ans,(p[mi].x-p[ma+1].x)*(p[mi].x-p[ma+1].x)+(p[mi].y-p[ma+1].y)*(p[mi].y-p[ma+1].y));
ma++;
}
if(a<0){
ans=max(ans,(p[mi+1].x-p[ma].x)*(p[mi+1].x-p[ma].x)+(p[mi+1].y-p[ma].y)*(p[mi+1].y-p[ma].y));
ans=max(ans,(p[mi+1].x-p[ma+1].x)*(p[mi+1].x-p[ma+1].x)+(p[mi+1].y-p[ma+1].y)*(p[mi+1].y-p[ma+1].y));
break;
}
}
}
return ans;
}
int HalfplaneIntersection(Line* l,int n,Point* poly){
sort(l,l+n);
int first,last;
Point *p=new Point
;Line *q=new Line
;
q[first=last=0]=l[0];
for(int i=1;i<n;i++){
while(first<last&&!OnLeft(l[i],p[last-1]))last--;
while(first<last&&!OnLeft(l[i],p[first]))first++;
q[++last]=l[i];
if(fabs(Cross(q[last].v,q[last-1].v))<eps){
last--;
if(OnLeft(q[last],l[i].p))q[last]=l[i];
}
if(first<last)p[last-1]=GetLineIntersection(q[last-1],q[last]);
}
while(first<last&&!OnLeft(q[first],p[last-1]))last--;
if(last-first<=1)return 0;
p[last]=GetLineIntersection(q[last],q[first]);
int m=0;
for(int i=first;i<=last;i++)poly[m++]=p[i];
return m;
}
void FindBorder(Point* p,int n){
//寻找多边形边界
int num;
Point ans[10005],z[205];
double x[205];
int a,min1=0;
for(int q=0;q<n;q++)if(p[q]<p[min1])min1=q;
p
=p[0];ans[0]=p[min1];
double s=-pi;num=1;
for(int idx=0;;idx++){
a=0;
for(int q=0;q<n;q++)
if(p[q]==ans[num-1]){
z[a]=p[(q+1)%n];
x[a]=atan2(-p[q].y+z[a].y,-p[q].x+z[a].x)-s-2*pi;
for(;dcmp(x[a])<=0;)x[a]+=2*pi;
a++;z[a]=p[(q+n-1)%n];
x[a]=atan2(-p[q].y+z[a].y,-p[q].x+z[a].x)-s-2*pi;
for(;dcmp(x[a])<=0;)x[a]+=2*pi;
a++;
}
for(int q=0;q<n;q++)
if(OnSegment(ans[num-1],p[q],p[q+1])){
z[a]=p[(q+1)];
x[a]=atan2(-ans[num-1].y+z[a].y,-ans[num-1].x+z[a].x)-s-2*pi;
for(;dcmp(x[a])<=0;)x[a]+=2*pi;
a++;
z[a]=p[q];
x[a]=atan2(-ans[num-1].y+z[a].y,-ans[num-1].x+z[a].x)-s-2*pi;
for(;dcmp(x[a])<=0;)x[a]+=2*pi;
a++;
}
int mi=0;
for(int q=0;q<a;q++)if(x[q]<x[mi])mi=q;
double ma=101000000;
Point f=z[mi],g;
for(int q=0;q<n;q++)
if(SegmentProperIntersection(z[mi],ans[num-1],p[q],p[q+1])){
g=GetLineIntersection(z[mi],z[mi]-ans[num-1],p[q],p[q]-p[q+1]);
if(Length(g-ans[num-1])<ma)ma=Length(g-ans[num-1]),f=g;
}
ans[num]=f;
s=atan2(ans[num-1].y-ans[num].y,ans[num-1].x-ans[num].x);num++;
if(ans[num-1]==ans[0])break;
}
printf("%d\n",num-1);
for(int q=0;q<num-1;q++)
printf("%.4lf %.4lf\n",ans[q].x,ans[q].y);
}
Polygon cutPolygon(Polygon& po,Point a,Point s){
//切割多边形，返回直线左边的部分,要正反调用两次
Point x,c;Polygon ans;
int sum=0;
for(int q=0;q<po.n;q++){
x=po[q];
c=po[(q+1)%po.n];
if(OnLeft(Line(a,s-a),x)||OnSegment(x,a,s))ans.p[sum]=x,sum++;
if(SegmentProperIntersection(x,c,a,s))ans.p[sum]=GetLineIntersection(x,c-x,a,s-a),sum++;
}
ans.n=sum;
return ans;
}
struct Point3{
double x,y,z;
Point3(){
x=0;y=0;z=0;
}
Point3(double xx,double yy,double zz){
x=xx;y=xx;z=zz;
}
bool operator==(Point3 p) const{
return dcmp(x-p.x)==0&&dcmp(y-p.y)==0&&dcmp(z-p.z)==0;
}
};
Point3 operator+(Point3 a,Point3 b){
return Point3(a.x+b.x,a.y+b.y,a.z+b.z);
}
Point3 operator-(Point3 a,Point3 b){
return Point3(a.x-b.x,a.y-b.y,a.z-b.z);
}
Point3 operator*(Point3 a,double p){
return Point3(a.x*p,a.y*p,a.z*p);
}
Point3 operator/(Point3 a,double p){
return Point3(a.x/p,a.y/p,a.z/p);
}
double Dot(Point3 a,Point3 b){
return a.x*b.x+a.y*b.y+a.z*b.y;
}
double Length(Point3 a){
return sqrt(Dot(a,a));
}
double Angle(Point3 a,Point3 b){
return acos(Dot(a,b)/Length(a)/Length(b));
}
//平面用点法式表示(p0,n),p0是平面上一个点，n 是平面法向量
double DistanceToPlane(Point3 &p,Point3 &p0,Point3 &n){
//点p到平面 p0-n 的的距离，n 是单位向量
return fabs(Dot(p-p0,n));
//如果不取绝对值，得到的是有向距离
}
Point3 GetPlaneProjection(Point3 &p,Point3 &p0,Point3 &n){
//点 p 在平面 p0-n 上的投影，n 是单位向量
return p-n*Dot(p-p0,n);
}
Point3 LinePlaneIntersection(Point3 p1,Point3 p2,Point3 p0,Point3 n){
//直线 p1-p2 到平面 p0-n 的交点,分母为 0 对应直线和平面平行，直线在平面上
Point3 v=p2-p1;
double t=(Dot(n,p0-p1)/Dot(n,p2-p1));//判断分母是否为 0
return p1+v*t;//如果是线段，判断 t 是不是在 0 和 1 之间
}
Point3 Cross(Point3 a,Point3 b){
return Point3(a.y*b.z-a.z*b.y,a.z*b.x-a.x*b.z,a.x*b.y-a.y*b.x);
}
double Area2(Point3 a,Point3 b,Point3 c){
return Length(Cross(b-a,c-a));
}
bool PointInTri(Point3 p,Point3 p0,Point3 p1,Point3 p2){
//判断点在三角形内，先判断是否在三角形所在的平面，再用面积关系
double area1=Area2(p,p0,p1);
double area2=Area2(p,p2,p1);
double area3=Area2(p,p0,p2);
return dcmp(area1+area2+area3-Area2(p0,p1,p2))==0;
}
bool TriSegIntersection(Point3 p0,Point3 p1,Point3 p2,Point3 a,Point3 b,Point3 &p){
//判断线段和三角形是否相交,没有考虑线段在平面上的情况
Point3 n=Cross(p1-p0,p2-p0);
if(dcmp(Dot(n,b-a))==0)return false;
else{
double t=Dot(n,p0-a)/Dot(n,b-a);
if(dcmp(t)<0||dcmp(t-1)>0)return false;
p=a+(b-a)*t;//计算交点
return PointInTri(p,p0,p1,p2);//判断交点在三角形内
}
}
double DistanceToLine(Point3 p,Point3 a,Point3 b){
//点 p 到直线 ab 距离
Point3 v1=b-a,v2=p-a;
return Length(Cross(v1,v2))/Length(v1);
}
double DistanceToSegment(Point3 p,Point3 a,Point3 b){
//点 p 到线段 ab 距离
if(a==b)return Length(p-a);
Point3 v1=b-a,v2=p-a,v3=p-b;
if(dcmp(Dot(v1,v2))<0) return Length(v2);
else if(dcmp(Dot(v1,v3))>0) return Length(v3);
else return Length(Cross(v1,v2))/Length(v1);
}
double Volume6(Point3 a,Point3 b,Point3 c,Point3 d){
return Dot(d-a,Cross(b-a,c-a));
}
double TriCir(Point a,Point s,Point d,Circle c,vector<Point> v){
//三角形与圆面积交
//其中圆心是 a ,三角形 abc
double x1=Length(s-a),x2=Length(d-a);
double vv,bb;
if(dcmp(x1-c.r)<=0&&dcmp(x2-c.r)<=0)return Area2(a,s,d)/2;
if(dcmp(x1-c.r)>0&&dcmp(x2-c.r)>0){
if(DistanceToLine(c.c,s,d)>=c.r)return Angle(d-a,s-a)*c.r*c.r/2;
else{
v.clear();
getLineCircleIntersection(Line(s,d-s),c,vv,bb,v);vv=Angle(v[0]-a,v[1]-a);
return min(Angle(v[0]-a,d-a)+Angle(v[1]-a,s-a),Angle(v[1]-a,d-a)+Angle(v[0]-a,s-a))/2*c.r*c.r+sin(Angle(v[0]-a,v[1]-a))/2*c.r*c.r;
}
}
if(dcmp(x1-c.r)>0&&dcmp(x2-c.r)<=0){
v.clear();
getLineCircleIntersection(Line(s,d-s),c,vv,bb,v);
if(OnSegment(v[1],s,d)){
v[0]=v[1];printf("$$");
}
return Area2(a,d,v[0])/2+Angle(v[0]-a,s-a)*c.r*c.r/2;
}
if(dcmp(x1-c.r)<=0&&dcmp(x2-c.r)>0){
v.clear();
getLineCircleIntersection(Line(s,d-s),c,vv,bb,v);
if(OnSegment(v[1],s,d))v[0]=v[1];
return Area2(a,s,v[0])/2+Angle(v[0]-a,d-a)*c.r*c.r/2;
}
}
Point p[600],ch[600];
Line rot[600*500];
double fun(int k,int n,int m){
double ret=0;
int cnt=0;
for(int i=0;i<n;i++){
if(i==k)continue;
rot[cnt]=Line(p[k],p[i]-p[k]);
++cnt;
}
sort(rot,rot+cnt);
int j=0;
while(1){
if(dcmp(Cross(rot[0].v,ch[j]-p[k]))>=0&&dcmp(Cross(rot[0].v,ch[(j+m-1)%m]-p[k]))<=0)break;
j=(j+1)%m;
}
for(int i=0;i<cnt;i++){
Point tmp=rot[i].v+rot[i].p;
while(DistanceToLine(ch[j+1],p[k],tmp)>DistanceToLine(ch[j],p[k],tmp)||dcmp(DistanceToLine(ch[j],p[k],tmp))==0)j=(j+1)%m;
ret=max(ret,DistanceToLine(ch[j],p[k],tmp));
}
return ret;
}
int main(){
//freopen("data.in","r",stdin);
//freopen("data.out","w",stdout);
int n;
while(~scanf("%d",&n)){
for(int i=0;i<n;i++)scanf("%lf%lf",&p[i].x,&p[i].y);
int m=ConvexHull(p,n,ch);
double ans=0;
for(int i=0;i<n;i++)ans=max(ans,fun(i,n,m));
printf("%.10lf\n",ans);
}
return 0;
}
G题:暂时没做，有时间补。
H题：状压枚举
代码：#include<stdio.h>
#include<string.h>
#include<time.h>
#include<math.h>
#include<algorithm>
#include<iostream>
#include<vector>
#include<set>
#include<map>
using namespace std;
typedef long long lld;
#define X first
#define Y second
#define mp make_pair
#define pb push_back
struct E
{
int x,y;
lld s;
void in()
{
cin >> x >> y >> s;
x--;
y--;
}
}e[100010];
lld pp[100010];
int qq;
int main()
{
int n,m,k;
while(scanf("%d %d %d",&n,&m,&k)!=EOF)
{
for(int i=0;i<m;i++)
e[i].in();
lld ans=(1LL)<<60;
for(int mask=0;mask<(1<<n);mask++)
{
if(!(mask&(1<<0)))
continue;
if(mask&(1<<(n-1)))
continue;
qq=0;
for(int i=0;i<m;i++)
{
int x=e[i].x;
int y=e[i].y;
bool flag=false;
if((mask&(1<<x)) && !(mask&(1<<y)))
flag=true;
swap(x,y);
if((mask&(1<<x)) && !(mask&(1<<y)))
flag=true;
if(flag)
pp[qq++]=e[i].s;
}
sort(pp,pp+qq);
lld now=0;
for(int i=0;i<qq-k;i++)
now+=pp[i];
ans=min(ans,now);
}
cout << ans << endl;
}
return 0;
}
I题:splay裸题，自己没实现出来，贴一个baka巨巨的，慢慢学习。
代码：
#include <iostream>
#include<cstdio>
#include<cstring>
#include<string>
#include<algorithm>
#include<map>
#include<queue>
#include<stack>
#include<cmath>
#include<vector>
#define inf 0x3f3f3f3f
#define Inf 0x3FFFFFFFFFFFFFFFLL
#define eps 1e-9
#define pi acos(-1.0)
using namespace std;
typedef long long ll;
const int maxn=300000+10;
int ch[maxn][2],pre[maxn],size[maxn],val[maxn],st[maxn];
int gcdst[maxn][2],root,tot;
int num[200010],fs[200010];
int gcd(int a,int b)
{
if(a==-1) return b;
if(b==-1) return a;
int c;
while(b)
{
c=a%b;
a=b;
b=c;
}
return a;
}
void NewNode(int &rt,int fa,int v,int state)
{
rt=++tot;
ch[rt][0]=ch[rt][1]=0;
size[rt]=1;val[rt]=v;
pre[rt]=fa;
gcdst[rt][state]=v;gcdst[rt][state^1]=-1;
st[rt]=state;
}
void PushUp(int rt)
{
size[rt]=size[ch[rt][0]]+size[ch[rt][1]]+1;
gcdst[rt][0]=gcd(gcdst[ch[rt][0]][0],gcdst[ch[rt][1]][0]);
gcdst[rt][1]=gcd(gcdst[ch[rt][0]][1],gcdst[ch[rt][1]][1]);
gcdst[rt][st[rt]]=gcd(gcdst[rt][st[rt]],val[rt]);
}
void Rotate(int x,int kind)
{
int y=pre[x];
ch[y][kind^1]=ch[x][kind];
pre[ch[x][kind]]=y;
if(pre[y])
ch[pre[y]][ch[pre[y]][1]==y]=x;
pre[x]=pre[y];
ch[x][kind]=y;
pre[y]=x;
PushUp(y);
}
void Splay(int rt,int goal)
{
while(pre[rt]!=goal)
{
int y=pre[rt];
if(pre[y]==goal)
Rotate(rt,ch[y][0]==rt);
else
{
int kind=ch[pre[y]][0]==y;
if(ch[y][kind]==rt)
{
Rotate(rt,kind^1);
Rotate(rt,kind);
}
else
{
Rotate(y,kind);
Rotate(rt,kind);
}
}
}
PushUp(rt);
if(goal==0) root=rt;
}
void RotateTo(int k,int goal)
{
int rt=root;
while(size[ch[rt][0]]!=k)
{
if(size[ch[rt][0]]>k)
rt=ch[rt][0];
else
{
k-=(size[ch[rt][0]]+1);
rt=ch[rt][1];
}
}
Splay(rt,goal);
}
void build(int l,int r,int &rt,int fa)
{
if(l>r) return ;
int m=(l+r)>>1;
NewNode(rt,fa,num[m],fs[m]);
build(l,m-1,ch[rt][0],rt);
build(m+1,r,ch[rt][1],rt);
PushUp(rt);
}
int Query(int L,int R,int state)
{
RotateTo(L-1,0);
RotateTo(R+1,root);
return gcdst[ch[ch[root][1]][0]][state];
}
void Insert(int pos,int v,int state)
{
RotateTo(pos,0);
if(ch[root][0]==0)
{
NewNode(ch[root][0],root,v,state);
PushUp(root);
return ;
}
int rc=ch[root][0];
while(ch[rc][1])
rc=ch[rc][1];
Splay(rc,root);
NewNode(ch[rc][1],rc,v,state);
PushUp(rc);
PushUp(root);
}
void Delete(int pos)
{
RotateTo(pos,0);
if(ch[root][0]==0)
{
root=ch[root][1];
pre[root]=0;
return ;
}
int rc=ch[root][0];
while(ch[rc][1]) rc=ch[rc][1];
Splay(rc,root);
int rt=ch[root][1];
ch[rc][1]=rt;
pre[rt]=rc;
root=rc;
pre[rc]=0;
PushUp(root);
}
void Changes(int pos)
{
RotateTo(pos,0);
st[root]^=1;
PushUp(root);
}
void Modify(int pos,int v)
{
RotateTo(pos,0);
val[root]=v;
PushUp(root);
}
void Init(int n)
{
for(int i=0;i<n;++i)
scanf("%d%d",&num[i],&fs[i]);
ch[0][0]=ch[0][1]=0;
pre[0]=size[0]=0;val[0]=-1;
gcdst[0][0]=gcdst[0][1]=-1;
root=tot=0;
NewNode(root,0,-1,0);
NewNode(ch[root][1],root,-1,0);
build(0,n-1,ch[ch[root][1]][0],ch[root][1]);
PushUp(ch[root][1]);
PushUp(root);
}
int main()
{
//freopen("in.txt","r",stdin);
//freopen("out.txt","w",stdout);
int n,q;
char str[5];
while(~scanf("%d%d",&n,&q))
{
Init(n);
while(q--)
{
scanf("%s",str);
int L,R,x,s;
if(str[0]=='Q')
{
scanf("%d%d%d",&L,&R,&s);
printf("%d\n",Query(L,R,s));
}
else if(str[0]=='I')
{
scanf("%d%d%d",&L,&x,&s);
Insert(L+1,x,s);
}
else if(str[0]=='D')
{
scanf("%d",&L);
Delete(L);
}
else if(str[0]=='R')
{
scanf("%d",&L);
Changes(L);
}
else
{
scanf("%d%d",&L,&x);
Modify(L,x);
}
}
}
return 0;
}
J题：没做，有时间补。
B题：素数判定，可以暴力，也可以米勒罗宾法。
米勒罗宾法代码：
/* ***********************************************
Author :xianxingwuguan
Created Time :2014/3/2 12:17:25
File Name :4.cpp
************************************************ */
#pragma comment(linker, "/STACK:102400000,102400000")
#include <stdio.h>
#include <iostream>
#include <algorithm>
#include <sstream>
#include <stdlib.h>
#include <string.h>
#include <limits.h>
#include <string>
#include <time.h>
#include <math.h>
#include <queue>
#include <stack>
#include <set>
#include <map>
using namespace std;
#define INF 0x3f3f3f3f
#define eps 1e-8
#define pi acos(-1.0)
typedef unsigned long long ll;
ll change(ll a,ll b){
ll ans=1;
for(ll i=1;i<a;i++)
ans=ans*b+1;
return ans;
}
const int S=2;
ll mult_mod(ll a,ll b,ll c)
{
a%=c;
b%=c;
ll ret=0;
while(b)
{
if(b&1){ret+=a;ret%=c;}
a<<=1;
if(a>=c)a%=c;
b>>=1;
}
return ret;
}
ll pow_mod(ll x,ll n,ll mod)
{
if(n==1)return x%mod;
x%=mod;
ll tmp=x;
ll ret=1;
while(n)
{
if(n&1)ret=mult_mod(ret,tmp,mod);
tmp=mult_mod(tmp,tmp,mod);
n>>=1;
}
return ret;
}
ll check(ll a,ll n,ll x,ll t)
{
ll ret=pow_mod(a,x,n);
ll last=ret;
for(int i=1;i<=t;i++)
{
ret=mult_mod(ret,ret,n);
if(ret==1 && last!=1 &&last!=n-1)return true;
last=ret;
}
if(ret!=1)return true;
return false;
}
bool judge(ll n)
{
if(n<2)return false;
if(n==2)return true;
if((n&1)==0)return false;
ll x=n-1;
ll t=0;
while((x&1)==0){x>>=1;t++;}
for(int i=0;i<S;i++)
{
ll a=rand()%(n-1)+1;
if(check(a,n,x,t))
return false;
}
return true;
}
int main()
{
//freopen("data.in","r",stdin);
//freopen("data.out","w",stdout);
ll a,b;
while(cin>>a>>b){
//cout<<change(a,b)<<endl;
//cout<<change(a,b)<<endl;
if(judge(change(b,a)))puts("YES");
else puts("NO");
}
return 0;
}
C题:暂时没做，有时间补。
D题：奇偶建图最大流，最大独立点集，
代码：
/* ***********************************************
Author :xianxingwuguan
Created Time :2014/3/3 21:35:47
File Name :F.cpp
************************************************ */
#pragma comment(linker, "/STACK:102400000,102400000")
#include <stdio.h>
#include <iostream>
#include <algorithm>
#include <sstream>
#include <stdlib.h>
#include <string.h>
#include <limits.h>
#include <string>
#include <time.h>
#include <math.h>
#include <queue>
#include <stack>
#include <set>
#include <map>
using namespace std;
#define INF 0x3f3f3f3f
#define eps 1e-8
#define pi acos(-1.0)
typedef long long ll;
const ll maxn=1010;
const ll maxm=400010;
struct Edge{
ll to,next,cap,flow;
Edge(){};
Edge(ll _next,ll _to,ll _cap,ll _flow){
next=_next;to=_to;cap=_cap;flow=_flow;
}
}edge[maxm];
ll head[maxn],tol,gap[maxn],dep[maxn],cur[maxn];
void addedge(ll u,ll v,ll flow){
edge[tol]=Edge(head[u],v,flow,0);head[u]=tol++;
edge[tol]=Edge(head[v],u,0,0);head[v]=tol++;
}
ll Q[maxn];
void bfs(ll start,ll end){
memset(dep,-1,sizeof(dep));
memset(gap,0,sizeof(gap));
gap[0]++;ll front=0,rear=0;
dep[end]=0;Q[rear++]=end;
while(front!=rear){
ll u=Q[front++];
for(ll i=head[u];i!=-1;i=edge[i].next){
ll v=edge[i].to;if(dep[v]==-1&&edge[i].cap)
Q[rear++]=v,dep[v]=dep[u]+1,gap[dep[v]]++;
}
}
}
ll S[maxn];
ll sap(ll start,ll end,ll N){
bfs(start,end);
memcpy(cur,head,sizeof(head));
ll top=0,u=start,ans=0;
while(dep[start]<N){
if(u==end){
ll MIN=INF,id;
for(ll i=0;i<top;i++)
if(MIN>edge[S[i]].cap-edge[S[i]].flow)
MIN=edge[S[i]].cap-edge[S[i]].flow,id=i;
for(ll i=0;i<top;i++)
edge[S[i]].flow+=MIN,edge[S[i]^1].flow-=MIN;
ans+=MIN,top=id,u=edge[S[top]^1].to;
continue;
}
bool flag=0;ll v;
for(ll i=cur[u];i!=-1;i=edge[i].next){
v=edge[i].to;
if(edge[i].cap-edge[i].flow&&dep[v]+1==dep[u]){
flag=1;cur[u]=i;break;
}
}
if(flag){
S[top++]=cur[u];u=v;continue;
}
ll MIN=N;
for(ll i=head[u];i!=-1;i=edge[i].next)
if(edge[i].cap-edge[i].flow&&dep[edge[i].to]<MIN)
MIN=dep[edge[i].to],cur[u]=i;
if(--gap[dep[u]]==0)break;gap[dep[u]=MIN+1]++;
if(u!=start)u=edge[S[--top]^1].to;
}
return ans;
}
ll gcd(ll a,ll b){
if(a==0)return b;
return gcd(b%a,a);
}
ll a[1000],b[1000];
int main()
{
//freopen("data.in","r",stdin);
//freopen("data.out","w",stdout);
ll n,p;
while(~scanf("%lld%lld",&n,&p)){
memset(head,-1,sizeof(head));tol=0;
ll sum=0;
for(ll i=1;i<=n;i++){
ll j,k;
scanf("%lld%lld",&j,&k);
a[i]=j^k;
b[i]=j&k;
sum+=b[i];
}
for(ll i=1;i<=n;i++)
if(a[i]&1){
addedge(0,i,b[i]);
for(ll j=1;j<=n;j++)
if(gcd(a[i]^a[j],p)==1)addedge(i,j,INF);
}
else addedge(i,n+1,b[i]);
printf("%lld\n",sum-sap(0,n+1,n+3));
}
return 0;
}
E题：暂时没做，有时间补。
F题:两种办法，一种是枚举每一个点，找最长边，然后枚举，n^2效率。
另一种:要以每个点为中心先极角排序，然后枚举边
这个时候假设是有向地枚举边
然后另一个点始终在左侧范围内
也就是那个第三个点
这样的话如果当前点到底边的高度比前一个点大，就一直移动,
动不了的时候，再将底边往下移动一条
第一种效率较高，
方法一代码：
/* ***********************************************
Author :xianxingwuguan
Created Time :2014/2/27 10:15:39
File Name :\代码\acm程序\计算几何.cpp
************************************************ */
#pragma comment(linker, "/STACK:102400000,102400000")
#include <stdio.h>
#include <iostream>
#include <algorithm>
#include <sstream>
#include <stdlib.h>
#include <string.h>
#include <limits.h>
#include <string>
#include <time.h>
#include <math.h>
#include <queue>
#include <stack>
#include <set>
#include <map>
using namespace std;
#define INF 0x3f3f3f3f
#define eps 1e-8
#define pi acos(-1.0)
typedef long long ll;
int dcmp(double x){
if(fabs(x)<eps)return 0;
return x>0?1:-1;
}
struct Point{
double x,y;
Point(double _x=0,double _y=0){
x=_x;y=_y;
}
};
Point operator + (Point a,Point b){
return Point(a.x+b.x,a.y+b.y);
}
Point operator - (Point a,Point b){
return Point(a.x-b.x,a.y-b.y);
}
Point operator * (Point a,double p){
return Point(a.x*p,a.y*p);
}
Point operator / (Point a,double p){
return Point(a.x/p,a.y/p);
}
bool operator < (const Point &a,const Point &b){//重载大小排序，
return a.x<b.x||(a.x==b.x&&a.y<b.y);
}
bool operator == (const Point &a,const Point &b){//重载相等
return dcmp(a.x-b.x)==0&&dcmp(a.y-b.y)==0;
}
double Dot(Point a,Point b){//点积
return a.x*b.x+a.y*b.y;
}
double Cross(Point a,Point b){//叉积
return a.x*b.y-a.y*b.x;
}
double Length(Point a){//向量长度
return sqrt(Dot(a,a));
}
double Angle(Point a,Point b){//a->b顺时针为负，逆时针为正。
return acos(Dot(a,b)/Length(a)/Length(b));
}
Point Rotate(Point a,double rad){ //向量逆时针旋转，rad为弧度。
return Point(a.x*cos(rad)-a.y*sin(rad),a.x*sin(rad)+a.y*cos(rad));
}
double Area2(Point a,Point b,Point c){//三角形面积二倍。
return fabs(Cross(b-a,c-a));
}
Point Normal(Point a){//计算向量的单位法线，即左转90度，长度归一化。
double l=Length(a);//用前确保A不是零向量
return Point(-a.y/l,a.x/l);
}
Point GetLineIntersection(Point P,Point v,Point Q,Point w){
Point u=P-Q;//求解直线P+t*v Q+t*w 的交点
double t=Cross(w,u)/Cross(v,w);//用前确保两条直线有唯一交点
return P+v*t;//当且仅当cross(v,w)!=0
}
double DistanceToLine(Point p,Point a,Point b){
Point v1=b-a;//计算点P到直线AB距离
Point v2=p-a;//如果不取绝对值，得的是有向距离
return fabs(Cross(v1,v2))/(Length(v1));
}
double DistanceToSegment(Point p,Point a,Point b){
if(a==b)return Length(p-a);//点到线段AB的距离。
Point v1=b-a,v2=p-a,v3=p-b;
if(dcmp(Dot(v1,v2))<0)return Length(v2);
else if(dcmp(Dot(v1,v3))>0)return Length(v3);
else return fabs(Cross(v1,v2))/Length(v1);
}
Point GetLineProjection(Point p,Point a,Point b){
//求点在直线上投影
//AB直线:A+tv , 投影Q的参数为 t0
//由 Dot(v,P-(A+t0*v))==0 得出
Point v=b-a;
return a+v*(Dot(v,p-a)/Dot(v,v));
}
void getLineGenralEquation(Point p1,Point p2,double &a,double &b,double &c){
a=p2.y-p1.y;//直线两点式转一般式。
b=p1.x-p2.x;
c=-a*p1.x-b*p1.y;
}
bool SegmentProperIntersection(Point a1,Point a2,Point b1,Point b2){
//线段相交判定
double c1=Cross(a2-a1,b1-a1);
double c2=Cross(a2-a1,b2-a1);
double c3=Cross(b2-b1,a1-b1);
double c4=Cross(b2-b1,a2-b1);
return dcmp(c1)*dcmp(c2)<0&&dcmp(c3)*dcmp(c4)<0;
}
bool OnSegment(Point p,Point a1,Point a2){
//判断一个点在线段上
return dcmp(Cross(a1-p,a2-p))==0&&dcmp(Dot(a1-p,a2-p))<0;
}
double ConvexPolygonArea(Point* p,int n){
//从第一个顶点出发把凸多边形分成 n-2 个三角形
double area=0;
for(int i=1;i<n-1;i++)area+=Cross(p[i]-p[0],p[i+1]-p[0]);
return area/2;
}
struct Circle{
Point c;
double r;
Circle(){}
Circle(Point cc,double rr){
c=cc; r=rr;
}
Point point(double a){
//通过圆心角来求坐标
return Point(c.x+cos(a)*r,c.y+sin(a)*r);
}
};
struct Line{
//有向直线,左边就是对应的半平面
Point p,v;
double ang;
Line(){}
Line(Point pp,Point vv){
p=pp;v=vv;
ang=atan2(v.y,v.x);
}
Point point(double a){
return p+(v*a);
}
bool operator<(const Line& L) const{
return ang<L.ang;
}
};
bool OnLeft(Line l,Point pp){
//点Ｐ在有向直线的左边，线上不算
return Cross(l.v,pp-l.p)>0;
}
Point GetLineIntersection(Line a,Line b){
//二直线相交，假设交点唯一存在
Point u=a.p-b.p;
double t=Cross(b.v,u)/Cross(a.v,b.v);
return a.p+a.v*t;
}
int getLineCircleIntersection(Line l,Circle cc,double& t1,double& t2,vector<Point>& sol){
//直线和圆的交点，直线为AB,圆的圆心为C，半径为R，第一种方法是解方程组
//设交点为P＝A+t*(B-A),代入圆方程整理后先得到 (a*t+b)^2+(c*t+d)^2=r^2
//进一步整理后得到一元二次方程 e*t^2+f*t+g=0
//根据判别式的值可以区分3种情况，
//
//函数返回的是交点个数，参数sol存放的是交点本身
//没有清空，可以反复调用函数，把所有的交点放在一个sol里
//---
//另一种方法是几何法
//先求C在AB上的投影P,AB的单位向量V
//交点P－LV，P+LV，由勾股定理算出
double a=l.v.x;
double b=l.p.x-cc.c.x;
double c=l.v.y;
double d=l.p.y-cc.c.y;
double e=a*a+c*c;
double f=2*(a*b+c*d);
double g=b*b+d*d-cc.r*cc.r;
double delta=f*f-4*e*g;//判别式
if(dcmp(delta)<0)return 0;//相离
if(dcmp(delta)==0){
t1=t2=-f/(2*e);
sol.push_back(l.point(t1));
return 1;//相切
}
t1=(-f-sqrt(delta))/(2*e);sol.push_back(l.point(t1));
t2=(-f+sqrt(delta))/(2*e);sol.push_back(l.point(t2));
return 2;//相交
}
double angle(Point v){
//计算向量的极角
return atan2(v.y,v.x);
}
int getCircleCircleIntersection(Circle c1,Circle c2,vector<Point>& sol){ //计算两圆交点
double d=Length(c1.c-c2.c);
if(dcmp(d)==0){
if(dcmp(c1.r-c2.r)==0)return -1;//两圆重合
return 0;
}
if(dcmp(c1.r+c2.r-d)<0)return 0;
if(dcmp(fabs(c1.r-c2.r)-d)>0)return 0;
double a=angle(c2.c-c1.c);//向量C1C2的极角
double da=acos((c1.r*c1.r+d*d-c2.r*c2.r)/(2*c1.r*d));//计算C1C2到C1P1的角
Point p1=c1.point(a-da);Point p2=c1.point(a+da);
sol.push_back(p1);if(p1==p2)return 1;
sol.push_back(p2);
return 2;
}
int getTangents(Point p,Circle c,Point* v){
//过点p到圆C的切线，v[i]是第i条切线的向量，返回切线条数
Point u=c.c-p;
double dist=Length(u);
if(dist<c.r)return 0;
else if(dcmp(dist-c.r)==0){
v[0]=Rotate(u,pi/2);
return 1;
}
else{
double ang=asin(c.r/dist);
v[0]=Rotate(u,-ang);v[1]=Rotate(u,ang);
return 2;
}
}
int getTangents(Circle a,Circle b,Point* aa,Point* bb){
//返回切线条数，-1表示无穷条切线,a[i]和b[i]分别表示第i条切线在圆A和圆B上的切点
int cnt=0;
if(a.r<b.r)swap(a,b), swap(aa,bb);
int d2=(int)((a.c.x-b.c.x)*(a.c.x-b.c.x)+(a.c.y-b.c.y)*(a.c.y-b.c.y));
int rdiff=(int)(a.r-b.r),rsum=(int)(a.r+b.r);
if(d2<rdiff*rdiff)return 0;//内含
double base=atan2(b.c.y-a.c.y,b.c.x-a.c.x);
if(d2==0&&a.r==b.r) return -1;//无穷多条切线
if(d2==rdiff*rdiff){
aa[cnt]=a.point(base);
bb[cnt]=b.point(base);
cnt++;
return 1;//内切，1条
}
//有外公切线
double ang=acos((a.r-b.r)/sqrt(d2));
aa[cnt]=a.point(base+ang);
bb[cnt]=b.point(base+ang);
cnt++;
aa[cnt]=a.point(base-ang);
bb[cnt]=b.point(base-ang);
cnt++;
if(d2==rsum*rsum){
aa[cnt]=a.point(base);
bb[cnt]=b.point(base+pi);
cnt++;//一条内公切线
}
else
if(d2>rsum*rsum){
double ang=acos((a.r+b.r)/sqrt(d2));
aa[cnt]=a.point(base+ang);
bb[cnt]=b.point(pi+base+ang);
cnt++;
aa[cnt]=a.point(base-ang);
bb[cnt]=b.point(pi+base-ang);
cnt++;//两条内公切线
}
return cnt;
}
Circle CircumscribedCircle(Point p1,Point p2,Point p3){
//三角形外接圆
double bx=p2.x-p1.x,by=p2.y-p1.y;
double cx=p3.x-p1.x,cy=p3.y-p1.y;
double d=2*(bx*cy-cx*by);
double ccx=(cy*(bx*bx+by*by)-by*(cx*cx+cy*cy))/d+p1.x;
double ccy=(bx*(cx*cx+cy*cy)-cx*(bx*bx+by*by))/d+p1.y;
Point p=Point(ccx,ccy);
return Circle(p,Length(p1-p));
}
Circle InscribedCircle(Point p1,Point p2,Point p3){
//三角形内切圆
double a=Length(p2-p3),b=Length(p3-p1);
double c=Length(p1-p2);
Point p=(p1*a+p1*b+p3*c)/(a+b+c);
return Circle(p,DistanceToLine(p,p1,p2));
}
double torad(double deg){
//角度转弧度
return deg/180*pi;
}
void get_coord(double r,double lat,double lng,double& x,double& y,double& z){
//经纬度(角度)转化为空间坐标
lat=torad(lat);lng=torad(lng);
x=r*cos(lat)*cos(lng);y=r*cos(lat)*sin(lng);
z=r*sin(lat);
}
//求球面上两点最短距离，走大圆弧
//半径r，弦长d，圆心角 2arcsin(d/2r)
struct Polygon{
int n;
Point p[1000];
Polygon(){ n=0;}
Polygon(int nn,Point* pp){
n=nn;
for(int q=0;q<n;q++) p[q]=pp[q];
}
Point operator[](int a) const{
return p[a];
}
};
int isPointInPolygon(Point p,Polygon po){
//假想有一条向右的射线，统计多边形穿过这条射线正反多少次，把这个数记为绕数
//逆时穿过加1，顺时针穿过减1
int wn=0, n=po.n;
for(int i=0;i<n;i++){
if(OnSegment(p,po[i],po[(i+1)%n]))return -1;//在边界上
int k=dcmp(Cross(po[(i+1)%n]-po[i],p-po[i]));
int d1=dcmp(po[i].y-p.y);
int d2=dcmp(po[(i+1)%n].y-p.y);
if(k>0&&d1<=0&&d2>0)wn++;
if(k<0&&d2<=0&&d1>0)wn--;
}
if(wn!=0) return 1;//内部
return 0;//外部
}
//点在凸多边形内判定更简单，只要判断是否在所有边的左边
//各顶点按逆时针顺序
int ConvexHull(Point* p,int n,Point* ch)
{
//计算凸包，输入点数组p，个数为n，输出点数组ch，函数返回凸包顶点数
//输入不能有重复点，函数执行完后输入点的顺序被破坏
//如果不希望在凸包的边上有输入点，把两个<=改成<
//在精度要求高时建议用dcmp比较
sort(p,p+n);
//n=unique(p,p+n)-p;
int m=0;
for(int i=0;i<n;i++){
while(m>1&&Cross(ch[m-1]-ch[m-2],p[i]-ch[m-2])<=0)m--;
ch[m++]=p[i];
}
int k=m;
for(int i=n-2;i>=0;i--){
while(m>k&&Cross(ch[m-1]-ch[m-2],p[i]-ch[m-2])<=0)m--;
ch[m++]=p[i];
}
if(n>1)m--;
return m;
}
double diameter(Point* o,int n,Point *p){
//求点集直径(最大距离),旋转卡壳
int num=ConvexHull(o,n,p);
int mi=0,ma=0;
for(int q=0;q<num;q++){
if(p[q]<p[mi])mi=q;
if(p[ma]<p[q])ma=q;
}
int v=ma;
double ans=(p[mi].x-p[ma].x)*(p[mi].x-p[ma].x)+(p[mi].y-p[ma].y)*(p[mi].y-p[ma].y);
ans=max(ans,(p[num-1].x-p[ma].x)*(p[num-1].x-p[ma].x)+(p[num-1].y-p[ma].y)*(p[num-1].y-p[ma].y));
ans=max(ans,(p[mi].x-p[ma-1].x)*(p[mi].x-p[ma-1].x)+(p[mi].y-p[ma-1].y)*(p[mi].y-p[ma-1].y));
ans=max(ans,(p[num-1].x-p[ma-1].x)*(p[num-1].x-p[ma-1].x)+(p[num-1].y-p[ma-1].y)*(p[num-1].y-p[ma-1].y));
for(;mi<v;mi++){
for(;;){
int a=dcmp(Cross(p[mi+1]-p[mi],p[ma+1]-p[ma]));
if(a==0){
ans=max(ans,(p[mi].x-p[ma].x)*(p[mi].x-p[ma].x)+(p[mi].y-p[ma].y)*(p[mi].y-p[ma].y));
ans=max(ans,(p[mi+1].x-p[ma].x)*(p[mi+1].x-p[ma].x)+(p[mi+1].y-p[ma].y)*(p[mi+1].y-p[ma].y));
ans=max(ans,(p[mi].x-p[ma+1].x)*(p[mi].x-p[ma+1].x)+(p[mi].y-p[ma+1].y)*(p[mi].y-p[ma+1].y));
ans=max(ans,(p[mi+1].x-p[ma+1].x)*(p[mi+1].x-p[ma+1].x)+(p[mi+1].y-p[ma+1].y)*(p[mi+1].y-p[ma+1].y));
ma++;
break;
}
if(a>0){
ans=max(ans,(p[mi].x-p[ma].x)*(p[mi].x-p[ma].x)+(p[mi].y-p[ma].y)*(p[mi].y-p[ma].y));
ans=max(ans,(p[mi].x-p[ma+1].x)*(p[mi].x-p[ma+1].x)+(p[mi].y-p[ma+1].y)*(p[mi].y-p[ma+1].y));
ma++;
}
if(a<0){
ans=max(ans,(p[mi+1].x-p[ma].x)*(p[mi+1].x-p[ma].x)+(p[mi+1].y-p[ma].y)*(p[mi+1].y-p[ma].y));
ans=max(ans,(p[mi+1].x-p[ma+1].x)*(p[mi+1].x-p[ma+1].x)+(p[mi+1].y-p[ma+1].y)*(p[mi+1].y-p[ma+1].y));
break;
}
}
}
return ans;
}
int HalfplaneIntersection(Line* l,int n,Point* poly){
sort(l,l+n);
int first,last;
Point *p=new Point
;Line *q=new Line
;
q[first=last=0]=l[0];
for(int i=1;i<n;i++){
while(first<last&&!OnLeft(l[i],p[last-1]))last--;
while(first<last&&!OnLeft(l[i],p[first]))first++;
q[++last]=l[i];
if(fabs(Cross(q[last].v,q[last-1].v))<eps){
last--;
if(OnLeft(q[last],l[i].p))q[last]=l[i];
}
if(first<last)p[last-1]=GetLineIntersection(q[last-1],q[last]);
}
while(first<last&&!OnLeft(q[first],p[last-1]))last--;
if(last-first<=1)return 0;
p[last]=GetLineIntersection(q[last],q[first]);
int m=0;
for(int i=first;i<=last;i++)poly[m++]=p[i];
return m;
}
void FindBorder(Point* p,int n){
//寻找多边形边界
int num;
Point ans[10005],z[205];
double x[205];
int a,min1=0;
for(int q=0;q<n;q++)if(p[q]<p[min1])min1=q;
p
=p[0];ans[0]=p[min1];
double s=-pi;num=1;
for(int idx=0;;idx++){
a=0;
for(int q=0;q<n;q++)
if(p[q]==ans[num-1]){
z[a]=p[(q+1)%n];
x[a]=atan2(-p[q].y+z[a].y,-p[q].x+z[a].x)-s-2*pi;
for(;dcmp(x[a])<=0;)x[a]+=2*pi;
a++;z[a]=p[(q+n-1)%n];
x[a]=atan2(-p[q].y+z[a].y,-p[q].x+z[a].x)-s-2*pi;
for(;dcmp(x[a])<=0;)x[a]+=2*pi;
a++;
}
for(int q=0;q<n;q++)
if(OnSegment(ans[num-1],p[q],p[q+1])){
z[a]=p[(q+1)];
x[a]=atan2(-ans[num-1].y+z[a].y,-ans[num-1].x+z[a].x)-s-2*pi;
for(;dcmp(x[a])<=0;)x[a]+=2*pi;
a++;
z[a]=p[q];
x[a]=atan2(-ans[num-1].y+z[a].y,-ans[num-1].x+z[a].x)-s-2*pi;
for(;dcmp(x[a])<=0;)x[a]+=2*pi;
a++;
}
int mi=0;
for(int q=0;q<a;q++)if(x[q]<x[mi])mi=q;
double ma=101000000;
Point f=z[mi],g;
for(int q=0;q<n;q++)
if(SegmentProperIntersection(z[mi],ans[num-1],p[q],p[q+1])){
g=GetLineIntersection(z[mi],z[mi]-ans[num-1],p[q],p[q]-p[q+1]);
if(Length(g-ans[num-1])<ma)ma=Length(g-ans[num-1]),f=g;
}
ans[num]=f;
s=atan2(ans[num-1].y-ans[num].y,ans[num-1].x-ans[num].x);num++;
if(ans[num-1]==ans[0])break;
}
printf("%d\n",num-1);
for(int q=0;q<num-1;q++)
printf("%.4lf %.4lf\n",ans[q].x,ans[q].y);
}
Polygon cutPolygon(Polygon& po,Point a,Point s){
//切割多边形，返回直线左边的部分,要正反调用两次
Point x,c;Polygon ans;
int sum=0;
for(int q=0;q<po.n;q++){
x=po[q];
c=po[(q+1)%po.n];
if(OnLeft(Line(a,s-a),x)||OnSegment(x,a,s))ans.p[sum]=x,sum++;
if(SegmentProperIntersection(x,c,a,s))ans.p[sum]=GetLineIntersection(x,c-x,a,s-a),sum++;
}
ans.n=sum;
return ans;
}
struct Point3{
double x,y,z;
Point3(){
x=0;y=0;z=0;
}
Point3(double xx,double yy,double zz){
x=xx;y=xx;z=zz;
}
bool operator==(Point3 p) const{
return dcmp(x-p.x)==0&&dcmp(y-p.y)==0&&dcmp(z-p.z)==0;
}
};
Point3 operator+(Point3 a,Point3 b){
return Point3(a.x+b.x,a.y+b.y,a.z+b.z);
}
Point3 operator-(Point3 a,Point3 b){
return Point3(a.x-b.x,a.y-b.y,a.z-b.z);
}
Point3 operator*(Point3 a,double p){
return Point3(a.x*p,a.y*p,a.z*p);
}
Point3 operator/(Point3 a,double p){
return Point3(a.x/p,a.y/p,a.z/p);
}
double Dot(Point3 a,Point3 b){
return a.x*b.x+a.y*b.y+a.z*b.y;
}
double Length(Point3 a){
return sqrt(Dot(a,a));
}
double Angle(Point3 a,Point3 b){
return acos(Dot(a,b)/Length(a)/Length(b));
}
//平面用点法式表示(p0,n),p0是平面上一个点，n 是平面法向量
double DistanceToPlane(Point3 &p,Point3 &p0,Point3 &n){
//点p到平面 p0-n 的的距离，n 是单位向量
return fabs(Dot(p-p0,n));
//如果不取绝对值，得到的是有向距离
}
Point3 GetPlaneProjection(Point3 &p,Point3 &p0,Point3 &n){
//点 p 在平面 p0-n 上的投影，n 是单位向量
return p-n*Dot(p-p0,n);
}
Point3 LinePlaneIntersection(Point3 p1,Point3 p2,Point3 p0,Point3 n){
//直线 p1-p2 到平面 p0-n 的交点,分母为 0 对应直线和平面平行，直线在平面上
Point3 v=p2-p1;
double t=(Dot(n,p0-p1)/Dot(n,p2-p1));//判断分母是否为 0
return p1+v*t;//如果是线段，判断 t 是不是在 0 和 1 之间
}
Point3 Cross(Point3 a,Point3 b){
return Point3(a.y*b.z-a.z*b.y,a.z*b.x-a.x*b.z,a.x*b.y-a.y*b.x);
}
double Area2(Point3 a,Point3 b,Point3 c){
return Length(Cross(b-a,c-a));
}
bool PointInTri(Point3 p,Point3 p0,Point3 p1,Point3 p2){
//判断点在三角形内，先判断是否在三角形所在的平面，再用面积关系
double area1=Area2(p,p0,p1);
double area2=Area2(p,p2,p1);
double area3=Area2(p,p0,p2);
return dcmp(area1+area2+area3-Area2(p0,p1,p2))==0;
}
bool TriSegIntersection(Point3 p0,Point3 p1,Point3 p2,Point3 a,Point3 b,Point3 &p){
//判断线段和三角形是否相交,没有考虑线段在平面上的情况
Point3 n=Cross(p1-p0,p2-p0);
if(dcmp(Dot(n,b-a))==0)return false;
else{
double t=Dot(n,p0-a)/Dot(n,b-a);
if(dcmp(t)<0||dcmp(t-1)>0)return false;
p=a+(b-a)*t;//计算交点
return PointInTri(p,p0,p1,p2);//判断交点在三角形内
}
}
double DistanceToLine(Point3 p,Point3 a,Point3 b){
//点 p 到直线 ab 距离
Point3 v1=b-a,v2=p-a;
return Length(Cross(v1,v2))/Length(v1);
}
double DistanceToSegment(Point3 p,Point3 a,Point3 b){
//点 p 到线段 ab 距离
if(a==b)return Length(p-a);
Point3 v1=b-a,v2=p-a,v3=p-b;
if(dcmp(Dot(v1,v2))<0) return Length(v2);
else if(dcmp(Dot(v1,v3))>0) return Length(v3);
else return Length(Cross(v1,v2))/Length(v1);
}
double Volume6(Point3 a,Point3 b,Point3 c,Point3 d){
return Dot(d-a,Cross(b-a,c-a));
}
double TriCir(Point a,Point s,Point d,Circle c,vector<Point> v){
//三角形与圆面积交
//其中圆心是 a ,三角形 abc
double x1=Length(s-a),x2=Length(d-a);
double vv,bb;
if(dcmp(x1-c.r)<=0&&dcmp(x2-c.r)<=0)return Area2(a,s,d)/2;
if(dcmp(x1-c.r)>0&&dcmp(x2-c.r)>0){
if(DistanceToLine(c.c,s,d)>=c.r)return Angle(d-a,s-a)*c.r*c.r/2;
else{
v.clear();
getLineCircleIntersection(Line(s,d-s),c,vv,bb,v);vv=Angle(v[0]-a,v[1]-a);
return min(Angle(v[0]-a,d-a)+Angle(v[1]-a,s-a),Angle(v[1]-a,d-a)+Angle(v[0]-a,s-a))/2*c.r*c.r+sin(Angle(v[0]-a,v[1]-a))/2*c.r*c.r;
}
}
if(dcmp(x1-c.r)>0&&dcmp(x2-c.r)<=0){
v.clear();
getLineCircleIntersection(Line(s,d-s),c,vv,bb,v);
if(OnSegment(v[1],s,d)){
v[0]=v[1];printf("$$");
}
return Area2(a,d,v[0])/2+Angle(v[0]-a,s-a)*c.r*c.r/2;
}
if(dcmp(x1-c.r)<=0&&dcmp(x2-c.r)>0){
v.clear();
getLineCircleIntersection(Line(s,d-s),c,vv,bb,v);
if(OnSegment(v[1],s,d))v[0]=v[1];
return Area2(a,s,v[0])/2+Angle(v[0]-a,d-a)*c.r*c.r/2;
}
}
Point p[600];
double Max(double a,double b){
if(a<=b)return b;
return a;
}
int main(){
//freopen("data.in","r",stdin);
//freopen("data.out","w",stdout);
int n;
while(~scanf("%d",&n)){
for(int i=0;i<n;i++)scanf("%lf%lf",&p[i].x,&p[i].y);
double ans=0.0;
for(int i=0;i<n;i++){
int id=0;
double han=Length(p[i]-p[0]);
for(int j=1;j<n;j++){
if(i==j)continue;
double hh=Length(p[i]-p[j]);
if(hh>han){
han=hh;
id=j;
}
}
for(int j=0;j<n;j++){
if(j==i||j==id)continue;
ans=max(ans,DistanceToLine(p[i],p[id],p[j]));
ans=max(ans,DistanceToLine(p[id],p[i],p[j]));
ans=max(ans,DistanceToLine(p[j],p[i],p[id]));
}
}
printf("%.10lf\n",ans);
}
return 0;
}
方法二代码：
/* ***********************************************
Author :xianxingwuguan
Created Time :2014/2/27 10:15:39
File Name :\代码\acm程序\计算几何.cpp
************************************************ */
#pragma comment(linker, "/STACK:102400000,102400000")
#include <stdio.h>
#include <iostream>
#include <algorithm>
#include <sstream>
#include <stdlib.h>
#include <string.h>
#include <limits.h>
#include <string>
#include <time.h>
#include <math.h>
#include <queue>
#include <stack>
#include <set>
#include <map>
using namespace std;
#define INF 0x3f3f3f3f
#define eps 1e-8
#define pi acos(-1.0)
typedef long long ll;
int dcmp(double x){
if(fabs(x)<eps)return 0;
return x>0?1:-1;
}
struct Point{
double x,y;
Point(double _x=0,double _y=0){
x=_x;y=_y;
}
};
Point operator + (Point a,Point b){
return Point(a.x+b.x,a.y+b.y);
}
Point operator - (Point a,Point b){
return Point(a.x-b.x,a.y-b.y);
}
Point operator * (Point a,double p){
return Point(a.x*p,a.y*p);
}
Point operator / (Point a,double p){
return Point(a.x/p,a.y/p);
}
bool operator < (const Point &a,const Point &b){//重载大小排序，
return a.x<b.x||(a.x==b.x&&a.y<b.y);
}
bool operator == (const Point &a,const Point &b){//重载相等
return dcmp(a.x-b.x)==0&&dcmp(a.y-b.y)==0;
}
double Dot(Point a,Point b){//点积
return a.x*b.x+a.y*b.y;
}
double Cross(Point a,Point b){//叉积
return a.x*b.y-a.y*b.x;
}
double Length(Point a){//向量长度
return sqrt(Dot(a,a));
}
double Angle(Point a,Point b){//a->b顺时针为负，逆时针为正。
return acos(Dot(a,b)/Length(a)/Length(b));
}
Point Rotate(Point a,double rad){ //向量逆时针旋转，rad为弧度。
return Point(a.x*cos(rad)-a.y*sin(rad),a.x*sin(rad)+a.y*cos(rad));
}
double Area2(Point a,Point b,Point c){//三角形面积二倍。
return fabs(Cross(b-a,c-a));
}
Point Normal(Point a){//计算向量的单位法线，即左转90度，长度归一化。
double l=Length(a);//用前确保A不是零向量
return Point(-a.y/l,a.x/l);
}
Point GetLineIntersection(Point P,Point v,Point Q,Point w){
Point u=P-Q;//求解直线P+t*v Q+t*w 的交点
double t=Cross(w,u)/Cross(v,w);//用前确保两条直线有唯一交点
return P+v*t;//当且仅当cross(v,w)!=0
}
double DistanceToLine(Point p,Point a,Point b){
Point v1=b-a;//计算点P到直线AB距离
Point v2=p-a;//如果不取绝对值，得的是有向距离
return fabs(Cross(v1,v2))/(Length(v1));
}
double DistanceToSegment(Point p,Point a,Point b){
if(a==b)return Length(p-a);//点到线段AB的距离。
Point v1=b-a,v2=p-a,v3=p-b;
if(dcmp(Dot(v1,v2))<0)return Length(v2);
else if(dcmp(Dot(v1,v3))>0)return Length(v3);
else return fabs(Cross(v1,v2))/Length(v1);
}
Point GetLineProjection(Point p,Point a,Point b){
//求点在直线上投影
//AB直线:A+tv , 投影Q的参数为 t0
//由 Dot(v,P-(A+t0*v))==0 得出
Point v=b-a;
return a+v*(Dot(v,p-a)/Dot(v,v));
}
void getLineGenralEquation(Point p1,Point p2,double &a,double &b,double &c){
a=p2.y-p1.y;//直线两点式转一般式。
b=p1.x-p2.x;
c=-a*p1.x-b*p1.y;
}
bool SegmentProperIntersection(Point a1,Point a2,Point b1,Point b2){
//线段相交判定
double c1=Cross(a2-a1,b1-a1);
double c2=Cross(a2-a1,b2-a1);
double c3=Cross(b2-b1,a1-b1);
double c4=Cross(b2-b1,a2-b1);
return dcmp(c1)*dcmp(c2)<0&&dcmp(c3)*dcmp(c4)<0;
}
bool OnSegment(Point p,Point a1,Point a2){
//判断一个点在线段上
return dcmp(Cross(a1-p,a2-p))==0&&dcmp(Dot(a1-p,a2-p))<0;
}
double ConvexPolygonArea(Point* p,int n){
//从第一个顶点出发把凸多边形分成 n-2 个三角形
double area=0;
for(int i=1;i<n-1;i++)area+=Cross(p[i]-p[0],p[i+1]-p[0]);
return area/2;
}
struct Circle{
Point c;
double r;
Circle(){}
Circle(Point cc,double rr){
c=cc; r=rr;
}
Point point(double a){
//通过圆心角来求坐标
return Point(c.x+cos(a)*r,c.y+sin(a)*r);
}
};
struct Line{
//有向直线,左边就是对应的半平面
Point p,v;
double ang;
Line(){}
Line(Point pp,Point vv){
p=pp;v=vv;
ang=atan2(v.y,v.x);
}
Point point(double a){
return p+(v*a);
}
bool operator<(const Line& L) const{
return ang<L.ang;
}
};
bool OnLeft(Line l,Point pp){
//点Ｐ在有向直线的左边，线上不算
return Cross(l.v,pp-l.p)>0;
}
Point GetLineIntersection(Line a,Line b){
//二直线相交，假设交点唯一存在
Point u=a.p-b.p;
double t=Cross(b.v,u)/Cross(a.v,b.v);
return a.p+a.v*t;
}
int getLineCircleIntersection(Line l,Circle cc,double& t1,double& t2,vector<Point>& sol){
//直线和圆的交点，直线为AB,圆的圆心为C，半径为R，第一种方法是解方程组
//设交点为P＝A+t*(B-A),代入圆方程整理后先得到 (a*t+b)^2+(c*t+d)^2=r^2
//进一步整理后得到一元二次方程 e*t^2+f*t+g=0
//根据判别式的值可以区分3种情况，
//
//函数返回的是交点个数，参数sol存放的是交点本身
//没有清空，可以反复调用函数，把所有的交点放在一个sol里
//---
//另一种方法是几何法
//先求C在AB上的投影P,AB的单位向量V
//交点P－LV，P+LV，由勾股定理算出
double a=l.v.x;
double b=l.p.x-cc.c.x;
double c=l.v.y;
double d=l.p.y-cc.c.y;
double e=a*a+c*c;
double f=2*(a*b+c*d);
double g=b*b+d*d-cc.r*cc.r;
double delta=f*f-4*e*g;//判别式
if(dcmp(delta)<0)return 0;//相离
if(dcmp(delta)==0){
t1=t2=-f/(2*e);
sol.push_back(l.point(t1));
return 1;//相切
}
t1=(-f-sqrt(delta))/(2*e);sol.push_back(l.point(t1));
t2=(-f+sqrt(delta))/(2*e);sol.push_back(l.point(t2));
return 2;//相交
}
double angle(Point v){
//计算向量的极角
return atan2(v.y,v.x);
}
int getCircleCircleIntersection(Circle c1,Circle c2,vector<Point>& sol){ //计算两圆交点
double d=Length(c1.c-c2.c);
if(dcmp(d)==0){
if(dcmp(c1.r-c2.r)==0)return -1;//两圆重合
return 0;
}
if(dcmp(c1.r+c2.r-d)<0)return 0;
if(dcmp(fabs(c1.r-c2.r)-d)>0)return 0;
double a=angle(c2.c-c1.c);//向量C1C2的极角
double da=acos((c1.r*c1.r+d*d-c2.r*c2.r)/(2*c1.r*d));//计算C1C2到C1P1的角
Point p1=c1.point(a-da);Point p2=c1.point(a+da);
sol.push_back(p1);if(p1==p2)return 1;
sol.push_back(p2);
return 2;
}
int getTangents(Point p,Circle c,Point* v){
//过点p到圆C的切线，v[i]是第i条切线的向量，返回切线条数
Point u=c.c-p;
double dist=Length(u);
if(dist<c.r)return 0;
else if(dcmp(dist-c.r)==0){
v[0]=Rotate(u,pi/2);
return 1;
}
else{
double ang=asin(c.r/dist);
v[0]=Rotate(u,-ang);v[1]=Rotate(u,ang);
return 2;
}
}
int getTangents(Circle a,Circle b,Point* aa,Point* bb){
//返回切线条数，-1表示无穷条切线,a[i]和b[i]分别表示第i条切线在圆A和圆B上的切点
int cnt=0;
if(a.r<b.r)swap(a,b), swap(aa,bb);
int d2=(int)((a.c.x-b.c.x)*(a.c.x-b.c.x)+(a.c.y-b.c.y)*(a.c.y-b.c.y));
int rdiff=(int)(a.r-b.r),rsum=(int)(a.r+b.r);
if(d2<rdiff*rdiff)return 0;//内含
double base=atan2(b.c.y-a.c.y,b.c.x-a.c.x);
if(d2==0&&a.r==b.r) return -1;//无穷多条切线
if(d2==rdiff*rdiff){
aa[cnt]=a.point(base);
bb[cnt]=b.point(base);
cnt++;
return 1;//内切，1条
}
//有外公切线
double ang=acos((a.r-b.r)/sqrt(d2));
aa[cnt]=a.point(base+ang);
bb[cnt]=b.point(base+ang);
cnt++;
aa[cnt]=a.point(base-ang);
bb[cnt]=b.point(base-ang);
cnt++;
if(d2==rsum*rsum){
aa[cnt]=a.point(base);
bb[cnt]=b.point(base+pi);
cnt++;//一条内公切线
}
else
if(d2>rsum*rsum){
double ang=acos((a.r+b.r)/sqrt(d2));
aa[cnt]=a.point(base+ang);
bb[cnt]=b.point(pi+base+ang);
cnt++;
aa[cnt]=a.point(base-ang);
bb[cnt]=b.point(pi+base-ang);
cnt++;//两条内公切线
}
return cnt;
}
Circle CircumscribedCircle(Point p1,Point p2,Point p3){
//三角形外接圆
double bx=p2.x-p1.x,by=p2.y-p1.y;
double cx=p3.x-p1.x,cy=p3.y-p1.y;
double d=2*(bx*cy-cx*by);
double ccx=(cy*(bx*bx+by*by)-by*(cx*cx+cy*cy))/d+p1.x;
double ccy=(bx*(cx*cx+cy*cy)-cx*(bx*bx+by*by))/d+p1.y;
Point p=Point(ccx,ccy);
return Circle(p,Length(p1-p));
}
Circle InscribedCircle(Point p1,Point p2,Point p3){
//三角形内切圆
double a=Length(p2-p3),b=Length(p3-p1);
double c=Length(p1-p2);
Point p=(p1*a+p1*b+p3*c)/(a+b+c);
return Circle(p,DistanceToLine(p,p1,p2));
}
double torad(double deg){
//角度转弧度
return deg/180*pi;
}
void get_coord(double r,double lat,double lng,double& x,double& y,double& z){
//经纬度(角度)转化为空间坐标
lat=torad(lat);lng=torad(lng);
x=r*cos(lat)*cos(lng);y=r*cos(lat)*sin(lng);
z=r*sin(lat);
}
//求球面上两点最短距离，走大圆弧
//半径r，弦长d，圆心角 2arcsin(d/2r)
struct Polygon{
int n;
Point p[1000];
Polygon(){ n=0;}
Polygon(int nn,Point* pp){
n=nn;
for(int q=0;q<n;q++) p[q]=pp[q];
}
Point operator[](int a) const{
return p[a];
}
};
int isPointInPolygon(Point p,Polygon po){
//假想有一条向右的射线，统计多边形穿过这条射线正反多少次，把这个数记为绕数
//逆时穿过加1，顺时针穿过减1
int wn=0, n=po.n;
for(int i=0;i<n;i++){
if(OnSegment(p,po[i],po[(i+1)%n]))return -1;//在边界上
int k=dcmp(Cross(po[(i+1)%n]-po[i],p-po[i]));
int d1=dcmp(po[i].y-p.y);
int d2=dcmp(po[(i+1)%n].y-p.y);
if(k>0&&d1<=0&&d2>0)wn++;
if(k<0&&d2<=0&&d1>0)wn--;
}
if(wn!=0) return 1;//内部
return 0;//外部
}
//点在凸多边形内判定更简单，只要判断是否在所有边的左边
//各顶点按逆时针顺序
int ConvexHull(Point* p,int n,Point* ch)
{
//计算凸包，输入点数组p，个数为n，输出点数组ch，函数返回凸包顶点数
//输入不能有重复点，函数执行完后输入点的顺序被破坏
//如果不希望在凸包的边上有输入点，把两个<=改成<
//在精度要求高时建议用dcmp比较
sort(p,p+n);
//n=unique(p,p+n)-p;
int m=0;
for(int i=0;i<n;i++){
while(m>1&&Cross(ch[m-1]-ch[m-2],p[i]-ch[m-2])<=0)m--;
ch[m++]=p[i];
}
int k=m;
for(int i=n-2;i>=0;i--){
while(m>k&&Cross(ch[m-1]-ch[m-2],p[i]-ch[m-2])<=0)m--;
ch[m++]=p[i];
}
if(n>1)m--;
return m;
}
double diameter(Point* o,int n,Point *p){
//求点集直径(最大距离),旋转卡壳
int num=ConvexHull(o,n,p);
int mi=0,ma=0;
for(int q=0;q<num;q++){
if(p[q]<p[mi])mi=q;
if(p[ma]<p[q])ma=q;
}
int v=ma;
double ans=(p[mi].x-p[ma].x)*(p[mi].x-p[ma].x)+(p[mi].y-p[ma].y)*(p[mi].y-p[ma].y);
ans=max(ans,(p[num-1].x-p[ma].x)*(p[num-1].x-p[ma].x)+(p[num-1].y-p[ma].y)*(p[num-1].y-p[ma].y));
ans=max(ans,(p[mi].x-p[ma-1].x)*(p[mi].x-p[ma-1].x)+(p[mi].y-p[ma-1].y)*(p[mi].y-p[ma-1].y));
ans=max(ans,(p[num-1].x-p[ma-1].x)*(p[num-1].x-p[ma-1].x)+(p[num-1].y-p[ma-1].y)*(p[num-1].y-p[ma-1].y));
for(;mi<v;mi++){
for(;;){
int a=dcmp(Cross(p[mi+1]-p[mi],p[ma+1]-p[ma]));
if(a==0){
ans=max(ans,(p[mi].x-p[ma].x)*(p[mi].x-p[ma].x)+(p[mi].y-p[ma].y)*(p[mi].y-p[ma].y));
ans=max(ans,(p[mi+1].x-p[ma].x)*(p[mi+1].x-p[ma].x)+(p[mi+1].y-p[ma].y)*(p[mi+1].y-p[ma].y));
ans=max(ans,(p[mi].x-p[ma+1].x)*(p[mi].x-p[ma+1].x)+(p[mi].y-p[ma+1].y)*(p[mi].y-p[ma+1].y));
ans=max(ans,(p[mi+1].x-p[ma+1].x)*(p[mi+1].x-p[ma+1].x)+(p[mi+1].y-p[ma+1].y)*(p[mi+1].y-p[ma+1].y));
ma++;
break;
}
if(a>0){
ans=max(ans,(p[mi].x-p[ma].x)*(p[mi].x-p[ma].x)+(p[mi].y-p[ma].y)*(p[mi].y-p[ma].y));
ans=max(ans,(p[mi].x-p[ma+1].x)*(p[mi].x-p[ma+1].x)+(p[mi].y-p[ma+1].y)*(p[mi].y-p[ma+1].y));
ma++;
}
if(a<0){
ans=max(ans,(p[mi+1].x-p[ma].x)*(p[mi+1].x-p[ma].x)+(p[mi+1].y-p[ma].y)*(p[mi+1].y-p[ma].y));
ans=max(ans,(p[mi+1].x-p[ma+1].x)*(p[mi+1].x-p[ma+1].x)+(p[mi+1].y-p[ma+1].y)*(p[mi+1].y-p[ma+1].y));
break;
}
}
}
return ans;
}
int HalfplaneIntersection(Line* l,int n,Point* poly){
sort(l,l+n);
int first,last;
Point *p=new Point
;Line *q=new Line
;
q[first=last=0]=l[0];
for(int i=1;i<n;i++){
while(first<last&&!OnLeft(l[i],p[last-1]))last--;
while(first<last&&!OnLeft(l[i],p[first]))first++;
q[++last]=l[i];
if(fabs(Cross(q[last].v,q[last-1].v))<eps){
last--;
if(OnLeft(q[last],l[i].p))q[last]=l[i];
}
if(first<last)p[last-1]=GetLineIntersection(q[last-1],q[last]);
}
while(first<last&&!OnLeft(q[first],p[last-1]))last--;
if(last-first<=1)return 0;
p[last]=GetLineIntersection(q[last],q[first]);
int m=0;
for(int i=first;i<=last;i++)poly[m++]=p[i];
return m;
}
void FindBorder(Point* p,int n){
//寻找多边形边界
int num;
Point ans[10005],z[205];
double x[205];
int a,min1=0;
for(int q=0;q<n;q++)if(p[q]<p[min1])min1=q;
p
=p[0];ans[0]=p[min1];
double s=-pi;num=1;
for(int idx=0;;idx++){
a=0;
for(int q=0;q<n;q++)
if(p[q]==ans[num-1]){
z[a]=p[(q+1)%n];
x[a]=atan2(-p[q].y+z[a].y,-p[q].x+z[a].x)-s-2*pi;
for(;dcmp(x[a])<=0;)x[a]+=2*pi;
a++;z[a]=p[(q+n-1)%n];
x[a]=atan2(-p[q].y+z[a].y,-p[q].x+z[a].x)-s-2*pi;
for(;dcmp(x[a])<=0;)x[a]+=2*pi;
a++;
}
for(int q=0;q<n;q++)
if(OnSegment(ans[num-1],p[q],p[q+1])){
z[a]=p[(q+1)];
x[a]=atan2(-ans[num-1].y+z[a].y,-ans[num-1].x+z[a].x)-s-2*pi;
for(;dcmp(x[a])<=0;)x[a]+=2*pi;
a++;
z[a]=p[q];
x[a]=atan2(-ans[num-1].y+z[a].y,-ans[num-1].x+z[a].x)-s-2*pi;
for(;dcmp(x[a])<=0;)x[a]+=2*pi;
a++;
}
int mi=0;
for(int q=0;q<a;q++)if(x[q]<x[mi])mi=q;
double ma=101000000;
Point f=z[mi],g;
for(int q=0;q<n;q++)
if(SegmentProperIntersection(z[mi],ans[num-1],p[q],p[q+1])){
g=GetLineIntersection(z[mi],z[mi]-ans[num-1],p[q],p[q]-p[q+1]);
if(Length(g-ans[num-1])<ma)ma=Length(g-ans[num-1]),f=g;
}
ans[num]=f;
s=atan2(ans[num-1].y-ans[num].y,ans[num-1].x-ans[num].x);num++;
if(ans[num-1]==ans[0])break;
}
printf("%d\n",num-1);
for(int q=0;q<num-1;q++)
printf("%.4lf %.4lf\n",ans[q].x,ans[q].y);
}
Polygon cutPolygon(Polygon& po,Point a,Point s){
//切割多边形，返回直线左边的部分,要正反调用两次
Point x,c;Polygon ans;
int sum=0;
for(int q=0;q<po.n;q++){
x=po[q];
c=po[(q+1)%po.n];
if(OnLeft(Line(a,s-a),x)||OnSegment(x,a,s))ans.p[sum]=x,sum++;
if(SegmentProperIntersection(x,c,a,s))ans.p[sum]=GetLineIntersection(x,c-x,a,s-a),sum++;
}
ans.n=sum;
return ans;
}
struct Point3{
double x,y,z;
Point3(){
x=0;y=0;z=0;
}
Point3(double xx,double yy,double zz){
x=xx;y=xx;z=zz;
}
bool operator==(Point3 p) const{
return dcmp(x-p.x)==0&&dcmp(y-p.y)==0&&dcmp(z-p.z)==0;
}
};
Point3 operator+(Point3 a,Point3 b){
return Point3(a.x+b.x,a.y+b.y,a.z+b.z);
}
Point3 operator-(Point3 a,Point3 b){
return Point3(a.x-b.x,a.y-b.y,a.z-b.z);
}
Point3 operator*(Point3 a,double p){
return Point3(a.x*p,a.y*p,a.z*p);
}
Point3 operator/(Point3 a,double p){
return Point3(a.x/p,a.y/p,a.z/p);
}
double Dot(Point3 a,Point3 b){
return a.x*b.x+a.y*b.y+a.z*b.y;
}
double Length(Point3 a){
return sqrt(Dot(a,a));
}
double Angle(Point3 a,Point3 b){
return acos(Dot(a,b)/Length(a)/Length(b));
}
//平面用点法式表示(p0,n),p0是平面上一个点，n 是平面法向量
double DistanceToPlane(Point3 &p,Point3 &p0,Point3 &n){
//点p到平面 p0-n 的的距离，n 是单位向量
return fabs(Dot(p-p0,n));
//如果不取绝对值，得到的是有向距离
}
Point3 GetPlaneProjection(Point3 &p,Point3 &p0,Point3 &n){
//点 p 在平面 p0-n 上的投影，n 是单位向量
return p-n*Dot(p-p0,n);
}
Point3 LinePlaneIntersection(Point3 p1,Point3 p2,Point3 p0,Point3 n){
//直线 p1-p2 到平面 p0-n 的交点,分母为 0 对应直线和平面平行，直线在平面上
Point3 v=p2-p1;
double t=(Dot(n,p0-p1)/Dot(n,p2-p1));//判断分母是否为 0
return p1+v*t;//如果是线段，判断 t 是不是在 0 和 1 之间
}
Point3 Cross(Point3 a,Point3 b){
return Point3(a.y*b.z-a.z*b.y,a.z*b.x-a.x*b.z,a.x*b.y-a.y*b.x);
}
double Area2(Point3 a,Point3 b,Point3 c){
return Length(Cross(b-a,c-a));
}
bool PointInTri(Point3 p,Point3 p0,Point3 p1,Point3 p2){
//判断点在三角形内，先判断是否在三角形所在的平面，再用面积关系
double area1=Area2(p,p0,p1);
double area2=Area2(p,p2,p1);
double area3=Area2(p,p0,p2);
return dcmp(area1+area2+area3-Area2(p0,p1,p2))==0;
}
bool TriSegIntersection(Point3 p0,Point3 p1,Point3 p2,Point3 a,Point3 b,Point3 &p){
//判断线段和三角形是否相交,没有考虑线段在平面上的情况
Point3 n=Cross(p1-p0,p2-p0);
if(dcmp(Dot(n,b-a))==0)return false;
else{
double t=Dot(n,p0-a)/Dot(n,b-a);
if(dcmp(t)<0||dcmp(t-1)>0)return false;
p=a+(b-a)*t;//计算交点
return PointInTri(p,p0,p1,p2);//判断交点在三角形内
}
}
double DistanceToLine(Point3 p,Point3 a,Point3 b){
//点 p 到直线 ab 距离
Point3 v1=b-a,v2=p-a;
return Length(Cross(v1,v2))/Length(v1);
}
double DistanceToSegment(Point3 p,Point3 a,Point3 b){
//点 p 到线段 ab 距离
if(a==b)return Length(p-a);
Point3 v1=b-a,v2=p-a,v3=p-b;
if(dcmp(Dot(v1,v2))<0) return Length(v2);
else if(dcmp(Dot(v1,v3))>0) return Length(v3);
else return Length(Cross(v1,v2))/Length(v1);
}
double Volume6(Point3 a,Point3 b,Point3 c,Point3 d){
return Dot(d-a,Cross(b-a,c-a));
}
double TriCir(Point a,Point s,Point d,Circle c,vector<Point> v){
//三角形与圆面积交
//其中圆心是 a ,三角形 abc
double x1=Length(s-a),x2=Length(d-a);
double vv,bb;
if(dcmp(x1-c.r)<=0&&dcmp(x2-c.r)<=0)return Area2(a,s,d)/2;
if(dcmp(x1-c.r)>0&&dcmp(x2-c.r)>0){
if(DistanceToLine(c.c,s,d)>=c.r)return Angle(d-a,s-a)*c.r*c.r/2;
else{
v.clear();
getLineCircleIntersection(Line(s,d-s),c,vv,bb,v);vv=Angle(v[0]-a,v[1]-a);
return min(Angle(v[0]-a,d-a)+Angle(v[1]-a,s-a),Angle(v[1]-a,d-a)+Angle(v[0]-a,s-a))/2*c.r*c.r+sin(Angle(v[0]-a,v[1]-a))/2*c.r*c.r;
}
}
if(dcmp(x1-c.r)>0&&dcmp(x2-c.r)<=0){
v.clear();
getLineCircleIntersection(Line(s,d-s),c,vv,bb,v);
if(OnSegment(v[1],s,d)){
v[0]=v[1];printf("$$");
}
return Area2(a,d,v[0])/2+Angle(v[0]-a,s-a)*c.r*c.r/2;
}
if(dcmp(x1-c.r)<=0&&dcmp(x2-c.r)>0){
v.clear();
getLineCircleIntersection(Line(s,d-s),c,vv,bb,v);
if(OnSegment(v[1],s,d))v[0]=v[1];
return Area2(a,s,v[0])/2+Angle(v[0]-a,d-a)*c.r*c.r/2;
}
}
Point p[600],ch[600];
Line rot[600*500];
double fun(int k,int n,int m){
double ret=0;
int cnt=0;
for(int i=0;i<n;i++){
if(i==k)continue;
rot[cnt]=Line(p[k],p[i]-p[k]);
++cnt;
}
sort(rot,rot+cnt);
int j=0;
while(1){
if(dcmp(Cross(rot[0].v,ch[j]-p[k]))>=0&&dcmp(Cross(rot[0].v,ch[(j+m-1)%m]-p[k]))<=0)break;
j=(j+1)%m;
}
for(int i=0;i<cnt;i++){
Point tmp=rot[i].v+rot[i].p;
while(DistanceToLine(ch[j+1],p[k],tmp)>DistanceToLine(ch[j],p[k],tmp)||dcmp(DistanceToLine(ch[j],p[k],tmp))==0)j=(j+1)%m;
ret=max(ret,DistanceToLine(ch[j],p[k],tmp));
}
return ret;
}
int main(){
//freopen("data.in","r",stdin);
//freopen("data.out","w",stdout);
int n;
while(~scanf("%d",&n)){
for(int i=0;i<n;i++)scanf("%lf%lf",&p[i].x,&p[i].y);
int m=ConvexHull(p,n,ch);
double ans=0;
for(int i=0;i<n;i++)ans=max(ans,fun(i,n,m));
printf("%.10lf\n",ans);
}
return 0;
}
G题:暂时没做，有时间补。
H题：状压枚举
代码：#include<stdio.h>
#include<string.h>
#include<time.h>
#include<math.h>
#include<algorithm>
#include<iostream>
#include<vector>
#include<set>
#include<map>
using namespace std;
typedef long long lld;
#define X first
#define Y second
#define mp make_pair
#define pb push_back
struct E
{
int x,y;
lld s;
void in()
{
cin >> x >> y >> s;
x--;
y--;
}
}e[100010];
lld pp[100010];
int qq;
int main()
{
int n,m,k;
while(scanf("%d %d %d",&n,&m,&k)!=EOF)
{
for(int i=0;i<m;i++)
e[i].in();
lld ans=(1LL)<<60;
for(int mask=0;mask<(1<<n);mask++)
{
if(!(mask&(1<<0)))
continue;
if(mask&(1<<(n-1)))
continue;
qq=0;
for(int i=0;i<m;i++)
{
int x=e[i].x;
int y=e[i].y;
bool flag=false;
if((mask&(1<<x)) && !(mask&(1<<y)))
flag=true;
swap(x,y);
if((mask&(1<<x)) && !(mask&(1<<y)))
flag=true;
if(flag)
pp[qq++]=e[i].s;
}
sort(pp,pp+qq);
lld now=0;
for(int i=0;i<qq-k;i++)
now+=pp[i];
ans=min(ans,now);
}
cout << ans << endl;
}
return 0;
}
I题:splay裸题，自己没实现出来，贴一个baka巨巨的，慢慢学习。
代码：
#include <iostream>
#include<cstdio>
#include<cstring>
#include<string>
#include<algorithm>
#include<map>
#include<queue>
#include<stack>
#include<cmath>
#include<vector>
#define inf 0x3f3f3f3f
#define Inf 0x3FFFFFFFFFFFFFFFLL
#define eps 1e-9
#define pi acos(-1.0)
using namespace std;
typedef long long ll;
const int maxn=300000+10;
int ch[maxn][2],pre[maxn],size[maxn],val[maxn],st[maxn];
int gcdst[maxn][2],root,tot;
int num[200010],fs[200010];
int gcd(int a,int b)
{
if(a==-1) return b;
if(b==-1) return a;
int c;
while(b)
{
c=a%b;
a=b;
b=c;
}
return a;
}
void NewNode(int &rt,int fa,int v,int state)
{
rt=++tot;
ch[rt][0]=ch[rt][1]=0;
size[rt]=1;val[rt]=v;
pre[rt]=fa;
gcdst[rt][state]=v;gcdst[rt][state^1]=-1;
st[rt]=state;
}
void PushUp(int rt)
{
size[rt]=size[ch[rt][0]]+size[ch[rt][1]]+1;
gcdst[rt][0]=gcd(gcdst[ch[rt][0]][0],gcdst[ch[rt][1]][0]);
gcdst[rt][1]=gcd(gcdst[ch[rt][0]][1],gcdst[ch[rt][1]][1]);
gcdst[rt][st[rt]]=gcd(gcdst[rt][st[rt]],val[rt]);
}
void Rotate(int x,int kind)
{
int y=pre[x];
ch[y][kind^1]=ch[x][kind];
pre[ch[x][kind]]=y;
if(pre[y])
ch[pre[y]][ch[pre[y]][1]==y]=x;
pre[x]=pre[y];
ch[x][kind]=y;
pre[y]=x;
PushUp(y);
}
void Splay(int rt,int goal)
{
while(pre[rt]!=goal)
{
int y=pre[rt];
if(pre[y]==goal)
Rotate(rt,ch[y][0]==rt);
else
{
int kind=ch[pre[y]][0]==y;
if(ch[y][kind]==rt)
{
Rotate(rt,kind^1);
Rotate(rt,kind);
}
else
{
Rotate(y,kind);
Rotate(rt,kind);
}
}
}
PushUp(rt);
if(goal==0) root=rt;
}
void RotateTo(int k,int goal)
{
int rt=root;
while(size[ch[rt][0]]!=k)
{
if(size[ch[rt][0]]>k)
rt=ch[rt][0];
else
{
k-=(size[ch[rt][0]]+1);
rt=ch[rt][1];
}
}
Splay(rt,goal);
}
void build(int l,int r,int &rt,int fa)
{
if(l>r) return ;
int m=(l+r)>>1;
NewNode(rt,fa,num[m],fs[m]);
build(l,m-1,ch[rt][0],rt);
build(m+1,r,ch[rt][1],rt);
PushUp(rt);
}
int Query(int L,int R,int state)
{
RotateTo(L-1,0);
RotateTo(R+1,root);
return gcdst[ch[ch[root][1]][0]][state];
}
void Insert(int pos,int v,int state)
{
RotateTo(pos,0);
if(ch[root][0]==0)
{
NewNode(ch[root][0],root,v,state);
PushUp(root);
return ;
}
int rc=ch[root][0];
while(ch[rc][1])
rc=ch[rc][1];
Splay(rc,root);
NewNode(ch[rc][1],rc,v,state);
PushUp(rc);
PushUp(root);
}
void Delete(int pos)
{
RotateTo(pos,0);
if(ch[root][0]==0)
{
root=ch[root][1];
pre[root]=0;
return ;
}
int rc=ch[root][0];
while(ch[rc][1]) rc=ch[rc][1];
Splay(rc,root);
int rt=ch[root][1];
ch[rc][1]=rt;
pre[rt]=rc;
root=rc;
pre[rc]=0;
PushUp(root);
}
void Changes(int pos)
{
RotateTo(pos,0);
st[root]^=1;
PushUp(root);
}
void Modify(int pos,int v)
{
RotateTo(pos,0);
val[root]=v;
PushUp(root);
}
void Init(int n)
{
for(int i=0;i<n;++i)
scanf("%d%d",&num[i],&fs[i]);
ch[0][0]=ch[0][1]=0;
pre[0]=size[0]=0;val[0]=-1;
gcdst[0][0]=gcdst[0][1]=-1;
root=tot=0;
NewNode(root,0,-1,0);
NewNode(ch[root][1],root,-1,0);
build(0,n-1,ch[ch[root][1]][0],ch[root][1]);
PushUp(ch[root][1]);
PushUp(root);
}
int main()
{
//freopen("in.txt","r",stdin);
//freopen("out.txt","w",stdout);
int n,q;
char str[5];
while(~scanf("%d%d",&n,&q))
{
Init(n);
while(q--)
{
scanf("%s",str);
int L,R,x,s;
if(str[0]=='Q')
{
scanf("%d%d%d",&L,&R,&s);
printf("%d\n",Query(L,R,s));
}
else if(str[0]=='I')
{
scanf("%d%d%d",&L,&x,&s);
Insert(L+1,x,s);
}
else if(str[0]=='D')
{
scanf("%d",&L);
Delete(L);
}
else if(str[0]=='R')
{
scanf("%d",&L);
Changes(L);
}
else
{
scanf("%d%d",&L,&x);
Modify(L,x);
}
}
}
return 0;
}
J题：没做，有时间补。
内容来自用户分享和网络整理，不保证内容的准确性，如有侵权内容，可联系管理员处理 
相关文章推荐
- ZOJ 3723 (浙大月赛)状压DP
- zoj 3802 Easy 2048 Again (14.8 浙大月赛 E)
- 【ZOJ3921 2016年浙大2月月赛G】【爆搜+智商贪心剪枝 好题】Guan Dan 2副牌取27张最少出光牌次数
- ZOJ 3723 (浙大月赛)状压DP
- 【ZOJ3921 2016年浙大2月月赛I】【物理题 底乘高积分式思维】In the Rain 圆柱体人行走 吸收雨体积
- HDU/HDOJ 1395 ACM浙大月赛 2^x mod n = 1
- 记次浙大月赛 134 - ZOJ Monthly, June 2014
- 【ZOJ3921 2016年浙大2月月赛J】【题目观察 暴力 大组合数】Musical Notes n个数，数值为1 2 4 8 16，和恰为m方案数
- ZOJ Monthly, July 2012浙大月赛解题报告
- 【ZOJ3921 2016年浙大2月月赛L】【随机化 答案猜解】Parity Modulo P 大数mod大数mod2
- HDU/HDOJ 1396 浙大月赛2003年
- Zoj 3647 Gao the Grid （ 2012年浙大9月月赛）组合数学
- 浙大2014 3月月赛I题(伸展树)
- 浙大月赛C题（2012/8）Cinema in Akiba（线段树）
- 浙大五月月赛E：P-norm（输入处理）
- ZOJ 3890 Wumpus （2015年7月浙大月赛、广搜）
- 浙大月赛 ZOJ Monthly, March 2014（简单题的题解）
- zoj_3882_7/27浙大月赛B
- zoj_3888_7/27浙大月赛H
- HDU/HDOJ 1395 ACM浙大月赛 2^x mod n = 1 可暴力 可用欧拉定理