计算几何模板
(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)
(33)求长方体上到原点距离的点
(34)Computational_Geometry
(1)凸包
/*凸包cug_1038*/
#include <stdio.h>
#include <stdlib.h>
struct point
{
int x, y;
}pp;
point p[100005];
int stack[100005], top;
int dis(point a, point b)
{
return ((a.x - b.x) * (a.x - b.x) + (a.y - b.y) * (a.y - b.y));
}
int multi(point b, point c, point a)
{
return (b.x - a.x) * (c.y - a.y) - (b.y - a.y) * (c.x - a.x);
}
void swap(point p[], int s, int t)
{
point tmp;
tmp = p[s];
p[s] = p[t];
p[t] = tmp;
}
int cmp(const void *a, const void *b)
{
point *c = (point *)a;
point *d = (point *)b;
double k = multi(*c, *d, pp);
if(k < 0) return 1;
else if(k == 0 && dis(*c, pp) >= dis(*d, pp)) return 1;
else return -1;
}
void Graham(point p[], int n, int stack[], int &top)
{
int i, u;
u = 0;
for(i = 1;i < n;i++){
if(p[i].y == p[u].y && p[i].x < p[u].x) u = i;
else if(p[i].y < p[u].y) u = i;
}
swap(p, 0, u);
pp = p[0];
qsort(p + 1, n - 1, sizeof(p[0]), cmp);
stack[0] = 0;
stack[1] = 1;
top = 1;
for(i = 2;i < n;i++){
while(multi(p[i], p[stack[top]], p[stack[top - 1]]) >= 0){
if(top == 0) break;
top--;
}
top++;
stack[top] = i;
}
}
int main()
{
int ca, i, j, n;
int area;
scanf("%d", &ca);
for(i = 1;i <= ca;i++){
scanf("%d", &n);
for(j = 0;j < n;j++){
scanf("%d%d", &p[j].x, &p[j].y);
}
Graham(p, n, stack, top);
area = 0;
for(j = 1;j <= top - 1;j++){
area += abs(multi(p[stack[0]], p[stack[j]], p[stack[j + 1]]));
}
printf("%.1lf\n", (double)area / 2);
}
return 0;
}
---------------------------------------------------------------------------------------------------------------------
(2)判断两条线段是否相交(平行，不平行)
bool isIntersected(TPoint s1, TPoint e1, TPoint s2, TPoint e2)
{
//判断线段是否相交
//1.快速排斥试验判断以两条线段为对角线的两个矩形是否相交
//2.跨立试验
if(
(max(s1.x, e1.x) >= min(s2.x, e2.x)) &&
(max(s2.x, e2.x) >= min(s1.x, e1.x)) &&
(max(s1.y, e1.y) >= min(s2.y, e2.y)) &&
(max(s2.y, e2.y) >= min(s1.y, e1.y)) &&
(multi(s2, e1, s1) * multi(e1, e2, s1) >= 0) &&
(multi(s1, e2, s2) * multi(e2, e1, s2) >= 0)
) return true;
return false;
}
(3)三角形的外接圆(已知不在同一直线上的三点求经过三点的圆)
/*三角形的外接圆pku_1329*/
#include <stdio.h>
#include <math.h>
const double eps = 1e-6;
typedef struct TPoint
{
double x;
double y;
}TPoint;
typedef struct TTriangle
{
TPoint t[3];
}TTriangle;
typedef struct TCircle
{
TPoint centre;
double r;
}TCircle;
double distance(TPoint p1, TPoint p2)
{
//计算平面上两个点之间的距离
return sqrt((p1.x - p2.x) * (p1.x - p2.x) + (p1.y - p2.y) * (p1.y - p2.y));
}
double triangleArea(TTriangle t)
{
//已知三角形三个顶点的坐标，求三角形的面积
return fabs(t.t[0].x * t.t[1].y + t.t[1].x * t.t[2].y + t.t[2].x * t.t[0].y
- t.t[1].x * t.t[0].y - t.t[2].x * t.t[1].y - t.t[0].x * t.t[2].y) / 2;
}
TCircle circumcircleOfTriangle(TTriangle t)
{
//三角形的外接圆
TCircle tmp;
double a, b, c, c1, c2;
double xA, yA, xB, yB, xC, yC;
a = distance(t.t[0], t.t[1]);
b = distance(t.t[1], t.t[2]);
c = distance(t.t[2], t.t[0]);
//根据S = a * b * c / R / 4;求半径R
tmp.r = a * b * c / triangleArea(t) / 4;
xA = t.t[0].x; yA = t.t[0].y;
xB = t.t[1].x; yB = t.t[1].y;
xC = t.t[2].x; yC = t.t[2].y;
c1 = (xA * xA + yA * yA - xB * xB - yB * yB) / 2;
c2 = (xA * xA + yA * yA - xC * xC - yC * yC) / 2;
tmp.centre.x = -(c1 * (yA - yC) - c2 * (yA - yB)) /
((xA - xB) * (yA - yC) - (xA - xC) * (yA - yB));
tmp.centre.y = -(c1 * (xA - xC) - c2 * (xA - xB)) /
((yA - yB) * (xA - xC) - (yA - yC) * (xA - xB));
return tmp;
}
int main()
{
TTriangle t;
TCircle circle;
double c, d, e;
while(scanf("%lf%lf%lf%lf%lf%lf", &t.t[0].x, &t.t[0].y,
&t.t[1].x, &t.t[1].y, &t.t[2].x, &t.t[2].y) != EOF){
circle = circumcircleOfTriangle(t);
// printf("%lf %lf %lf\n", circle.centre.x, circle.centre.y, circle.r);
if(fabs(circle.centre.x) < eps) printf("x^2");
else if(circle.centre.x < 0) printf("(x - %.3lf)^2 + ", -circle.centre.x);
else printf("(x + %.3lf)^2 + ", circle.centre.x);
if(fabs(circle.centre.y) < eps) printf("y^2 = ");
else if(circle.centre.y < 0) printf("(y - %.3lf)^2 = ", -circle.centre.y);
else printf("(y + %.3lf)^2 = ", circle.centre.y);
printf("%.3lf^2\n", circle.r);
c = 2 * circle.centre.x;
d = 2 * circle.centre.y;
e = circle.centre.x * circle.centre.x +
circle.centre.y * circle.centre.y - circle.r * circle.r;
printf("x^2 + y^2 ");
//if(fabs(c) < eps)
if(c < 0) printf("- %.3lfx ", -c);
else printf("+ %.3lfx ", c);
if(d < 0) printf("- %.3lfy ", -d);
else printf("+ %.3lfy ", d);
if(e < 0) printf("- %.3lf = 0\n", -e);
else printf("+ %.3lf = 0\n", e);
printf("\n");
}
return 0;
}
(4)三角形的垂心内心重心中垂线
/*cug_1011_垂心内心重心中垂线.cpp*/
#include <iostream>
#include <cmath>
using namespace std;
const double eps = 1e-6;
struct point
{
double x, y;
};
void K()
{
//到三边距离和最短
}
void L(double a, double b, double c, double A, double B, double C)
{ //垂线的交点
double t1, t2, t3;
t1 = c * cos(A) / cos(M_PI / 2 - C);
t2 = c * cos(B) / cos(M_PI / 2 - C);
t3 = a * cos(C) / cos(M_PI / 2 - A);
t1 += t2 + t3;
printf("%.3lf\n", t1);
}
struct TLine
{
double a, b, c;
};
TLine lineFromSegment(point p1, point p2)
{
//线段所在直线,返回直线方程的三个系统
TLine tmp;
tmp.a = p2.y - p1.y;
tmp.b = p1.x - p2.x;
tmp.c = p2.x * p1.y - p1.x * p2.y;
return tmp;
}
point LineInter(TLine l1, TLine l2)
{
//求两直线得交点坐标
point tmp;
if(fabs(l1.b) < eps){
tmp.x = -l1.c / l1.a;
tmp.y = (-l2.c - l2.a * tmp.x) / l2.b;
}
else{
tmp.x = (l1.c * l2.b - l1.b * l2.c) / (l1.b * l2.a - l2.b * l1.a);
tmp.y = (-l1.c - l1.a * tmp.x) / l1.b;
}
return tmp;
}
double dis(point a, point b)
{
return sqrt((a.x - b.x) * (a.x - b.x) + (a.y - b.y) * (a.y - b.y));
}
void F(double a, double b, double c, double A, double B, double C)
{
//到三顶点的距离和最短, 费马点
/*当三角形最大的顶角小于120度的时候，三角形内一点到
三顶点之间的距离最小是与三顶点夹角都成120度的点P当
最到顶点大于等于120度，该顶点取最小值
补充一下，当三角形的最大角小于120度时，费尔码点在三
角形内，作法有多种，可以从任二办向外作等边三角形，联
接正三角形的顶点和原三角形的对角，两者的联线即所求。
当三角形的最大角大于等于120度时，
费尔码点在三角形的钝角上。*/
if(A - 2 * M_PI / 3 > -eps){
printf("%.3lf ", b + c);
return ;
}
else if(B - 2 * M_PI / 3 > -eps){
printf("%.3lf ", a + c);
return ;
}
else if(C - 2 * M_PI / 3 > -eps){
printf("%.3lf ", a + b);
return ;
}
point pa, pb, pc, pc1, pa1;
pa.x = 0, pa.y = 0;
pb.x = c, pb.y = 0;
pc.x = b * cos(A);
pc.y = b * sin(A);
pc1.x = c * cos(-M_PI / 3);
pc1.y = c * sin(-M_PI / 3);
pa1.x = a * cos(2 * M_PI / 3 - B) + c;
pa1.y = a * sin(2 * M_PI / 3 - B);
TLine l1, l2;
l1 = lineFromSegment(pa, pa1);
l2 = lineFromSegment(pc, pc1);
point f = LineInter(l1, l2);
printf("%.3lf ", dis(pa, f) + dis(pb, f) + dis(pc, f));
}
void I(double a, double b, double c, double A, double B, double C)
{
//角平分线的交点到三顶点的距离和
double t, ans;
t = (a + b - c) / 2;
ans = t / cos(C / 2) + (b - t) / cos(A / 2) + (a - t) / cos(B / 2);
printf("%.3lf ", ans);
}
void G(double a, double b, double c, double A, double B, double C)
{
//中线的交点
double t1, t2, t3;
t1 = sqrt((b / 2) * (b / 2) + a * a - 2 * a * b / 2 * cos(C));
t2 = sqrt((a / 2) * (a / 2) + c * c - 2 * a * c / 2 * cos(B));
t3 = sqrt((c / 2) * (c / 2) + b * b - 2 * b * c / 2 * cos(A));
t1 += t2 + t3;
printf("%.3lf ", t1 * 2 / 3);
}
void O(double a, double b, double c, double A, double B, double C)
{ //垂线的交点
double t = (A + C - B) / 2;
printf("%.3lf\n", 3 * b / 2 / cos(t));
}
int main()
{
int i, ca;
double a, b, c;
double A, B, C;
cin >> ca;
for(i = 1;i <= ca;i++){
cin >> a >> b >> c;
A = (b * b + c * c - a * a) / 2 / b / c;
B = (a * a + c * c - b * b) / 2 / a / c;
C = (a * a + b * b - c * c) / 2 / a / b;
A = acos(A), B = acos(B), C = acos(C);
F(a, b, c, A, B, C);
I(a, b, c, A, B, C);
G(a, b, c, A, B, C);
O(a, b, c, A, B, C);
}
return 0;
}
============================================================================================--------------------------------------------------------------------------------------------
(5)求直线的交点
/*求直线的交点，注意平形的情况无解，避免RE*/
TPoint LineInter(TLine l1, TLine l2)
{
//求两直线得交点坐标
TPoint tmp;
double a1 = l1.a;
double b1 = l1.b;
double c1 = l1.c;
double a2 = l2.a;
double b2 = l2.b;
double c2 = l2.c;
//注意这里b1 = 0
if(fabs(b1) < eps){
tmp.x = -c1 / a1;
tmp.y = (-c2 - a2 * tmp.x) / b2;
}
else{
tmp.x = (c1 * b2 - b1 * c2) / (b1 * a2 - b2 * a1);
tmp.y = (-c1 - a1 * tmp.x) / b1;
}
//cout << "交点坐标" << endl;
//cout << a1 * tmp.x + b1 * tmp.y + c1 << endl;
//cout << a2 * tmp.x + b2 * tmp.y + c2 << endl;
return tmp;
}
(6)根据线段两端点的坐标求垂直平分线上除中点外的另一点
TPoint GetOtherPoint(TPoint pre, TPoint tmp)
{
/*根据线段两端点的坐标求垂直平分线上除中点外的另一点*/
double kx, ky;
TPoint other, mid;
mid.x = (pre.x + tmp.x) / 2;
mid.y = (pre.y + tmp.y) / 2;
kx = pre.x - tmp.x;
ky = pre.y - tmp.y;
if(fabs(kx) < eps){
other.y = mid.y;
other.x = 1.0;
if(fabs(other.x - mid.x) < eps) other.x = 2.0;
}
else if(fabs(ky) < eps){
other.x = mid.x;
other.y = 1.0;
if(fabs(other.y - mid.y) < eps) other.y = 2.0;
}
else {
double k = -kx / ky;
other.x = 1.0;
if(fabs(other.x - mid.x) < eps) other.x = 2.0;
other.y = mid.y + k * (other.x - mid.x);
}
return other;
}
(7)根据两点坐标求直线方程
TLine lineFromSegment(TPoint p1, TPoint p2)
{
//线段所在直线,返回直线方程的三个系统
TLine tmp;
tmp.a = p2.y - p1.y;
tmp.b = p1.x - p2.x;
tmp.c = p2.x * p1.y - p1.x * p2.y;
return tmp;
}
(8)差积的应用
double multi(TPoint p1, TPoint p2, TPoint p0)
{
//求矢量[p0, p1], [p0, p2]的叉积
//p0是顶点
return (p1.x - p0.x) * (p2.y - p0.y) - (p2.x - p0.x) * (p1.y - p0.y);
//若结果等于0，则这三点共线
//若结果大于0，则p0p2在p0p1的逆时针方向
//若结果小于0，则p0p2在p0p1的顺时针方向
}
/*
折线的拐向的判断（从p0向p1看过去的左边）
若 (p2 - p1) 叉乘 (p1 - p0) < 0 ,则p0p1在p1点拐向左侧后得到p1p2
若 (p2 - p1) 叉乘 (p1 - p0) = 0 ,则 p0, p1, p2三点共线
若 (p2 - p1) 叉乘 (p1 - p0) > 0 ,则p0p1在p1点拐向右侧后得到p1p2
*/
(9)三角形的面积公式
//角形面积的计算
// S = ah / 2
// S = absinC / 2
// S = sqrt(p * (p - a) * (p - b) * (p - c)), p = (a + b + c) / 2
// S = abc / R / 4
double triangleArea(TPoint t[])
{
return fabs((t[1].x - t[0].x) * (t[2].y - t[0].y)
- (t[2].x - t[0].x) * (t[1].y - t[0].y)) / 2;
}
(10)三角形的内接圆(未检验正确性)
TCircle incircleOfTriangle(TTriangle t)
{
//三角形的内接圆(未检验正确性)
TCircle tmp;
double a, b, c, angleA, angleB, angleC, p, p2, p3, f1, f2;
double xA, yA, xB, yB, xC, yC;
a = distance(t.t[0], t.t[1]);
b = distance(t.t[1], t.t[2]);
c = distance(t.t[2], t.t[0]);
/*
S = p * r
p = (a + b + c) / 2;
r = S / P = 2 * S / (a + b + c)
*/
tmp.r = 2 * triangleArea(t) / (a + b +c);
angleA = acos((b * b + c * c - a * a) / (2 * b * c));
angleB = acos((a * a + c * c - b * b) / (2 * a * c));
angleC = acos((a * a + b * b - c * c) / (2 * a * b));
p = sin(angleA / 2);
p2 = sin(angleB / 2);
p3 = sin(angleC / 2);
xA = t.t[0].x; yA = t.t[0].y;
xB = t.t[1].x; yB = t.t[1].y;
xC = t.t[2].x; yC = t.t[2].y;
f1 = ((tmp.r / p2) * (tmp.r / p2) - (tmp.r / p) * (tmp.r / p) +
xA * xA - xB * xB + yA * yA - yB * yB) / 2;
f2 = ((tmp.r / p3) * (tmp.r / p3) - (tmp.r / p) * (tmp.r / p) +
xA * xA - xC * xC + yA * yA - yC * yC) / 2;
tmp.centre.x = (f1 * (yA - yC) - f2 * (yA - yB)) /
((xA - xB) * (yA - yC) - (xA - xC) * (yA - yB));
tmp.centre.y = (f1 * (xA - xC) - f2 * (xA - xB)) /
((yA - yB) * (xA - xC) - (yA - yC) * (xA - xB));
return tmp;
}
(11)多边形的面积(适合凹多边形)
double polygonArea(TPolygon p, int n)
{
//已知多边形各顶点的坐标，求其面积
double area;
int i;
area = 0;
for(i = 1;i <= n;i++){
area += (p.p[i - 1].x * p.p[i % n].y - p.p[i % n].x * p.p[i - 1].y);
}
return area / 2;
}
(12)判断点是否在线段上
bool isPointOnSegment(TPoint p, TPoint p1, TPoint p2)
{
// 判断p点是否在线段p1p2上
//1.p是否在直线p1p2上
//2.p是否在以p1p2为对角线的矩形中
if(multi(p1, p2, p) != 0) return false ;
if((p.x > p1.x && p.x > p2.x) || (p.x < p1.x && p.x < p2.x)) return false;
if((p.y > p1.y && p.y > p2.y) || (p.y < p1.y && p.y < p2.y)) return false;
return true;
}
(13)平面上两个点之间的距离
double distance(TPoint p1, TPoint p2)
{
//计算平面上两个点之间的距离
return sqrt((p1.x - p2.x) * (p1.x - p2.x) + (p1.y - p2.y) * (p1.y - p2.y));
}
(14)p点关于直线L的对称点
TPoint symmetricalPointofLine(TPoint p, TLine L)
{
//p点关于直线L的对称点
TPoint p2;
double d;
d = L.a * L.a + L.b * L.b;
p2.x = (L.b * L.b * p.x - L.a * L.a * p.x -
2 * L.a * L.b * p.y - 2 * L.a * L.c) / d;
p2.y = (L.a * L.a * p.y - L.b * L.b * p.y -
2 * L.a * L.b * p.x - 2 * L.b * L.c) / d;
return p2;
}
//点关于线段的对称点
//首先可以根据线段的两个端点求出线段所在的直线L，然后再来求指定
//点关于直线L的对称点
(15)判断一个矩形是否在另一个矩形中
/*
判断一个矩形是否在另一个矩形中
*/
#include <stdio.h>
#include <math.h>
#define pi acos(-1)
void swap(double &a, double &b)
{
double tmp;
tmp = a;
a = b;
b = tmp;
}
int main()
{
double a, b, x, y;
while(scanf("%lf%lf%lf%lf", &a, &b, &x, &y)){
if(a == 0 && b == 0 && x== 0 && y== 0) break;
if(a * b <= x * y){
printf("Box cannot be dropped.\n");
continue;
}
if(a < b) swap(a, b);
if(x < y) swap(x, y);
if((a >= x && b >= y) || (x <= b)){
printf("Escape is possible.\n");
continue;
}
int tag = 0;
for(double i = 0;i <= 90;i += 0.2){
double angle = i * pi / 180;
if(x * cos(angle) + y * sin(angle) < a
&& y * cos(angle) + x * sin(angle ) < b){
printf("Escape is possible.\n");
tag = 1;
break;
}
}
if(tag == 0) printf("Box cannot be dropped.\n");
}
return 0;
}
(16)直线和圆的交点+点关于线的对称点+点到线的距离+直线方程(fzu_1035)
/*
fzu_1035
1.直线和圆的交点
2.点关于线的对称点
3.点到线的距离
4.直线方程
*/
#include <iostream>
#include <cmath>
using namespace std;
#define INF 999999999
const double eps = 1e-6;
int up;
typedef struct TPoint
{
double x;
double y;
}TPoint;
typedef struct TCircle
{
TPoint center;
double r;
}TCircle;
typedef struct TLine
{
//直线标准式中的系数
double a, b, c;
}TLine;
void SloveLine(TLine &line, TPoint start, TPoint dir)
{
//根据直线上一点和直线的方向求直线的方程
if(dir.x == 0){
line.a = 1;
line.b = 0;
line.c = start.x;
}
else {
double k = dir.y / dir.x;
line.a = k;
line.b = -1;
line.c = start.y - k * start.x;
}
}
TLine lineFromSegment(TPoint p1, TPoint p2)
{
//线段所在直线,返回直线方程的三个系统
TLine tmp;
tmp.a = p2.y - p1.y;
tmp.b = p1.x - p2.x;
tmp.c = p2.x * p1.y - p1.x * p2.y;
return tmp;
}
TPoint symmetricalPointofLine(TPoint p, TLine L)
{
//p点关于直线L的对称点
TPoint p2;
double d;
d = L.a * L.a + L.b * L.b;
p2.x = (L.b * L.b * p.x - L.a * L.a * p.x -
2 * L.a * L.b * p.y - 2 * L.a * L.c) / d;
p2.y = (L.a * L.a * p.y - L.b * L.b * p.y -
2 * L.a * L.b * p.x - 2 * L.b * L.c) / d;
return p2;
}
double distanc(TPoint p1, TPoint p2)
{
//计算平面上两个点之间的距离
return sqrt((p1.x - p2.x) * (p1.x - p2.x) + (p1.y - p2.y) * (p1.y - p2.y));
}
bool samedir(TPoint dir, TPoint start, TPoint point)
{
//判断方向
TPoint tmp;
tmp.x = point.x - start.x;
tmp.y = point.y - start.y;
if(tmp.x != 0 && dir.x != 0){
if(tmp.x / dir.x > 0) return true;
else return false;
}
else if(tmp.y != 0 && dir.y != 0){
if(tmp.y / dir.y > 0) return true;
else return false;
}
return true;
}
bool Intersected(TPoint &point, TLine line, const TCircle circle[],
TPoint start, TPoint dir, int which)
{
//如果圆与直线有(有效交点)交点就存放在变量point中
double a = line.a, b = line.b, c = line.c;
double x0 = circle[which].center.x, y0 = circle[which].center.y;
double r = circle[which].r;
//有交点，求交点
double x2front = b * b + a * a;
double x1front = -2 * x0 * b * b + 2 * a * b * y0 + 2 * a * c;
double front = x0 * x0 * b * b + y0 * y0 * b * b
+ c * c + 2 * c * y0 * b - b * b * r * r;
double d = x1front * x1front - 4 * x2front * front;
TPoint p1, p2;
bool k1, k2;
if(fabs(d) < eps){
//x2front不可能等于零
point.x = -x1front / x2front / 2;
point.y = (-c - a * point.x) / b;
//判断方向
if(samedir(dir, start, point)) return true;
else return false;
}
else if(d < 0) return false;
else {
p1.x = (-x1front + sqrt(d)) / 2 / x2front;
p1.y = (-c - a * p1.x) / b;
p2.x = (-x1front - sqrt(d)) / 2 / x2front;
p2.y = (-c - a * p2.x) / b;
k1 = samedir(dir, start, p1);
k2 = samedir(dir, start, p2);
if(k1 == false && k2 == false) return false;
if(k1 == true && k2 == true){
double dis1 = distanc(p1, start);
double dis2 = distanc(p2, start);
if(dis1 < dis2) point = p1;
else point = p2;
return true;
}
else if(k1 == true) point = p1;
else point = p2;
return true;
}
}
void Reflect(int &num, TCircle circle[], TPoint start, TPoint dir, int n)
{
//反复反射
int i;
TLine line;
TPoint interpoint, newstart;
int u;
SloveLine(line, start, dir);
int tag = 0;
double mindis = INF;
for(i = 1;i <= n;i++){
if(i != up && Intersected(interpoint, line, circle, start, dir, i)){
double dis = distanc(start, interpoint);
if(dis < mindis){
tag = 1;
u = i;
mindis = dis;
newstart = interpoint;
}
}
}
if(tag == 0){
cout << "inf" << endl;
return ;
}
else {
if(num == 10){
cout << "..." << endl;
return ;
}
cout << u << " ";
num++;
//新的方向
TLine line1;
TPoint p;
line1 = lineFromSegment(newstart, circle[u].center);
if(fabs(line1.a * start.x + line1.b * start.y +line1.c) <= eps){
dir.x = -dir.x;
dir.y = -dir.y;
}
else {
p = symmetricalPointofLine(start, line1);//start的对称点
dir.x = p.x - newstart.x;
dir.y = p.y - newstart.y;
}
start = newstart;
up = u;
Reflect(num, circle, start, dir, n);
}
}
int main()
{
//freopen("fzu_1035.in", "r", stdin);
//freopen("fzu_1035.out", "w", stdout);
int n, i, j, num, test = 1;
TCircle circle[30];
TPoint start, dir;
while(cin >> n && n){
for(i = 1;i <= n;i++){
cin >> circle[i].center.x >> circle[i].center.y >> circle[i].r;
}
cin >> start.x >> start.y >> dir.x >> dir.y;
cout << "Scene " << test++ << endl;
num = 0;
up = -1;
Reflect(num, circle, start, dir, n);
cout << endl;
}
return 0;
}
(17)判断点是否在多边形内
#include <stdio.h>
#include <math.h>
#define Maxn 1000
#define eps 1e-6
#define pi acos(-1)
#define inf 9999999
struct TPoint
{
double x, y;
};
double dispoint(TPoint p1, TPoint p2)
{
return sqrt((p1.x - p2.x) * (p1.x - p2.x) + (p1.y - p2.y) * (p1.y - p2.y));
}
int pall(TPoint p1, TPoint p2, TPoint p3, TPoint p4)
{
if(fabs((p1.y - p2.y) * (p3.x - p4.y)
- (p1.x - p2.x) * (p3.y - p4.y)))
return 1;
else return 0;
}
int IsOnSegment(TPoint p, TPoint p1, TPoint p2)
{
double d1 = dispoint(p1, p);
double d2 = dispoint(p, p2);
double d3 = dispoint(p1, p2);
if(fabs(d1 + d2 - d3) < eps) return 1;
else return 0;
}
double max(double x, double y)
{
//比较两个数的大小，返回大的数
if(x > y) return x;
else return y;
}
double min(double x, double y)
{
//比较两个数的大小，返回小的数
if(x < y) return x;
else return y;
}
double multi(TPoint p1, TPoint p2, TPoint p0)
{
//求矢量[p0, p1], [p0, p2]的叉积
//p0是顶点
return (p1.x - p0.x) * (p2.y - p0.y) - (p2.x - p0.x) * (p1.y - p0.y);
//若结果等于0，则这三点共线
//若结果大于0，则p0p2在p0p1的逆时针方向
//若结果小于0，则p0p2在p0p1的顺时针方向
}
bool isIntersected(TPoint s1, TPoint e1, TPoint s2, TPoint e2)
{
//判断线段是否相交
//1.快速排斥试验判断以两条线段为对角线的两个矩形是否相交
//2.跨立试验
if(
(max(s1.x, e1.x) >= min(s2.x, e2.x)) &&
(max(s2.x, e2.x) >= min(s1.x, e1.x)) &&
(max(s1.y, e1.y) >= min(s2.y, e2.y)) &&
(max(s2.y, e2.y) >= min(s1.y, e1.y)) &&
(multi(s2, e1, s1) * multi(e1, e2, s1) >= 0) &&
(multi(s1, e2, s2) * multi(e2, e1, s2) >= 0)
) return true;
return false;
}
int IsInPolygon(TPoint p0, TPoint p[], int n)
{
double a, k;
TPoint p2;
int i, j, tag;
for(a = 1.0;a < pi / 2; a += 1.0){
k = tan(a);
p2.x = inf;
p2.y = p0.y + k * (p2.x - p0.x);
tag = 0;
for(i = 0;i < n;i++){
if(IsOnSegment(p[i], p0, p2)){
tag = 1;
break;
}
}
if(tag == 1) continue;
for(i = 0;i < n;i++){
if(pall(p0, p2, p[i], p[i + 1])){
tag = 1;
break;
}
}
if(tag == 0) break;
}
int num = 0;
for(i = 0;i < n;i++){
if(isIntersected(p0, p2, p[i], p[i + 1])){
num++;
}
}
if(num % 2 == 0) return 0;
else return 1;
}
int main()
{
//freopen("in.in", "r", stdin);
//freopen("out.out", "w", stdout);
TPoint p[Maxn], p0;
int n, m, i, j, tag, test = 1;
while(scanf("%d", &n) && n){
if(test != 1) printf("\n");
printf("Problem %d:\n", test++);
scanf("%d", &m);
for(i = 0;i < n;i++){
scanf("%lf%lf", &p[i].x, &p[i].y);
}
p[n] = p[0];
for(i = 1;i <= m;i++){
scanf("%lf%lf", &p0.x, &p0.y);
tag = 0;
for(j = 0;j < n;j++){
if(IsOnSegment(p0, p[j], p[j + 1])){
printf("Within\n");
tag = 1;
break;
}
}
if(tag == 1) continue;
if(IsInPolygon(p0, p, n)) printf("Within\n");
else printf("Outside\n");
}
}
return 0;
}
/*
判断点q是否在多边形内的一种方法，过q作水平射线L，
如果L与多边形P的边界不相交，则q在P的外部。否则，
L和P的边界相交，具体地说，交点个数为奇（偶）数
时，点q在P的内（外）部。 */
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#define MaxNode 50
#define INF 999999999
typedef struct TPoint
{
double x;
double y;
}TPoiont;
typedef struct TSegment
{
TPoint p1;
TPoint p2;
}TSegment;
typedef struct TPolygon
{
TPoint point[MaxNode];
int n;
}TPolygon;
double multi(TPoint p1, TPoint p2, TPoint p0)
{
//求矢量[p0, p1], [p0, p2]的叉积
//p0是顶点
return (p1.x - p0.x) * (p2.y - p0.y) - (p2.x - p0.x) * (p1.y - p0.y);
//若结果等于0，则这三点共线
//若结果大于0，则p0p2在p0p1的逆时针方向
//若结果小于0，则p0p2在p0p1的顺时针方向
}
double max(double x, double y)
{
//比较两个数的大小，返回大的数
if(x > y) return x;
else return y;
}
double min(double x, double y)
{
//比较两个数的大小，返回小的数
if(x < y) return x;
else return y;
}
bool Intersect(TSegment L1, TSegment L2)
{
//线段l1与l2相交而且不在端点上时，返回true
//判断线段是否相交
//1.快速排斥试验判断以两条线段为对角线的两个矩形是否相交
TPoint s1 = L1.p1;
TPoint e1 = L1.p2;
TPoint s2 = L2.p1;
TPoint e2 = L2.p2;
//2.跨立试验
if(
(max(s1.x, e1.x) > min(s2.x, e2.x)) &&
(max(s2.x, e2.x) > min(s1.x, e1.x)) &&
(max(s1.y, e1.y) > min(s2.y, e2.y)) &&
(max(s2.y, e2.y) > min(s1.y, e1.y)) &&
(multi(s2, e1, s1) * multi(e1, e2, s1) > 0) &&
(multi(s1, e2, s2) * multi(e2, e1, s2) > 0)
) return true;
return false;
}
bool Online(TSegment L, TPoint p)
{
//p在L上(不在端点)时返回true
//1.在L所在的直线上
//2.在L为对角线的矩形中
double dx, dy, dx1, dy1;
dx = L.p2.x - L.p1.x;
dy = L.p2.y - L.p1.y;
dx1 = p.x - L.p1.x;
dy1 = p.y - L.p1.y;
if(dx * dy1 - dy * dx1 != 0) return false; //叉积
if(dx1 * (dx1 - dx) < 0 || dy1 * (dy1 - dy) < 0) return true;
return false;
}
bool same1(TSegment L, TPoint p1, TPoint p2)
{
//判断p1, p2是否在L的同侧
if(multi(p1, L.p2, L.p1) * multi(L.p2, p2, L.p1)< 0) return true;
return false;
}
bool Inside(TPoint q, TPolygon polygon)
{
int c, i;
TSegment L1, L2;
c = 0;
L1.p1 = q;
L1.p2 = q;
L1.p2.x = INF;
/*
（1）相交
1.p[i]和p[i+1]在L的两侧
2.p[i]和p[i+2]在L的同侧
3.p[u]和p[i+3]在L的同侧或异侧
*/
for(i = 0;i <= polygon.n - 1;i++){
L2.p1 = polygon.point[i];
L2.p2 = polygon.point[(i + 1) % polygon.n];
if(Intersect(L1, L2)){
c++;
continue;
}
if(!Online(L1, polygon.point[(i + 1) % polygon.n])) continue;
if(!Online(L1, polygon.point[(i + 2) % polygon.n]) &&
!same1(L1, polygon.point[i], polygon.point[(i + 2) % polygon.n])){
c++;
continue;
}
if(Online(L1, polygon.point[(i + 2) % polygon.n]) &&
!same1(L1, polygon.point[i], polygon.point[(i + 3) % polygon.n]))
c++;
}
if(c % 2 == 0) return false;
else return true;
}
int main()
{
int i, test, k;
int primp, primq;
TPoint p;
p.x = 0;
p.y = 0;
test = 1;
TPolygon polygon;
while(scanf("%d", &polygon.n) != EOF && polygon.n){
printf("Pilot %d\n", test++);
for(i = 0;i <= polygon.n - 1;i++){
scanf("%lf%lf", &polygon.point[i].x, &polygon.point[i].y);
}
scanf("%d%d", &primp, &primq);
if(Inside(p, polygon)){
printf("The pilot is in danger!\n");
k = (primp - 1) * (primq - 1) / 2;
printf("The secret number is %d.\n", k);
}
else printf("The pilot is safe.\n");
printf("\n");
}
return 0;
}
(18)N点中三个点组成三角形面积最大
//Rotating Calipers algorithm
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#define MaxNode 50005
int stack[MaxNode];
int top;
double max;
typedef struct TPoint
{
int x;
int y;
}TPoint;
TPoint point[MaxNode];
void swap(TPoint point[], int i, int j)
{
TPoint tmp;
tmp = point[i];
point[i] = point[j];
point[j] = tmp;
}
double multi(TPoint p1, TPoint p2, TPoint p0)
{
return (p1.x - p0.x) * (p2.y - p0.y) - (p2.x - p0.x) * (p1.y - p0.y);
}
double distance(TPoint p1, TPoint p2)
{
return (p1.x - p2.x) * (p1.x - p2.x) + (p1.y - p2.y) * (p1.y - p2.y);
}
int cmp(const void *a, const void *b)
{
TPoint *c = (TPoint *)a;
TPoint *d = (TPoint *)b;
double k = multi(*c, *d, point[0]);
if(k< 0) return 1;
else if(k == 0 && distance(*c, point[0]) >= distance(*d, point[0]))
return 1;
else return -1;
}
void grahamScan(int n)
{
//Graham扫描求凸包
int i, u;
//将最左下的点调整到p[0]的位置
u = 0;
for(i = 1;i <= n - 1;i++){
if((point[i].y < point[u].y) ||
(point[i].y == point[u].y && point[i].x < point[u].x))
u = i;
}
swap(point, 0, u);
//将平p[1]到p[n - 1]按按极角排序，可采用快速排序
qsort(point + 1, n - 1, sizeof(point[0]), cmp);
for(i = 0;i <= 2;i++) stack[i] = i;
top = 2;
for(i = 3;i <= n - 1;i++){
while(multi(point[i], point[stack[top]], point[stack[top - 1]]) >= 0){
top--;
if(top == 0) break;
}
top++;
stack[top] = i;
}
}
int main()
{
double triangleArea(int i, int j, int k);
void PloygonTriangle();
int i, n;
while(scanf("%d", &n) && n != -1){
for(i = 0;i < n;i++)
scanf("%d%d", &point[i].x, &point[i].y);
if(n <= 2){
printf("0.00\n");
continue;
}
if(n == 3){
printf("%.2lf\n", triangleArea(0, 1, 2));
continue;
}
grahamScan(n);
PloygonTriangle();
printf("%.2lf\n", max);
}
return 0;
}
void PloygonTriangle()
{
double triangleArea(int i, int j, int k);
int i, j , k;
double area, area1;
max = -1;
for(i = 0;i <= top - 2;i++){
k = -1;
for(j = i + 1; j <= top - 1;j++){
if(k <= j) k= j + 1;
area = triangleArea(stack[i], stack[j], stack[k]);
if(area > max) max = area;
while(k + 1 <= top){
area1= triangleArea(stack[i], stack[j], stack[k + 1]);
if(area1 < area) break;
if(area1 > max) max = area1;
area = area1;
k++;
}
}
}
}
double triangleArea(int i, int j, int k)
{
//已知三角形三个顶点的坐标，求三角形的面积
double l = fabs(point[i].x * point[j].y + point[j].x * point[k].y
+ point[k].x * point[i].y - point[j].x * point[i].y
- point[k].x * point[j].y - point[i].x * point[k].y) / 2;
return l;
}
(19)扇形的重心
//Xc = 2*R*sinA/3/A
//A为圆心角的一半
#include <stdio.h>
#include <math.h>
int main()
{
double r, angle;
while(scanf("%lf%lf", &r, &angle) != EOF){
angle /= 2;
printf("%.6lf\n", 2 * r * sin(angle) / 3 / angle);
}
return 0;
}
(20)多边形的重心
/*
题目描述：
有一个密度均匀的平面N多边形(3 <= N <= 1000000)，可能凹也可能凸，但没有边相交叉，
另外已知N个有序(顺时针或逆时针)顶点的坐标值，第j个顶点坐标为（Xi ， Yi ），且满
足 (|Xi|, |Yi| <= 20000)，求这个平面多边形的重心。
解题过程：
从第1个顶点出发，分别连接第i, i+1个顶点组成三角形Ti,1 < i < n,
一共n-2个三角形正好是多连形的一个划分，分别求出每个三角形的面积Si，
总面积为各个面积相加
根据物理学知识得：n个点(xi,yi)每个重量是mi,则重心是
X = (x1*M1+x2*M2+...+xn*Mn) / (M1+M2+....+Mn)
Y = (y1*M1+y2*M2+...+yn*Mn) / (M1+M2+....+Mn)
另个需要用的知识有：
已知3点求三角形的面积，设三点分别为p[0].x, p[0].y p[1].x, p[1].y p[1].x, p[1].y
面积s =[ p[0].x*p[1].y-p[1].x*p[0].y + p[1].x*p[2].y-p[2].x*p[1].y + p[2].x*p[0].y-p[0].x*p[2].y ] / 2 ,
这是这3个点是逆时针的值，顺时针取负。
已知3点求重心，x = (p[0].x+p[1].x+p[2].x)/3.0 , y = (p[0].y+p[1].y+p[2].y)/3.0
另外在求解的过程中，不需要考虑点的输入顺序是顺时针还是逆时针，相除后就抵消了，
还要注意 一点是不必在求每个小三角形的重心时都除以3，可以在最后除一下
*/
/*fzu_1132*/
#include <stdio.h>
#include <math.h>
typedef struct TPoint
{
double x;
double y;
}TPoint;
double triangleArea(TPoint p0, TPoint p1, TPoint p2)
{
//已知三角形三个顶点的坐标，求三角形的面积
double k = p0.x * p1.y + p1.x * p2.y
+ p2.x * p0.y - p1.x * p0.y
- p2.x * p1.y - p0.x * p2.y;
//if(k >= 0) return k / 2;
// else return -k / 2;
return k / 2;
}
int main()
{
int i, n, test;
TPoint p0, p1, p2, center;
double area, sumarea, sumx, sumy;
scanf("%d", &test);
while(test--){
scanf("%d", &n);
scanf("%lf%lf", &p0.x, &p0.y);
scanf("%lf%lf", &p1.x, &p1.y);
sumx = 0;
sumy = 0;
sumarea = 0;
for(i = 2;i < n;i++){
scanf("%lf%lf", &p2.x, &p2.y);
center.x = p0.x + p1.x + p2.x;
center.y = p0.y + p1.y + p2.y;
area = triangleArea(p0, p1, p2);
sumarea += area;
sumx += center.x * area;
sumy += center.y * area;
p1 = p2;
}
printf("%.2lf %.2lf\n", sumx / sumarea / 3, sumy / sumarea / 3);
}
return 0;
}
(21)判断N点是否共面
/*fzu_1393*/
#include <iostream>
#include <cmath>
using namespace std;
typedef struct TPoint
{
double x;
double y;
double z;
}TPoint;
const double eps = 1e-6;
int main()
{
TPoint a, b, c, d, ab, ac, abc;
int n, i;
cin >> n;
while(n--){
cin >> a.x >> a.y >> a.z
>> b.x >> b.y >> b.z
>> c.x >> c.y >> c.z
>> d.x >> d.y >> d.z;
ab.x = b.x - a.x;
ab.y = b.y - a.y;
ab.z = b.z - a.z;
ac.x = c.x - a.x;
ac.y = c.y - a.y;
ac.z = c.z - a.z;
abc.x = ab.y * ac.z - ab.z * ac.y;
abc.y = -(ab.x * ac.z - ab.z * ac.x);
abc.z = ab.x * ac.y - ab.y * ac.x;
//abc.x * (x - a.x) + abc.y * (y - a.y) + abc.z * (z - a.z) = 0;
if(fabs(abc.x * (d.x - a.x) + abc.y * (d.y - a.y)
+ abc.z * (d.z - a.z)) < eps)
printf("Yes\n");
else printf("No\n");
}
return 0;
}
(22)求共线的点最多为多少
/*
2617120 chenhaifeng 1118 Accepted 512K 1890MS C++ 977B 2007-09-04 18:43:26
直接O(n^3)超时，用一个标记数组，标记i,j所做直线已经查找过，可以跳过
大牛的思想
朴素做法是 O(n3) 的，超时。我的做法是枚举每个点，
然后求其它点和它连线的斜率，再排序。这样就得到经过
该点的直线最多能经过几个点。求个最大值就行了。复
杂度是 O(n2logn) 的。把排序换成 hash，
可以优化到 O(n2)。
2617134 chenhaifeng 1118 Accepted 276K 312MS G++ 1394B 2007-09-04 18:49:08
*/
/*pku_1118_共线最多的点的个数.cpp*/
#include <stdio.h>
#include <math.h>
bool f[705][705];
int a[705];
int main()
{
int n, i, j, s, num, maxn;
int x[705], y[705];
int t, m;
while(scanf("%d", &n) != EOF && n){
for(i = 0;i <= n - 1;i++){
scanf("%d%d", &x[i], &y[i]);
}
maxn = -1;
for(i = 0;i <= n - 1;i++){
for(j = i;j <= n - 1;j++){
f[i][j] = false;
}
}
for(i = 0;i <= n - 1;i++){
for(j = i + 1;j <= n - 1;j++){
if(f[i][j] == true) continue;
if(n - j < maxn) break;
num = 2;
t = 2;
a[0] = i;
a[1] = j;
f[i][j] = true;
for(s = j + 1;s <= n - 1;s++){
if(f[i][s] == true || f[j][s] == true) continue;
if((y[i] - y[s]) * (x[j] - x[s]) == (x[i] - x[s]) * (y[j] - y[s])){
num++;
a[t] = s;
for(m = 0;m <= t - 1;m++){
f[m][s] = true;
}
t++;
}
}
if(num > maxn) maxn = num;
}
}
printf("%d\n", maxn);
}
return 0;
}
(23)N个矩形的相交的面积
/*
大牛的思想
题目给出 n 个矩形，要求它们的面积并。具体做法是离散化。
先把 2n 个 x 坐标排序去重，然后再把所有水平线段（
要记录是矩形上边还是下边）按 y 坐标排序。
最后对于每一小段区间 (x[i], x[i + 1]) 扫描所有的水平线段，
求出这些水平线段在小区间内覆盖的面积。总的时间复杂度是 O(n^2)。
利用线段树，可以优化到 O(nlogn)。
pku_1151
*/
#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#define up 1
#define down -1
typedef struct TSeg
{
double l, r;
double y;
int UpOrDown;
}TSeg;
TSeg seg[210];
int segn;
double x[210];
int xn;
int cmp1(const void *a, const void *b)
{
if(*(double *)a < *(double *)b) return -1;
else return 1;
}
int cmp2(const void *a, const void *b)
{
TSeg *c = (TSeg *)a;
TSeg *d = (TSeg *)b;
if(c->y < d->y) return -1;
else return 1;
}
void movex(int t, int &xn)
{
int i;
for(i = t;i <= xn - 1;i++){
x[i] = x[i + 1];
}
xn--;
}
int main()
{
//freopen("in.in", "r", stdin);
//freopen("out.out", "w", stdout);
int n, i, j, cnt, test = 1;
double x1, y1, x2, y2, ylow, area;
while(scanf("%d", &n) != EOF && n){
xn = 0;
segn = 0;
for(i = 0;i < n;i++){
scanf("%lf%lf%lf%lf", &x1, &y1, &x2, &y2);
x[xn++] = x1;
x[xn++] = x2;
seg[segn].l = x1;
seg[segn].r = x2;
seg[segn].y = y1;
seg[segn++].UpOrDown = up;
seg[segn].l = x1;
seg[segn].r = x2;
seg[segn].y = y2;
seg[segn++].UpOrDown = down;
}
qsort(x, xn, sizeof(x[0]), cmp1);
/*除掉重复的x*/
for(i = 1;i < xn;){
if(x[i] == x[i - 1]) movex(i, xn);
else i++;
}
qsort(seg, segn, sizeof(seg[0]), cmp2);
area = 0.0;
for(i = 0;i < xn - 1;i++){
cnt = 0;
for(j = 0;j < segn;j++){
if(seg[j].l <= x[i] && seg[j].r >= x[i + 1]){
if(cnt == 0) ylow = seg[j].y;
if(seg[j].UpOrDown == down) cnt++;
else cnt--;
if(cnt == 0) area += (x[i + 1] - x[i]) * (seg[j].y - ylow);
}
}
}
printf("Test case #%d\n", test++);
printf("Total explored area: %.2lf\n", area);
}
return 0;
}
(24)三角形外接圆+圆的参数方程(pku_1266)
//pku_1266_三角形外接圆+圆的参数方程.cpp
#include <stdio.h>
#include <math.h>
const double eps = 1e-6;
const double pi = acos(-1);
const double inf = 9999999.0;
typedef struct TPoint
{
double x, y;
}TPoint;
typedef struct TCircle
{
TPoint c;
double r;
}TCircle;
typedef struct TTriangle
{
TPoint t[3];
}TTriangle;
double max(double x, double y)
{
if(x > y) return x;
else return y;
}
double min(double x, double y)
{
if(x < y) return x;
else return y;
}
double distance(TPoint p1, TPoint p2)
{
return sqrt((p1.x - p2.x) * (p1.x - p2.x)
+ (p1.y - p2.y) * (p1.y - p2.y));
}
double triangleArea(TTriangle t)
{
//已知三角形三个顶点的坐标，求三角形的面积
return fabs(t.t[0].x * t.t[1].y + t.t[1].x * t.t[2].y
+ t.t[2].x * t.t[0].y - t.t[1].x * t.t[0].y
- t.t[2].x * t.t[1].y - t.t[0].x * t.t[2].y) / 2;
}
TCircle circumcircleOfTriangle(TTriangle t)
{
//三角形的外接圆
TCircle tmp;
double a, b, c, c1, c2;
double xA, yA, xB, yB, xC, yC;
a = distance(t.t[0], t.t[1]);
b = distance(t.t[1], t.t[2]);
c = distance(t.t[2], t.t[0]);
//根据S = a * b * c / R / 4;求半径R
tmp.r = a * b * c / triangleArea(t) / 4;
xA = t.t[0].x; yA = t.t[0].y;
xB = t.t[1].x; yB = t.t[1].y;
xC = t.t[2].x; yC = t.t[2].y;
c1 = (xA * xA + yA * yA - xB * xB - yB * yB) / 2;
c2 = (xA * xA + yA * yA - xC * xC - yC * yC) / 2;
tmp.c.x = (c1 * (yA - yC) - c2 * (yA - yB)) /
((xA - xB) * (yA - yC) - (xA - xC) * (yA - yB));
tmp.c.y = (c1 * (xA - xC) - c2 * (xA - xB)) /
((yA - yB) * (xA - xC) - (yA - yC) * (xA - xB));
return tmp;
}
double Get_angle(TPoint p, TCircle c)
{
double ans;
ans = acos((p.x - c.c.x) / c.r);
if(p.y < c.c.y){
if(ans >= M_PI / 2){
ans = (M_PI - ans) + M_PI;
}
else {
ans = 2 * M_PI - ans;
}
}
return ans;
}
double swap(double &a, double &b)
{
double c;
c = a;
a = b;
b = c;
}
int UP(double a)
{
int a0;
if(a >= 0){
a0 = int(a);
if(fabs(a0 - a) < eps) return a0;
else return a0 + 1;
}
else {
a0 = (int)a;
if(fabs(a0 - 1 - a) < eps) return a0 - 1;
return a0;
}
}
int DOWN(double a)
{
int a0;
if(a >= 0){
a0 = (int)a;
if(fabs(a - (a0 + 1)) < eps) return a0 + 1;
return a0;
}
else {
a0 = (int)a;
if(fabs(a0 - a) < eps) return a0;
else return a0 - 1;
}
}
int main()
{
TCircle c;
TTriangle t;
TPoint p1, p2, p3;
double a1, a2, a3, tmp;
int max_x, min_x, max_y, min_y;
scanf("%lf%lf", &p1.x, &p1.y);
scanf("%lf%lf%lf%lf", &p2.x, &p2.y, &p3.x, &p3.y);
t.t[0] = p1;
t.t[1] = p2;
t.t[2] = p3;
c = circumcircleOfTriangle(t);
a1 = Get_angle(p1, c);
a2 = Get_angle(p2, c);
a3 = Get_angle(p3, c);
if(a1 > a2) swap(a1, a2);
//讨论
if(a3 >= a1 && a3 <= a2){
tmp = 0;
if(tmp >= a1 && tmp <= a2){
max_x = UP(c.c.x + c.r);
}
else max_x = UP(max(p1.x, p2.x));
tmp = pi / 2;
if(tmp >= a1 && tmp <= a2){
max_y = UP(c.c.y + c.r);
}
else max_y = UP(max(p1.y, p2.y));
tmp = pi;
if(tmp >= a1 && tmp <= a2){
min_x = DOWN(c.c.x - c.r);
}
else min_x = DOWN(min(p1.x, p2.x));
tmp = pi * 3 / 2;
if(tmp >= a1 && tmp <= a2){
min_y = DOWN(c.c.y - c.r);
}
else min_y = DOWN(min(p1.y, p2.y));
}
else {
tmp = 0;
if(tmp <= a1 || tmp >= a2){
max_x = UP(c.c.x + c.r);
}
else max_x = UP(max(p1.x, p2.x));
tmp = pi / 2;
if(tmp <= a1 || tmp >= a2){
max_y = UP(c.c.y + c.r);
}
else max_y = UP(max(p1.y, p2.y));
tmp = pi;
if(tmp <= a1 || tmp >= a2){
min_x = DOWN(c.c.x - c.r);
}
else min_x = DOWN(min(p1.x, p2.x));
tmp = pi * 3 / 2;
if(tmp <= a1 || tmp >= a2){
min_y = DOWN(c.c.y - c.r);
}
else min_y = DOWN(min(p1.y, p2.y));
}
printf("%d\n", (max_x - min_x) * (max_y - min_y));
return 0;
}
(25)判断线段是否有交点并求交点(cug_1078)
/*cug_1078判断线段是否有交点并求交点*/
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
const double eps = 1e-6;
typedef struct TPodouble
{
double x;
double y;
}TPodouble;
typedef struct TLine
{
double a, b, c;
}TLine;
double max(double x, double y)
{
if(x > y) return x;
else return y;
}
double min(double x, double y)
{
if(x < y) return x;
else return y;
}
double multi(TPodouble p1, TPodouble p2, TPodouble p0)
{
return (p1.x - p0.x) * (p2.y - p0.y) - (p2.x - p0.x) * (p1.y - p0.y);
}
TLine lineFromSegment(TPodouble p1, TPodouble p2)
{
TLine tmp;
tmp.a = p2.y - p1.y;
tmp.b = p1.x - p2.x;
tmp.c = p2.x * p1.y - p1.x * p2.y;
return tmp;
}
void Linedoubleer(TLine l1, TLine l2)
{
double x, y;
double a1 = l1.a;
double b1 = l1.b;
double c1 = l1.c;
double a2 = l2.a;
double b2 = l2.b;
double c2 = l2.c;
if(b1 == 0){
x = -c1 / (double)a1;
y = (-c2 - a2 * x) / (double)b2;
}
else{
x = (double)(c1 * b2 - b1 * c2) / (double)(b1 * a2 - b2 * a1);
y = (double)(-c1 - a1 * x) /(double)b1;
}
printf("%.3lf %.3lf\n", x, y);
}
bool isIntersected(TPodouble s1, TPodouble e1, TPodouble s2, TPodouble e2){
if(
(max(s1.x, e1.x) >= min(s2.x, e2.x)) &&
(max(s2.x, e2.x) >= min(s1.x, e1.x)) &&
(max(s1.y, e1.y) >= min(s2.y, e2.y)) &&
(max(s2.y, e2.y) >= min(s1.y, e1.y)) &&
(multi(s2, e1, s1) * multi(e1, e2, s1) >= 0) &&
(multi(s1, e2, s2) * multi(e2, e1, s2) >= 0)
) return true;
return false;
}
double dis(TPodouble a, TPodouble b)
{
return sqrt((a.x - b.x) * (a.x - b.x) + (a.y - b.y) * (a.y - b.y));
}
int ISSAME(TPodouble a, TPodouble b)
{
if(fabs(a.x - b.x) > eps) return 0;
if(fabs(a.y - b.y) > eps) return 0;
return 1;
}
void between(TPodouble s1, TPodouble e1, TPodouble s2, TPodouble e2)
{
double d1 = dis(s1, s2);
double d2 = dis(s1, e2);
double d3 = dis(e1, s2);
double d4 = dis(e1, e2);
double d5 = dis(s1, e1);
double d6 = dis(s2, e2);
if(ISSAME(s1, s2) && fabs(d5 + d6 - d4) < eps) printf("%.3lf %.3lf\n", s1.x, s1.y);
else if(ISSAME(s1, e2) && fabs(d3 - d5 - d6) < eps) printf("%.3lf %.3lf\n", s1.x, s1.y);
else if(ISSAME(e1, s2) && fabs(d2 - d5 - d6) < eps) printf("%.3lf %.3lf\n", s2.x, s2.y);
else if(ISSAME(e1, e2) && fabs(d1 - d5 - d6) < eps) printf("%.3lf %.3lf\n", e2.x, e2.y);
else printf("Too much points in common\n");
}
int main()
{
TPodouble s1, e1, s2, e2;
TLine l1, l2;
while(scanf("%lf%lf%lf%lf%lf%lf%lf%lf", &s1.x, &s1.y, &e1.x, &e1.y, &s2.x, &s2.y, &e2.x, &e2.y) != EOF){
if(fabs((s1.y - e1.y) * (s2.x - e2.x) - (s1.x - e1.x) * (s2.y - e2.y)) < eps){
l1 = lineFromSegment(s1, e1);
l2.c = -(l1.a * s2.x + l1.b * s2.y);
if(fabs((l1.c - l2.c) / sqrt((double)l1.a * l1.a + l1.b * l1.b)) > 0){
printf("No point in common\n");
}
else {
if(!isIntersected(s1, e1, s2, e2)) printf("No point in common\n");
else between(s1, e1, s2, e2);
}
}
else {
if(isIntersected(s1, e1, s2, e2)){
l1 = lineFromSegment(s1, e1);
l2 = lineFromSegment(s2, e2);
Linedoubleer(l1, l2);
}
else printf("No point in common\n");
}
}
return 0;
}
(26)简单多边形的核
/*多边形的核*/
#include <stdio.h>
#include <math.h>
#define Maxn 3005
const double eps = 1e-10;
typedef struct TPodouble
{
double x;
double y;
}TPoint;
typedef struct TPolygon
{
TPoint p[Maxn];
int n;
};
typedef struct TLine
{
double a, b, c;
}TLine;
bool same(TPoint p1, TPoint p2)
{
if(p1.x != p2.x) return false;
if(p1.y != p2.y) return false;
return true;
}
double multi(TPoint p1, TPoint p2, TPoint p0)
{
//求矢量[p0, p1], [p0, p2]的叉积
//p0是顶点
return (p1.x - p0.x) * (p2.y - p0.y) - (p2.x - p0.x) * (p1.y - p0.y);
//若结果等于0，则这三点共线
//若结果大于0，则p0p2在p0p1的逆时针方向
//若结果小于0，则p0p2在p0p1的顺时针方向
}
TLine lineFromSegment(TPoint p1, TPoint p2)
{
//线段所在直线,返回直线方程的三个系统
TLine tmp;
tmp.a = p2.y - p1.y;
tmp.b = p1.x - p2.x;
tmp.c = p2.x * p1.y - p1.x * p2.y;
return tmp;
}
TPoint LineInter(TLine l1, TLine l2)
{
//求两直线得交点坐标
TPoint tmp;
double a1 = l1.a;
double b1 = l1.b;
double c1 = l1.c;
double a2 = l2.a;
double b2 = l2.b;
double c2 = l2.c;
//注意这里b1 = 0
if(fabs(b1) < eps){
tmp.x = -c1 / a1;
tmp.y = (-c2 - a2 * tmp.x) / b2;
}
else{
tmp.x = (c1 * b2 - b1 * c2) / (b1 * a2 - b2 * a1);
tmp.y = (-c1 - a1 * tmp.x) / b1;
}
return tmp;
}
TPolygon Cut_polygon(TPoint p1, TPoint p2, TPolygon polygon)
{
TPolygon new_polygon;
TPoint interp;
TLine l1, l2;
int i, j;
double t1, t2;
new_polygon.n = 0;
for(i = 0;i <= polygon.n - 1;i++){
t1 = multi(p2, polygon.p[i], p1);
t2 = multi(p2, polygon.p[i + 1], p1);
if(fabs(t1) < eps || fabs(t2) < eps){
if(fabs(t1) < eps) new_polygon.p[new_polygon.n++] = polygon.p[i];
if(fabs(t2) < eps) new_polygon.p[new_polygon.n++] = polygon.p[i + 1];
}
else if(t1 < 0 && t2 < 0){
new_polygon.p[new_polygon.n++] = polygon.p[i];
new_polygon.p[new_polygon.n++] = polygon.p[i + 1];
}
else if(t1 * t2 < 0){
l1 = lineFromSegment(p1, p2);
l2 = lineFromSegment(polygon.p[i], polygon.p[i + 1]);
interp = LineInter(l1, l2);
if(t1 < 0) {
new_polygon.p[new_polygon.n++] = polygon.p[i];
new_polygon.p[new_polygon.n++] = interp;
}
else {
new_polygon.p[new_polygon.n++] = interp;
new_polygon.p[new_polygon.n++] = polygon.p[i + 1];
}
}
}
polygon.n = 0;
if(new_polygon.n == 0) return polygon;
polygon.p[polygon.n++] = new_polygon.p[0];
for(i = 1;i < new_polygon.n;i++){
if(!same(new_polygon.p[i], new_polygon.p[i - 1])){
polygon.p[polygon.n++] = new_polygon.p[i];
}
}
if(polygon.n != 1 && same(polygon.p[polygon.n - 1], polygon.p[0])) polygon.n--;
polygon.p[polygon.n] = polygon.p[0];
return polygon;
}
double polygonArea(TPolygon p)
{
//已知多边形各顶点的坐标，求其面积
int i, n;
double area;
n = p.n;
area = 0;
for(i = 1;i <= n;i++){
area += (p.p[i - 1].x * p.p[i % n].y - p.p[i % n].x * p.p[i - 1].y);
}
return area / 2;
}
void ChangeClockwise(TPolygon &polygon)
{
TPoint tmp;
int i;
for(i = 0;i <= (polygon.n - 1) / 2;i++){
tmp = polygon.p[i];
polygon.p[i] = polygon.p[polygon.n - 1 - i];
polygon.p[polygon.n - 1 - i] = tmp;
}
}
int main()
{
int test, i, j;
double area;
TPolygon polygon, new_polygon;
scanf("%d", &test);
while(test--){
scanf("%d", &polygon.n);
for(i = 0;i <= polygon.n - 1;i++){
scanf("%lf%lf", &polygon.p[i].x, &polygon.p[i].y);
}
/*若是逆时针转化为顺时针*/
if(polygonArea(polygon) > 0) ChangeClockwise(polygon);
polygon.p[polygon.n] = polygon.p[0];
new_polygon = polygon;
for(i = 0;i <= polygon.n - 1;i++){
new_polygon = Cut_polygon(polygon.p[i], polygon.p[i + 1], new_polygon);
}
area = polygonArea(new_polygon);
if(area < 0) printf("%.2lf\n", -area);
else printf("%.2lf\n", area);
}
return 0;
}
(27)线段重叠+投影(pku_1375)
#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#define eps 1e-8
#define pi acos(-1)
typedef struct SEG
{
double l, r;
}SEG;
typedef struct TPoint
{
double x, y;
}TPoint;
typedef struct TLine
{
double a, b, c;
}TLine;
typedef struct Circle
{
TPoint c;
double r;
}Circle;
double distance(double x1, double y1, double x2, double y2)
{
return sqrt((x1 - x2) * (x1 - x2) + (y1 - y2) * (y1 - y2));
}
TPoint LineInter(TLine l1, TLine l2)
{
//求两直线得交点坐标
TPoint tmp;
double a1 = l1.a;
double b1 = l1.b;
double c1 = l1.c;
double a2 = l2.a;
double b2 = l2.b;
double c2 = l2.c;
//注意这里b1 = 0
if(fabs(b1) < eps){
tmp.x = -c1 / a1;
tmp.y = (-c2 - a2 * tmp.x) / b2;
}
else{
tmp.x = (c1 * b2 - b1 * c2) / (b1 * a2 - b2 * a1);
tmp.y = (-c1 - a1 * tmp.x) / b1;
}
return tmp;
}
TLine lineFromSegment(TPoint p1, TPoint p2)
{
//线段所在直线,返回直线方程的三个系统
TLine tmp;
tmp.a = p2.y - p1.y;
tmp.b = p1.x - p2.x;
tmp.c = p2.x * p1.y - p1.x * p2.y;
return tmp;
}
int cmp(const void *a, const void *b)
{
SEG *c = (SEG *)a;
SEG *d = (SEG *)b;
if(c->l < d->l) return -1;
else return 1;
}
int main()
{
int n, i, j;
TPoint b, inter;
double lbc, lcinter, a, bangle;
double start, end;
TLine l1;
Circle circle[505];
SEG seg[505];
while(scanf("%d", &n) && n){
scanf("%lf%lf", &b.x, &b.y);
for(i = 0;i < n;i++){
scanf("%lf%lf%lf", &circle[i].c.x, &circle[i].c.y, &circle[i].r);
if(fabs(b.x - circle[i].c.x) < eps) a = pi / 2;
else a = atan((circle[i].c.y - b.y) / (circle[i].c.x - b.x));
if(a < 0) a += pi;
lbc = distance(b.x, b.y, circle[i].c.x, circle[i].c.y);
bangle = asin(circle[i].r / lbc);
l1 = lineFromSegment(b, circle[i].c);
inter.x = -l1.c / l1.a;
inter.y = 0.0;
lcinter = distance(inter.x, inter.y, b.x, b.y);
seg[i].l = inter.x - lcinter * circle[i].r / lbc / sin(a - bangle);
seg[i].r = inter.x + lcinter * circle[i].r / lbc / sin(pi - a - bangle);
}
qsort(seg, n, sizeof(seg[0]), cmp);
start = seg[0].l;
end = seg[0].r;
for(i = 1;i < n;i++){
if(seg[i].l > end){
printf("%.2lf %.2lf\n", start, end);
start = seg[i].l;
end = seg[i].r;
}
else {
if(seg[i].r > end) end = seg[i].r;
}
}
printf("%.2lf %.2lf\n\n", start, end);
}
return 0;
}
(28)二分+圆的参数方程(pku_2600)
/*二分+圆的参数方程(pku_2600)*/
#include <stdio.h>
#include <math.h>
const double eps = 1e-4;
const double pi = acos(-1.0);
struct TPoint
{
double x, y;
}p[60], a[60];
double angle[60];
double multi(TPoint p1, TPoint p2, TPoint p0)
{
return (p1.x - p0.x) * (p2.y - p0.y) - (p2.x - p0.x) * (p1.y - p0.y);
}
TPoint fine_a2(TPoint a1, TPoint m, double angle1)
{
TPoint a2;
double r, angle2, angle3;
r = sqrt((a1.x - m.x) * (a1.x - m.x) + (a1.y - m.y) * (a1.y - m.y));
angle2 = acos((a1.x - m.x) / r);
if(a1.y < m.y) {
if(angle2 <= pi / 2) angle2 = -angle2;
if(angle2 > pi / 2) angle2 = 3 * pi / 2 - (angle2 - pi / 2);
}
angle3 = angle2 - angle1;
a2.x = m.x + r * cos(angle3);
a2.y = m.y + r * sin(angle3);
if(multi(m, a2, a1) < 0) return a2;
angle3 = angle2 + angle1;
a2.x = m.x + r * cos(angle3);
a2.y = m.y + r * sin(angle3);
if(multi(m, a2, a1) < 0) return a2;
}
int main()
{
int n, i, j;
while(scanf("%d", &n) != EOF){
for(i = 0;i < n;i++){
scanf("%lf%lf", &p[i].x, &p[i].y);
}
for(i = 0;i < n;i++){
scanf("%lf", &angle[i]);
angle[i] = angle[i] * pi / 180;
}
a[0].x = 0;
a[0].y = 0;
while(1){
for(i = 1;i <= n;i++){
a[i] = fine_a2(a[i - 1], p[i - 1], angle[i - 1]);
}
if(fabs(a[n].x - a[0].x) <= eps
&& fabs(a[n].y - a[0].y) <= eps) break;
else {
a[0].x = (a[0].x + a[n].x) / 2;
a[0].y = (a[0].y + a[n].y) / 2;
}
}
for(i = 0;i < n;i++){
printf("%.0lf %.0lf\n", a[i].x, a[i].y);
}
}
return 0;
}
(29)Pick公式
//A = b / 2 + i -1 其中 b 与 i 分別表示在边界上及內部的格子点之个數
//http://www.hwp.idv.tw/bbs1/htm/%A6V%B6q%B7L%BFn%A4%C0/%A6V%B6q%B7L%BFn%A4%C0.htm
// http://acm.pku.edu.cn/JudgeOnline/problem?id=2954
#include <stdio.h>
#include <stdlib.h>
typedef struct TPoint
{
int x;
int y;
}TPoint;
typedef struct TLine
{
int a, b, c;
}TLine;
int triangleArea(TPoint p1, TPoint p2, TPoint p3)
{
//已知三角形三个顶点的坐标，求三角形的面积
int k = p1.x * p2.y + p2.x * p3.y + p3.x * p1.y
- p2.x * p1.y - p3.x * p2.y - p1.x * p3.y;
if(k < 0) return -k;
else return k;
}
TLine lineFromSegment(TPoint p1, TPoint p2)
{
//线段所在直线,返回直线方程的三个系统
TLine tmp;
tmp.a = p2.y - p1.y;
tmp.b = p1.x - p2.x;
tmp.c = p2.x * p1.y - p1.x * p2.y;
return tmp;
}
void swap(int &a, int &b)
{
int t;
t = a;
a = b;
b = t;
}
int Count(TPoint p1, TPoint p2)
{
int i, sum = 0, y;
TLine l1 = lineFromSegment(p1, p2);
if(l1.b == 0) return abs(p2.y - p1.y) + 1;
if(p1.x > p2.x) swap(p1.x, p2.x); //这里没有交换WA两次
for(i = p1.x;i <= p2.x;i++){
y = -l1.c - l1.a * i;
if(y % l1.b == 0) sum++;
}
return sum;
}
int main()
{
//freopen("in.in", "r", stdin);
//freopen("OUT.out", "w", stdout);
TPoint p1, p2, p3;
while(scanf("%d%d%d%d%d%d", &p1.x, &p1.y, &p2.x, &p2.y, &p3.x, &p3.y) != EOF){
if(p1.x == 0 && p1.y == 0 && p2.x == 0 && p2.y == 0 && p3.x == 0 && p3.y == 0) break;
int A = triangleArea(p1, p2, p3);//A为面积的两倍
int b = 0;
int i;
b = Count(p1, p2) + Count(p1, p3) + Count(p3, p2) - 3;//3个顶点多各多加了一次
//i = A / 2- b / 2 + 1;
i = (A - b) / 2 + 1;
printf("%d\n", i);
}
return 0;
}
(30)根据经度纬度求球面距离
/*
假设地球是球体，
设地球上某点的经度为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;
}
(31)两圆切线的交点
/*zju_1199_切线交点.cpp*/
#include <stdio.h>
#include <math.h>
#define eps 1e-6
#define pi acos(-1.0)
struct TPoint
{
double x, y;
};
struct TLine
{
double a, b, c;
};
struct TCircle
{
TPoint c;
double r;
};
double distanc(TPoint p1, TPoint p2)
{
return sqrt((p1.x - p2.x) * (p1.x - p2.x) + (p1.y - p2.y) * (p1.y - p2.y));
}
void swap(TCircle &c1, TCircle &c2)
{
TCircle c;
c = c1;
c1 = c2;
c2 = c;
}
int IsOnSegment(TPoint p, TPoint p1, TPoint p2)
{
double d1 = distanc(p1, p);
double d2 = distanc(p, p2);
double d3 = distanc(p1, p2);
if(fabs(d1 + d2 - d3) < eps) return 1;
else return 0;
}
int main()
{
//freopen("in.in", "r", stdin);
//freopen("out.out", "w", stdout);
int test;
TLine l1;
TPoint p1;
TCircle c1, c2;
double k, a;
scanf("%d", &test);
while(test--){
scanf("%lf%lf%lf", &c1.c.x, &c1.c.y, &c1.r);
scanf("%lf%lf%lf", &c2.c.x, &c2.c.y, &c2.r);
if(fabs(c1.r - c2.r) < eps){
printf("Impossible.\n");
continue;
}
double d = distanc(c1.c, c2.c);
//if(fabs(c2.r- c1.r) - d > 0) 这样就wa
if(fabs(c2.r- c1.r) - d > -eps){
printf("Impossible.\n");
continue;
}
if(c2.r < c1.r) swap(c1, c2);
double x = c1.r *d / (c2.r - c1.r);
if(fabs(c2.c.x - c1.c.x) < eps){
a = pi / 2;
}
k = (c2.c.y - c1.c.y) / (c2.c.x - c1.c.x);
a = atan(k);
p1.x = c1.c.x + x * cos(a);
p1.y = c1.c.y + x * sin(a);
if(!IsOnSegment(p1, c1.c, c2.c)){
printf("%.2lf %.2lf\n", p1.x, p1.y);
}
else {
a += pi;
p1.x = c1.c.x + x * cos(a);
p1.y = c1.c.y + x * sin(a);
if(!IsOnSegment(p1, c1.c, c2.c)){
printf("%.2lf %.2lf\n", p1.x, p1.y);
}
}
}
return 0;
}
(31)两圆切线的交点(32)线段与三角形的交(usaco_3.4)
/*usaco3_4*/
#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#define Maxn 205
#define eps 1e-6
typedef struct TPoint
{
double x, y;
}TPoint;
typedef struct SEG
{
TPoint p1, p2;
int p1num, p2num;
}SEG;
typedef struct ANGLE
{
double l, r;
}ANGLE;
typedef struct TLine
{
double a, b, c;
}TLine;
struct ANS
{
int left, right;
};
double max(double x, double y)
{
//比较两个数的大小，返回大的数
if(x > y) return x;
else return y;
}
double min(double x, double y)
{
//比较两个数的大小，返回小的数
if(x < y) return x;
else return y;
}
int IsSame(TPoint p1, TPoint p2)
{
if(fabs(p1.x - p2.x) > eps) return 0;
if(fabs(p1.y - p2.y) > eps) return 0;
return 1;
}
double distance(TPoint p1, TPoint p2)
{
//计算平面上两个点之间的距离
return sqrt((p1.x - p2.x) * (p1.x - p2.x) + (p1.y - p2.y) * (p1.y - p2.y));
}
double multi(TPoint p1, TPoint p2, TPoint p0)
{
//求矢量[p0, p1], [p0, p2]的叉积
//p0是顶点
return (p1.x - p0.x) * (p2.y - p0.y) - (p2.x - p0.x) * (p1.y - p0.y);
//若结果等于0，则这三点共线
//若结果大于0，则p0p2在p0p1的逆时针方向
//若结果小于0，则p0p2在p0p1的顺时针方向
}
int parallel(TPoint s1, TPoint e1, TPoint s2, TPoint e2)
{
if(fabs((s1.y - e1.y) * (s2.x - e2.x)
- (s2.y - e2.y) * (s1.x - e1.x)) < eps) return 1;
else return 0;
}
int cmp(const void *a, const void *b)
{
ANGLE *c = (ANGLE *)a;
ANGLE *d = (ANGLE *)b;
if(c->l < d->l) return -1;
else return 1;
}
int cmp1(const void *a, const void *b)
{
struct ANS *c = (struct ANS *)a;
struct ANS *d = (struct ANS *)b;
if(c->right < d->right) return -1;
else if(c->right > d->right) return 1;
else {
if(c->left < d->left) return -1;
else return 1;
}
}
int isIntersected(TPoint s1, TPoint e1, TPoint s2, TPoint e2)
{
//判断线段是否相交
//1.快速排斥试验判断以两条线段为对角线的两个矩形是否相交
//2.跨立试验
if(
(max(s1.x, e1.x) >= min(s2.x, e2.x)) &&
(max(s2.x, e2.x) >= min(s1.x, e1.x)) &&
(max(s1.y, e1.y) >= min(s2.y, e2.y)) &&
(max(s2.y, e2.y) >= min(s1.y, e1.y)) &&
(multi(s2, e1, s1) * multi(e1, e2, s1) >= 0) &&
(multi(s1, e2, s2) * multi(e2, e1, s2) >= 0)
) return 1;
return 0;
}
TPoint LineInter(TLine l1, TLine l2)
{
//求两直线得交点坐标
TPoint tmp;
double a1 = l1.a;
double b1 = l1.b;
double c1 = l1.c;
double a2 = l2.a;
double b2 = l2.b;
double c2 = l2.c;
//注意这里b1 = 0
if(fabs(b1) < eps){
tmp.x = -c1 / a1;
tmp.y = (-c2 - a2 * tmp.x) / b2;
}
else{
tmp.x = (c1 * b2 - b1 * c2) / (b1 * a2 - b2 * a1);
tmp.y = (-c1 - a1 * tmp.x) / b1;
}
return tmp;
}
TLine lineFromSegment(TPoint p1, TPoint p2)
{
//线段所在直线,返回直线方程的三个系统
TLine tmp;
tmp.a = p2.y - p1.y;
tmp.b = p1.x - p2.x;
tmp.c = p2.x * p1.y - p1.x * p2.y;
return tmp;
}
double AreaOfTriangle(TPoint p1, TPoint p2, TPoint p3)
{
return fabs((p1.x - p3.x) * (p2.y - p3.y)
- (p1.y - p3.y) * (p2.x - p3.x)) / 2;
}
int isPointInTriangle(TPoint p, TPoint p1, TPoint p2, TPoint p3)
{
double area1 = 0, area2 = 0;
area1 = AreaOfTriangle(p1, p2, p)
+ AreaOfTriangle(p1, p3, p)
+ AreaOfTriangle(p2, p3, p);
area2 = AreaOfTriangle(p1, p2, p3);
if(fabs(area1 - area2) < eps) return 1;
else return 0;
}
double get_angle(TPoint p1, TPoint p2, TPoint point)
{
//这里可以用余弦定理做
TPoint p3, p4;
double tmpa, tmpb, tmpc, tmpd;
p3.x = p1.x - point.x;
p3.y = p1.y - point.y;
p4.x = p2.x - point.x;
p4.y = p2.y - point.y;
tmpa = p3.x * p4.x + p3.y * p4.y;
tmpb = distance(p1, point);
tmpc = distance(p2, point);
tmpd = tmpa / tmpb / tmpc;
if(tmpd > 1) tmpd = 1;
if(tmpd < -1) tmpd = -1;
return acos(tmpd);
}
int find_point_2(TPoint point, SEG seg1, SEG seg2, SEG &segtmp)
{
/*找出线段seg2在point，和线段seg1组成的三角形的两个交点*/
TLine l1, l2;
TPoint p[3];
int n = 0;
if(parallel(seg2.p1, seg2.p2, point, seg1.p1)
|| parallel(seg2.p1, seg2.p2, point, seg1.p2)){ return 0;}
if(isPointInTriangle(seg2.p1, point, seg1.p1, seg1.p2)
&& isPointInTriangle(seg2.p2, point, seg1.p1, seg1.p2)){
segtmp.p1 = seg2.p1;
segtmp.p2 = seg2.p2;
return 1;
}
if(isIntersected(point, seg1.p1, seg2.p1, seg2.p2)){
l1 = lineFromSegment(point, seg1.p1);
l2 = lineFromSegment(seg2.p1, seg2.p2);
p[n] = LineInter(l1, l2);
n++;
}
if(isIntersected(point, seg1.p2, seg2.p1, seg2.p2)){
l1 = lineFromSegment(point, seg1.p2);
l2 = lineFromSegment(seg2.p1, seg2.p2);
p[n] = LineInter(l1, l2);
n++;
}
if(isIntersected(seg1.p1, seg1.p2, seg2.p1, seg2.p2)){
l1 = lineFromSegment(seg1.p1, seg1.p2);
l2 = lineFromSegment(seg2.p1, seg2.p2);
p[n] = LineInter(l1, l2);
n++;
}
if(n == 3){
if(!IsSame(p[0], p[1])){
segtmp.p1 = p[0];
segtmp.p2 = p[1];
}
else {
segtmp.p1 = p[0];
segtmp.p2 = p[2];
}
return 1;
}
else if(n == 0) return 0;
else if(n == 2){
if(!IsSame(p[0], p[1])){
segtmp.p1 = p[0];
segtmp.p2 = p[1];
return 1;
}
else {
if(isPointInTriangle(seg2.p1, point, seg1.p1, seg1.p2)){
segtmp.p1 = p[0];
segtmp.p2 = seg2.p1;
if(IsSame(segtmp.p2, segtmp.p1)) return 0;
return 1;
}
if(isPointInTriangle(seg2.p2, point, seg1.p1, seg1.p2)){
segtmp.p1 = p[0];
segtmp.p2 = seg2.p2;
if(IsSame(segtmp.p2, segtmp.p1)) return 0;
return 1;
}
}
return 0;
}
else if(n == 1){
segtmp.p1 = p[0];
if(isPointInTriangle(seg2.p1, point, seg1.p1, seg1.p2)
&& !IsSame(seg2.p1, segtmp.p1)){
segtmp.p2 = seg2.p1;
return 1;
}
else if(isPointInTriangle(seg2.p2, point, seg1.p1, seg1.p2)
&& !IsSame(seg2.p2, segtmp.p1)){
segtmp.p2 = seg2.p2;
return 1;
}
return 0;
}
}
int main()
{
freopen("fence4.in", "r", stdin);
freopen("fence4.out", "w", stdout);
int i, j, pn, segn, anglen, ansn;
int flag1, flag2, flag3, tag;
double max_angle, angle1, angle2, end;
ANGLE angle[Maxn];
TPoint p[Maxn];
ANS ans[Maxn];
TPoint point, mid;
SEG seg[Maxn], segtmp;
scanf("%d", &pn);
segn = 0;
scanf("%lf%lf", &point.x, &point.y);
scanf("%lf%lf", &p[0].x, &p[0].y);
for(i = 1;i < pn;i++){
scanf("%lf%lf", &p[i].x, &p[i].y);
seg[segn].p1 = p[i - 1];
seg[segn].p2 = p[i];
seg[segn].p1num = i - 1;
seg[segn++].p2num = i;
}
seg[segn].p1 = p[pn - 1];
seg[segn].p2 = p[0];
seg[segn].p1num = pn - 1;
seg[segn++].p2num = 0;
ansn = 0;
for(i = 0;i <= segn - 1;i++){
max_angle = get_angle(seg[i].p1, seg[i].p2, point);
if(fabs(max_angle) < eps) continue;
anglen = 0;
for(j = 0;j <= segn - 1;j++){
if(i == j) continue;
if(find_point_2(point, seg[i], seg[j], segtmp)){
if(IsSame(segtmp.p1, segtmp.p2))continue;
angle1 = get_angle(segtmp.p1, seg[i].p2, point);
angle2 = get_angle(segtmp.p2, seg[i].p2, point);
if(angle1 < angle2){
angle[anglen].l = angle1;
angle[anglen++].r = angle2;
}
else {
angle[anglen].l = angle2;
angle[anglen++].r = angle1;
}
}
}
qsort(angle, anglen, sizeof(angle[0]), cmp);
tag = 0;
if(anglen == 0 || angle[0].l > eps) tag = 1;
else {
end = angle[0].r;
for(j = 1;j < anglen;j++){
if(angle[j].l - end > eps){
tag = 1;
break;
}
else if(angle[j].r - end > eps) end = angle[j].r;
}
if(max_angle- end > eps) tag = 1;
}
if(tag == 1){
if(seg[i].p1num < seg[i].p2num) {
ans[ansn].left = seg[i].p1num;
ans[ansn++].right = seg[i].p2num;
}
else {
ans[ansn].left = seg[i].p2num;
ans[ansn++].right = seg[i].p1num;
}
}
}
if(ansn == 0){
printf("NOFENCE\n");
return 0;
}
qsort(ans, ansn, sizeof(ans[0]), cmp1);
printf("%d\n", ansn);
for(i = 0;i < ansn;i++){
printf("%.0lf %.0lf %.0lf %.0lf\n", p[ans[i].left].x,
p[ans[i].left].y, p[ans[i].right].x, p[ans[i].right].y);
}
return 0;
}
/*
* Computational_Geometry
* summarized by ZouYu
* at 2006.9.24
*/
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
//****************************************常量定义*****************************************//
const int MAX = 10000; //注意数据范围！！！
const double EPS = 1e-10;
const double DOUBLE_INF = 1e300;
const double PI = 3.1415926535897932384626433;
//***************************************数据类型定义***************************************//
typedef struct{ //点
double x, y;
}Point, Vector;
typedef struct{ //多边形 (n个顶点)
Point vertex[MAX+1];
int num;
}PSet, Polygon;
typedef struct{ //三角
Point vertex[4];
}Triangle;
typedef struct{ //圆
double r;
Point center;
}Circle;
typedef struct{ //直线
double A, B, C;
}Line;
//****************************************数值运算*****************************************//
//判断某数是否等于，亦用于判断两数是否相等
bool isZero(double a){
return a<EPS && a>-EPS;
}
//某数为-0.0时返回.0
double CheckZero(double zero){
if(zero<EPS && zero>-EPS)
return 0.0;
else return zero;
}
//最大值
double Max(double a, double b){
return a>b?a:b;
}
//最小值
double Min(double a, double b){
return a<b?a:b;
}
bool Equ(double a, double b){
return isZero(a-b);
}
//****************************************点的算法*****************************************//
//判断两点是否相同
bool SamePoint(Point a, Point b){
return isZero(a.x-b.x) && isZero(a.y-b.y);
}
//平面上两点之间距离
double Distance(Point a, Point b){
return sqrt((a.x-b.x)*(a.x-b.x)+((a.y-b.y)*(a.y-b.y)));
}
//****************************************向量算法*****************************************//
//计算向量的点积v1 * v2
double DotProduct(Vector v1, Vector v2){
return v1.x*v2.x+v1.y*v2.y;
}
//计算向量的叉积v1 x v2
double CrossProduct(Vector v1, Vector v2){
return v1.x*v2.y-v1.y*v2.x;
}
//矢量点乘
double DotProduct(Point p_begin, Point p_end, Point q_begin, Point q_end){
Vector v1, v2;
v1.x = p_end.x-p_begin.x;
v1.y = p_end.y-p_begin.y;
v2.x = q_end.x-q_begin.x;
v2.y = q_end.y-q_begin.y;
return DotProduct(v1, v2);
}
//矢量叉乘
double CrossProduct(Point p_begin, Point p_end, Point q_begin, Point q_end){
Vector v1, v2;
v1.x = p_end.x-p_begin.x;
v1.y = p_end.y-p_begin.y;
v2.x = q_end.x-q_begin.x;
v2.y = q_end.y-q_begin.y;
return CrossProduct(v1, v2);
}
//****************************************点与线*******************************************//
//点p到线段ab所在直线的距离
double DistanceToLine(Point p, Point a, Point b){
Vector v1,v2;
v1.x = a.x-p.x; v1.y = a.y-p.y;
v2.x = b.x-a.x; v2.y = b.y-a.y;
return fabs(CrossProduct(v1,v2)/Distance(a,b));
}
//判断点是否在线段上(包括端点)
bool isPointOnSegment(Point p, Point p1, Point p2){
Vector v1,v2;
v1.x = p.x-p1.x; v1.y = p.y-p1.y;
v2.x = p2.x-p.x; v2.y = p2.y-p.y;
if(isZero(CrossProduct(v1, v2))==false) return false;
if(((p.x>p1.x)&&(p.x>p2.x))||((p.x<p1.x)&&(p.x<p2.x))) return false;
if(((p.y>p1.y)&&(p.y>p2.y))||((p.y<p1.y)&&(p.y<p2.y))) return false;
return true;
}
//点到直线的垂足
void GetFoot(Point p, Point a, Point b, Point &foot){
double t;
if(SamePoint(a, b)) foot=a;
else{
t = (-(a.x-p.x)*(b.x-a.x)-(a.y-p.y)*(b.y-a.y))/
((b.x-a.x)*(b.x-a.x)+(b.y-a.y)*(b.y-a.y));
foot.x = a.x+t*(b.x-a.x);
foot.y = a.y+t*(b.y-a.y);
}
}
//点到线段的距离
double DistanceToSegment(Point p, Point a, Point b){
Point foot;
double tresa, tresb;
GetFoot(p, a, b, foot);
if(isPointOnSegment(foot, a, b)) return Distance(p, foot);
else{
tresa = Distance(p, a);
tresb = Distance(p, b);
if(tresa>tresb) return tresb;
else return tresa;
}
}
//点关于直线的对称点 o:原始点 t:对称点
void symmetricalPointofLine(Point o, Line l, Point* t){
double d;
d = l.A*l.A+l.B*l.B;
t->x = (l.B*l.B*o.x - l.A*l.A*o.x - 2*l.A*l.B*o.y - 2*l.A*l.C)/d;
t->y = (l.A*l.A*o.y - l.B*l.B*o.y - 2*l.A*l.B*o.x - 2*l.B*l.C)/d;
}
//****************************************线的算法*****************************************//
//计算两直线ab,cd的交点,结果存于p中
bool LinesIntersect(Point *p, Point a, Point b, Point c, Point d){
double mul, t;
Vector v1, v2;
v1.x = b.x-a.x; v1.y = b.y-a.y;
v2.x = d.x-c.x; v2.y = d.y-c.y;
mul = CrossProduct(v1, v2);
if(Equ(mul, 0)) return false;
else{
t = ((c.x-a.x)*v2.y+(a.y-c.y)*v2.x)/mul;
p->x = a.x+v1.x*t;
p->y = a.y+v1.y*t;
return true;
}
}
//线段所在直线(保证A>=0)
void LineFromSegment(Line *l, Point a, Point b)
{
int sign = 1;
l->A = a.y-b.y;
if(l->A < 0){
sign = -1;
l->A *= sign;
}
l->B = sign*(b.x-a.x);
l->C = sign*(a.x*b.y-a.y*b.x);
}
//判断线段是否相交(端点接触不算相交;跨立实验中'>'换成'>='时，端点接触也算)
bool isSegmentIntersect(Point s1, Point e1, Point s2, Point e2){
return (Max(s1.x, e1.x)>=Min(s2.x, e2.x)) &&
(Max(s2.x, e2.x)>=Min(s1.x, e1.x)) &&
(Max(s1.y, e1.y)>=Min(s2.y, e2.y)) &&
(Max(s2.y, e2.y)>=Min(s1.y, e1.y)) && //快速排斥试验(点稀疏时会加快速度)
(CrossProduct(s1, s2, s1, e1)*CrossProduct(s1, e1, s1, e2)>0) &&
(CrossProduct(s2, s1, s2, e2)*CrossProduct(s2, e2, s2, e1)>0);//跨立实验
}
//判断线段是否与直线相交:线段(s, e)直线(lf, lt)
bool isIntercectWithLine(Point s, Point e, Point lf, Point lt){
return CrossProduct(lf, s, lf, lt)*CrossProduct(lf, lt, lf, e)>0;
}
//线段到线段的距离
double SegmentToSegment(Point a, Point b, Point c, Point d){
double res;
if(isSegmentIntersect(a, b, c, d)) return 0;
else{
res = DOUBLE_INF;
res = Min(res, DistanceToSegment(a, c, d));
res = Min(res, DistanceToSegment(b, c, d));
res = Min(res, DistanceToSegment(c, a, b));
res = Min(res, DistanceToSegment(d, a, b));
return res;
}
}
// 根据直线解析方程返回直线的斜率k,水平线返回,竖直线返回DOUBLE_INF
double slope(Line l)
{
if(fabs(l.A) < EPS) return 0;
if(fabs(l.B) < EPS) return DOUBLE_INF;
return -(l.A/l.B);
}
// 返回直线的倾斜角alpha ( 0 - pi)
double alpha(Line l)
{
if(fabs(l.A)< EPS)return 0;
if(fabs(l.B)< EPS)return PI/2;
double k=slope(l);
if(k>0) return atan(k);
else return PI+atan(k);
}
//判断线段和圆相交
bool SegmentIntersectCircle(Point a, Point b, Circle c){
double l1, l2, l3;
double cosa, h;
l1 = Distance(a, c.center);
l2 = Distance(b, c.center);
l3 = Distance(a, b);
if(l1<c.r && l2<c.r) return false;
if((l1-c.r)*(l2-c.r) < EPS) return true;
cosa = (l1*l1+l3*l3-l2*l2)/2/l1/l3;
h = l1*sqrt(1-cosa*cosa);
if(h-c.r>EPS) return false;
if(l1<l2) cosa = l1*l1+l3*l3-l2*l2;
else cosa = l3*l3+l2*l2-l1*l1;
if(cosa<-EPS) return false;
else return true;
}
//****************************************三角形*******************************************//
//海伦公式求三角形面积:已知三角形三边长度
double HeroArea(double a, double b, double c){
double s = (a+b+c)/2.0;
return sqrt(s*(s-a)*(s-b)*(s-c));
}
//已知三角形三边长度以及外接圆半径
double RArea(double a, double b, double c, double R){
return a*b*c/(4.0*R);
}
//已知三角形三顶点坐标求面积
double TriangleArea(Triangle t){
double area = 0.0;
area += t.vertex[0].x*t.vertex[1].y - t.vertex[0].y*t.vertex[1].x;
area += t.vertex[1].x*t.vertex[2].y - t.vertex[1].y*t.vertex[2].x;
area += t.vertex[2].x*t.vertex[0].y - t.vertex[2].y*t.vertex[0].x;
return 0.5*fabs(area);
}
//三角形的外接圆
void CircumCircle(Circle *circle, Triangle t){
double a = Distance(t.vertex[0], t.vertex[1]);
double b = Distance(t.vertex[1], t.vertex[2]);
double c = Distance(t.vertex[2], t.vertex[0]);
circle->r = a*b*c/(4.0*TriangleArea(t));
double ax = t.vertex[0].x, ay = t.vertex[0].y;
double bx = t.vertex[1].x, by = t.vertex[1].y;
double cx = t.vertex[2].x, cy = t.vertex[2].y;
double c1 = (ax*ax+ay*ay-bx*bx-by*by)/2.0;
double c2 = (ax*ax+ay*ay-cx*cx-cy*cy)/2.0;
circle->center.x = (c1*(ay-cy)-c2*(ay-by))/((ax-bx)*(ay-cy)-(ax-cx)*(ay-by));
circle->center.y = (c1*(ax-cx)-c2*(ax-bx))/((ay-by)*(ax-cx)-(ay-cy)*(ax-bx));
}
//三角形的内切圆
void InCircleOfTriangle(Circle *circle, Triangle t){
double a = Distance(t.vertex[0], t.vertex[1]);
double b = Distance(t.vertex[1], t.vertex[2]);
double c = Distance(t.vertex[2], t.vertex[0]);
circle->r = (2.0*TriangleArea(t))/(a+b+c);
double angleA = acos((b*b+c*c-a*a)/(2*b*c));
double angleB = acos((a*a+c*c-b*b)/(2*a*c));
double angleC = acos((b*b+a*a-c*c)/(2*b*a));
double p1 = sin(angleA/2), p2 = sin(angleB/2), p3 = sin(angleC/2);
double ax = t.vertex[0].x, ay = t.vertex[0].y;
double bx = t.vertex[1].x, by = t.vertex[1].y;
double cx = t.vertex[2].x, cy = t.vertex[2].y;
double f1 = 0.5*((circle->r/p2)*(circle->r/p2)-(circle->r/p1)*(circle->r/p1)+ax*ax-bx*bx+ay*ay-by*by);
double f2 = 0.5*((circle->r/p3)*(circle->r/p3)-(circle->r/p1)*(circle->r/p1)+ax*ax-cx*cx+ay*ay-cy*cy);
circle->center.x = (f1*(ay-cy)-f2*(ay-by))/((ax-bx)*(ay-cy)-(ax-cx)*(ay-by));
circle->center.y = (f1*(ax-cx)-f2*(ax-bx))/((ay-by)*(ax-cx)-(ay-cy)*(ax-bx));
}
//用面积判断点是否在三角形内(边界也算在内)
bool isPointInTriangle(Point p, Triangle t){
double area = 0;
Point temp;
int i;
for(i = 0; i < 3; i++)
{
temp.x = t.vertex[i].x; temp.y = t.vertex[i].y;
t.vertex[i].x = p.x; t.vertex[i].y = p.y;
area += TriangleArea(t);
t.vertex[i].x = temp.x; t.vertex[i].y = temp.y;
}
return isZero(area-TriangleArea(t));
}
//用线段拐向判断点是否在三角形内(不包括边界,亦可以修改成包括边界)
bool isPointIntriangle2(Point p, Triangle t){
double k1 = CrossProduct(p, t.vertex[0], p, t.vertex[1]);
double k2 = CrossProduct(p, t.vertex[1], p, t.vertex[2]);
double k3 = CrossProduct(p, t.vertex[2], p, t.vertex[0]);
return k1*k2>0 && k2*k3>0;
}
//三角形垂心
void VerticalPointOfTriangle(Point* vertical, Point A, Point B, Point C){
double xac = A.x-C.x, xcb = C.x-B.x, xba = B.x-A.x;
double yac = A.y-C.y, ycb = C.y-B.y, yba = B.y-A.y;
vertical->x = CheckZero((A.y*A.x*xcb+B.y*B.x*xac+C.y*C.x*xba-yac*ycb*yba)/(A.y*xcb+B.y*xac+C.y*xba));
vertical->y = CheckZero((A.x*A.y*ycb+B.x*B.y*yac+C.x*C.y*yba-xac*xcb*xba)/(A.x*ycb+B.x*yac+C.x*yba));
}
//*****************************************矩形********************************************//
//判断矩形X*Y能否放入矩形A*B中(X,A为各自矩形较长的一边)
bool isRectangleIn(double A, double B, double X, double Y){
double la, lb, xl1, xl2;
if (A >= X && B >= Y) return 1;
if (Y >= B) return 0;
xl1 = sqrt(A * A + B * B);
xl2 = sqrt(X * X + Y * Y);
if (xl2 >= xl1) return 0;
la = (A - sqrt(xl2 * xl2 - B * B)) / 2;
lb = (B - sqrt(xl2 * xl2 - A * A)) / 2;
if (la * la + lb * lb >= Y * Y) return 1;
return 0;
}
//判断矩形X*Y能否放入矩形A*B中(X,A为各自矩形较长的一边)
bool TryPut ( double a , double b , double x , double y ){
double l = sqrt ( x * x + y * y ) , alfa1 , alfa2 , l2;
alfa1 = acos ( b / l );
alfa2 = atan ( y / x ) + alfa1;
if ( alfa2 > PI / 2 ) return false;
l2 = x * sin ( alfa2 ) + y * cos ( alfa2 );
return l2 < a + 1e-10;
}
bool isRectangleIn2 ( double A , double B , double X , double Y ){
if ( A >= X && B >= Y ) return true;
if ( A < X && B < Y ) return false;
return TryPut ( A , B , X , Y );
}
//****************************************多边形*******************************************//
//简单多边形面积,n为顶点个数
double PolygonArea(Polygon poly){
double area = 0;
int n = poly.num;
poly.vertex[n].x = poly.vertex[0].x, poly.vertex[n].y = poly.vertex[0].y;
for(int i = 1; i <= n; i++)
area += poly.vertex[i-1].x*poly.vertex[i].y-poly.vertex[i].x*poly.vertex[i-1].y;
return fabs(area)/2.0;
}
//判断点a是否在简单多边形p内(不考虑落在边界上的情况的话，删掉第一个for循环)
int PointInPolygon(Point a, Polygon p){
for(int onBoard = 0; onBoard < p.num; onBoard++)
if(isPointOnSegment(a,p.vertex[onBoard],p.vertex[(onBoard+1)%p.num]))
return 1;
int i, flag, times;
Point k;
k = p.vertex[0];
for (i = 1; i < p.num; i++) {
if (p.vertex[i].x < k.x) k.x = p.vertex[i].x;
if (p.vertex[i].y > k.y) k.y = p.vertex[i].y;
}
k.y++;
do{
k.x--;
flag = 1;
for (i=0; i<p.num; i++) if (isPointOnSegment(p.vertex[i], a, k))
{
flag = 0; break;
}
}while (!flag);
times = 0;
for(i = 0; i < p.num; i++)
if(isSegmentIntersect(a, k, p.vertex[i], p.vertex[(i+1)%p.num]))
times++;
return times%2 == 1;
}
//用面积判断点是否在简单多边形内(边界也算在内)
bool isPointInPolygon(Point p, Polygon poly){
double area = 0;
int n = poly.num;
poly.vertex[n].x = poly.vertex[0].x, poly.vertex[n].y = poly.vertex[0].y;
Triangle t;
t.vertex[2].x = p.x; t.vertex[2].y = p.y;
for(int i = 1; i <= n; i++) {
t.vertex[0].x = poly.vertex[i-1].x; t.vertex[0].y = poly.vertex[i-1].y;
t.vertex[1].x = poly.vertex[i].x; t.vertex[1].y = poly.vertex[i].y;
area += TriangleArea(t);
}
return isZero(area-TriangleArea(t));
}
//用线段拐向判断点是否在简单多边形内(不包括边界,亦可以修改成包括边界)
bool isPointIntriangle2(Point p, Polygon poly){
int n = poly.num;
double k = CrossProduct(p, poly.vertex[n-1], p, poly.vertex[0]);
double kk;
for(int i = 1; i < n; i++)
{
kk = CrossProduct(p, poly.vertex[i-1], p, poly.vertex[i]);
if(kk*k<=0)
return false;
}
return true;
}
//求多边形p的重心(BaryCenter), 结果存入c(坐标),s(面积)中, 2004.12.1
void Calc(int n, Point vertex[], Point center[], double area[]){
center[n].x = (vertex[0].x+vertex[1].x+vertex[2].x)/3.0;
center[n].y = (vertex[0].y+vertex[1].y+vertex[2].y)/3.0;
area[n] = (vertex[2].y-vertex[1].y)*(vertex[1].x-vertex[0].x)
- (vertex[2].x-vertex[1].x)*(vertex[1].y-vertex[0].y);
area[n] /= 2;
}
bool ADD(int a, int b, Point center[], double area[]){
if(fabs(area[a]+area[b]) < EPS) return false;
center[a].x += (center[b].x-center[a].x)*area[b]/(area[a]+area[b]);
center[a].y += (center[b].y-center[a].y)*area[b]/(area[a]+area[b]);
area[a] += area[b];
return true;
}
void BaryCenter(Point *c, double *s, const Polygon *p){
int i;
Point vertex[3], center[3];
double area[3];
for (i=0; i<3; i++) vertex[i] = p->vertex[i];
Calc(0, vertex, center, area);
i = 3;
while(i < p->num){
vertex[1] = vertex[2];
vertex[2] = p->vertex[i];
i++;
Calc(1, vertex, center, area);
while(!ADD(0, 1, center, area)){
vertex[1] = vertex[2];
vertex[2] = p->vertex[i];
i++;
Calc(2, vertex, center, area);
if (!ADD(0, 2, center, area))
ADD(1, 2, center, area);
}
}
*c = center[0];
*s = area[0];
}
//******************************************凸包*******************************************//
//Graham扫描法在点集s中寻找凸包,结果存入p中
//O(N^2) SelectSort
void GrahamScan(Polygon *p, PSet s){
int i, j, k = 0;
double mul;
Point t;
for(i = 1; i < s.num; i++)
if (s.vertex[i].y<s.vertex[k].y ||
(Equ(s.vertex[i].y,s.vertex[k].y) &&
s.vertex[i].x<s.vertex[k].x))
k = i;
t = s.vertex[0];
s.vertex[0] = s.vertex[k];
s.vertex[k] = t;
for(i = 1; i < s.num-1; i++){
k = i;
for(j = i+1; j < s.num; j++){
mul = CrossProduct(s.vertex[0], s.vertex[j], s.vertex[0], s.vertex[k]);
if (mul>0 || (Equ(mul,0) && Distance(s.vertex[0],s.vertex[j])
< Distance(s.vertex[0],s.vertex[k])))
k = j;
}
t = s.vertex[i];
s.vertex[i] = s.vertex[k];
s.vertex[k] = t;
}
*p = s;
for(i = 0; (i+1)%p->num != 0;){
mul = CrossProduct(p->vertex[i], p->vertex[(i+1)%p->num],
p->vertex[(i+1)%p->num], p->vertex[(i+2)%p->num]);
if (mul >= 0) i = (i+1)%p->num;
else{
for(j = (i+1)%p->num+1; j < p->num; j++)
p->vertex[j-1] = p->vertex[j];
p->num--;
i = (i-1+p->num)%p->num;
}
}
}
//Graham扫描法在点集s中寻找凸包结果存入p中
//O(NlogN) QuickSort and Stack, 2004.11.30
Point workarr[MAX];
int stk[MAX], top;
double Multi(Point p1, Point p2, Point p0){
return (p1.x-p0.x)*(p2.y-p0.y) - (p2.x-p0.x)*(p1.y-p0.y);
}
int Comp(const void *e1, const void *e2){
Point p1 = *((Point *)e1), p2 = *((Point *)e2);
Point p0 = workarr[0];
double mul = Multi(p1, p2, p0);
if (mul > 0 || (Equ(mul, 0) && Distance(p0, p1)< Distance(p0, p2)))
return -1;
else return 1;
}
void GrahamScan2(Polygon *p, const PSet *s){
int i, k;
Point t;
if(s->num < 3){
p->num = s->num;
for(i = 0; i < s->num; i++)
p->vertex[i] = s->vertex[i];
return;
}
for(i = 0; i < s->num; i++) workarr[i] = s->vertex[i];
k = 0;
for(i = 1; i < s->num; i++)
if (workarr[i].y < workarr[k].y ||
(Equ(workarr[i].y, workarr[k].y) &&
workarr[i].x<workarr[k].x))
k = i;
t = workarr[0];
workarr[0] = workarr[k];
workarr[k] = t;
qsort(workarr+1, s->num-1, sizeof(Point), Comp);
for(i = 1; i <= 3; i++) stk[i] = i-1;
top = 3;
for(i = 3; i < s->num; i++){
while (Multi(workarr[i], workarr[stk[top]], workarr[stk[top-1]]) > 0) top--;
top++;
stk[top] = i;
}
p->num = top;
for(i = 1; i <= top; i++) p->vertex[i-1] = workarr[stk[i]];
}