标签： 模板
[模板]圆与多边形交的面积
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 |
#include<stdio.h> #include<string.h> #include<stdlib.h> #include<math.h> #include<algorithm> const double eps = 1e-8; const double pi = acos(-1.0); int dcmp(double x) { if(x > eps) return 1; return x < -eps ? -1 : 0; } struct Point { double x, y; Point(){x = y = 0;} Point(double a, double b) {x = a, y = b;} inline void read() {scanf("%lf%lf", &x, &y);} inline Point operator-(const Point &b)const {return Point(x - b.x, y - b.y);} inline Point operator+(const Point &b)const {return Point(x + b.x, y + b.y);} inline Point operator*(const double &b)const {return Point(x * b, y * b);} inline double dot(const Point &b)const {return x * b.x + y * b.y;} inline double cross(const Point &b, const Point &c)const {return (b.x - x) * (c.y - y) - (c.x - x) * (b.y - y);} inline double Dis(const Point &b)const {return sqrt((*this - b).dot(*this - b));} inline bool InLine(const Point &b, const Point &c)const//三点共线 {return !dcmp(cross(b, c));} inline bool OnSeg(const Point &b, const Point &c)const//点在线段上，包括端点 {return InLine(b, c) && (*this - c).dot(*this - b) < eps;} }; inline double min(double a, double b) {return a < b ? a : b;} inline double max(double a, double b) {return a > b ? a : b;} inline double Sqr(double x) {return x * x;} inline double Sqr(const Point &p) {return p.dot(p);} Point LineCross(const Point &a, const Point &b, const Point &c, const Point &d) { double u = a.cross(b, c), v = b.cross(a, d); return Point((c.x * v + d.x * u) / (u + v), (c.y * v + d.y * u) / (u + v)); } double LineCrossCircle(const Point &a, const Point &b, const Point &r, double R, Point &p1, Point &p2) { Point fp = LineCross(r, Point(r.x + a.y - b.y, r.y + b.x - a.x), a, b); double rtol = r.Dis(fp); double rtos = fp.OnSeg(a, b) ? rtol : min(r.Dis(a), r.Dis(b)); double atob = a.Dis(b); double fptoe = sqrt(R * R - rtol * rtol) / atob; if(rtos > R - eps) return rtos; p1 = fp + (a - b) * fptoe; p2 = fp + (b - a) * fptoe; return rtos; } double SectorArea(const Point &r, const Point &a, const Point &b, double R) //不大于180度扇形面积，r->a->b逆时针 { double A2 = Sqr(r - a), B2 = Sqr(r - b), C2 = Sqr(a - b); return R * R * acos((A2 + B2 - C2) * 0.5 / sqrt(A2) / sqrt(B2)) * 0.5; } double TACIA(const Point &r, const Point &a, const Point &b, double R) //TriangleAndCircleIntersectArea，逆时针，r为圆心 { double adis = r.Dis(a), bdis = r.Dis(b); if(adis < R + eps && bdis < R + eps) return r.cross(a, b) * 0.5; Point ta, tb; if(r.InLine(a, b)) return 0.0; double rtos = LineCrossCircle(a, b, r, R, ta, tb); if(rtos > R - eps) return SectorArea(r, a, b, R); if(adis < R + eps) return r.cross(a, tb) * 0.5 + SectorArea(r, tb, b, R); if(bdis < R + eps) return r.cross(ta, b) * 0.5 + SectorArea(r, a, ta, R); return r.cross(ta, tb) * 0.5 + SectorArea(r, a, ta, R) + SectorArea(r, tb, b, R); } const int N = 505; Point p[N]; double SPICA(int n, Point r, double R)//SimplePolygonIntersectCircleArea { int i; double res = 0, if_clock_t; for(i = 0; i < n; ++ i) { if_clock_t = dcmp(r.cross(p[i], p[(i + 1) % n])); if(if_clock_t < 0) res -= TACIA(r, p[(i + 1) % n], p[i], R); else res += TACIA(r, p[i], p[(i + 1) % n], R); } return fabs(res); } double r; int n; int main() { while (~scanf("%lf", &r)) { scanf("%d", &n); for (int i = 0; i < n; i++) p[i].read(); printf("%.2f\n", SPICA(n, Point(0, 0), r)); } return 0; } |
Read More »
[模板]旋转卡壳
旋转卡壳就是对凸包用两条平行线夹着，然后绕着凸包的定点和边旋转，在旋转的过程中我们可以得到凸包的宽，直径等特征值，当然还可以计算其他东西。
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 |
#include <bits/stdc++.h> using namespace std; const int maxn=100+10; const double eps=1e-8; int dcmp(double x) { if (fabs(x)<eps) return 0; if (x>0) return 1; return -1; } struct point{ double x,y; point(){} point(double x,double y):x(x),y(y){} point operator + (const point & o) const{ return point(x+o.x,y+o.y); } point operator - (const point & o) const{ return point(x-o.x,y-o.y); } double operator ^ (const point & o) const { return x*o.y-y*o.x; } bool operator < (const point & o) const{ if (y==o.y) return x<o.x; return y>o.y; } void read() { scanf("%lf%lf",&x,&y); } double dis(const point & o) { return sqrt((x-o.x)*(x-o.x)+(y-o.y)*(y-o.y)); } } p[maxn]; double dis(point a,point b,point c) { double res= ((a-b)^(a-c))/b.dis(c); cerr<<res<<endl; return res; } bool cmp_x(const point &p,const point &q) { if (p.x!=q.x) return p.x<q.x; return p.y<q.y; } vector<point> convex_hull(point * ps,int n) { sort(ps,ps+n,cmp_x); int k=0; vector<point> qs(n*2); for (int i=0;i<n;i++) { while(k>1&&dcmp((qs[k-1]-qs[k-2])^(ps[i]-qs[k-1]))<=0) k--; qs[k++]=ps[i]; } int t=k; for (int i=n-2;i>=0;i--) { while(k>t&&dcmp((qs[k-1]-qs[k-2])^(ps[i]-qs[k-1]))<=0) k--; qs[k++]=ps[i]; } qs.resize(k-1); return qs; } int main() { int n; while(scanf("%d",&n)!=EOF) { for (int i=0;i<n;i++) { p[i].read(); } vector<point> pt=convex_hull(p,n); int l,r; int maxy,miny; maxy=-100005; miny=100005; for (int i=0;i<pt.size();i++) { if (pt[i].y>maxy) { maxy=pt[i].y; r=i; } if (pt[i].y<miny) { miny=pt[i].y; l=i; } } n=pt.size(); double ans=4*100000.0; int sl=l,sr=r; while(l!=sr||r!=sl) { int nl=(sl+1)%n; int nr=(sr+1)%n; if (dcmp((pt[nl]-pt[sl])^(pt[nr]-pt[sr]))<0) { ans=min(ans,dis(pt[sr],pt[sl],pt[nl])); sl=nl; } else { ans=min(ans,dis(pt[sl],pt[sr],pt[nr])); sr=nr; } } printf("%.10lf\n",ans); } return 0; } |
Read More »
[模板]ntt和fft
FFT：
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 |
struct Complex { double real, image; Complex( double _real, double _image) { real = _real; image = _image; } Complex(){} }; Complex operator + (const Complex &c1, const Complex &c2) { return Complex(c1.real + c2.real, c1.image + c2.image); } Complex operator - (const Complex &c1, const Complex &c2) { return Complex(c1.real - c2.real, c1.image - c2.image); } Complex operator * (const Complex &c1, const Complex &c2) { return Complex(c1.real*c2.real - c1.image*c2.image, c1.real*c2.image + c1.image*c2.real); } int rev(int id, int len) { int ret = 0; for(int i = 0; (1 << i) < len; i++) { ret <<= 1; if(id & (1 << i)) ret |= 1; } return ret; } Complex A[1 << 19]; void FFT(Complex *a, int len, int DFT) { for(int i = 0; i < len; i++) A[rev(i, len)] = a[i]; for(int s = 1; (1 << s) <= len; s++) { int m = (1 << s); Complex wm = Complex(cos(DFT*2*PI/m), sin(DFT*2*PI/m)); for(int k = 0; k < len; k += m) { Complex w = Complex(1, 0); for(int j = 0; j < (m >> 1); j++) { Complex t = w*A[k + j + (m >> 1)]; Complex u = A[k + j]; A[k + j] = u + t; A[k + j + (m >> 1)] = u - t; w = w*wm; } } } if(DFT == -1) for(int i = 0; i < len; i++) A[i].real /= len, A[i].image /= len; for(int i = 0; i < len; i++) a[i] = A[i]; return; } |
NTT
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 |
#include <iostream> #include <cstdio> #include <cstring> #include <algorithm> #include <cmath> using namespace std; const int maxn=300020; const int P=880803841; const int G=26; int wn[25]; long long mul(long long x,long long y) { return (x*y-(long long)(x/(long double)P*y+1e-3)*P+P) % P; } int qpow(int x,int k,int p) { int ret=1; while(k) { if (k&1) ret=1LL*ret*x%p; k>>=1; x=1LL*x*x%p; } return ret; } void getwn() { for (int i=1;i<=21;i++) { int t=1<<i; wn[i]=qpow(G,(P-1)/t,P); } } void change(int *y,int len) { for (int i=1,j=len/2;i<len-1;i++) { if (i<j) swap(y[i],y[j]); int k=len/2; while(j>=k) { j-=k; k/=2; } j+=k; } } void NTT(int *y,int len,int on) { change(y,len); int id=0; for (int h=2;h<=len;h<<=1) { id++; for (int j=0;j<len;j+=h) { int w=1; for (int k=j;k<j+h/2;k++) { int u=y[k]; int t=1LL*y[k+h/2]*w%P; y[k]=u+t; if (y[k]>=P) y[k]-=P; y[k+h/2]=u-t+P; if (y[k+h/2]>=P) y[k+h/2]-=P; w=1LL*w*wn[id]%P; } } } if (on==-1) { for (int i=1;i<len/2;i++) { swap(y[i],y[len-i]); } int inv=qpow(len,P-2,P); for (int i=1;i<len;i++) { y[i]=1LL*y[i]*inv%P; } } } int main() { getwn(); return 0; } |
[模板]半平面交
给定一系列半平面，求其交集的面积 下面代码p为给定点集，用于得到半平面l
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 |
#include <iostream> #include <cstdio> #include <cstring> #include <algorithm> #include <cmath> using namespace std; const int maxn=1500+10; const double eps=1e-8; int n,pn,dq[maxn],top,bot; struct point{ double x,y; point(){} point(double x,double y):x(x),y(y){} void read() { scanf("%lf%lf",&x,&y); } void print() { printf("(%lf,%lf)\n",x,y); } point operator +(const point &o)const { return point(x+o.x,y+o.y); } point operator -(const point &o)const { return point(x-o.x,y-o.y); } point operator *(double k) const { return point(x*k,y*k); } double operator *(const point &o)const { return x*o.x+y*o.y; } double operator ^(const point &o)const { return x*o.y-y*o.x; } } p[maxn]; struct line{ point a,b; double angle; line(){} line(point oa,point ob) { a=oa; b=ob; angle=atan2(ob.y-oa.y,ob.x-oa.x); } line & operator =(line l) { a.x=l.a.x; a.y=l.a.y; b.x=l.b.x; b.y=l.b.y; angle=l.angle; return *this; } } l[maxn]; int dblcmp(double k) { if (fabs(k)<eps) return 0; return k>0?1:-1; } double multi(point p0,point p1,point p2) { return (p1-p0)^(p2-p0); } bool cmp(const line &l1,const line &l2) { int d=dblcmp(l1.angle-l2.angle); if (!d) { return dblcmp(multi(l1.a,l2.a,l2.b))<0; } return d<0; } void getIntersect(line l1,line l2,point &p) { double A1=l1.b.y-l1.a.y; double B1=l1.a.x-l1.b.x; double C1=(l1.b.x-l1.a.x)*l1.a.y-(l1.b.y-l1.a.y)*l1.a.x; double A2=l2.b.y-l2.a.y; double B2=l2.a.x-l2.b.x; double C2=(l2.b.x-l2.a.x)*l2.a.y-(l2.b.y-l2.a.y)*l2.a.x; p.x=(C2*B1-C1*B2)/(A1*B2-A2*B1); p.y=(C1*A2-C2*A1)/(A1*B2-A2*B1); } bool judge(line l0,line l1,line l2) { point p; getIntersect(l1,l2,p); return dblcmp(multi(p,l0.a,l0.b))>0; } void HalfPlaneIntersect() { int i,j; sort(l,l+n,cmp); for (i=0,j=0;i<n;i++) { if (dblcmp(l[i].angle-l[j].angle)>0) { l[++j]=l[i]; } } n=j+1; dq[0]=0; dq[1]=1; top=1; bot=0; for (int i=2;i<n;i++) { while(top>bot&&judge(l[i],l[dq[top]],l[dq[top-1]])) top--; while(top>bot&&judge(l[i],l[dq[bot]],l[dq[bot+1]])) bot++; dq[++top]=i; } while(top>bot&&judge(l[dq[bot]],l[dq[top]],l[dq[top-1]])) top--; while(top>bot&&judge(l[dq[top]],l[dq[bot]],l[dq[bot+1]])) bot++; dq[++top]=dq[bot]; for (pn=0,i=bot;i<top;i++,pn++) { getIntersect(l[dq[i+1]],l[dq[i]],p[pn]); } } double getArea() { if (pn<3) return 0; double area=0; for (int i=1;i<pn-1;i++) { //p[i].print(); area+=multi(p[0],p[i],p[i+1]); } if (area<0) area=-area; return area/2; } int main() { int t,i; scanf("%d",&t); while(t--) { scanf("%d",&n); for (i=0;i<n;i++) { p[i].read(); } for (i=0;i<n-1;i++) { l[i]=line(p[i],p[i+1]); } l[i]=line(p[i],p[0]); HalfPlaneIntersect(); printf("%.2lf\n",getArea()); } return 0; } |
Read More »
[模板]最大空凸包
最大空凸包指的是以给定点集的子集为定点的凸包，且该凸包内部不包含其他点，求其中面积最大的那个。 下面代码中的dot数组表示输入的点集，n为点集大小。
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 |
#include <iostream> #include <cstdio> #include <cstring> #include <cmath> #include <algorithm> using namespace std; const int maxn=55; const double zero=1e-8; struct Vector{ double x,y; }; inline Vector operator - (Vector a, Vector b) { Vector c; c.x=a.x-b.x; c.y=a.y-b.y; return c; } inline double Sqr(double a) { return a*a; } inline int Sign(double a) { if (fabs(a)<=zero) return 0; return a<0?-1:1; } inline bool operator < (Vector a,Vector b) { return Sign(b.y-a.y)>0||Sign(b.y-a.y)==0&&Sign(b.x-a.x)>0; } inline double Max(double a,double b) { return a>b ? a : b; } inline double Length(Vector a) { return sqrt(Sqr(a.x)+Sqr(a.y)); } inline double Cross(Vector a,Vector b) { return a.x*b.y-a.y*b.x; } Vector dot[maxn],List[maxn]; double opt[maxn][maxn]; int seq[maxn]; int n,len; double ans; bool Compare(Vector a,Vector b) { int temp=Sign(Cross(a,b)); if(temp!=0) return temp>0; temp=Sign(Length(b)-Length(a)); return temp>0; } void solve(int vv) { int t,i,j,_len; for (i=len=0;i<n;i++) { if (dot[vv]<dot[i]) List[len++]=dot[i]-dot[vv]; } for (i=0;i<len;i++) { for (j=0;j<len;j++) { opt[i][j]=0; } } sort(List,List+len,Compare); double v; for (t=1;t<len;t++) { _len=0; for (i=t-1;i>=0&&Sign(Cross(List[t],List[i]))==0;i--); while(i>=0) { v=Cross(List[i],List[t])/2; seq[_len++]=i; for (j=i-1;j>=0&&Sign(Cross(List[i]-List[t],List[j]-List[t]))>0;j--); if (j>=0) v+=opt[i][j]; ans=Max(ans,v); opt[t][i]=v; i=j; } for (i=_len-2;i>=0;i--) { opt[t][seq[i]]=Max(opt[t][seq[i]],opt[t][seq[i+1]]); } } } int main() { int T; scanf("%d",&T); while(T--) { scanf("%d",&n); for (int i=0;i<n;i++) { scanf("%lf%lf",&dot[i].x,&dot[i].y); } ans=0; for (int i=0;i<n;i++) { solve(i); } printf("%.1lf\n",ans); } return 0; } |
Read More »
插值？插值！
这是一个神奇的东西，mark！
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 |
#include <cstdio> #include <cstring> #include <cmath> #include <algorithm> #include <vector> #include <string> #include <map> #include <set> #include <cassert> using namespace std; #define rep(i,a,n) for (int i=a;i<n;i++) #define per(i,a,n) for (int i=n-1;i>=a;i--) #define pb push_back #define mp make_pair #define all(x) (x).begin(),(x).end() #define fi first #define se second #define SZ(x) ((int)(x).size()) typedef vector<int> VI; typedef long long ll; typedef pair<int,int> PII; const ll mod=1000000007; ll powmod(ll a,ll b) {ll res=1;a%=mod; assert(b>=0); for(;b;b>>=1){if(b&1)res=res*a%mod;a=a*a%mod;}return res;} // head int _; ll n; namespace linear_seq { const int N=10010; ll res[N],base[N],_c[N],_md[N]; vector<ll> Md; void mul(ll *a,ll *b,ll k) { rep(i,0,k+k) _c[i]=0; rep(i,0,k) if (a[i]) rep(j,0,k) _c[i+j]=(_c[i+j]+a[i]*b[j])%mod; for (int i=k+k-1;i>=k;i--) if (_c[i]) rep(j,0,SZ(Md)) _c[i-k+Md[j]]=(_c[i-k+Md[j]]-_c[i]*_md[Md[j]])%mod; rep(i,0,k) a[i]=_c[i]; } int solve(ll n,VI a,VI b) { ll ans=0,pnt=0; ll k=SZ(a); assert(SZ(a)==SZ(b)); rep(i,0,k) _md[k-1-i]=-a[i];_md[k]=1; Md.clear(); rep(i,0,k) if (_md[i]!=0) Md.push_back(i); rep(i,0,k) res[i]=base[i]=0; res[0]=1; while ((1ll<<pnt)<=n) pnt++; for (int p=pnt;p>=0;p--) { mul(res,res,k); if ((n>>p)&1) { for (int i=k-1;i>=0;i--) res[i+1]=res[i];res[0]=0; rep(j,0,SZ(Md)) res[Md[j]]=(res[Md[j]]-res[k]*_md[Md[j]])%mod; } } rep(i,0,k) ans=(ans+res[i]*b[i])%mod; if (ans<0) ans+=mod; return ans; } VI BM(VI s) { VI C(1,1),B(1,1); int L=0,m=1,b=1; rep(n,0,SZ(s)) { ll d=0; rep(i,0,L+1) d=(d+(ll)C[i]*s[n-i])%mod; if (d==0) ++m; else if (2*L<=n) { VI T=C; ll c=mod-d*powmod(b,mod-2)%mod; while (SZ(C)<SZ(B)+m) C.pb(0); rep(i,0,SZ(B)) C[i+m]=(C[i+m]+c*B[i])%mod; L=n+1-L; B=T; b=d; m=1; } else { ll c=mod-d*powmod(b,mod-2)%mod; while (SZ(C)<SZ(B)+m) C.pb(0); rep(i,0,SZ(B)) C[i+m]=(C[i+m]+c*B[i])%mod; ++m; } } return C; } int gao(VI a,ll n) { VI c=BM(a); c.erase(c.begin()); rep(i,0,SZ(c)) c[i]=(mod-c[i])%mod; return solve(n,c,VI(a.begin(),a.begin()+SZ(c))); } }; int main() { while(~scanf("%lld",&n)){ printf("%d\n",linear_seq::gao(VI{1,1, 5, 11, 36, 95, 281, 781, 2245, 6336, 18061, 51205, 145601, 413351, 1174500, 3335651, 9475901, 26915305, 76455961, 217172736},n)); } } |
Read More »
[HDU-3487]Play with Chain
F – Play with Chain Time Limit:2000MS Memory Limit:32768KB 64bit IO Format:%I64d & %I64u Submit Status Description YaoYao is fon…
Read More »