A polygon surface area calculation problem.
The formula for the regular n-gon’s area is:
And since the interior angle of a polygon is
C++ implementation:
#include <iostream>
#include <cmath>
using namespace std;
typedef long double gtype;
const gtype pi = M_PI;
#define polA(n) (((n) - 2) * pi / (n))
int main() {
gtype r, n;
cout.precision(3);
while(cin >> r >> n)
cout << fixed << (r * r * sin(polA(n)) * n * 0.5) << endl;
return 0;
}
A circumcircle problem.
Given a triangle, find the circumference of its circumcircle. There’s an easy way and a hard way.
The hard way is to find the circle’s centre then calculate the radius. I used this solution for a previous problem. Code:
#include <iostream>
#include <complex>
#include <limits>
using namespace std;
typedef long double gtype;
#define equ(a, b, e) ( abs((a) - (b)) <= (e) )
const gtype len_before_point_digits = 1E6;
const gtype epsLen = numeric_limits<gtype>::epsilon() * len_before_point_digits;
const gtype pi = M_PI;
typedef complex<gtype> point;
#define x real()
#define y imag()
#define crs(a, b) ( (conj(a) * (b)).y )
#define perp(p) ( point((p).y, -(p).x) )
#define intrN(a, b, p, q) ( crs((p)-(a),(b)-(a)) * (q) - crs((q)-(a),(b)-(a)) * (p) )
#define intrD(a, b, p, q) ( crs((p)-(a),(b)-(a)) - crs((q)-(a),(b)-(a)) )
#define circC(r) (2 * pi * abs(r))
#define plus(x) ((x) > 0 ? " + " : " - ")
int main() {
point p[3], v[2][2], c;
gtype r;
cout.precision(2);
while (cin >> p[0].x >> p[0].y >> p[1].x >> p[1].y >> p[2].x >> p[2].y) {
// perp 1
v[0][0] = 0.5l * (p[0] + p[1]);
v[0][1] = perp(p[1] - p[0]) + v[0][0];
// perp 2
v[1][0] = 0.5l * (p[1] + p[2]);
v[1][1] = perp(p[2] - p[1]) + v[1][0];
// intersect
gtype d = intrD(v[0][0], v[0][1], v[1][0], v[1][1]);
if (!equ(d, 0, epsLen)) {
c = intrN(v[0][0], v[0][1], v[1][0], v[1][1]) / d;
r = abs(c - p[0]);
} else
r = 0;
// out
cout << fixed << circC(r) << endl;
}
return 0;
}
The easy way is to find the circumcircle’s radius through:
And since:
This method is less prone to error. Implementation:
#include <iostream>
#include <complex>
#include <limits>
using namespace std;
typedef long double gtype;
const gtype pi = M_PI;
typedef complex<gtype> point;
#define x real()
#define y imag()
#define crs(a, b) ( (conj(a) * (b)).y )
#define circC(r) (2 * pi * abs(r))
#define area(a, b, c) ( (gtype) 0.5 * crs((b) - (a), (c) - (a)) )
#define plus(x) ((x) > 0 ? " + " : " - ")
int main() {
point p[3];
gtype r;
cout.precision(2);
while (cin >> p[0].x >> p[0].y >> p[1].x >> p[1].y >> p[2].x >> p[2].y) {
r = abs(p[0] - p[1]) * abs(p[1] - p[2]) * abs(p[2] - p[0]) / area(p[0], p[1], p[2]) * 0.25;
cout << fixed << circC(r) << endl;
}
return 0;
}
A geometry problem.
Given three points
Given triangle
We need to transform this process into actual computation with vectors. A mid-point is easily determined. A perpendicular to a vector can be found by swapping the coordinates and negating only one of them. Intersection is also a simple task.
C++ implementation using some geometry tools:
#include <algorithm>
#include <iostream>
#include <complex>
#include <cstdlib>
#include <limits>
using namespace std;
typedef double gtype;
// const
const gtype len_before_point_digits = 1E6;
const gtype epsLen = numeric_limits<gtype>::epsilon() * len_before_point_digits;
// point
typedef complex<gtype> point;
#define x real()
#define y imag()
// vector
#define crs(a, b) ( (conj(a) * (b)).y )
#define perp(p) ( point((p).y, -(p).x) )
// line
#define intrN(a, b, p, q) ( crs((p)-(a),(b)-(a)) * (q) - crs((q)-(a),(b)-(a)) * (p) )
#define intrD(a, b, p, q) ( crs((p)-(a),(b)-(a)) - crs((q)-(a),(b)-(a)) )
// I/O
#define plus(x) ((x) > 0 ? " + " : " - ")
int main() {
point p[3], v[2][2], c;
gtype r;
cout.precision(3);
while (cin >> p[0].x >> p[0].y >> p[1].x >> p[1].y >> p[2].x >> p[2].y) {
// perp 1
v[0][0] = 0.5 * (p[0] + p[1]);
v[0][1] = perp(p[1] - p[0]) + v[0][0];
// perp 2
v[1][0] = 0.5 * (p[1] + p[2]);
v[1][1] = perp(p[2] - p[1]) + v[1][0];
// intersect
c = intrN(v[0][0], v[0][1], v[1][0], v[1][1]) / intrD(v[0][0], v[0][1], v[1][0], v[1][1]);
r = abs(c - p[0]);
// out
gtype h = c.x, k = c.y;
cout << fixed << "(x" << plus(-h) << abs(h) << ")^2 + (y" <<
plus(-k) << abs(k) << ")^2 = " << r << "^2" << endl;
gtype a = 2 * c.x, b = 2 * c.y, d = c.x * c.x + c.y * c.y - r * r;
cout << fixed << "x^2 + y^2" << plus(-a) << abs(a) << "x" <<
plus(-b) << abs(b) << "y" << plus(d) << abs(d) << " = 0" << endl;
cout << endl;
}
return 0;
}
A well-tested C++ library for computational geometry problems. The code has bad style and conventions, but it’ll suffice for hackathons.
The numerical type involved in computation varies depending on the problem. So, for flexibility we’ll name this type gtype that we’ll use throughout the code.
The macros equ and inRange should be self-explanatory, however, note that the comparison they make use some tolerance bound e.
After that, follow some constants that are of type gtype. len_before_point_digits and ang_before_point_digits are also problem dependent and they’re VERY critical. They define how many decimal digits before the decimal point we’re going to use. For example, if we know from our problem that distances will never exceed 100,000, len_before_point_digits should be 1E6. The same applies for ang_before_point_digits.
epsLen, epsAng and epsArea are the tolerance bounds we’re likely to use for distances, radian angles and areas. Yep! It’s that simple.
typedef long double gtype; #define equ(a, b, e) ( abs((a) - (b)) <= (e) ) #define inRange(a, b, p, e) ( (p) >= min(a, b) - (e) && (p) <= max(a, b) + (e)) // CONST const gtype len_digits_before_point = 1E6; const gtype epsLen = numeric_limits<gtype>::epsilon() * len_digits_before_point; const gtype ang_digits_before_point = 1E2; const gtype epsAng = numeric_limits<gtype>::epsilon() * ang_digits_before_point; const gtype epsArea = epsLen * epsLen; const gtype pi = M_PI;
Let’s get down to business. Let’s define a point. What’s a point? It’s an ordered pair of some type. So, instead of writing it from scratch, we’ll use C++‘s std::complex which is also an ordered pair that can be given a numerical type. std::complex spares us a lot of effort in overriding operators and has some very useful methods. We define the type std::complex <gtype> as our point.
To get a point‘s x and y, we get std::complex‘s real and imaginary parts respectively.
Another good feature of std::complex is that it allows both Cartesian and polar constructors. But, the polar constructor must be given the correct type for the magnitude component, that’s why we defined the macro.
The last two macros do simple angle conversion.
typedef complex<gtype> point;
#define x real()
#define y imag()
#define polar(r, t) polar((gtype) (r), (gtype) (t))
bool cmp_point(point const& p1, point const& p2) {
return p1.x == p2.x ? (p1.y < p2.y) : (p1.x < p2.x);
}
#define deg(a) ( (a) * 180 /pi )
#define rad(a) ( (a) * pi /180 )
Vectors are the building blocks of any geometrical structure. Rotations, products and transformations are very important in any problem. Fortunately, we don’t have to implement all of that. norm, abs, conj, arg are provided by std::complex. We’ll use some of these functions and the rest we’ll code on our own. I’ve chosen the fastest and most precise functions I could lay my hands on:
// VECTOR //#define norm(a) dot(a,a) // fast for [long double] //#define abs(a) sqrt(norm(a)) // fast for [int] #define arg(a) ( atan2((a).y, (a).x) ) #define unit(v) ( (v) / abs(v) ) #define dot(a, b) ( (conj(a) * (b)).x ) #define crs(a, b) ( (conj(a) * (b)).y ) #define rot(v, t) ( (v) * polar(1, t) ) #define rot2(v, t, o) ( rot((v) - (o), (t)) + (o) ) #define refelct(v, m) ( conj((v) / (m)) * (m) ) #define perp(v) ( point((v).y, -(v).x) ) #define vecVecProj(u, v) ( dot(u, v) * (v) / norm(v) )
Lines-point / line-line intersections are very frequently used. So, here are some useful macros:
// LINE #define pntLinProj(a, b, p) ( (a) + vecVecProj((p) - (a), (b) - (a)) ) #define pntLinDist(a, b, p) ( abs(crs((b)-(a), (p)-(a)) / abs((b)-(a))) ) #define linIntrN(a, b, p, q) ( crs((p)-(a),(b)-(a)) * (q) - crs((q)-(a),(b)-(a)) * (p) ) #define linIntrD(a, b, p, q) ( crs((p)-(a),(b)-(a)) - crs((q)-(a),(b)-(a)) ) #define linOnLin(a, b, p) ( abs(crs((b) - (a), (p) - (a))) <= epsLen ) #define linOnRay(a, b, p) ( linOnLin(a, b, p) && dot((b) - (a), (p) - (a)) >= -epsLen ) #define linOnSeg(a, b, p) ( linOnRay(b, a, p) && dot((b) - (a), (p) - (a)) >= -epsLen ) #define linInRect(a, b, p, e) ( inRange((a).x, (b).x, (p).x, e) && inRange((a).y, (b).y, (p).y, e) )
Perimeter, surface, and cosine law:
// TRIANGLE #define triPerim(a, b, c) ( (a) + (b) + (c) ) #define triSurface(a, b, c) ( (gtype) 0.5 * crs((b) - (a), (c) - (a)) ) #define triSurfaceH(s, a, b, c) ( sqrt((s) * ((s) - (a)) * ((s) - (b)) * ((s) - (c))) ) #define cosLaw(a, b, c) ( ((a) * (a) + (b) * (b) - (c) * (c)) / (2 * (a) * (b)) )
Arc length, chord length, sector surface area, and circle circle intersection:
// CIRCLE
#define cirArcLen(r, a) ((a) * abs(r))
#define cirSectorSurface(r, a) ((a) * norm(r) / 2)
#define cirChordLen(r, a) sqrt(2 * norm(r) * (1 - cos(a)))
int cirCircInter(point & c0, point & c1, gtype r0, gtype r1, vector<point> & out) {
point cT, v;
gtype rT, phi;
// init
if (r1 > r0)
rT = r0, r0 = r1, r1 = rT, cT = c0, c0 = point(c1), c1 = point(cT);
out.clear(), v = c1 - c0, phi = acos(cosLaw(r0, abs(v), r1));
// same center
if (equ(abs(v), 0, epsLen))
equ(r0, r1, epsLen) ? (equ(r0, 0, epsLen) ? out.push_back(c0) : out.resize(3)) : out.clear();
// none
else if (abs(abs(v) - r0) > r1 + epsLen)
out.clear();
// one point
else if (equ(abs(abs(v) - r0), r1, epsLen))
out.push_back(unit(v) * r0 + c0);
// two points
else
out.push_back(rot(unit(v) * r0, -phi) + c0), out.push_back(rot(unit(v) * r0, phi) + c0);
return out.size();
}
Area, convexity check, and rotating calipers (min-width). Also includes convex hull creation:
// POLYGON (must be closed, start & end vertex overlap)
#define polAngle(n) (((n) - 2) * pi / (n))
long double polArea(vector<point> & p) {
long double area = 0;
for (size_t i = 0; i < p.size() - 1; ++i)
area += crs(p[i], p[i + 1]);
return abs(area * 0.5l);
}
bool polIsConvex(vector<point> & p) {
int a[2] = { 0, 0 }, n = p.size() - 1;
for (size_t i = 1; i < p.size(); ++i) {
long double cross = crs(p[i % n] - p[i - 1], p[(i + 1) % n] - p[i - 1]);
if (cross != 0)
++a[cross > 0];
}
return !(a[0] && a[1]);
}
// CONVEX HULL
// O(n.h) - jarvis's march
// TODO
// O(n.log(n)) - monotone chain
vector<point> mcH;
void monotoneChain(vector<point> &ps) {
vector<point> p(ps.begin(), ps.end() - 1);
int n = p.size(), k = 0;
mcH = vector<point>(2 * n);
sort(p.begin(), p.end(), cmp_point);
for (int i = 0; i < n; i++) {
while (k >= 2 && crs(mcH[k - 1] - mcH[k - 2], p[i] - mcH[k - 2]) <= 0)
k--;
mcH[k++] = p[i];
}
for (int i = n - 2, t = k + 1; i >= 0; i--) {
while (k >= t && crs(mcH[k - 1] - mcH[k - 2], p[i] - mcH[k - 2]) <= 0)
k--;
mcH[k++] = p[i];
}
mcH.resize(k);
}
// O(n) - rotating calipers (works on a ccw closed convex hull)
gtype rotatingCalipers(vector<point> &ps) {
int aI = 0, bI = 0;
for (size_t i = 1; i < ps.size(); ++i)
aI = (ps[i].y < ps[aI].y ? i : aI), bI = (ps[i].y > ps[bI].y ? i : bI);
gtype minWidth = ps[bI].y - ps[aI].y, aAng, bAng;
point aV = point(1, 0), bV = point(-1, 0);
for (gtype ang = 0; ang < pi; ang += min(aAng, bAng)) {
aAng = acos(dot(ps[aI + 1] - ps[aI], aV) / abs(aV) / abs(ps[aI + 1] - ps[aI]));
bAng = acos(dot(ps[bI + 1] - ps[bI], bV) / abs(bV) / abs(ps[bI + 1] - ps[bI]));
aV = rot(aV, min(aAng, bAng)), bV = rot(bV, min(aAng, bAng));
if (aAng < bAng)
minWidth = min(minWidth, pntLinDist(ps[aI], ps[aI] + aV, ps[bI])), aI = (aI + 1) % (ps.size() - 1);
else
minWidth = min(minWidth, pntLinDist(ps[bI], ps[bI] + bV, ps[aI])), bI = (bI + 1) % (ps.size() - 1);
}
return minWidth;
}
This is some useful C++ code that does some tests on three points:
#include <algorithm>
#include <iostream>
#include <complex>
#include <cstdlib>
#include <limits>
using namespace std;
typedef long double gtype;
#define before_point_digits ( 1E6 )
#define epsLen ( (gtype) numeric_limits<gtype>::epsilon() * before_point_digits )
#define epsAng ( (gtype) numeric_limits<gtype>::epsilon() * 1E2)
typedef complex<gtype> point;
#define x real()
#define y imag()
#define equ(a, b, e) ( abs(a - b) <= e )
#define dot(a, b) ( (conj(a) * (b)).x )
#define crs(a, b) ( (conj(a) * (b)).y )
int main() {
int nCase;
gtype xP, yP, l[3];
double a[3];
point p[3];
cin >> nCase;
for (int caseI = 0; caseI < nCase; ++caseI) {
// input
for (int i = 0; i < 3; ++i) {
cin >> xP >> yP;
p[i].x = xP;
p[i].y = yP;
}
// calc
l[0] = abs(p[0] - p[1]);
l[1] = abs(p[1] - p[2]);
l[2] = abs(p[2] - p[0]);
a[0] = dot(p[1] - p[0], p[2] - p[0]);
a[1] = dot(p[0] - p[1], p[2] - p[1]);
a[2] = dot(p[0] - p[2], p[1] - p[2]);
// test
if (!equ(0.5 * crs(p[1] - p[0], p[2] - p[0]), 0, epsLen)) {
// sides
if (equ(l[0], l[1], epsLen) && equ(l[1], l[2], epsLen))
cout << "equilateral ";
else if (equ(l[0], l[1], epsLen) || equ(l[1], l[2], epsLen) || equ(l[2], l[0], epsLen))
cout << "isosceles ";
else
cout << "scalene ";
// angles
if (equ(a[0], 0, epsAng) || equ(a[1], 0, epsAng) || equ(a[2], 0, epsAng))
cout << "right ";
else if (a[0] < epsAng || a[1] < epsAng || a[2] < epsAng)
cout << "obtuse ";
else
cout << "acute ";
} else
cout << "not a ";
cout << "triangle" << endl;
}
return 0;
}
Sample input:
7 0 0 0 4 1 2 1 1 1 4 3 2 0 0 0 0 0 0 3 3 3 4 5 3 4 4 4 5 5 6 5 5 5 6 6 5 6 6 6 7 6 8
Sample output:
isosceles obtuse triangle scalene acute triangle not a triangle scalene right triangle scalene obtuse triangle isosceles right triangle not a triangle
gsourcecode.wordpress.com 