#pragma once #include "vde/foundation/exact_predicates.h" #ifdef VDE_USE_GMP #include namespace vde::foundation::exact { /** * @brief GMP 精确算术谓词后端 * * 基于 GNU MP(多精度有理数)的精确几何谓词。 * 使用高精度有理数运算,零舍入误差,适用于对鲁棒性要求极高的场景。 * * @note 仅在定义了 VDE_USE_GMP 宏且链接了 GMP 库时可用 * @warning GMP 版本比自适应双精度版本慢约 10-50 倍,仅在双精度不足以确定结果时使用 * * @ingroup foundation */ struct GmpExact { /// GMP 有理数类型 using Rational = mpq_class; /** * @brief 精确有理数二维定向测试 * * 使用 GMP 有理数精确计算三点定向,零误差。 * * @param a 第一个点 * @param b 第二个点 * @param c 第三个点 * @return 定向结果(保证正确,无舍入误差) * * @code{.cpp} * // 对近共线点仍能给出确定结果 * Point2D a(0, 0), b(1, 1e-15), c(2, 2e-15); * auto ori = GmpExact::orient_2d(a, b, c); * @endcode * * @see GmpExact::orient_3d, GmpExact::in_circle */ 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; } /** * @brief 精确有理数三维定向测试 * * 判断点 d 相对于平面 abc 的位置,使用有理数精确运算。 * * @param a,b,c 定义平面的三点 * @param d 待测试点 * @return 定向结果(零误差) * * @see GmpExact::orient_2d */ 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; } /** * @brief 精确有理数圆内测试 * * 判断点 d 是否在三角形 abc 的外接圆内,零误差。 * * @param a,b,c 三角形三顶点 * @param d 待测试点 * @return 圆测试结果(零误差) * * @see GmpExact::orient_2d */ 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