http://acm.hdu.edu.cn/showproblem.php?pid=5839
给定空间中n个不相同的点，求所有特殊四面体的个数。(n <= 200)
特殊四边形是满足以下条件的四面体：
1、至少有4条边相等。
2、如果只有四条边相等，则剩余两条边不相邻。
分析题意，可以知道我们要求的是这样的四面体：
图中四条红边必须相等，但黑边不一定相等。蓝色平面是i，j两点所连线段的垂直平分面，由此可知两条蓝边也必须相等。
由此可以很自然地得到做法：枚举i和j，对所有在ij中垂面上的点k，计算其到ij的距离，排序后相同距离的归为一组。若一组有m个，则有m*(m-1)/2种方法可以生成目标四面体。
为了避免计算重复并减小常数，限制 i < j,k，即对于某个四面体，取其四个顶点中编号最小的那个点 i，用 i 所在的黑边去给它计数。时间复杂度O(n^3*logn)
不过，麻烦的问题是正四面体和四点共面的存在。
假设不是正四面体的所求四面体数为x，正四面体数为y，满足上述统计方法中性质的四点共面有z种。那么直接按上述方法统计所得结果为 x+3y+z。注意正四面体会被与i相连的三条边统计三次。
如何去掉正四面体和四点共面的影响呢？只要统计一下他们两类的个数就行了。具体方法和多校 How many triangles 类似，对于同一类的，按照角度排序（图中theta），再扫一遍就行了。不过这是三维的，角度的处理比较复杂。我的具体方法是以k循环中第一个点k1为起始点（角度为0），用点积算出虚拟二维坐标系中的x，用叉积算出y（因为三维叉积是一个向量，要取其模为y的模，根据向量和ij向量是否同向确定其正负），一个atan2就ok了。
其它技巧，因为题目中所有数字均为正数而且 -2000<=x<=2000，所以除了极角排序用double外，其余全用int就行了。输入的所有点坐标都可以乘以二，这样线段 ij 的中点坐标也可以用整数表示了。
具体代码因为复制粘贴，比较长：
1 #include <iostream> 2 #include <cstdio> 3 #include <string> 4 #include <cstring> 5 #include <algorithm> 6 #include <cmath> 7 using namespace std; 8 9 const int N = 205; 10 const double pi = acos(-1.0); 11 12 struct point { 13 int x,y,z; 14 point() {} 15 point(int xx,int yy,int zz) { 16 x = xx,y = yy,z = zz; 17 } 18 void input() { 19 cin >> x >> y >> z; 20 x *= 2, y *= 2, z *= 2; 21 } 22 point operator -(const point &n2) const { 23 return point(x-n2.x, y-n2.y, z-n2.z); 24 } 25 point operator +(const point &n2) const { 26 return point(x+n2.x, y+n2.y, z+n2.z); 27 } 28 point operator /(const int n) const { 29 return point(x/n, y/n, z/n); 30 } 31 bool operator ==(const point &n2) const { 32 return x==n2.x && y==n2.y && z==n2.z; 33 } 34 }a[N]; 35 36 struct sta { 37 int x, id; 38 sta() {} 39 sta(int xx, int idd) { 40 x = xx; id = idd; 41 } 42 bool operator <(const sta &n2) const { 43 return x<n2.x; 44 } 45 }b[N*2]; 46 47 struct angle { 48 double x; int id; 49 angle() {} 50 angle(double xx, int idd) { 51 x = xx; id = idd; 52 } 53 bool operator <(const angle &n2) const { 54 return x<n2.x; 55 } 56 void output() { 57 cout << "angle " << x <<' ' << id << endl; 58 } 59 }c[N*2]; 60 61 point det(const point &p1, const point &p2) 62 { 63 return point(p1.y*p2.z-p2.y*p1.z, -p1.x*p2.z+p2.x*p1.z, p1.x*p2.y-p2.x*p1.y); 64 } 65 int len_sqr(const point &p) 66 { 67 return p.x*p.x + p.y*p.y + p.z*p.z; 68 } 69 int dot(const point &p1, const point &p2) 70 { 71 return p1.x*p2.x + p1.y*p2.y + p1.z*p2.z; 72 } 73 int sig(int x) 74 { 75 return x==0? 0: (x>0? 1:-1); 76 } 77 78 int main() 79 { 80 int test,n,i,j,k; 81 cin >> test; 82 for(int t = 1; t<=test; t++) 83 { 84 cin >> n; 85 for(i = 1; i<=n; i++) a[i].input(); 86 long long ans = 0;// x+3y+z 87 long long ans2 = 0;// 3y 88 long long ans3 = 0;//z 89 for(i = 1; i<=n; i++) 90 for(j = i+1; j<=n; j++) 91 { 92 point mid = (a[i]+a[j])/2, z = a[i]-a[j]; 93 int l2 = len_sqr(z); 94 int tot = 0; 95 for(k = i+1; k<=n; k++) 96 { 97 if(dot(a[k]-mid, z)==0) 98 b[++tot] = sta(len_sqr(mid-a[k]), k); 99 } 100 sort(b+1, b+1+tot); 101 102 int now; 103 for(k = 1; k<=tot; k++) 104 { 105 now = 1; 106 while(k<tot && b[k].x==b[k+1].x) { 107 now++; k++; 108 } 109 if(now>1) 110 ans += now*(now-1)/2; 111 } 112 113 for(k = 1; k<=tot; k++) 114 { 115 116 int from = k; 117 while(k<tot && b[k].x==b[k+1].x) { 118 k++; 119 } 120 int to = k; 121 if(b[k].x*4==l2*3) { 122 123 int m = 0; 124 for(k = from; k<=to; k++) 125 { 126 point temp = det(a[b[from].id]-mid, a[b[k].id]-mid); 127 int x = dot(a[b[from].id]-mid, a[b[k].id]-mid); 128 if(temp==point(0,0,0)) 129 { 130 if(x>0) 131 c[++m] = angle(.0, b[k].id); 132 else 133 c[++m] = angle(pi, b[k].id); 134 } 135 else { 136 double y = sqrt(len_sqr(temp)),theta; 137 if(sig(temp.x)==sig(z.x) && sig(temp.y)==sig(z.y) && sig(temp.z)==sig(z.z)) 138 theta = atan2(y, x); 139 else 140 theta = atan2(-y, x); 141 if(theta<0) theta += 2.0*pi; 142 c[++m] = angle(theta, b[k].id); 143 } 144 } 145 sort(c+1, c+1+m); 146 for(k = m+1; k<=m*2; k++) 147 { 148 c[k] = c[k-m]; 149 c[k].x += 2*pi; 150 } 151 int p = 1; 152 for(k = 1; k<=m; k++) 153 { 154 while(c[p+1].x-c[k].x<pi && dot(a[c[p+1].id]-mid, a[c[k].id]-mid)*4>=l2) 155 p++; 156 if(dot(a[c[p].id]-mid, a[c[k].id]-mid)*4==l2) ans2++; 157 } 158 } 159 // 160 int m = 0; 161 for(k = from; k<=to; k++) 162 { 163 point temp = det(a[b[from].id]-mid, a[b[k].id]-mid); 164 int x = dot(a[b[from].id]-mid, a[b[k].id]-mid); 165 if(temp==point(0,0,0)) 166 { 167 if(x>0) 168 c[++m] = angle(.0, b[k].id); 169 else 170 c[++m] = angle(pi, b[k].id); 171 } 172 else { 173 double y = sqrt(len_sqr(temp)),theta; 174 if(sig(temp.x)==sig(z.x) && sig(temp.y)==sig(z.y) && sig(temp.z)==sig(z.z)) 175 theta = atan2(y, x); 176 else 177 theta = atan2(-y, x); 178 if(theta<0) theta += 2.0*pi; 179 c[++m] = angle(theta, b[k].id); 180 } 181 } 182 sort(c+1, c+1+m); 183 for(k = m+1; k<=m*2; k++) 184 { 185 c[k] = c[k-m]; 186 c[k].x += 2*pi; 187 } 188 int p = 1; 189 for(k = 1; k<=m; k++) 190 { 191 while(1) { 192 if(p==2*m) break; 193 point temp = det(a[c[k].id]-mid, a[c[p+1].id]-mid); 194 if(c[p+1].x-c[k].x<4 &&( sig(temp.x)==sig(z.x) && sig(temp.y)==sig(z.y) && sig(temp.z)==sig(z.z) || temp==point(0,0,0))) 195 p++; 196 else break; 197 } 198 199 if(c[p].x-c[k].x > 1 && det(a[c[p].id]-mid, a[c[k].id]-mid)==point(0,0,0)) ans3++; 200 } 201 k = to; 202 } 203 204 } 205 206 207 ans = ans-ans2-ans3/2+ans2/3; 208 209 printf("Case #%d: %I64d ", t, ans); 210 } 211 return 0; 212 }