四维的01背包，注意一下所有代价都是0时的方案输出即可
#include<queue> #include<cstdio> #include<cstdlib> #include<cstring> #include<algorithm> #define maxn 37 using namespace std; int p[maxn],a[maxn],c[maxn],m[maxn],g[maxn]; int dp[maxn][maxn][maxn][maxn],use[maxn]; int n,P,A,C,M; vector <int> mem; int main(){ scanf("%d",&n); int cnt=0; for(int i=1;i<=n;i++) { scanf("%d%d%d%d%d",&p[i],&a[i],&c[i],&m[i],&g[i]); if(p[i]==0&&a[i]==0&&c[i]==0&&m[i]==0) { mem.push_back(i); use[i]=1; } } scanf("%d%d%d%d",&P,&A,&C,&M); dp[0][0][0][0]=0; for(int i=1;i<=n;i++){ if(use[i]==1) continue; for(int j=P;j>=0;j--) for(int k=A;k>=0;k--) for(int t=C;t>=0;t--) for(int l=M;l>=0;l--) if(j>=p[i] && k>=a[i] && t>=c[i] && l>=m[i]){ if(dp[j][k][t][l]<dp[j-p[i]][k-a[i]][t-c[i]][l-m[i]]+g[i]){ dp[j][k][t][l]=dp[j-p[i]][k-a[i]][t-c[i]][l-m[i]]+g[i]; } } } int ans=0,l1=P,l2=A,l3=C,l4=M; for(int j=P;j>=0;j--) for(int k=A;k>=0;k--) for(int t=C;t>=0;t--) for(int l=M;l>=0;l--) if(dp[j][k][t][l]>=ans){ l1=j; l2=k; l3=t; l4=l; ans=dp[j][k][t][l]; } while(1){ int t=233; for(int i=1;i<=n;i++) if(!(p[i]==0 && a[i]==0 && c[i]==0 && m[i]==0) && !use[i] && l1>=p[i] && l2>=a[i] && l3>=c[i] && l4>=m[i] && dp[l1][l2][l3][l4]==dp[l1-p[i]][l2-a[i]][l3-c[i]][l4-m[i]]+g[i]){ use[i]=1; t=i; break; } if(t==233) break; mem.push_back(t); l1-=p[t]; l2-=a[t]; l3-=c[t]; l4-=m[t]; } printf("%d\n",mem.size()); for(int i=0;i<mem.size();i++) printf("%d ",mem[i]-1); return 0; }
Splay基本操作，可以把目标区间[l,r]删掉，再提到开头，也可以在[l,r],[1,r-l+1],[r-l+2,r]打三次翻转标记
1 #include<cstdio> 2 #include<cstdlib> 3 #include<cstring> 4 #include<algorithm> 5 6 #define maxn 100000+5 7 8 using namespace std; 9 10 struct Splay_Tree{ 11 int ch[2]; 12 int sz,key; 13 }tr[maxn<<1]; 14 15 int a[maxn]; 16 int n,m,cnt,root; 17 18 inline int in(){ 19 int x=0,f=1; 20 char ch=getchar(); 21 while(ch<'0' || ch>'9'){ if(ch=='-') f=-1; ch=getchar(); } 22 while(ch>='0' && ch<='9') x=x*10+ch-'0',ch=getchar(); 23 return x*f; 24 } 25 26 inline void Update(int k){ 27 tr[k].sz=tr[tr[k].ch[0]].sz+tr[tr[k].ch[1]].sz+1; 28 } 29 30 inline void Rotate(int &k,int d){ 31 int t=tr[k].ch[d^1]; 32 tr[k].ch[d^1]=tr[t].ch[d]; tr[t].ch[d]=k; 33 Update(k); Update(t); k=t; 34 } 35 36 void Splay(int &k,int x){ 37 int d1=(tr[tr[k].ch[0]].sz<x?1:0),t=tr[k].ch[d1]; 38 if(d1==1) x-=tr[tr[k].ch[0]].sz+1; 39 if(x){ 40 int d2=(tr[tr[t].ch[0]].sz<x?1:0); 41 if(d2==1) x-=tr[tr[t].ch[0]].sz+1; 42 if(x){ 43 Splay(tr[t].ch[d2],x); 44 if(d1==d2) Rotate(k,d1^1); 45 else Rotate(tr[k].ch[d1],d1); 46 } 47 Rotate(k,d1^1); 48 } 49 Update(k); 50 } 51 52 void dfs(int x){ 53 if(tr[x].ch[0]) dfs(tr[x].ch[0]); 54 if(tr[x].key>=1 && tr[x].key<=n) printf("%d ",tr[x].key); 55 if(tr[x].ch[1]) dfs(tr[x].ch[1]); 56 } 57 58 void Insert(int x,char val){ 59 Splay(root,x+1); Splay(tr[root].ch[1],x+1-tr[tr[root].ch[0]].sz); 60 //本应为Splay(tr[root].r,x+2-(tr[tr[root].ch[0]].sz+1)); 61 int k=++cnt; 62 tr[k].key=val; tr[k].ch[0]=tr[k].ch[1]=0; 63 tr[tr[root].ch[1]].ch[0]=k; 64 Update(k); Update(tr[root].ch[1]); Update(root); 65 } 66 67 void Delete(int l,int r){ 68 Splay(root,l); 69 Splay(tr[root].ch[1],1+r-tr[tr[root].ch[0]].sz); 70 int k=tr[tr[root].ch[1]].ch[0]; 71 tr[tr[root].ch[1]].ch[0]=0; 72 Update(tr[root].ch[1]); Update(root); 73 Splay(root,1); Splay(tr[root].ch[1],1-tr[tr[root].ch[0]].sz); 74 tr[tr[root].ch[1]].ch[0]=k; 75 Update(k); Update(tr[root].ch[1]); Update(root); 76 77 78 } 79 80 void Build(int l,int r,int &k){ 81 if(l>r) return; 82 if(l==r){ 83 k=++cnt; tr[k].ch[0]=tr[k].ch[1]=0; 84 tr[k].key=a[l]; 85 tr[k].sz=1; 86 return; 87 } 88 int mid=(l+r)>>1; k=++cnt; 89 Build(l,mid-1,tr[k].ch[0]); 90 Build(mid+1,r,tr[k].ch[1]); 91 tr[k].key=a[mid]; 92 Update(k); 93 } 94 95 int main(){ 96 n=in(); m=in(); 97 for(int i=1;i<=n;i++) a[i]=i; 98 Build(0,n+1,root); 99 // dfs(root); puts(""); 100 for(int i=1;i<=m;i++){ 101 int l,r,x; 102 l=in(); x=in(); 103 r=l+x-1; 104 Delete(l,r); 105 } 106 dfs(root); 107 return 0; 108 }
先把串拼成二倍长，预处理整串hash值。然后遍历取n个字符的hash值，map一下看看有没有重复并且存到vector中，最后直接输出
然鹅被卡常了还是getchar快一点
#include <cstdio> #include <cstdlib> #include <iostream> #include <cstring> #include <algorithm> #include <cmath> #include <queue> #include <vector> #include <string> #include <map> #include <unordered_map> #define LSON l,m,x<<1 #define RSON m+1,r,x<<1|1 using namespace std; const int MAX=2e6+5; const unsigned long long base = 163; char s[MAX]; unsigned long long p[MAX],hh[MAX]; unordered_map<unsigned long long,int> mp; int cnt; vector<int> L[MAX]; void init(){//处理hh值 p[0] = 1; hh[0] = 0; int n = strlen(s + 1); for(int i = 1; i < MAX/2; i ++) p[i] =p[i-1] * base; for(int i = 1; i <= n; i ++) hh[i] = hh[i - 1] * base + (s[i] - 'a'); } inline unsigned long long get(int l, int r){//取出g里l - r里面的字符串的hh值 return hh[r] - hh[l - 1] * p[r - l + 1]; } inline string read()//inline继续加快速度 { char ch=getchar(); string st1=""; while (!((ch>='a')&&(ch<='z')))//把前面没用的东西去掉,当然,ch在什么范围内可以依据需要修改 ch=getchar(); while ((ch>='a')&&(ch<='z')) st1+=ch,ch=getchar(); return st1;//返回 }//在主程序内可以写st=read(),这样子要读的字符串就到了st内 int main(){ int i,j; char ch; i=0; while((ch=getchar())!=EOF){ s[++i]=ch; } if(s[i]=='\n') i--; int n=i; for(i=n+1;i<=2*n;i++) s[i]=s[i-n]; init(); for(i=1;i<=n;i++){ unsigned long long h=get(i,i+n-1); if(!mp[h]){ mp[h]=++cnt; L[cnt].push_back(i-1); } else L[mp[h]].push_back(i-1); } printf("%d\n",cnt); for(i=1;i<=cnt;i++){ int size=L[i].size(); printf("%d ",size); for(j=0;j<size;j++){ printf("%d%c",L[i][j],j==size-1?'\n':' '); } } return 0; }
求十六进制串的SOD，发现有规律只和数字的和有关，需要统计区间内所有子串的SOD，遂用线段树维护每个区间内SOD=0-15方案的个数，向上更新的时候卷积就好。
#include <cstdio> #include <cstdlib> #include <iostream> #include <cstring> #include <algorithm> #include <cmath> #include <queue> #include <vector> #include <string> #include <map> #include <ctime> #define LSON l,m,x<<1 #define RSON m+1,r,x<<1|1 using namespace std; const int MAX=1e5+5; const int D=16; const int mod=1e9+7; struct node{ long long v[D]; }xds[MAX*4]; int n,q; char s[MAX]; int R[D][D]; long long num[D]; node operator + (const node &a,const node &b){ node ret; for(int i=0;i<D;i++) ret.v[i]=0; for(int i=0;i<D;i++){ for(int j=0;j<D;j++){ if(a.v[i]&&b.v[j]){ int val=(a.v[i]*b.v[j])%mod; int to=R[i][j]; ret.v[to]=(ret.v[to]+val)%mod; } } } return ret; } int trans(char c){ if(isdigit(c)) return c-'0'; return c-'A'+10; } void pushup(int x){ xds[x]=xds[x<<1]+xds[x<<1|1]; } void js(int l,int r,int x) { int m=(l+r)/2; if(l==r){ xds[x].v[0]++; xds[x].v[trans(s[l])]++; return; } js(LSON); js(RSON); pushup(x); } void xg(int pos,int val,int l,int r,int x){ int m=(l+r)/2; if(l==r&&pos==l){ for(int i=0;i<D;i++) xds[x].v[i]=0; xds[x].v[0]++; xds[x].v[val]++; return; } if(pos<=m) xg(pos,val,LSON); else xg(pos,val,RSON); pushup(x); } node cx(int L,int R,int l,int r,int x){ int m = (l + r) / 2; if(L==l&&R==r) return xds[x]; if(R<=m) return cx(L,R,LSON); if(L>m) return cx(L,R,RSON); return cx(L,m,LSON)+cx(m+1,R,RSON); } long long query(int l,int r){ long long ans=0; node tmp=cx(l,r,1,n,1); for(int i=0;i<D;i++) ans=(ans+(tmp.v[i]*num[i])%mod)%mod; return ans-1; } void init(void){ for(int i=0;i<D;i++) for(int j=0;j<D;j++){ int to=(i+j)%(D-1); if((i+j)!=0&&to==0) to=15; R[i][j]=to; } num[0]=1; for(int i=1;i<D;i++) num[i]=(num[i-1]*1021)%mod; } int main(){ int i,j; init(); scanf("%d%d",&n,&q); scanf("%s",s+1); js(1,n,1); while(q--){ int op,l,r,pos; char c; scanf("%d",&op); if(op==1){ scanf("%d %c",&pos,&c); xg(pos,trans(c),1,n,1); } else{ scanf("%d%d",&l,&r); printf("%lld\n",query(l,r)); } } return 0; }
题目要求两个相同的染色最近距离为恰好为D，设函数f(D)为对某点距离D以内没有相同的颜色的点的方案数，即两个相同的染色最近距离大于等于D+1，则答案为f(D-1)-f(D)。
故按照bfs顺序计数，发现bfs过的节点里和当前节点距离小于D的点两两之间距离一定小于D。故每处理一个点x，在已经bfs过的点中搜索距离这个点D以内的点的个数，计为cnt，则x可染色的方案数为k-cnt，然后利用乘法原理计数。
#include<queue> #include<cstdio> #include<cstdlib> #include<cstring> #include<algorithm> #define maxn 5000+5 using namespace std; typedef long long LL; const int mod=1e9+7; struct Edge { int u,v,nxt; Edge() { } Edge(int a,int b,int c) { u=a,v=b,nxt=c; } } e[maxn<<1]; int head[maxn],p[maxn],out[maxn]; int ind,n,k,d; void addedge(int x,int y) { e[ind]=Edge(x,y,head[x]),head[x]=ind++; e[ind]=Edge(y,x,head[y]),head[y]=ind++; } int qpow(int x,int b) { int base=x,res=1; while(b) { if(b&1) res=(LL)res*base%mod; base=(LL)base*base%mod; b>>=1; } return res; } int bfs(int x,int pa,int dep) { if(dep==0) return 1; int res=1; for(int i=head[x];i!=-1;i=e[i].nxt) if(e[i].v!=pa && out[e[i].v]){ res+=bfs(e[i].v,x,dep-1); } return res; } int f(int dis) { if(dis==0) return qpow(k,n); memset(p,0,sizeof(p)); memset(out,0,sizeof(out)); queue <int> q; q.push(1); p[1]=1; int res=1; while(!q.empty()) { int x=q.front(); q.pop(); out[x]=1; int cnt=bfs(x,0,dis)-1; if(cnt>=k) return 0; res=(LL)res*(k-cnt)%mod; for(int i=head[x];i!=-1;i=e[i].nxt) if(!p[e[i].v]){ p[e[i].v]=1; q.push(e[i].v); } } return res; } int main() { memset(head,-1,sizeof(head)); scanf("%d%d%d",&n,&k,&d); for(int i=1; i<n; i++) { int x,y; scanf("%d%d",&x,&y); addedge(x,y); } printf("%d",((f(d-1)-f(d))%mod+mod)%mod); return 0; }
随便推推公式就好，发现和小于n/num的质数数量有关，求一求质数数量的前缀和
#include <cstdio> #include <cstdlib> #include <iostream> #include <cstring> #include <algorithm> #include <cmath> #include <queue> #include <vector> #include <string> #include <map> using namespace std; int pri[6000010], cnt; bool flag[10000010]; long long sum[10000010]; void init() { cnt=0; for(int i = 2; i <= 10000000; i ++) { if(!flag[i]) { pri[cnt ++] = i; flag[i] = 1; } for(int j = 0; j < cnt && i * pri[j] <= 10000000; j ++) { flag[i * pri[j]] = 1; if(i % pri[j] == 0) { break; } } } } int main() { int i; init(); for(i=0;i<cnt;i++) sum[pri[i]]++; for(i=0;i<=10000000;i++) sum[i]+=sum[i-1]; int n; scanf("%d",&n); long long ans=0; for(i=1;i<=n;i++) { ans+=sum[n/i]*(sum[n/i]-1); } printf("%lld\n",ans); return 0; }
I：Expected Size of Random Convex Hull
本地用随机数模拟打表
#include <cstdio> #include <cstdlib> #include <iostream> #include <cstring> #include <algorithm> #include <cmath> #include <queue> #include <vector> #include <string> #include <map> #include <ctime> #define LSON l,m,x<<1 #define RSON m+1,r,x<<1|1 #define RAND_MAX 80000 using namespace std; const int MAX=30000+5; //1.1 Point定义 const double eps=1e-8; const double PI=acos(-1.0); int sgn(double x) { if(fabs(x)<eps)return 0; if(x < 0)return -1; else return 1; } struct Point { double x,y; Point(){} //默认构造器（），如果没有自定义构造器，则编译器会自动生成默认构造器 Point(double _x,double _y) //自定义构造器 { x=_x;y=_y; } /* Point(double _x,double _y) : x(_x), y(_y) //自定义构造器版本2，效果同上，但效率更高 { } */ Point operator - (const Point &b)const { return Point(x-b.x,y-b.y); } double operator ^ (const Point &b)const//叉积 { return x*b.y-y*b.x; } double operator * (const Point &b)const//点积 { return x*b.x+y*b.y; } bool operator < (const Point &b)const { if(x==b.x) return y<b.y; return x<b.x; } void transXY(double B)//绕原点旋转角度B（弧度值），后x,y的变化 { double tx=x,ty=y; x=tx*cos(B)-ty*sin(B); y=tx*sin(B)+ty*cos(B); } }; //输入点数组p，个数为n，输出点数组ch。返回凸包顶点数 //输入不能有重复点。函数执行完之后输入点的顺序被破坏 //如果不希望在凸包的边上有输入点，把两个<=改为< //在精度要求高时建议用cmp比较 int ConvexHull(Point* p,int n,Point *ch) { sort(p,p+n); //先比较x坐标，再比较y坐标 int m=0; for(int i=0;i<n;i++) { while(m > 1 && ((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 && ((ch[m-1]-ch[m-2])^(p[i]-ch[m-2])) <= 0) m--; ch[m++] = p[i]; } if(n > 1) m--; return m; } int m=1e7,n; Point p[MAX],ch[MAX]; double ans; long long sum; int Rand(int l,int r){ return (((long long)rand())*rand()%(r-l+1)+l); } int main(){ int i,j; srand(time(0)); printf("%d\n",RAND_MAX); for(n=4;n<=10;n++){ sum=0; for(j=0;j<m;j++){ //cout<<j<<endl; for(i=0;i<n;i++){ p[i].x=rand(); p[i].y=rand(); if(p[i].x>p[i].y) swap(p[i].x,p[i].y); } sum+=ConvexHull(p,n,ch); } printf("%lf,",(double)sum/m); } return 0; } //double ans[100]={0,0,0,3.000000,3.000000,4.000000,4.000000,5.000000,6.000000,6.000000,6.000000}; //int main(){ // int n; // int a,b,c,x,y,z; // scanf("%d%d%d%d%d%d%d",&a,&b,&c,&x,&y,&z,&n); // printf("%lf\n",ans[n]); // // return 0; //}
二分圆的大小，用多边形与圆的面积交的模板进行检查
1 #include <cstdio> 2 #include <cstring> 3 #include <cmath> 4 #include <vector> 5 #include <algorithm> 6 7 using namespace std; 8 const double pi = 4 * atan(1); 9 const double eps = 1e-10; 10 11 inline int dcmp (double x) { if (fabs(x) < eps) return 0; else return x < 0 ? -1 : 1; } 12 inline double getDistance (double x, double y) { return sqrt(x * x + y * y); } 13 14 struct Point { 15 double x, y; 16 Point (double x = 0, double y = 0): x(x), y(y) {} 17 void read () { scanf("%lf%lf", &x, &y); } 18 void write () { printf("%lf %lf", x, y); } 19 20 bool operator == (const Point& u) const { return dcmp(x - u.x) == 0 && dcmp(y - u.y) == 0; } 21 bool operator != (const Point& u) const { return !(*this == u); } 22 bool operator < (const Point& u) const { return x < u.x || (x == u.x && y < u.y); } 23 bool operator > (const Point& u) const { return u < *this; } 24 bool operator <= (const Point& u) const { return *this < u || *this == u; } 25 bool operator >= (const Point& u) const { return *this > u || *this == u; } 26 Point operator + (const Point& u) { return Point(x + u.x, y + u.y); } 27 Point operator - (const Point& u) { return Point(x - u.x, y - u.y); } 28 Point operator * (const double u) { return Point(x * u, y * u); } 29 Point operator / (const double u) { return Point(x / u, y / u); } 30 double operator * (const Point& u) { return x*u.y - y*u.x; } 31 }; 32 33 typedef Point Vector; 34 35 struct Line { 36 double a, b, c; 37 Line (double a = 0, double b = 0, double c = 0): a(a), b(b), c(c) {} 38 }; 39 40 struct Circle { 41 Point o; 42 double r; 43 Circle () {} 44 Circle (Point o, double r = 0): o(o), r(r) {} 45 void read () { o.read(), scanf("%lf", &r); } 46 Point point(double rad) { return Point(o.x + cos(rad)*r, o.y + sin(rad)*r); } 47 double getArea (double rad) { return rad * r * r / 2; } 48 }; 49 50 51 namespace Punctual { 52 double getDistance (Point a, Point b) { double x=a.x-b.x, y=a.y-b.y; return sqrt(x*x + y*y); } 53 }; 54 55 namespace Vectorial { 56 /* 点积: 两向量长度的乘积再乘上它们夹角的余弦, 夹角大于90度时点积为负 */ 57 double getDot (Vector a, Vector b) { return a.x * b.x + a.y * b.y; } 58 59 /* 叉积: 叉积等于两向量组成的三角形有向面积的两倍, cross(v, w) = -cross(w, v) */ 60 double getCross (Vector a, Vector b) { return a.x * b.y - a.y * b.x; } 61 62 double getLength (Vector a) { return sqrt(getDot(a, a)); } 63 double getPLength (Vector a) { return getDot(a, a); } 64 double getAngle (Vector u) { return atan2(u.y, u.x); } 65 double getAngle (Vector a, Vector b) { return acos(getDot(a, b) / getLength(a) / getLength(b)); } 66 Vector rotate (Vector a, double rad) { return Vector(a.x*cos(rad)-a.y*sin(rad), a.x*sin(rad)+a.y*cos(rad)); } 67 /* 单位法线 */ 68 Vector getNormal (Vector a) { double l = getLength(a); return Vector(-a.y/l, a.x/l); } 69 }; 70 71 namespace Linear { 72 using namespace Vectorial; 73 74 Line getLine (double x1, double y1, double x2, double y2) { return Line(y2-y1, x1-x2, y1*(x2-x1)-x1*(y2-y1)); } 75 Line getLine (double a, double b, Point u) { return Line(a, -b, u.y * b - u.x * a); } 76 77 bool getIntersection (Line p, Line q, Point& o) { 78 if (fabs(p.a * q.b - q.a * p.b) < eps) 79 return false; 80 o.x = (q.c * p.b - p.c * q.b) / (p.a * q.b - q.a * p.b); 81 o.y = (q.c * p.a - p.c * q.a) / (p.b * q.a - q.b * p.a); 82 return true; 83 } 84 85 /* 直线pv和直线qw的交点 */ 86 bool getIntersection (Point p, Vector v, Point q, Vector w, Point& o) { 87 if (dcmp(getCross(v, w)) == 0) return false; 88 Vector u = p - q; 89 double k = getCross(w, u) / getCross(v, w); 90 o = p + v * k; 91 return true; 92 } 93 94 /* 点p到直线ab的距离 */ 95 double getDistanceToLine (Point p, Point a, Point b) { return fabs(getCross(b-a, p-a) / getLength(b-a)); } 96 double getDistanceToSegment (Point p, Point a, Point b) { 97 if (a == b) return getLength(p-a); 98 Vector v1 = b - a, v2 = p - a, v3 = p - b; 99 if (dcmp(getDot(v1, v2)) < 0) return getLength(v2); 100 else if (dcmp(getDot(v1, v3)) > 0) return getLength(v3); 101 else return fabs(getCross(v1, v2) / getLength(v1)); 102 } 103 104 /* 点p在直线ab上的投影 */ 105 Point getPointToLine (Point p, Point a, Point b) { Vector v = b-a; return a+v*(getDot(v, p-a) / getDot(v,v)); } 106 107 /* 判断线段是否存在交点 */ 108 bool haveIntersection (Point a1, Point a2, Point b1, Point b2) { 109 double c1=getCross(a2-a1, b1-a1), c2=getCross(a2-a1, b2-a1), c3=getCross(b2-b1, a1-b1), c4=getCross(b2-b1,a2-b1); 110 return dcmp(c1)*dcmp(c2) < 0 && dcmp(c3)*dcmp(c4) < 0; 111 } 112 113 /* 判断点是否在线段上 */ 114 bool onSegment (Point p, Point a, Point b) { return dcmp(getCross(a-p, b-p)) == 0 && dcmp(getDot(a-p, b-p)) < 0; } 115 } 116 117 namespace Triangular { 118 using namespace Vectorial; 119 120 double getAngle (double a, double b, double c) { return acos((a*a+b*b-c*c) / (2*a*b)); } 121 double getArea (double a, double b, double c) { double s =(a+b+c)/2; return sqrt(s*(s-a)*(s-b)*(s-c)); } 122 double getArea (double a, double h) { return a * h / 2; } 123 double getArea (Point a, Point b, Point c) { return fabs(getCross(b - a, c - a)) / 2; } 124 double getDirArea (Point a, Point b, Point c) { return getCross(b - a, c - a) / 2; } 125 }; 126 127 namespace Polygonal { 128 using namespace Vectorial; 129 using namespace Linear; 130 131 double getArea (Point* p, int n) { 132 double ret = 0; 133 for (int i = 1; i < n-1; i++) 134 ret += getCross(p[i]-p[0], p[i+1]-p[0]); 135 return fabs(ret)/2; 136 } 137 138 /* 凸包 */ 139 int getConvexHull (Point* p, int n, Point* ch) { 140 sort(p, p + n); 141 int m = 0; 142 for (int i = 0; i < n; i++) { 143 /* 可共线 */ 144 //while (m > 1 && dcmp(getCross(ch[m-1]-ch[m-2], p[i]-ch[m-1])) < 0) m--; 145 while (m > 1 && dcmp(getCross(ch[m-1]-ch[m-2], p[i]-ch[m-1])) <= 0) m--; 146 ch[m++] = p[i]; 147 } 148 int k = m; 149 for (int i = n-2; i >= 0; i--) { 150 /* 可共线 */ 151 //while (m > k && dcmp(getCross(ch[m-1]-ch[m-2], p[i]-ch[m-2])) < 0) m--; 152 while (m > k && dcmp(getCross(ch[m-1]-ch[m-2], p[i]-ch[m-2])) <= 0) m--; 153 ch[m++] = p[i]; 154 } 155 if (n > 1) m--; 156 return m; 157 } 158 159 int isPointInPolygon (Point o, Point* p, int n) { 160 int wn = 0; 161 for (int i = 0; i < n; i++) { 162 int j = (i + 1) % n; 163 if (onSegment(o, p[i], p[j])) return 0; // 边界上 164 int k = dcmp(getCross(p[j] - p[i], o-p[i])); 165 int d1 = dcmp(p[i].y - o.y); 166 int d2 = dcmp(p[j].y - o.y); 167 if (k > 0 && d1 <= 0 && d2 > 0) wn++; 168 if (k < 0 && d2 <= 0 && d1 > 0) wn--; 169 } 170 return wn ? -1 : 1; 171 } 172 }; 173 174 namespace Circular { 175 using namespace Linear; 176 using namespace Vectorial; 177 using namespace Triangular; 178 179 /* 直线和原的交点 */ 180 int getLineCircleIntersection (Point p, Point q, Circle O, double& t1, double& t2, vector<Point>& sol) { 181 Vector v = q - p; 182 /* 使用前需清空sol */ 183 //sol.clear(); 184 double a = v.x, b = p.x - O.o.x, c = v.y, d = p.y - O.o.y; 185 double e = a*a+c*c, f = 2*(a*b+c*d), g = b*b+d*d-O.r*O.r; 186 double delta = f*f - 4*e*g; 187 if (dcmp(delta) < 0) return 0; 188 if (dcmp(delta) == 0) { 189 t1 = t2 = -f / (2 * e); 190 sol.push_back(p + v * t1); 191 return 1; 192 } 193 194 t1 = (-f - sqrt(delta)) / (2 * e); sol.push_back(p + v * t1); 195 t2 = (-f + sqrt(delta)) / (2 * e); sol.push_back(p + v * t2); 196 return 2; 197 } 198 199 /* 圆和圆的交点 */ 200 int getCircleCircleIntersection (Circle o1, Circle o2, vector<Point>& sol) { 201 double d = getLength(o1.o - o2.o); 202 203 if (dcmp(d) == 0) { 204 if (dcmp(o1.r - o2.r) == 0) return -1; 205 return 0; 206 } 207 208 if (dcmp(o1.r + o2.r - d) < 0) return 0; 209 if (dcmp(fabs(o1.r-o2.r) - d) > 0) return 0; 210 211 double a = getAngle(o2.o - o1.o); 212 double da = acos((o1.r*o1.r + d*d - o2.r*o2.r) / (2*o1.r*d)); 213 214 Point p1 = o1.point(a-da), p2 = o1.point(a+da); 215 216 sol.push_back(p1); 217 if (p1 == p2) return 1; 218 sol.push_back(p2); 219 return 2; 220 } 221 222 /* 过定点作圆的切线 */ 223 int getTangents (Point p, Circle o, Vector* v) { 224 Vector u = o.o - p; 225 double d = getLength(u); 226 if (d < o.r) return 0; 227 else if (dcmp(d - o.r) == 0) { 228 v[0] = rotate(u, pi / 2); 229 return 1; 230 } else { 231 double ang = asin(o.r / d); 232 v[0] = rotate(u, -ang); 233 v[1] = rotate(u, ang); 234 return 2; 235 } 236 } 237 238 /* a[i] 和 b[i] 分别是第i条切线在O1和O2上的切点 */ 239 int getTangents (Circle o1, Circle o2, Point* a, Point* b) { 240 int cnt = 0; 241 if (o1.r < o2.r) { swap(o1, o2); swap(a, b); } 242 double d2 = getLength(o1.o - o2.o); d2 = d2 * d2; 243 double rdif = o1.r - o2.r, rsum = o1.r + o2.r; 244 if (d2 < rdif * rdif) return 0; 245 if (dcmp(d2) == 0 && dcmp(o1.r - o2.r) == 0) return -1; 246 247 double base = getAngle(o2.o - o1.o); 248 if (dcmp(d2 - rdif * rdif) == 0) { 249 a[cnt] = o1.point(base); b[cnt] = o2.point(base); cnt++; 250 return cnt; 251 } 252 253 double ang = acos( (o1.r - o2.r) / sqrt(d2) ); 254 a[cnt] = o1.point(base+ang); b[cnt] = o2.point(base+ang); cnt++; 255 a[cnt] = o1.point(base-ang); b[cnt] = o2.point(base-ang); cnt++; 256 257 if (dcmp(d2 - rsum * rsum) == 0) { 258 a[cnt] = o1.point(base); b[cnt] = o2.point(base); cnt++; 259 } else if (d2 > rsum * rsum) { 260 double ang = acos( (o1.r + o2.r) / sqrt(d2) ); 261 a[cnt] = o1.point(base+ang); b[cnt] = o2.point(base+ang); cnt++; 262 a[cnt] = o1.point(base-ang); b[cnt] = o2.point(base-ang); cnt++; 263 } 264 return cnt; 265 } 266 267 /* 三点确定外切圆 */ 268 Circle CircumscribedCircle(Point p1, Point p2, Point p3) { 269 double Bx = p2.x - p1.x, By = p2.y - p1.y; 270 double Cx = p3.x - p1.x, Cy = p3.y - p1.y; 271 double D = 2 * (Bx * Cy - By * Cx); 272 double cx = (Cy * (Bx * Bx + By * By) - By * (Cx * Cx + Cy * Cy)) / D + p1.x; 273 double cy = (Bx * (Cx * Cx + Cy * Cy) - Cx * (Bx * Bx + By * By)) / D + p1.y; 274 Point p = Point(cx, cy); 275 return Circle(p, getLength(p1 - p)); 276 } 277 278 /* 三点确定内切圆 */ 279 Circle InscribedCircle(Point p1, Point p2, Point p3) { 280 double a = getLength(p2 - p3); 281 double b = getLength(p3 - p1); 282 double c = getLength(p1 - p2); 283 Point p = (p1 * a + p2 * b + p3 * c) / (a + b + c); 284 return Circle(p, getDistanceToLine(p, p1, p2)); 285 } 286 287 /* 三角形一顶点为圆心 */ 288 double getPublicAreaToTriangle (Circle O, Point a, Point b) { 289 if (dcmp((a-O.o)*(b-O.o)) == 0) return 0; 290 int sig = 1; 291 double da = getPLength(O.o-a), db = getPLength(O.o-b); 292 if (dcmp(da-db) > 0) { 293 swap(da, db); 294 swap(a, b); 295 sig = -1; 296 } 297 298 double t1, t2; 299 vector<Point> sol; 300 int n = getLineCircleIntersection(a, b, O, t1, t2, sol); 301 302 if (dcmp(da-O.r*O.r) <= 0) { 303 if (dcmp(db-O.r*O.r) <= 0) return getDirArea(O.o, a, b) * sig; 304 305 int k = 0; 306 if (getPLength(sol[0]-b) > getPLength(sol[1]-b)) k = 1; 307 308 double ret = getArea(O.o, a, sol[k]) + O.getArea(getAngle(sol[k]-O.o, b-O.o)); 309 double tmp = (a-O.o)*(b-O.o); 310 return ret * sig * dcmp(tmp); 311 } 312 313 double d = getDistanceToSegment(O.o, a, b); 314 if (dcmp(d-O.r) >= 0) { 315 double ret = O.getArea(getAngle(a-O.o, b-O.o)); 316 double tmp = (a-O.o)*(b-O.o); 317 return ret * sig * dcmp(tmp); 318 } 319 320 321 double k1 = O.r / getLength(a - O.o), k2 = O.r / getLength(b - O.o); 322 Point p = O.o + (a - O.o) * k1, q = O.o + (b - O.o) * k2; 323 double ret1 = O.getArea(getAngle(p-O.o, q-O.o)); 324 double ret2 = O.getArea(getAngle(sol[0]-O.o, sol[1]-O.o)) - getArea(O.o, sol[0], sol[1]); 325 double ret = (ret1 - ret2), tmp = (a-O.o)*(b-O.o); 326 return ret * sig * dcmp(tmp); 327 } 328 329 double getPublicAreaToPolygon (Circle O, Point* p, int n) { 330 if (dcmp(O.r) == 0) return 0; 331 double area = 0; 332 for (int i = 0; i < n; i++) { 333 int u = (i + 1) % n; 334 area += getPublicAreaToTriangle(O, p[i], p[u]); 335 } 336 return fabs(area); 337 } 338 }; 339 340 using namespace Vectorial; 341 using namespace Circular; 342 using namespace Polygonal; 343 344 const int maxn = 5005; 345 const double inf = 1e4; 346 347 int N,M; 348 Point P[maxn], C; 349 350 int main () { 351 352 double x, y, p, q; 353 scanf("%d", &N); 354 for (int i = 0; i < N; i++) P[i].read(); 355 scanf("%d",&M); 356 for(int i=1;i<=M;i++){ 357 358 scanf("%lf%lf%lf%lf", &x, &y, &p, &q); 359 double S = getArea(P, N); 360 361 double l = 0, r = inf; 362 while (dcmp(r - l) > 0) { 363 double mid = (l + r) / 2; 364 double area = getPublicAreaToPolygon(Circle(Point(x, y), mid), P, N); 365 if (dcmp(q*area-q*S+p*S) >= 0) 366 r = mid; 367 else 368 l = mid; 369 } 370 double ans = (l + r) / 2; 371 printf("%.10lf\n", ans); 372 } 373 return 0; 374 }
扫一扫
1551
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
