More servicesWindows Live
HomeHotmailSpacesOneCare
 
MSN
Sign in
 
 
Spaces home  因为选择,所以喜欢!PhotosProfileFriendsMore Tools Explore the Spaces community

因为选择,所以喜欢!

ACM算法博客 http://chenhaifeng.blog.edu.cn

[计算几何]两个圆相交的面积

#include<math.h>  
#include<stdio.h>  

#define pi acos(-1.0)
#define sqr(x) ((x) * (x))
#define dis2(a, b) (sqrt(sqr(a.x - b.x) + sqr(a.y - b.y)))

struct point
{
	double x, y;
};

struct circle
{
	point c;
	double r;
};

double IntersectionArea(circle cir1, circle cir2)
{  
	double A1, A2, d, ang1, ang2;
	circle ctmp;
	if (cir1.r < cir2.r)  
    {  
		ctmp = cir1;
		cir1 = cir2;
		cir2 = ctmp;	
    }
    double ret = 0;  
    A1 = pi * sqr(cir1.r);  
    A2 = pi * sqr(cir2.r);  
    d = dis2(cir1.c, cir2.c);
      
    if (d >= cir1.r + cir2.r) return 0;  
    if (d <= cir1.r - cir2.r) return A2;  
       
    ang1 = acos((sqr(cir1.r)+ sqr(d) - sqr(cir2.r)) / 2 / cir1.r/ d);  
    ang2 = acos((sqr(cir2.r)+ sqr(d) - sqr(cir1.r)) / 2 / cir2.r/ d);  
    
    ret -= sqr(cir1.r) * ang1 - sqr(cir1.r) * sin(ang1) * cos(ang1);
    ret -= sqr(cir2.r) * ang2 - sqr(cir2.r) * sin(ang2) * cos(ang2); 
    return -ret;
} 

int main()
{
	int test;
	circle cir1, cir2;
	scanf("%d", &test);
	while(test--)
	{
		scanf("%lf%lf%lf", &cir1.c.x, &cir1.c.y, &cir1.r);
		scanf("%lf%lf%lf", &cir2.c.x, &cir2.c.y, &cir2.r);
        printf("%.3lf\n", IntersectionArea(cir1, cir2));
	}
	return 0;
}

[计算几何]判断两个多边形是否相似(pku1931)

http://acm.pku.edu.cn/JudgeOnline/problem?id=1931 
/*
http://acm.pku.edu.cn/JudgeOnline/problem?id=1931

Biometrics

Description

Recently, the term Biometrics been used to refer to the emerging field of technology devoted to identification of individuals using biological traits, such as those based on retinal or iris scanning, fingerprints, or face recognition. 
A simple biometric system translates a human image into a polygon by considering certain features (eyes, nose, ears, etc.) to be vertices and connecting them with line segments. The polygon has distinct vertices but may be degenerate in that the line segments could intersect. Because these polygons are generally created from remote images, there is some uncertainty as to their scale and rotation. Your job is to determine whether or not two polygons are similar; that is, can they be made equal by repositioning, rotating and magnifying them?
Input

Input consists of several test cases. Each test case consists of three lines containing: 

f, the number of features 

f coordinate pairs giving the vertices of the first polygon 

f coordinate pairs giving the vertices of the second polygon 

The vertices for both polygons correspond to the same set of features in the same order; for example, right ear tip, chin cleft, right eye, nose, left eye, left ear tip, space between front teeth. Each polygon has f vertices; each vertex is given as an x and y coordinate pair. There are at least three and no more than ten features. Coordinates are integers between -1000 and 1000. A line containing 0 follows the last test case. 
Output

For each case, output a line "similar" or "dissimilar" as appropriate. The two polygons are similar if, after some combination of translation, rotation, and scaling (but not reflection) both vertices corresponding to each feature are in the same position. 
Sample Input

4
0 0 0 1 1 1 1 0
0 1 1 0 0 -1 -1 0
3
0 0 10 0 10 10
0 0 -10 0 -10 10
3
0 0 10 10 20 20
0 0 11 11 22 22
3
0 0 10 10 20 20
0 0 11 11 20 20
0

Sample Output

similar
dissimilar
similar
dissimilar

Source
Waterloo local 2003.09.27
*/
/*
判断两多边形相似,只要对应边成比例,且对应的角相当,
由于这题多边形可能存在边的自交,所以还要加上每个角
的转向相同
*/

#include <math.h>
#include <stdio.h>

#define eps 1e-6
#define sqr(x) ((x) * (x))
#define same(a, b) (fabs((a) - (b)) < eps) 
#define dis2(a, b) (sqrt(sqr(a.x - b.x) + sqr(a.y - b.y)))

struct point 
{
	double x, y;
	point operator-(point &b)
	{
		point c;
		c.x = x - b.x;
		c.y = y - b.y;
		return c;
	}
};

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 get_anlge(point p1, point p2, point p3)
{
	//不需要求反余弦,点击相当,反余弦相等,避免使用反三角函数
	return dot(p1 - p2, p3 - p2) / dis2(p1, p2) / dis2(p2, p3);
}

int get_dir(point p1, point p2, point p3)
{
	//根据三点得到转向
	double t1 = cross(p2 - p1, p3 - p2);
	if(fabs(t1) < eps) return 1;
	if(t1 < 0) return 2;
	else return 3;	
}

int slove(double ang1[], double ang2[], double len1[], double len2[], int dir1[], int dir2[], int n)
{
	/*由于题目已经告诉了对应点的匹配顺序,所以只需要从
	0,开始匹配,如果没有告诉对应的匹配顺序,就还有枚举
	匹配对应边*/
	//int k;
	/*for(int i = 0;i < n;i++)
	{
		for(int j = 0;j < n;j++)
		{
			//第i条边和第j条边相对应*/
			int s, t, k;
			s = 0;
			t = 0;
			for(k = 0;k < n;k++)
			{
				if(!same(ang1[s], ang2[t]) || !same(len1[s], len2[t]) || dir1[s] != dir2[t]) break;
				s++;
				t++;
				s %= n;
				t %= n;
			} 
			if(k == n) return 1;
		/*}
	}*/
	return 0;
}

