2016 Multi-University Training Contest 1
A - Abandoned country
两个问题的组合，一个是求最小生成树，第二个是树形dp
// // main.cpp // multi2016.1.A // // Created by New_Life on 16/7/25. // Copyright © 2016年 chenhuan001. All rights reserved. // #include <iostream> #include <stdio.h> #include <string.h> #include <vector> #include <algorithm> using namespace std; #define N 100100 int n,m; int bin[N]; struct node { int to,next,w; }edge[2*N]; int pre[N],cnt; void add_edge(int u,int v,int w) { edge[cnt].to = v; edge[cnt].w = w; edge[cnt].next = pre[u]; pre[u] = cnt++; } int findfa(int x) { if(x==bin[x]) return x; return bin[x] = findfa(bin[x]); } struct EDGE { int x,y,z; }E[1001000]; int cmpE(EDGE e1,EDGE e2) { return e1.z<e2.z; } long long all = 0; int dfs(int s,int fa) { int tcnt = 0; for(int p =pre[s];p!=-1;p=edge[p].next) { int v = edge[p].to; if(v == fa) continue; int w = edge[p].w; int tmp = dfs(v,s); all += (long long)w*tmp*(n-tmp); tcnt += tmp; } return tcnt+1; } int main() { int T; cin>>T; while(T--) { cnt = 0; memset(pre,-1,sizeof(pre)); scanf("%d%d",&n,&m); for(int i=1;i<=n;i++) bin[i] = i; for(int i=0;i<m;i++) { int x,y,z; scanf("%d%d%d",&x,&y,&z); E[i].x = x; E[i].y = y; E[i].z = z; } sort(E, E+m, cmpE); long long ans=0; for(int i=0;i<m;i++) { int a = E[i].x; int b = E[i].y; int w = E[i].z; int fa = findfa(a); int fb = findfa(b); if(fa != fb) { bin[fa] = fb; ans += w; add_edge(a,b,w); add_edge(b,a,w); } } //然后搞一次dfs,e_e 自己写吧 all = 0; dfs(1,-1); double ansp = (double)all/((double)n*(n-1)/2); printf("%lld %.2lf\n",ans,ansp); } return 0; }
B - Chess
经典的多游戏博弈，了解sg函数即可秒之。
// // main.cpp // multi2016.1.b // // Created by New_Life on 16/7/25. // Copyright © 2016年 chenhuan001. All rights reserved. // #include <iostream> #include <stdio.h> #include <string.h> #include <algorithm> using namespace std; //求下sg函数 int mark[1<<21]; void dfs(int s) { int pre = -1; int next[22]; memset(next,0,sizeof(next)); for(int i=19;i>=0;i--) { if( ((1<<i)&s) != 0) { if(pre !=-1) { //跳到这里 int nexts = (s-(1<<i))|(1<<pre); if(mark[nexts] == -1) dfs(nexts); next[ mark[nexts] ] = 1; } } else pre = i; } for(int i=0;i<=20;i++) if(next[i] == 0) { mark[s] = i; return ; } } int main() { //满满地套路 memset(mark,-1,sizeof(mark)); for(int i=0;i<(1<<20);i++) { if(mark[i] != -1) continue; dfs(i); } int T; cin>>T; while(T--) { int n; scanf("%d",&n); int ans = 0; for(int i=0;i<n;i++) { int m; scanf("%d",&m); int tmp=0; for(int j=0;j<m;j++) { int x; scanf("%d",&x); tmp |= (1<<(x-1)); } ans ^= mark[tmp]; } if(ans == 0) printf("NO\n"); else printf("YES\n"); } return 0; }
D - GCD
用普通的st算法预处理出区间的gcd，然后二分即可
// // main.cpp // multi2016.1.first // // Created by New_Life on 16/7/25. // Copyright © 2016年 chenhuan001. All rights reserved. // #include <iostream> #include <stdio.h> #include <string.h> #include <map> #include <algorithm> using namespace std; #define N 100100 int dp[N][20]; int n; int s[N]; map<int,int>cao; int _gcd(int b,int d) { if(b == 0) return d; return _gcd(d%b,b); } void RMQ_init() { for(int i=1; i<=n; i++) dp[i][0]=s[i]; for(int j=1; (1<<j)<=n; j++) for(int i=1;i+(1<<j)-1<=n;i++) dp[i][j] = _gcd(dp[i][j-1],dp[i+(1<<(j-1))][j-1]); } int RMQ(int L,int R) { int k=0; while((1<<(k+1))<=R-L+1) k++; return _gcd(dp[L][k],dp[R-(1<<k)+1][k]); } long long save[N]; //long long save[N]; int ans[N]; int main() { int T; int tt = 1; cin>>T; while(T--) { cao.clear(); scanf("%d",&n); for(int i=1;i<=n;i++) scanf("%d",s+i); RMQ_init(); printf("Case #%d:\n",tt++); int q; scanf("%d",&q); int id = 1; for(int i=0;i<q;i++) { int b,d; scanf("%d%d",&b,&d); //printf("%d %I64d\n",RMQ(b,d),save[ RMQ(b, d) ]); int gcd=RMQ(b,d); //printf("%d %I64d\n",RMQ(b,d),0LL); if(cao[gcd] == 0) { cao[gcd] = id++; } ans[i] = gcd; } memset(save, 0, sizeof(save)); for(int i=1;i<=n;i++) { int tmp = s[i]; int b = i, d = n; do { int preb = b; d = n; tmp = RMQ(i, b);// while(b<d) { int mid = (b+d+1)/2; if(RMQ(i,mid)<tmp) d = mid-1; else b = mid; } int tid = cao[tmp]; if(tid != 0) { save[ tid ] += (b-preb+1);//搞错这个了 } b++; }while(b<=n); } for(int i=0;i<q;i++) { //printf("%d %I64d\n",ans[i],save[ cao[ans[i]] ]); printf("%d %lld\n",ans[i],save[ cao[ans[i]] ]); } } return 0; }
E - Necklace
暴力枚举阴珠子的圆排列情况，然后跑匈牙利就行了
// // main.cpp // multi2016.1.E // // Created by New_Life on 16/7/26. // Copyright © 2016年 chenhuan001. All rights reserved. // #include <iostream> #include <stdio.h> #include <stdlib.h> #include <string.h> #include <algorithm> using namespace std; #define N 11 int n,m; int nolink[N][N]; int savex[N]; int mark[N]; int mark1[N]; int mat[N][N]; int ans; int pre[N]; int dfs1(int s) { for(int i=1;i<=n;i++) { if(mat[s][i] == 0 || mark[i] == 1) continue; mark[i] = 1; if(pre[i] == -1 || dfs1(pre[i])) { pre[i] = s; return 1; } } return 0; } void gao() { //建图 memset(mat,0,sizeof(mat)); for(int i=1;i<=n;i++) { for(int j=1;j<=n;j++) { if(nolink[ savex[i] ][j]==0 && nolink[ savex[i==n?1:i+1] ][j]==0) { mat[j][i] = 1; } } } memset(pre, -1,sizeof(pre)); memset(mark,0,sizeof(mark)); int tmp=0; for(int i=1;i<=n;i++) { memset(mark,0,sizeof(mark)); tmp += dfs1(i); } ans = max(ans,tmp); } void dfs(int s) { if(s>n) { gao(); return ; } for(int i=2;i<=n;i++) { if(mark1[i] == 1) continue; mark1[i] = 1; savex[s] = i; dfs(s+1); mark1[i] = 0; } } int main(int argc, const char * argv[]) { while(scanf("%d%d",&n,&m)!=EOF) { memset(nolink,0,sizeof(nolink)); //换一换 for(int i=0;i<m;i++) { int x,y; scanf("%d%d",&x,&y); nolink[y][x] = 1; } //暴力枚举出X边的所有情况。 memset(mark1,0,sizeof(mark1)); savex[1] = 1; mark[1] = 1; ans = 0; dfs(2); printf("%d\n",n-ans); } return 0; } /* 3 5 1 2 1 3 2 1 2 2 3 1 */
F - PowMod
这题又是两个问题的叠加，第一个问题比较巧妙。把问题由某个素因子分解成子问题，妙啊，妙啊
// // main.cpp // multi2016.1.F // // Created by New_Life on 16/7/26. // Copyright © 2016年 chenhuan001. All rights reserved. // #include <iostream> #include <string.h> #include <stdlib.h> #include <algorithm> #include <stdlib.h> using namespace std; #define MOD 1000000007 int phi[10000100]; long long sum[10000100]; void getphis(int maxn) { phi[1]=1; phi[2]=1; for(int i=2;i<=maxn;i++) phi[i]=i; for(int i=2;i<=maxn;i+=2) phi[i]/=2; for(int i=3;i<=maxn;i+=2) { if(phi[i]==i)//为素数 { for(int j=i;j<=maxn;j+=i) { phi[j]=phi[j]-phi[j]/i; } } } for(int i=1;i<=maxn;i++) sum[i] = (sum[i-1]+phi[i])%MOD; } int save[10]; int cnt=0; long long gao(int n,int m,int sign) { int tmp = -1; int id = 0; for(int i=cnt-1;i>=0;i--) { if( ((1<<i)&sign) != 0) { id = i; tmp = save[i]; break; } } if(tmp == -1) { return sum[m]; } if(m==1) { return phi[n]; } if(m==0) return 0; return (gao(n,m/tmp,sign)+(long long)phi[tmp]*gao(n/tmp,m,sign-(1<<id)))%MOD; } /************** 快速幂模板 调用:Quk_Mul(a,b,mod) 返回:a^b%mod 复杂度:当mod>10^9，log(mod)*log(b),否则log(b) ***************/ long long Mod_Mul(long long a,long long b,long long mod) { long long msum=0; while(b) { if(b&1) msum = (msum+a)%mod; b>>=1; a = (a+a)%mod; } return msum; } long long Quk_Mul(long long a,long long b,long long mod) { bool qmflag=mod>1e9?1:0; long long qsum=1; while(b) { if(b&1) qsum = (qmflag==1) ? Mod_Mul(qsum,a,mod) : (qsum*a)%mod; b>>=1; a = (qmflag==1) ? Mod_Mul(a,a,mod) : (a*a)%mod; } return qsum; } long long fun(long long x,int p) { if(p==1) { return 0; } return Quk_Mul(x, fun(x,phi[p]) + phi[p], p); } int main() { getphis(10000000); int n,m,p; while(scanf("%d%d%d",&n,&m,&p)!=EOF) { int tn = n; cnt=0; for(int i=2;i*i<=tn;i++) { if(tn%i == 0) { tn/=i; save[cnt++] = i; } } if(tn!=1) save[cnt++] =tn; long long k = gao(n,m,(1<<cnt)-1); cout<<fun(k,p)<<endl; } return 0; } /* 1 2 6 1 100 9 */
G - Rigid Frameworks
题意很恐怖，想了挺久一直往找特殊规律，然后状压dp搞。 结果竟然可以看成，二分完全图的联通子图数。 仔细想想确实是可以这样转换。
然后就是一个经典的dp了。
我的做法稍有点不同，但是本质是一样的。
// // main.cpp // hdu5729 // // Created by New_Life on 16/7/27. // Copyright © 2016年 chenhuan001. All rights reserved. // #include <iostream> #include <stdio.h> #include <math.h> #include <string.h> #include <stdlib.h> #include <algorithm> using namespace std; #define N 110 #define MOD 1000000007 //组合数打表模板,适用于N<=3000 //c[i][j]表示从i个中选j个的选法。 long long C[N][N]; long long pow2[N]; void get_C(int maxn) { C[0][0] = 1; for(int i=1;i<=maxn;i++) { C[i][0] = 1; for(int j=1;j<=i;j++) // C[i][j] = C[i-1][j]+C[i-1][j-1]; C[i][j] = (C[i-1][j]+C[i-1][j-1])%MOD; } } long long dp[11][11][N]; long long dfs(int n,int m,int k) { if(dp[n][m][k]!=-1) return dp[n][m][k]; if(k<n+m-1) return dp[n][m][k] = 0; long long all = C[n*m][k]; all -= C[(n-1)*m][k]; long long sum = 0; for(int i=1;i<=n;i++) for(int j=1;j<=m;j++) for(int p=i+j-1;p<=k;p++) { if(i+j==n+m) continue; long long tmp = (dfs(i,j,p)*C[n-1][i-1])%MOD; tmp = tmp*C[m][j]%MOD; tmp = tmp*C[(n-i)*(m-j)][k-p]%MOD; sum += tmp; sum %= MOD; } return dp[n][m][k]=((all-sum)%MOD+MOD)%MOD; } int main() { for(int i=0;i<=100;i++) { long long tmp = 1; for(int j=1;j<=i;j++) tmp = (tmp*2)%MOD; pow2[i] = tmp; } get_C(100); memset(dp,-1,sizeof(dp)); int n,m; while(cin>>n>>m) { long long ans = 0; for(int i=1;i<=n*m;i++) { ans = ans + (dfs(n,m,i)*pow2[i])%MOD; ans %= MOD; } ans = (ans+MOD)%MOD; cout<<ans<<endl; } return 0; }
而且这题数据最多只有10*10,所以完全可以打表啊。
写了一个打表程序，10s之内可以跑出，但是比标准解法还要复杂。--!
// // main.cpp // multi2016.1.G // // Created by New_Life on 16/7/27. // Copyright © 2016年 chenhuan001. All rights reserved. // #include <iostream> #include <stdio.h> #include <string> #include <stdlib.h> #include <algorithm> using namespace std; #define MOD 1000000007 typedef long long ll; ll save[1100]; ll save1[11][11]; ll pow2[11]; ll mark[1<<10][1<<10][10][2]; ll dfs(int L,int R,int flag,ll cnt) { if(mark[L][R][cnt][flag]!=-1) return mark[L][R][cnt][flag]; if(L==0 && R==0) return 1; if(L==0) return mark[L][R][cnt][flag]=save1[save[R]][cnt]; if(R==0) return mark[L][R][cnt][flag]=save1[save[L]][cnt]; ll tcnt = 0; if(flag == 0)//左边选择 { int tr = R; for(;tr!=0;tr = ((tr-1)&R) )//必定选择一个以上 { tcnt += save1[ save[tr] ][cnt]*dfs(L, R-tr, flag^1, save[tr]); tcnt %= MOD; } } else//右边选择 { int tl = L; for(;tl!=0;tl = ((tl-1)&L) )//必定选择一个以上 { tcnt += save1[ save[tl] ][cnt]*dfs(L-tl, R, flag^1, save[tl]); tcnt %= MOD; } } return mark[L][R][cnt][flag]=tcnt; } ll C[12][12]; void get_C() { C[0][0] = 1; for(int i=1;i<=10;i++) { C[i][0] = 1; for(int j=1;j<=i;j++) C[i][j] = C[i-1][j]+C[i-1][j-1]; } } void get_fff() { for(int n=1;n<=10;n++) for(int m=1;m<=10;m++) { ll add[110]; for(int j=1;j<=m;j++) add[j] = C[m][j]; ll tmp[2][110]; memset(tmp,0,sizeof(tmp)); int a = 0; tmp[a][0] = 1; for(int i=1;i<=n;i++) { memset(tmp[a^1],0,sizeof(tmp[a^1])); for(int j=100;j>=1;j--) { for(int k=1;k<=m;k++) { if(j-k < 0) continue; tmp[a^1][j] = (tmp[a^1][j] + tmp[a][j-k]*add[k])%MOD; } } a = a^1; } ll sum = 0; for(int i=1;i<=100;i++) { sum = (sum+pow2[i]*tmp[a][i])%MOD; } save1[n][m] = sum; } } long long ans[10][10]={2LL,4LL,8LL,16LL,32LL,64LL,128LL,256LL,512LL,1024LL,4LL,48LL,448LL,3840LL,31744LL,258048LL,2080768LL,16711680LL,133955584LL,72693241LL,8LL,448LL,15008LL,429568LL,11596928LL,306009088LL,2369992LL,530422352LL,523796386LL,215597769LL,16LL,3840LL,429568LL,38529024LL,207602155LL,200364260LL,160231949LL,554488410LL,190587140LL,810639015LL,32LL,31744LL,11596928LL,207602155LL,519052755LL,105908051LL,72277468LL,23284327LL,286851052LL,865424583LL,64LL,258048LL,306009088LL,200364260LL,105908051LL,126514477LL,312460128LL,711245882LL,987939142LL,374773317LL,128LL,2080768LL,2369992LL,160231949LL,72277468LL,312460128LL,25696077LL,433482528LL,357533852LL,103280608LL,256LL,16711680LL,530422352LL,554488410LL,23284327LL,711245882LL,433482528LL,162005329LL,771914613LL,811634190LL,512LL,133955584LL,523796386LL,190587140LL,286851052LL,987939142LL,357533852LL,771914613LL,585919048LL,573168128LL,1024LL,72693241LL,215597769LL,810639015LL,865424583LL,374773317LL,103280608LL,811634190LL,573168128LL,935300639LL}; int main(int argc, const char * argv[]) { for(int i=0;i<(1<<10);i++) { int tmp = 0; for(int j=0;j<10;j++) { if( (i&(1<<j))!=0 ) tmp++; } save[i] = tmp; } for(int i=0;i<=100;i++) { int tmp =1 ; for(int j=0;j<i;j++) { tmp *= 2; tmp %= MOD; } pow2[i] = tmp; } //预处理3 get_C(); get_fff(); memset(mark,-1,sizeof(mark)); for(int n=1;n<=10;n++) { for(int m=1;m<=10;m++) { //cout<<dfs((1<<n)-1-1,(1<<m)-1,0,1)<<"LL,"; } } int n,m; while(cin>>n>>m) { cout<<ans[n-1][m-1]<<endl; } return 0; }
K - tetrahedron
计算几何，要是知道求四面体内切球心的公式当然可以瞬间秒，但是本弱完全不知，只能和它肛道理
找出面面之间的平分面，交点即使球心坐标。
如何面面的平分面呢。
1.直接用面面交板子，求出交线。
2.找到一个垂直于交线，且顶点在分别在两个平面上的三角形，然后求得一个面的点。
#include <iostream> #include <cmath> #include <stdio.h> #include <vector> #include <string.h> #include <stdlib.h> #include <algorithm> using namespace std; #define MAX_N 110 /*------------------常量区-------------------*/ const double INF = 1e10; // 无穷大 const double EPS = 1e-5; // 计算精度 const double PI = acos(-1.0);// PI const int PINXING = 0; // 平行 const int XIANGJIAO = 1; // 相交 const int XIANGLI = 0; // 相离 const int GONGXIAN = 2; // 共线 const int CHONGDIE = -1; // 重叠 const int INSIDE = 1; // 点在图形内部 const int OUTSIDE = 0; // 点在图形外部 const int BORDER = 2; // 点在图形边界 /*-----------------类型定义区----------------*/ struct Point { // 二维点或矢量 double x, y; //double angle, dis; Point() {} Point(double x0, double y0): x(x0), y(y0) {} void read() { scanf("%lf%lf",&x,&y); } }; struct Line { // 二维的直线或线段 Point p1, p2; Line() {} Line(Point p10, Point p20): p1(p10), p2(p20) {} void read() { scanf("%lf%lf%lf%lf",&p1.x,&p1.y,&p2.x,&p2.y); } }; struct Rect { // 用长宽表示矩形的方法 w, h分别表示宽度和高度 double w, h; Rect() {} Rect(double _w,double _h) : w(_w),h(_h) {} }; struct Rect_2 { // 表示矩形，左下角坐标是(xl, yl)，右上角坐标是(xh, yh) double xl, yl, xh, yh; Rect_2() {} Rect_2(double _xl,double _yl,double _xh,double _yh) : xl(_xl),yl(_yl),xh(_xh),yh(_yh) {} }; struct Circle { //圆 Point c; double r; Circle() {} Circle(Point _c,double _r) :c(_c),r(_r) {} }; typedef vector<Point> Polygon; // 二维多边形 typedef vector<Point> Points; // 二维点集 /*-------------------基本函数区---------------------*/ inline double max(double x,double y) { return x > y ? x : y; } inline double min(double x, double y) { return x > y ? y : x; } inline bool ZERO(double x) // x == 0 { return (fabs(x) < EPS); } inline bool ZERO(Point p) // p == 0 { return (ZERO(p.x) && ZERO(p.y)); } inline bool EQ(double x, double y) // eqaul, x == y { return (fabs(x - y) < EPS); } inline bool NEQ(double x, double y) // not equal, x != y { return (fabs(x - y) >= EPS); } inline bool LT(double x, double y) // less than, x < y { return ( NEQ(x, y) && (x < y) ); } inline bool GT(double x, double y) // greater than, x > y { return ( NEQ(x, y) && (x > y) ); } inline bool LEQ(double x, double y) // less equal, x <= y { return ( EQ(x, y) || (x < y) ); } inline bool GEQ(double x, double y) // greater equal, x >= y { return ( EQ(x, y) || (x > y) ); } // 输出浮点数前，防止输出-0.00调用该函数进行修正！ inline double FIX(double x) { return (fabs(x) < EPS) ? 0 : x; } /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ //-------------------3D 区域----------------------------// struct Point3D { //三维点或矢量 double x, y, z; Point3D() {} Point3D(double x0, double y0, double z0): x(x0), y(y0), z(z0) {} void read() { scanf("%lf%lf%lf",&x,&y,&z); } }; struct Line3D { // 三维的直线或线段 Point3D p1, p2; Line3D() {} Line3D(Point3D p10, Point3D p20): p1(p10), p2(p20) {} void read() { scanf("%lf%lf%lf%lf%lf%lf",&p1.x,&p1.y,&p1.z,&p2.x,&p2.y,&p2.z); } }; struct Area3D{ Point3D p1,p2,p3; Area3D(){} Area3D(Point3D p10, Point3D p20,Point3D p30): p1(p10), p2(p20), p3(p30){} void read() { p1.read(); p2.read(); p3.read(); //scanf("%lf%lf%lf%lf%lf%lf%lf%lf%lf",&p1.x,&p1.y,&p1.z,&p2.x,&p2.y,&p2.z,&p3.x,&p3.y,&p3.z); } }; inline bool ZERO(Point3D p) // p == 0 { return (ZERO(p.x) && ZERO(p.y) && ZERO(p.z)); } ////////////////////////////////////////////////////////////////////////////////////// //三维矢量运算 bool operator==(Point3D p1, Point3D p2) { return ( EQ(p1.x, p2.x) && EQ(p1.y, p2.y) && EQ(p1.z, p2.z) ); } bool operator<(Point3D p1, Point3D p2) { if (NEQ(p1.x, p2.x)) { return (p1.x < p2.x); } else if (NEQ(p1.y, p2.y)) { return (p1.y < p2.y); } else { return (p1.z < p2.z); } } Point3D operator+(Point3D p1, Point3D p2) { return Point3D(p1.x + p2.x, p1.y + p2.y, p1.z + p2.z); } Point3D operator-(Point3D p1, Point3D p2) { return Point3D(p1.x - p2.x, p1.y - p2.y, p1.z - p2.z); } Point3D operator * (const Point3D& A, double p) { return Point3D(A.x*p, A.y*p, A.z*p); } Point3D operator / (const Point3D& A, double p) { return Point3D(A.x/p, A.y/p, A.z/p); } Point3D operator*(Point3D p1, Point3D p2) // 计算叉乘 p1 x p2 { return Point3D(p1.y * p2.z - p1.z * p2.y, p1.z * p2.x - p1.x * p2.z, p1.x * p2.y - p1.y * p2.x ); } double operator&(Point3D p1, Point3D p2) { // 计算点积 p1·p2 return (p1.x * p2.x + p1.y * p2.y + p1.z * p2.z); } double Norm(Point3D p) // 计算矢量p的模 { return sqrt(p.x * p.x + p.y * p.y + p.z * p.z); } //取平面法向量 Point3D GetV(Area3D s){ return (s.p1-s.p2)*(s.p2-s.p3); } //判三点共线 int PointOnLine(Point3D p1,Point3D p2,Point3D p3){ return ZERO( (p1-p2)*(p2-p3) ); } //判四点共面 int PointOnArea(Point3D a,Point3D b,Point3D c,Point3D d){ return ZERO( GetV(Area3D(a, b, c))&(d-a) ); } //求三维空间中两点间的距离 double Dis(Point3D p1, Point3D p2) { return sqrt((p1.x-p2.x)*(p1.x-p2.x)+(p1.y-p2.y)*(p1.y-p2.y)+(p1.z-p2.z)*(p1.z-p2.z)); } // 求三维空间中点到直线的距离 double Dis(Point3D p, Line3D L) { if(L.p1==L.p2) return Dis(p, L.p1); return Norm((p - L.p1) * (L.p2 - L.p1)) / Norm(L.p2 - L.p1); } bool OnLine(Point3D p, Line3D L) // 判断三维空间中点p是否在直线L上 { if(L.p1==L.p2 && p==L.p1) return true;//共点时，返回true return ZERO( (p - L.p1) * (L.p2 - L.p1) ); } bool OnLineSeg(Point3D p, Line3D L) // 判断三维空间中点p是否在线段l上 { return ( ZERO((L.p1 - p) * (L.p2 - p)) && EQ( Norm(p - L.p1) + Norm(p - L.p2), Norm(L.p2 - L.p1)) ); } // 点p到平面Ap-Al的距离。 double Dis(Point3D p, Point3D Ap, Point3D Al) { return fabs((p-Ap)&Al)/Norm(Al); // 如果不取绝对值，得到的是有向距离 } // 点p在平面Ap-Al上的投影。 Point3D PointToArea(Point3D p,Point3D Ap, Point3D Al) { Al=Al/(Norm(Al));//把Al变成法向量。 return p-Al*((p-Ap)&Al); } //得到点p到直线L的距离，并返回p到直直线L的最近点rep double PointToLine(Point3D p,Line3D L,Point3D& rep) { if(L.p1==L.p2) { rep=L.p1; return Dis(p,L.p1); } Point3D a,b; a = L.p2-L.p1; b = p-L.p1; double dis12 = Dis(L.p1,L.p2); double dis = ( Norm(a*b) )/dis12; double k = (a&b)/(Norm(a)*dis12) ; rep = L.p1+(L.p2-L.p1)*k; return dis; } //求两条直线之间的关系（三维） //输入：两条不为点的直线 //输出：相交返回XIANGJIAO和交点p,平行返回PINGXING,共线返回GONGXIAN int LineAndLine(Line3D L1,Line3D L2,Point3D &p) { Point3D px,py; px = L1.p1 - L1.p2; py = L2.p1 - L2.p2; if( ZERO(px*py) )//平行或者共线 { if( ZERO( (L2.p1-L1.p1)*py ) ) //共线 { return GONGXIAN; } return PINXING; } //判断是否共面 Point3D tp=(L1.p1-L2.p1)*py; if( !ZERO(tp&px) ) return XIANGLI;//XIANGLI与平行相同 p = L1.p1; Point3D tp1=(L2.p1-L1.p1)*(L2.p1-L2.p2); Point3D tp2=(L1.p2-L1.p1)*(L2.p1-L2.p2); double _t = Norm(tp1)/Norm(tp2); //tp1和tp2肯定是共线的，如果反向则_t 为负 if( LT( (tp1&tp2),0 ) ) _t*=-1; p.x += (L1.p2.x-L1.p1.x)*_t; p.y += (L1.p2.y-L1.p1.y)*_t; p.z += (L1.p2.z-L1.p1.z)*_t; return XIANGJIAO; } //空间两直线最近点对。直线不能平行，直线不能为点. //ans1为直线a1,b1上的最近点 Point3D ans1,ans2; double SegSegDistance(Point3D a1, Point3D b1, Point3D a2, Point3D b2) { Point3D v1 = (a1-b1), v2 = (a2-b2); Point3D N = v1*v2; Point3D ab = (a1-a2); double ans = (N&ab) / Norm(N); Point3D p1 = a1, p2 = a2; Point3D d1 = b1-a1, d2 = b2-a2; double t1, t2; t1 = ((p2-p1)*d2 )&(d1*d2); t2 = ((p2-p1)*d1 )&(d1*d2); double dd = Norm( (d1*d2) ); t1 /= dd*dd; t2 /= dd*dd; ans1=a1+(b1-a1)*t1; ans2=a2+(b2-a2)*t2; return fabs(ans); } //直线与平面交 int LineAndArea(Line3D l1,Point3D Ap,Point3D Al,Point3D &p) { //输入直线，和平面的点法式 //第一步，判断直线与平面是否平行。 if( ZERO((l1.p2-l1.p1)&Al) ) return 0;//直线与平面平行。 double _t =( (Ap-l1.p1)&Al ) / ((l1.p1-l1.p2)&Al); p = l1.p1+(l1.p1-l1.p2)*_t; return 1; } void dfs(int x,double &len) { len++; dfs(x-1,len); dfs(x-2,len); } //空间两直线最近点对 //注意：直线不能平行 double LineAndLine(Line3D l1,Line3D l2,Point3D &p1,Point3D &p2) { //先求出法向量 Point3D v1,v2; v1 = l1.p2-l1.p1; v2 = l2.p2-l2.p1; Point3D vt=v1*v2; //然后先把l2投影到 l1所在的平面上 double len = ((l2.p1-l1.p1)&vt)/Norm(vt); double normvt = -len/Norm(vt); vt.x = vt.x*normvt; vt.y = vt.y*normvt; vt.z = vt.z*normvt; Line3D tl2; tl2.p1 = l2.p1+vt; tl2.p2 = l2.p2+vt; int sign=LineAndLine(l1, tl2, p1); /* //测试用 if(sign!=XIANGJIAO) { int x=0; printf("%lf\n",len/x); dfs(100000000,len); } */ return fabs(len); } //已知空间四面体6条边，求体积 double P( double a,double b,double c,double d,double e ){ return a*(b*c-d*e); } double Get4V(int OA,int OB,int OC,int AB,int CA,int BC) { OA*=OA;OB*=OB;OC*=OC;AB*=AB;CA*=CA;BC*=BC; double ans=0; ans+=P( OA,OB,OC,(OB+OC-BC)/2.0,(OB+OC-BC)/2.0 ); ans-=P( (OA+OB-AB)/2.0,(OA+OB-AB)/2.0,OC,(OA+OC-CA)/2.0,(OB+OC-BC)/2.0 ); ans+=P( (OA+OC-CA)/2.0,(OA+OB-AB)/2.0,(OB+OC-BC)/2.0,OB,(OA+OC-CA)/2.0); return sqrt(ans/36); } //求两面相交,平行或共面返回PINGXING,否则返回XIANGJIAO和直线rel int AreaAndArea(Area3D a1,Area3D a2,Line3D &rel) { Point3D va1 = GetV(a1),va2 = GetV(a2); Point3D lv = va1*va2;//相交直线的方向向量 if( ZERO(lv) )//平行 { return PINXING; } //然后得到某一个交点 Point3D p1; if(LineAndArea(Line3D(a1.p1,a1.p2), a2.p1, va2, p1) == 0) if(LineAndArea(Line3D(a1.p1,a1.p3), a2.p1, va2, p1) == 0) LineAndArea(Line3D(a1.p2,a1.p3), a2.p1, va2, p1); rel.p1 = p1; rel.p2 = p1 + (lv*10); return XIANGJIAO; } ////////////////////////////////////////////////////////////////////////////////////// /*---------------------代码区---------------------------*/ int equal3(Point3D p1,Point3D p2) { return p1.x==p2.x&&p1.y==p2.y&&p1.z==p2.z; } Area3D gettwoplanein(Area3D p1,Area3D p2,Point3D a,Point3D b) { Line3D l; l.p1 = a; l.p2 = b; Point3D pf1;//p1 的法向量 pf1 = GetV(p1); //找一个不在交线上的点 Point3D p1d,p2d; if( (!equal3(p1.p1,a)) && (!equal3(p1.p1, b)) ) p1d=p1.p1; if( (!equal3(p1.p2,a)) && (!equal3(p1.p2, b)) ) p1d=p1.p2; if( (!equal3(p1.p3,a)) && (!equal3(p1.p3, b)) ) p1d=p1.p3; if( (!equal3(p2.p1,a)) && (!equal3(p2.p1, b)) ) p2d=p2.p1; if( (!equal3(p2.p2,a)) && (!equal3(p2.p2, b)) ) p2d=p2.p2; if( (!equal3(p2.p3,a)) && (!equal3(p2.p3, b)) ) p2d=p2.p3; //然后对应到l上去 Point3D p1e,p2e; PointToLine(p1d, l, p1e); PointToLine(p2d, l, p2e); p2d = p2d + (p1e-p2e); //然后就是得到面上的一点 double l1 = Dis(p1d, p1e); double l2 = Dis(p1e, p2d); double l3 = Dis(p1d, p2d); double shita = acos( (l1*l1+l2*l2-l3*l3)/(2*l1*l2) ); shita/=2; double beta = acos( (l1*l1+l3*l3-l2*l2)/(2*l1*l3) ); beta = PI -shita -beta; double x1 = l1*sin(shita)/sin(beta); double x2 = l3-x1; Point3D pp; pp.x = p1d.x + (p2d.x-p1d.x)*x1/(x1+x2); pp.y = p1d.y + (p2d.y-p1d.y)*x1/(x1+x2); pp.z = p1d.z + (p2d.z-p1d.z)*x1/(x1+x2); Area3D p3; p3.p1 = a; p3.p2 = b; p3.p3 = pp; return p3; } int main() { Point3D A,B,C,D; while(scanf("%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf",&A.x,&A.y,&A.z,&B.x,&B.y,&B.z,&C.x,&C.y,&C.z,&D.x,&D.y,&D.z)!=EOF) { Area3D ABC,ABD,BCD,ACD; ABC.p1 = A; ABC.p2 = B; ABC.p3 = C; ABD.p1 = A; ABD.p2 = B; ABD.p3 = D; BCD.p1 = B; BCD.p2 = C; BCD.p3 = D; ACD.p1 = A; ACD.p2 = C; ACD.p3 = D; //step one, 判断是存在四面体 if( PointOnArea(A,B,C,D) ) { printf("O O O O\n"); continue; } // 先得到三个角平分面 Area3D p1,p2,p3; p1 = gettwoplanein(ABD, BCD, B, D);//讲道理，现在已经求出 p2 = gettwoplanein(ABC, BCD, B, C); p3 = gettwoplanein(ABC, ACD, A, C); //然后求三个面的交点 Line3D ll; AreaAndArea(p1, p2, ll); Point3D ans; LineAndArea(ll, p3.p1, GetV(p3), ans); double dis = Dis(ans, ABC.p1,GetV(ABC)); printf("%.4lf %.4lf %.4lf %.4lf\n",ans.x,ans.y,ans.z,dis); } return 0; }
