74 lines
2.6 KiB
C++
74 lines
2.6 KiB
C++
|
#pragma once
|
||
|
|
#include "vde/foundation/exact_predicates.h"
|
||
|
|
|
||
|
|
#ifdef VDE_USE_GMP
|
||
|
|
#include <gmpxx.h>
|
||
|
|
|
||
|
|
namespace vde::foundation::exact {
|
||
|
|
|
||
|
|
/// GMP-based exact predicates — zero rounding error, arbitrary precision
|
||
|
|
struct GmpExact {
|
||
|
|
using Rational = mpq_class;
|
||
|
|
|
||
|
|
/// orient_2d with exact rational arithmetic
|
||
|
|
static Orientation orient_2d(const Point2D& a, const Point2D& b, const Point2D& c) {
|
||
|
|
Rational ax(a.x()), ay(a.y());
|
||
|
|
Rational bx(b.x()), by(b.y());
|
||
|
|
Rational cx(c.x()), cy(c.y());
|
||
|
|
Rational det = (bx - ax) * (cy - ay) - (by - ay) * (cx - ax);
|
||
|
|
int s = mpq_sgn(det.get_mpq_t());
|
||
|
|
return s > 0 ? Orientation::CounterClockwise :
|
||
|
|
s < 0 ? Orientation::Clockwise : Orientation::Collinear;
|
||
|
|
}
|
||
|
|
|
||
|
|
/// orient_3d with exact rational arithmetic
|
||
|
|
static Orientation orient_3d(const Point3D& a, const Point3D& b,
|
||
|
|
const Point3D& c, const Point3D& d) {
|
||
|
|
Rational ax(a.x()), ay(a.y()), az(a.z());
|
||
|
|
Rational bx(b.x()), by(b.y()), bz(b.z());
|
||
|
|
Rational cx(c.x()), cy(c.y()), cz(c.z());
|
||
|
|
Rational dx(d.x()), dy(d.y()), dz(d.z());
|
||
|
|
|
||
|
|
Rational adx = ax - dx, ady = ay - dy, adz = az - dz;
|
||
|
|
Rational bdx = bx - dx, bdy = by - dy, bdz = bz - dz;
|
||
|
|
Rational cdx = cx - dx, cdy = cy - dy, cdz = cz - dz;
|
||
|
|
|
||
|
|
Rational det = adx * (bdy * cdz - bdz * cdy)
|
||
|
|
+ bdx * (cdy * adz - cdz * ady)
|
||
|
|
+ cdx * (ady * bdz - adz * bdy);
|
||
|
|
|
||
|
|
int s = mpq_sgn(det.get_mpq_t());
|
||
|
|
return s > 0 ? Orientation::CounterClockwise :
|
||
|
|
s < 0 ? Orientation::Clockwise : Orientation::Collinear;
|
||
|
|
}
|
||
|
|
|
||
|
|
/// in_circle with exact rational arithmetic
|
||
|
|
static CircleTest in_circle(const Point2D& a, const Point2D& b,
|
||
|
|
const Point2D& c, const Point2D& d) {
|
||
|
|
Rational ax(a.x()), ay(a.y());
|
||
|
|
Rational bx(b.x()), by(b.y());
|
||
|
|
Rational cx(c.x()), cy(c.y());
|
||
|
|
Rational dx(d.x()), dy(d.y());
|
||
|
|
|
||
|
|
Rational adx = ax - dx, ady = ay - dy;
|
||
|
|
Rational bdx = bx - dx, bdy = by - dy;
|
||
|
|
Rational cdx = cx - dx, cdy = cy - dy;
|
||
|
|
|
||
|
|
Rational ab = adx * bdy - ady * bdx;
|
||
|
|
Rational bc = bdx * cdy - bdy * cdx;
|
||
|
|
Rational ca = cdx * ady - cdy * adx;
|
||
|
|
Rational al = adx * adx + ady * ady;
|
||
|
|
Rational bl = bdx * bdx + bdy * bdy;
|
||
|
|
Rational cl = cdx * cdx + cdy * cdy;
|
||
|
|
|
||
|
|
Rational det = al * bc + bl * ca + cl * ab;
|
||
|
|
int s = mpq_sgn(det.get_mpq_t());
|
||
|
|
return s > 0 ? CircleTest::Inside :
|
||
|
|
s < 0 ? CircleTest::Outside : CircleTest::On;
|
||
|
|
}
|
||
|
|
};
|
||
|
|
|
||
|
|
} // namespace vde::foundation::exact
|
||
|
|
|
||
|
|
#endif // VDE_USE_GMP
|