double polygonArea(point p[], int n)
{
    //已知多边形各顶点的坐标,求其面积
    double area = 0.0;
    for(int i = 1;i <= n;i++)
        area += (p[i - 1].x * p[i % n].y - p[i % n].x * p[i - 1].y);  
    return area;   
}

int main()
{
	int n, similar;
	point p1[20], p2[20];
	double max1, max2;
	double ang1[20], ang2[20], len1[20], len2[20];
	int dir1[20], dir2[20];
	while(scanf("%d", &n) && n)
	{
		for(int i = 0;i < n;i++)
		{
			scanf("%lf%lf", &p1[i].x, &p1[i].y);
		}
		for(int i = 0;i < n;i++)
		{
			scanf("%lf%lf", &p2[i].x, &p2[i].y);
		}
		ang1[0] = get_anlge(p1[n - 1], p1[0], p1[1]);
		ang2[0] = get_anlge(p2[n - 1], p2[0], p2[1]);
		dir1[0] = get_dir(p1[n - 1], p1[0], p1[1]);
		dir2[0] = get_dir(p2[n - 1], p2[0], p2[1]);
		for(int i = 1;i < n;i++)
		{
			ang1[i] = get_anlge(p1[i - 1], p1[i], p1[(i + 1) % n]);
			ang2[i] = get_anlge(p2[i - 1], p2[i], p2[(i + 1) % n]);
			dir1[i] = get_dir(p1[i - 1], p1[i], p1[(i + 1) % n]);
			dir2[i] = get_dir(p2[i - 1], p2[i], p2[(i + 1) % n]);
		}
		max1 = -1, max2 = -1;
		for(int i = 0;i < n;i++)
		{
			len1[i] = dis2(p1[i], p1[(i + 1) % n]);
			if(len1[i] > max1) max1 = len1[i];
			len2[i] = dis2(p2[i], p2[(i + 1) % n]);
			if(len2[i] > max2) max2 = len2[i];
		}
		for(int i = 0;i < n;i++)
		{
			len1[i] /= max1;
			len2[i] /= max2;
		}
		if(slove(ang1, ang2, len1, len2, dir1, dir2, n)) printf("similar\n");
		else printf("dissimilar\n");		
	}
	return 0;
}

#include <stdio.h>
#include <math.h>

typedef struct {
  int x, y;
} Point;

/* counterclockwise, clockwise, or undefined */
enum {CCW, CW, CNEITHER};

int ccw(Point a, Point b, Point c)
{
  int dx1 = b.x - a.x;
  int dx2 = c.x - b.x;
  int dy1 = b.y - a.y;
  int dy2 = c.y - b.y;
  int t1 = dy2 * dx1;
  int t2 = dy1 * dx2;

  if (t1 == t2) {
    if (dx1 * dx2 < 0 || dy1 * dy2 < 0) {
      if (dx1*dx1 + dy1*dy1 >= dx2*dx2 + dy2*dy2) {
        return CNEITHER;
      } else {
        return CW;
      }
    } else {
      return CCW;
    }
  } else if (t1 > t2) {
    return CCW;
  } else {
    return CW;
  }
}

void read_poly(int n, Point *poly)
{
  int i;
  
  for (i = 0; i < n; i++) {
    scanf("%d %d", &(poly[i].x), &(poly[i].y));
  }
  poly[n] = poly[0];
  poly[n+1] = poly[1];
}

/* Note: we are really computing the cosine of the angle */
void compute(int n, Point *poly, double *angle, int *orient, double *dist)
{
  int dx1, dy1, dx2, dy2;
  int i;
  
  for (i = 0; i < n; i++) {
    /* i-th angle is angle formed by points i, i+1, i+2 */
    dx1 = poly[i].x - poly[i+1].x;  dy1 = poly[i].y - poly[i+1].y;
    dx2 = poly[i+2].x - poly[i+1].x; dy2 = poly[i+2].y - poly[i+1].y;
    angle[i] = (dx1*dx2+dy1*dy2)/sqrt(dx1*dx1+dy1*dy1)/sqrt(dx2*dx2+dy2*dy2);
    orient[i] = ccw(poly[i], poly[i+1], poly[i+2]);
    dist[i] = sqrt((poly[i+1].x-poly[i].x)*(poly[i+1].x-poly[i].x) +
      (poly[i+1].y-poly[i].y)*(poly[i+1].y-poly[i].y));
  }
}

int main(void)
{
  Point poly1[12], poly2[12];
  double angle1[10], angle2[10];
  int orient1[10], orient2[10];
  double dist1[10], dist2[10], ratio;
  int n, i;
  int similar;

  while (scanf("%d", &n) != EOF && n > 0) {
    read_poly(n, poly1);
    read_poly(n, poly2);

    compute(n, poly1, angle1, orient1, dist1);
    compute(n, poly2, angle2, orient2, dist2);
    ratio = dist2[0]/dist1[0];

    similar = 1;
    for (i = 0; i < n && similar; i++) {
      similar = (orient1[i] == orient2[i] && 
		 fabs(angle1[i]-angle2[i]) < 1e-8 &&
		 fabs(dist1[i]*ratio-dist2[i]) < 1e-8);
      
    }

    printf("%s\n", (similar) ? "similar" : "dissimilar");
  }
  
  return 0;
}

[计算几何]最短路径+角度旋转(pku3072)

http://acm.pku.edu.cn/JudgeOnline/problem?id=3072 
注意细节,没有啥好说的
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>

#define pi acos(-1.0)
#define eps 1e-6
#define inf 0x7fffffff
#define sqr(x) ((x) * (x))
#define dis2(a, b) sqrt(sqr((double)a.x - b.x) + sqr((double)a.y - b.y))

int f[30];
double ans;

struct point
{
	int x, y;
	void input()
	{
		scanf("%d%d", &x, &y);
	}
};
point p[30];

double get_dis(const double ang, point p1, point p2, const int R, double &new_ang)
{
	double d = dis2(p1, p2);
	if(d > R) return -1;
	double t = atan2(p2.y - p1.y, p2.x - p1.x);
	if(t < 0) t += 2 * pi;
	new_ang = t * 180 / pi;
	double a = fabs(ang - new_ang);
	if(a > 180) a = 360 - a;
	return d + a;
}

void dfs(int x, const int n, double dis, double ang, const int R)
{
	if(dis > ans) return;//这句话忘加,tle 
	if(x == n - 1) 
	{
		if(dis < ans) ans = dis;
		return ;
	}
	double d, new_ang;
	for(int i = 0;i < n;i++)
	{
		if(f[i] == 0)
		{
			d = get_dis(ang, p[x], p[i], R, new_ang);
			if(fabs(d + 1) < eps) continue;
			f[i] = 1;
			dfs(i, n, dis + d, new_ang, R);	
			f[i] = 0;
		}
	}
}

int main()
{
	int R, n;
	double d, ang;
	while(scanf("%d%d", &R, &n) != EOF)
	{
		if(R == -1 && n == -1) break;
		for(int i = 0;i < n;i++)
			p[i].input();
		memset(f, 0, sizeof(f));
		ans = 1e250;
		double t = atan2(p[n - 1].y - p[0].y, p[n - 1].x - p[0].x);
		if(t < 0) t += 2 * pi;
		ang = t * 180 / pi;
		f[0] = 1;
		dfs(0, n, 0, ang, R);
		f[0] = 0;
		if(ans > 1e200) printf("impossible\n");
		else printf("%.0lf\n", ans);	
	}
	return 0;	
}

[计算几何]旋转+三角形外接圆(pku2957)

http://acm.pku.edu.cn/JudgeOnline/problem?id=2957 
/*
将第一点旋转a1度(旋转半径为r1即p1到原点的距离)
再将第二个点旋转a2度(旋转半径为r2即p2到原点的距离) 
最后三点都在同一个圆上,这个圆即为M在p3时绕P旋转的
所在轨道(圆周运动),再通过这三点求外接圆,得到的圆心
即为在第三时刻p的位置坐标圆心到原点的距离即为p到s的距离
*/
#include <math.h>
#include <stdio.h>

#define pi acos(-1.0)
#define sqr(x) ((x) * (x))
#define triArea(a, b, c) (fabs(multi(a, b, c) / 2))
#define dis2(a, b) sqrt(sqr(a.x - b.x) + sqr(a.y - b.y))
#define multi(a, b, c) ((b.x - a.x) * (c.y - a.y) - (c.x - a.x) * (b.y - a.y))

struct point
{
	double x, y;
	void input(){scanf("%lf%lf", &x, &y);}
	double dis2_or(){return sqrt(sqr(x) + sqr(y));}
};

struct circle
{
	point center;
	double r;
};

circle circumcircleOfTriangle(point t0, point t1, point t2)
{
    //三角形的外接圆
    circle tmp;
    double a, b, c, c1, c2;
    double xA, yA, xB, yB, xC, yC;
    a = dis2(t0, t1);
    b = dis2(t1, t2);
    c = dis2(t2, t0);
    //根据S = a * b * c / R / 4;求半径R 
    tmp.r = a * b * c / triArea(t0, t1, t2) / 4;
    
    xA = t0.x; yA = t0.y;
    xB = t1.x; yB = t1.y;
    xC = t2.x; yC = t2.y;
    c1 = (xA * xA + yA * yA - xB * xB - yB * yB) / 2;
    c2 = (xA * xA + yA * yA - xC * xC - yC * yC) / 2;
    
    tmp.center.x = (c1 * (yA - yC) - c2 * (yA - yB)) / 
		((xA - xB) * (yA - yC) - (xA - xC) * (yA - yB)); 
    tmp.center.y = (c1 * (xA - xC) - c2 * (xA - xB)) / 
		((yA - yB) * (xA - xC) - (yA - yC) * (xA - xB)); 
    return tmp;     
}

int main()
{
	int T, k1, k2;
	circle ci;
	point p1, p2, p3, p1_, p2_;
	double a1, a2, r1, r2;
	while(scanf("%d%d%d", &T, &k1, &k2) != EOF)
	{
		if(T == 0 && k1 == 0 && k2 == 0) break;
		p1.input();
		p2.input();
		p3.input();
		a1 = atan2(p1.y, p1.x) + 2* pi * ((double)k1 + k2) / T;
		a2 = atan2(p2.y, p2.x) + 2 * pi * (double)k2 / T;
		r1 = p1.dis2_or();
		r2 = p2.dis2_or();
		p1_.x = r1 * cos(a1);
		p1_.y = r1 * sin(a1);
		p2_.x = r2 * cos(a2);
		p2_.y = r2 * sin(a2);
		ci = circumcircleOfTriangle(p1_, p2_, p3);
		printf("%.0lf\n", ci.center.dis2_or());		
	}
	return 0;
}

[计算几何]经度纬度的应用(pku2587,pku3407)

 

http://acm.pku.edu.cn/JudgeOnline/problem?id=2587
http://acm.pku.edu.cn/JudgeOnline/problem?id=3407

/*
Brookebond s'en va en guerre...
Description
Famous military leader marshal Brookebond who has never lost a battle
 (truth to be told, he has never won one either) is sincerely convinced 
 that all military operations can be planned on the globe. Let’s not reveal
  the depth of his misconceptions to the poor marshal. Instead, let’s help
   him by writing a program to compute the distance between two points on the
    surface of the Earth given in the geographic coordinates. An order from
	 the marshal Brookebond declares the Earth to be a perfect sphere having
	  radius of 6370 kilometers. And the orders, as we all know, are not to
	   be discussed…
Geographic latitude and longitude are measured in degrees and minutes with the
 accuracy to one minute (one degree containing 60 minutes, of course). The
  latitude is measured from 90 degrees of northern latitude (N) for the North
   Pole to 90 degrees of southern latitude (S) for the South Pole; the latitude
    of any point on the equator is 0 degrees (N or S, does not matter). The
	 longitude of is measured from 180 degrees of eastern longitude (E) to
	  180 degrees of western longitude (W); for point on meridians with 
	  longitude 180 or 0, E and W are equivalent.
Input
The first and second lines of the input contain the coordinates (latitude 
and longitude) of one point each. The designation of a latitude contains 
two integers (degrees and minutes) and a letter N or S. Similarly, the 
designation of a longitude contains two integers (degrees and minutes) 
and a letter E or W. Adjacent values on a line are separated by one or 
more spaces.
Output
The first and only line of the output must contain the distance
 (in kilometers) between the points from the input file with 
 the precision of 1 meter.
Sample Input

55 0 N 40 0 E
59 0 N 49 30 E
Sample Output

725.979

假设地球是球体, 
设地球上某点的经度为lambda,纬度为phi, 
则这点的空间坐标是 
x=cos(phi)*cos(lambda) 
y=cos(phi)*sin(lambda) 
z=sin(phi) 
设地球上两点的空间坐标分别为(x1,y1,z1),(x2,y2,z2) 
直线距离即为R*sqrt((x2-x1)*(x2-x1)+(y2-y1)*(y2-y1)+(z2-z1)*(z2-z1)),
则它们的夹角为 
A = acos(x1 * x2 + y1 * y2 + z1 * z2),
则两地距离为 A * R,其中R为地球平均半径6371 

这里坐标都要乘以半径R,但由于是求角度,所以统一都没有乘 
注意这里还要判断坐标的正负和经度纬度的规定有关 

pku_3407

#include <stdio.h>
#include <math.h>

const double pi = acos(-1.0);

struct TPoint
{
   double x, y, z;
};
 
int main()
{
    double w1, wm1, j1, jm1, wd1, wd2;
    double w2, wm2, j2, jm2, jd1, jd2;
    TPoint p1, p2;
    char chr1, chr2;
    while(scanf("%lf%lf ", &w1, &wm1) != EOF){
        scanf("%c ", &chr1);
        scanf("%lf %lf %c", &j1, &jm1, &chr2);
        wd1 = (w1 + wm1 / 60) * pi / 180;
        jd1 = (j1 + jm1 / 60) * pi / 180;
        if(chr1 == 'S') wd1 *= -1.0;
        if(chr2 == 'W') jd1 *= -1.0;
        p1.x = cos(wd1) * cos(jd1);
        p1.y = cos(wd1) * sin(jd1);
        p1.z = sin(wd1);
        scanf("%lf %lf %c %lf %lf %c", &w2, &wm2, &chr1, &j2, &jm2, &chr2);
        wd2 = (w2 + wm2 / 60) * pi / 180;
        jd2 = (j2 + jm2 / 60) * pi / 180;
        if(chr1 == 'S') wd2 *= -1.0;
        if(chr2 == 'W') jd2 *= -1.0;
        p2.x = cos(wd2) * cos(jd2);
        p2.y = cos(wd2) * sin(jd2);
        p2.z = sin(wd2);
        double a = acos(p1.x * p2.x + p1.y * p2.y + p1.z * p2.z);
        printf("%.3lf\n", a * 6370.0);   
    }
    return 0;
}
*/
//longitude  经度 
//latitude   纬度 
#include <math.h>
#include <stdio.h>

#define pi acos(-1.0)
#define sqr(a) ((a) * (a))

struct point 
{
	double x, y, z;
	double  latitude, longitude;
	double l1, l2;
	void input()
	{
		scanf("%lf%lf", &latitude, &longitude);
		l1 = latitude;
		l2 = longitude;
		latitude *= pi / 180;
		longitude *= pi / 180;
		x = cos(latitude) * cos(longitude);
        y = cos(latitude) * sin(longitude);
        z = sin(latitude);
	}
	double dis(point &b)
	{
		return acos(x * b.x + y * b.y + z * b.z);
	}
};

int main()
{
	int n, u;
	point p[1005];
	while(scanf("%d", &n) != EOF)
	{
		for(int i = 0;i < n;i++)
			p[i].input();
		u = 0;
		double mindis = 1e250, maxdis, tmp;
		for(int i = 0;i < n;i++)
		{
			maxdis = -2;
			for(int j = 0;j < n;j++)
			{
				tmp = p[i].dis(p[j]);
				if(tmp > maxdis) maxdis = tmp;
			}
			if(maxdis < mindis) 
			{
				mindis = maxdis;
				u = i;
			}
		}
		printf("%.2lf %.2lf\n", p[u].l1, p[u].l2);	
	}
	return 0;
}

[计算几何]计算几何模板下载

 

点击下载   计算几何_陈海丰.doc

(1)凸包
(2)判断两条线段是否相交(平行,不平行)
(3)三角形的外接圆(已知不在同一直线上的三点求经过三点的圆)
(4)三角形的垂心内心重心中垂线
(5)求直线的交点
(6)根据线段两端点的坐标求垂直平分线上除中点外的另一点
(7)根据两点坐标求直线方程
(8)差积的应用
(9)三角形的面积公式
(10)三角形的内接圆(未检验正确性)
(11)多边形的面积(适合凹多边形)
(12)判断点是否在线段上
(13)平面上两个点之间的距离
(14)p点关于直线L的对称点
(15)判断一个矩形是否在另一个矩形中
(16)直线和圆的交点+点关于线的对称点+点到线的距离+直线方程(fzu_1035)
(17)判断点是否在多边形内
(18)N点中三个点组成三角形面积最大
(19)扇形的重心
(20)多边形的重心
(21)判断N点是否共面
(22)求共线的点最多为多少
(23)N个矩形的相交的面积
(24)三角形外接圆+圆的参数方程(pku_1266)
(25)判断线段是否有交点并求交点(cug_1078)
(26)简单多边形的核
(27)线段重叠+投影(pku_1375)
(28)二分+圆的参数方程(pku_2600)
(29)Pick公式
(30)根据经度纬度求球面距离
(31)两圆切线的交点
(32)线段与三角形的交(usaco_3.4)

[C/C++函数]atan, atan2

atan, atan2

Calculates the arctangent of x (atan) or the arctangent of y/x (atan2).

double atan( double x );

double atan2( double y, double x );

Routine

Required Header

Compatibility

atan

<math.h>

ANSI, Win 95, Win NT

atan2

<math.h>

ANSI, Win 95, Win NT

 

For additional compatibility information, see Compatibility in the Introduction.

Libraries

LIBC.LIB

Single thread static library, retail version

LIBCMT.LIB

Multithread static library, retail version

MSVCRT.LIB

Import library for MSVCRT.DLL, retail version

 

Return Value

atan returns the arctangent of x. atan2 returns the arctangent of y/x. If x is 0, atan returns 0. If both parameters of atan2 are 0, the function returns 0. You can modify error handling by using the _matherr routine. atan returns a value in the range –π/2 to π/2 radians; atan2 returns a value in the range –π to π radians, using the signs of both parameters to determine the quadrant of the return value.

Parameters

x, y

Any numbers

Remarks

The atan function calculates the arctangent of x. atan2 calculates the arctangent of y/x. atan2 is well defined for every point other than the origin, even if x equals 0 and y does not equal 0.

Example

/* ATAN.C: This program calculates 
 * the arctangent of 1 and -1.
 */


void main( void ) 
{
   double x1, x2, y;

   printf( "Enter a real number: " );
   scanf( "%lf", &x1 );
   y = atan( x1 );
   printf( "Arctangent of %f: %f\n", x1, y );
   printf( "Enter a second real number: " );
   scanf( "%lf", &x2 );
   y = atan2( x1, x2 );
   printf( "Arctangent of %f / %f: %f\n", x1, x2, y );

Output

Enter a real number: -862.42
Arctangent of -862.420000: -1.569637
Enter a second real number: 78.5149
Arctangent of -862.420000 / 78.514900: -1.480006

Floating-Point Support Routines

See Also   acos, asin, cos, _matherr, sin, tan

 

[计算几何]几何公式

几何公式

【三角形】:

1. 半周长 P=(a+b+c)/2

2. 面积 S=aHa/2=absin(C)/2=sqrt(P(P-a)(P-b)(P-c))

3. 中线 Ma=sqrt(2(b^2+c^2)-a^2)/2=sqrt(b^2+c^2+2bccos(A))/2

4. 角平分线 Ta=sqrt(bc((b+c)^2-a^2))/(b+c)=2bccos(A/2)/(b+c)

5. 高线 Ha=bsin(C)=csin(B)=sqrt(b^2-((a^2+b^2-c^2)/(2a))^2)

6. 内切圆半径 r=S/P=asin(B/2)sin(C/2)/sin((B+C)/2)

               =4Rsin(A/2)sin(B/2)sin(C/2)=sqrt((P-a)(P-b)(P-c)/P)

               =Ptan(A/2)tan(B/2)tan(C/2)

7. 外接圆半径 R=abc/(4S)=a/(2sin(A))=b/(2sin(B))=c/(2sin(C))

【四边形】:

D1,D2为对角线,M对角线中点连线,A为对角线夹角

1. a^2+b^2+c^2+d^2=D1^2+D2^2+4M^2

2. S=D1D2sin(A)/2

(以下对圆的内接四边形)

3. ac+bd=D1D2

4. S=sqrt((P-a)(P-b)(P-c)(P-d)),P为半周长

【正n边形】:

R为外接圆半径,r为内切圆半径

1. 中心角 A=2PI/n

2. 内角 C=(n-2)PI/n

3. 边长 a=2sqrt(R^2-r^2)=2Rsin(A/2)=2rtan(A/2)

4. 面积 S=nar/2=nr^2tan(A/2)=nR^2sin(A)/2=na^2/(4tan(A/2))

【圆】:

1. 弧长 L=rA

2. 弦长 a=2sqrt(2hr-h^2)=2rsin(A/2)

3. 弓形高 h=r-sqrt(r^2-a^2/4)=r(1-cos(A/2))=atan(A/4)/2

4. 扇形面积 S1=rl/2=r^2A/2

5. 弓形面积 S2=(rl-a(r-h))/2=r^2(A-sin(A))/2

【棱柱】:

1. 体积 V=Ah,A为底面积,h为高

2. 侧面积 S=lp,l为棱长,p为直截面周长

3. 全面积 T=S+2A

【棱锥】:

1. 体积 V=Ah/3,A为底面积,h为高

(以下对正棱锥)

2. 侧面积 S=lp/2,l为斜高,p为底面周长

3. 全面积 T=S+A

【棱台】:

1. 体积 V=(A1+A2+sqrt(A1A2))h/3,A1.A2为上下底面积,h为高

(以下为正棱台)

2. 侧面积 S=(p1+p2)L/2,p1.p2为上下底面周长,l为斜高

3. 全面积 T=S+A1+A2

【圆柱】:

1. 侧面积 S=2PIrh

2. 全面积 T=2PIr(h+r)

3. 体积 V=PIr^2h

【圆锥】:

1. 母线 L=sqrt(h^2+r^2)

2. 侧面积 S=PIrl

3. 全面积 T=PIr(L+r)

4. 体积 V=PIr^2h/3

【圆台】:

1. 母线 L=sqrt(h^2+(r1-r2)^2)

2. 侧面积 S=PI(r1+r2)L

3. 全面积 T=PIr1(L+r1)+PIr2(L+r2)

4. 体积 V=PI(r1^2+r2^2+r1r2)h/3

【球】:

1. 全面积 T=4PIr^2

2. 体积 V=4PIr^3/3

【球台】:

1. 侧面积 S=2PIrh

2. 全面积 T=PI(2rh+r1^2+r2^2)

3. 体积 V=PIh(3(r1^2+r2^2)+h^2)/6

【球扇形】:

1. 全面积 T=PIr(2h+r0),h为球冠高,r0为球冠底面半径

2. 体积 V=2PIr^2h/3

clip_image002Euler的任意四面体体积公式(已知边长求体积)

clip_image004  

已知4点坐标求体积(其中四个点的坐标分别为(x1,y1,z1,x2,y2,z2,x3,y3,z3,x4,y4,z4clip_image006

 

 

 

注意事项:

1. 注意舍入方式(0.5的舍入方向);防止输出-0.

2. 几何题注意多测试不对称数据.

3. 整数几何注意xmultdmult是否会出界;

   符点几何注意eps的使用.

4. 避免使用斜率;注意除数是否会为0.

5. 公式一定要化简后再代入.

6. 判断同一个2*PI域内两角度差应该是

   abs(a1-a2)<beta||abs(a1-a2)>pi+pi-beta;

   相等应该是

   abs(a1-a2)<eps||abs(a1-a2)>pi+pi-eps;

7. 需要的话尽量使用atan2,注意:atan2(0,0)=0,

   atan2(1,0)=pi/2,atan2(-1,0)=-pi/2,atan2(0,1)=0,atan2(0,-1)=pi.

8. cross product = |u|*|v|*sin(a)

   dot product = |u|*|v|*cos(a)

9. (P1-P0)x(P2-P0)结果的意义:

   : <P0,P1><P0,P2>顺时针(0,pi)

   : <P0,P1><P0,P2>逆时针(0,pi)

   0 : <P0,P1>,<P0,P2>共线,夹角为0pi

 

 

[计算几何]长方体上一点到原点的距离(pku2977)

http://acm.pku.edu.cn/JudgeOnline/problem?id=2977

Box walking

Description

You are given a three-dimensional box of integral dimensions lx × ly × lz The edges of the box are axis-aligned, and one corner of the box is located at position (0, 0, 0). Given the coordinates (x, y, z) of some arbitrary position on the surface of the box, your goal is to return the square of the length of the shortest path along the box’s surface from (0, 0, 0) to (x, y, z).

If lx = 1, ly = 2, and lz = 1, then the shortest path from (0, 0, 0) to (1, 2, 1) is found by walking from (0, 0, 0) to (1, 1, 0), followed by walking from (1, 1, 0) to (1, 2, 1). The total path length is √8.

Input

The input test file will contain multiple test cases, each of which consists of six integers lx, ly, lz, x, y, z where 1 ≤ lx, ly, lz ≤ 1000. Note that the box may have zero volume, but the point (x, y, z) is always guaranteed to be on one of the six sides of the box. The end-of-file is marked by a test case with lx = ly = lz = x = y = z = 0 and should not be processed.

Output

For each test case, write a single line with a positive integer indicating the square of the shortest path length. (Note: The square of the path length is always an integer; thus, your answer is exact, not approximate.)

Sample Input

1 1 2 1 1 2
1 1 1 1 1 1
0 0 0 0 0 0

Sample Output

8
5

Source

Stanford Local 2005

/*
Box
---

Figuring out the shortest path on a box may seem trivial at first.
All we need to do is find the unfolding of the box that brings
the destination point closest.
If the final point is on one of the three faces touching the origin,
the shortest path is always a straight line directly to the destination;
if it's on one of the other three faces, it may seem obvious that we can
start on a face adjacent to the final face, and then just go around a
single corner.  It may seem obvious, that is, but it is wrong.

Consider a box that has a width of 100 units, and a height and depth
of 4 units, and let's say our destination is (100,3,3).  If we
use the face perpendicular to the depth and adjacent with the origin,
and then the destination face, we must traverse 103 units horizontally
and 3 units vertically for a total distance of 103 and some change:

   |--------~...~------------|----|
   |                         |   *|
   |                         |    |
   |                         |    |
   *--------~...~------------|----|

If on the other hand we use both the face perpendicular to the
depth and the one on top of that we can cut this down to 101
units horizontally and 7 units vertically for a total distance
of only 101 and some change:

   |--------~...~------------|----|
   |                         |*   |
   |                         |    |
   |                         |    |
   |--------~...~------------|----|
   |                         |
   |                         |
   |                         |
   *--------~...~------------|

So sometimes we will need to traverse three faces to get to
our final point the quickest way.  Do we ever need to traverse
four or more faces?  In this problem, probably not, but it
turns out we don't need to know that; we can easily enumerate
all single, two-face, three-face, four-face, and even more
using recursion.

Indeed, manually enumerating and coding all the cases is
terribly error-prone and very difficult to test; a single sign
error or coordinate swap can doom your submission.  So writing
the general case is much safer.

We'll start with a routine that enumerates the unfoldings
starting with a single face.  If the final destination is on
that face, the shortest path is always directly to that
point (every subpath of the final shortest path is itself a
straight line on the unfolded box).  Otherwise, we need to
traverse either over the right edge, or over the top edge,
to the next face.  As we cross each edge, we need to
accumulate how much horizontal and vertical distance we have
covered.

long best = Integer.MAX_VALUE ;
long sqadd(int x, int y) {
   return x * (long)x + y * (long)y ;
}
void rec(
   int x, int y,        // the distance we've traversed so far
   int w, int h, int d, // the width and height of the current face, and
                        // depth of the box in this orientation
   int dx, int dy, int dz, // the position of the destination
   int togo) {          // how many additional faces we allow ourselves
   if (dz == 0) { // the destination is on this face
      best = Math.min(best, sqadd(x+dx, y+dy)) ;
      return ;
   }
   if (togo == 0) return ; // no more faces allowed
   // now we consider going to the right edge; this rotates the box and point
   rec(x + w, y, d, h, w, dz, dy, w-dx, togo-1) ;
   // now we consider going to the top edge; this rotates the box and point
   rec(x, y + h, w, d, h, dx, dz, h-dy, togo-1) ;
}

Finally, in the main routine, we need to consider all three
possible starting faces.  We don't know how many faces a best path
might be, so we be generous and say 6.  This is only a maximum of
2^6 * 3 cases, so it runs very fast.

for (each case) {
   best = Integer.MAX_VALUE ;
   rec(0, 0, w, h, d, x, y, z, 6) ;
   rec(0, 0, h, d, w, y, z, x, 6) ;
   rec(0, 0, d, w, h, z, x, y, 6) ;
}

There's one other concern we may have.  We are considering all
unfoldings, even the unfoldings for which the straight line from
the origin to the (unfolded) destination goes off the actual
unfolded box.  In such a case, the path will always be longer
than the final best result, so considering it doesn't hurt
anything.  Essentially, all the space outside the unfolding
represents the edges of the box "opened up" or "stretched";
by opening those edges that way, you can only increase the
overall path length.  Consider the picture below:

         |----|
         |*   |
         /    |
        /|    |
   |---/-|----|
   |  /  |    |
   | /   |    |
   |/    |    |
   *-----|----|

Our straight line leaves the leftmost rectangle and reenters
the topmost rectangle.  But the two edges so represented
are actually joined in the real box, so a different unfolding
for which that line does not leave the unfolding always will
have a shorter distance:

   |-----|
   |     |
   |     |
   |*    |
   |-----|
   |     |
   |     |
   |     |
   *-----|

But if you were worried about this not being the case, it is
straightforward to extend the recursion to maintain a minimum
and maximum slope (represented by rationals) permitted by this
straight-line traversals, and thus eliminate all unfoldings for
which a straight line would leave the boundaries.
*/
#include <stdio.h>

#define sqr(x) ((x) * (x))
#define min(a, b) ((a) < (b) ? (a) : (b))

int length(int lx, int ly, int lz, int x, int y, int z)
{
	int d = sqr(lx + y) + sqr(z);
	d = min(d, sqr(lx + z) + sqr(y));
	d = min(d, sqr(lx + lz - z) + sqr(y + lz));
	return min(d, sqr(lx + ly - y) + sqr(z + ly));
}

int slove(int lx, int ly, int lz, int x, int y, int z)
{
	if(x == 0 || y == 0 || z == 0) return sqr(x) + sqr(y) + sqr(z);
	else if(x == lx) return length(lx, ly, lz, x, y, z);
	else if(y == ly) return length(ly, lx, lz, y, x, z);
	else if(z == lz) return length(lz, lx, ly, z, x, y);
}

int main()
{
	freopen("box.in", "r", stdin);
	freopen("out.txt", "w", stdout);
	int lx, ly, lz, x, y, z;
	while(scanf("%d%d%d%d%d%d", &lx, &ly, &lz, &x, &y, &z) != EOF)
	{
		if(!lx && !ly && !lz && !x && !y && !z) break;
		printf("%d\n", slove(lx, ly, lz, x, y, z));
	}
	return 0;
}

[计算几何]凸包的相关算法

R_Animation otating C_Animation alipers homepage

Some history:
In 1978, M.I. Shamos's Ph.D. thesis "Computational Geometry" marks the "birth" of the field within Computer Science. Among the results he presents is a very simple algorithm for finding the diameter of a convex polygon, i.e. the greatest distance determined a pair of points belonging to the polygon.
The diameter turns out to be determined by an anti-podal pair of points. Shamos provides a simple O(n) time algorithm for determining the anti-podal pairs of a convex n-gon. Since there are at most 3n/2 of these, the diameter can then be found in O(n) time.
As Toussaint later put it, Shamos' algorithm resembles rotating a pair of calipers around the polygon. Thus the term "Rotating Calipers". In 1983, Toussaint presents a paper in which a variety of problems are solved using the same technique. Since then, new algorithms based on the paradigm have been obtained, solving a multitude of problems.
They include:
Thesis
My Master's Thesis, Computational Geometry with the Rotating Calipers (in Postscript format) gathers the above problems, while providing details and proofs. The bibliography is available here for viewing.
Applet
To use the Rotating Calipers and solve many of the above problems yourself, go to the Rotating Calipers applet page.

Back to my homepage

November 26th, 1999


求凸包的算法如下-Graham扫描算法:(可参见算法导论)
GRAHAM-SCAN(Q)
1 let p0 be the point in Q with the minimum y-coordinate,
             or the leftmost such point in case of a tie
2 let 〈p1, p2, ..., pm〉 be the remaining points in Q,
             sorted by polar angle in counterclockwise order around p0
             (if more than one point has the same angle, remove all but
             the one that is farthest from p0)
3 PUSH(p0, S)
4 PUSH(p1, S)
5 PUSH(p2, S)
6 for i ← 3 to m
7       do while the angle formed by points NEXT-TO-TOP(S), TOP(S),
                     and pi makes a nonleft turn
8              do POP(S)
9          PUSH(pi, S)
10 return S

//测试http://lab.cug.edu.cn:8080/COJ/prepare.do?fun=viewProblem&pid=1038

#include <math.h>
#include <stdio.h>
#include <stdlib.h>

#define sqr(a) ((a) * (a))
#define dis2(a, b) sqrt(sqr(a.x - b.x) + sqr(a.y - b.y))
#define multi(a, b, c) ((b.x - a.x) * (c.y - a.y) - (b.y - a.y) * (c.x - a.x))

struct point
{
	int x, y;
	int operator<(point &b)
	{
		return ((y < b.y) || (y == b.y && x < b.x));
	}
}ptmp;

int cmp(const void *a, const void *b)
{
	point *c = (point *)a;
	point *d = (point *)b;
	point e, f;
	e = *c;
	f = *d;
	int t = multi(ptmp, e, f);
	if(t > 0) return -1;
	if(t == 0 && dis2(ptmp, e) <= dis2(ptmp, f)) return -1;
	return 1;
}

void GrahamScanGram(point p[], int stack[], int &top, int n)
{
	int u = 0, i;
	for(int i = 0;i < n;i++)
		if(p[i] < p[u]) u = i;
	ptmp = p[0], p[0] = p[u], p[u] = ptmp;
	ptmp = p[0];
	qsort(p + 1, n - 1, sizeof(p[0]), cmp);
	stack[0] = 0, stack[1] = 1, stack[2] = 2;
	top = 3;
	for(i = 3;i < n;i++)	
	{
		while(top >= 2 && multi(p[stack[top - 2]], p[stack[top - 1]], p[i]) <= 0) top--;
		stack[top++] = i;	
	}
}

int main()
{
	point p[100005];
	int stack[100005], top;
	int n, ca, i;
	scanf("%d", &ca);
	while(ca--)
	{
		scanf("%d", &n);
		for(i = 0;i < n;i++)
			scanf("%d%d", &p[i].x, &p[i].y);
		GrahamScanGram(p, stack, top, n);
		int area = 0;
		for(i = 1;i < top - 1;i++)
			area += abs(multi(p[stack[0]], p[stack[i]], p[stack[i + 1]]));
		printf("%.1lf\n", (double)area / 2);
	}
	return 0;	
}
计算几何学习网站:

http://chenhaifeng.blog.edu.cn/
Rotating Calipers http://cgm.cs.mcgill.ca/~orm/rotcal.frame.html
Computational Geometry http://dev.gameres.com/Program/Abstract/Geometry.htm
http://www.cs.uu.nl/geobook/
http://cgm.cs.mcgill.ca/~godfried/teaching/cg-web.html
http://learn.tsinghua.edu.cn:8080/2001315450/cg.html
Computational Geometry Code  http://compgeom.cs.uiuc.edu/~jeffe/compgeom/code.html

[算法分析]去掉K个分数,使分子和除以分母和最大(pku2976)

http://acm.pku.edu.cn/JudgeOnline/problem?id=2976
2 
/*
Drop
----

This problem was intended to be the hardest problem in the set, and indeed, nobody
solved it. When choosing to drop tests, there is a trade-off between tests with very
low percentage scores and tests with mediocre percentage scores but with a large number
of questions. Which one of these is better can depend on the other tests you are
keeping. Figuring out how to balance these requirements is the heart of this problem.

The first thing you have to realize is that a simple greedy algorithm is not going
to be correct. For example, it seems reasonable to drop tests one at a time, always
choosing to drop the test that will improve your cumulative average by as much as
possible. However, it turns out this is wrong. For example, consider dropping 2 tests
from a group of 100/100, 14/50, 14/50, 100/200, and 100/200. The greedy approach would
first drop a 14/50 test and then a 100/200 test to get a score of 214/350 = 61%, but
the optimal strategy is to drop both 100/200 tests for a score of 128/200 = 64%.

If you have done a lot of other programming competitions, you might also think of the
dynamic programming approach (see the stones write-up). By tracking the number of
questions and the number of tests, you can use this to solve the problem in
(number of questions)*(number of tests) time. Unfortunately, the maximum number of
questions over all the tests is a trillion, so that is no good either.

What do you do then? It turns out that the problem can be solved with a binary search.
We ask the following question: "Can you drop k tests to make your cumulative average
at least x?". It turns out that fixing x makes the problem substantially easier
because this is enough to determine which tests are better than others.

If we fix x, we need to choose tests so that
    (sum a_i) / (sum b_i) >= x
       <=> sum a_i >= sum (x*b_i)
       <=> sum (a_i - x*b_i) >= 0.
Thus, we compute c_i = a_i - x*b_i for each i. We now need to drop k of those values so
that their sum is at least 0. This is a much easier problem! Just sort the c_i's and
drop the k smallest values.

This reduces everything to a standard binary search problem. For each x, we can test
whether we can get an average of at least x, and we need to find the maximum average
that can be made. C++ code is shown bleow.

    double lb = 0, ub = 1;
    for (i = 0; i < 100; i++) {
      double x = (lb+ub)/2;

      for (j = 0; j < tests; j++)
	scores[j] = num[j] - x*den[j];
      sort(scores, scores+tests);
      double total = 0;
      for (j = k; j < tests; j++)
	total += scores[j];

      if (total >= 0)
	lb = x;
      else
	ub = x;
    }
    cout << int(100*lb + 0.5) << endl;
*/
#include <stdio.h>
#include <iostream>
#include <algorithm>

using namespace std;

#define eps 1e-6

int main()
{
	freopen("drop.in", "r", stdin);
	freopen("out.txt", "w", stdout);
    int a[1005], b[1005];
    int n, k;
    double d[1005], s;
    double l, r, mid;
    while(scanf("%d%d", &n, &k) != EOF)
    {
		if(n == k && n == 0) break;
		for(int i = 0;i < n;i++)
			scanf("%d", &a[i]);
		for(int i = 0;i < n;i++)
			scanf("%d", &b[i]);
		l = 0;
		r = 100;
		while(r > l + eps)
		{
			mid = (r + l) / 2;
			for(int i = 0;i < n;i++)
				d[i] = a[i] - mid * b[i];
			sort(d, d + n);
			s = 0;
			for(int i = k;i < n;i++)
				s += d[i];
			if(s >= 0) l = mid;
			else r = mid;
		}
		printf("%d\n", int(l * 100 + 0.5));
    }
    return 0;
}

[计算几何]凸多面体的表面积(点积叉积的应用)

http://acm.pku.edu.cn/JudgeOnline/problem?id=2974 

未命名

Source Code

Problem: 2974   User: chenhaifeng
Memory: 312K   Time: 32MS
Language: G++   Result: Accepted
  • Source Code
    //测试  http://acm.pku.edu.cn/JudgeOnline/problem?id=2974
    /***************************************************************************
    Area
    ----
    
    This problem was a test of your ability to do vector math.  Given a set of points in R^3,
    your goal is to find the surface area of the 3D convex hull.  You are guaranteed that
    all faces of the convex hull contain at most three point from your set.  At first glance,
    the problem sounds extremely difficult because it sounds like you need a really complex
    three-dimensional convex hull algorithm.
    
    However, the maximum allowed number of points in the set is only 25, so it turns out the
    following fairly simple solution suffices:
    
    1) Consider all possible n choose 3 subset