茂之钳
238 lines
6.8 KiB
C++
238 lines
6.8 KiB
C++
#pragma once
|
||
#include <cmath>
|
||
#include <limits>
|
||
#include <algorithm>
|
||
|
||
namespace vde::foundation {
|
||
|
||
/**
|
||
* @brief 双精度区间算术类
|
||
*
|
||
* 存储区间 [lo, hi] 的上下边界，支持区间算术运算和确定性比较。
|
||
* 主要用于验证浮点计算的可靠性：若区间结果不跨越零，则可确定符号；
|
||
* 若跨越零，则说明双精度不足以确定结果，需要更高精度或精确算术。
|
||
*
|
||
* 区间运算遵循保守原则（结果区间包含所有可能的真实值），
|
||
* 即 out = [min(所有组合), max(所有组合)]。
|
||
*
|
||
* @ingroup foundation
|
||
*/
|
||
class Interval {
|
||
public:
|
||
double lo, hi;
|
||
|
||
/// 默认构造：空区间 [0, 0]
|
||
Interval() : lo(0), hi(0) {}
|
||
|
||
/**
|
||
* @brief 从单值构造精确区间（退化为点）
|
||
* @param v 区间值，[v, v]
|
||
*/
|
||
Interval(double v) : lo(v), hi(v) {}
|
||
|
||
/**
|
||
* @brief 从上下界构造区间
|
||
*
|
||
* 自动交换以确保 lo ≤ hi。
|
||
*
|
||
* @param l 下界
|
||
* @param h 上界
|
||
*/
|
||
Interval(double l, double h) : lo(std::min(l,h)), hi(std::max(l,h)) {}
|
||
|
||
/**
|
||
* @brief 区间中点
|
||
* @return (lo + hi) / 2
|
||
*/
|
||
[[nodiscard]] double mid() const { return (lo + hi) * 0.5; }
|
||
|
||
/**
|
||
* @brief 区间宽度
|
||
* @return hi - lo
|
||
*/
|
||
[[nodiscard]] double width() const { return hi - lo; }
|
||
|
||
/**
|
||
* @brief 判断区间是否包含给定值
|
||
* @param v 测试值
|
||
* @return true 若 lo ≤ v ≤ hi
|
||
*/
|
||
[[nodiscard]] bool contains(double v) const { return lo <= v && v <= hi; }
|
||
|
||
/**
|
||
* @brief 判断区间是否足够窄（收敛）
|
||
* @param eps 宽度阈值（默认 1e-12）
|
||
* @return true 若区间宽度 < eps
|
||
*/
|
||
[[nodiscard]] bool is_exact(double eps = 1e-12) const { return width() < eps; }
|
||
|
||
/**
|
||
* @brief 判断两个区间是否有重叠
|
||
* @param o 另一个区间
|
||
* @return true 若存在公共点
|
||
*/
|
||
[[nodiscard]] bool overlaps(const Interval& o) const {
|
||
return !(hi < o.lo || o.hi < lo);
|
||
}
|
||
|
||
/**
|
||
* @brief 确定性区间比较
|
||
*
|
||
* 仅当两个区间无重叠时才能确定大小关系。
|
||
*
|
||
* @param o 另一个区间
|
||
* @return -1 若确定 this < o
|
||
* @return 1 若确定 this > o
|
||
* @return 0 若区间重叠，无法确定
|
||
*
|
||
* @code{.cpp}
|
||
* Interval a(1.0, 2.0), b(3.0, 4.0);
|
||
* int c = a.cmp(b); // c == -1 (确定 a < b)
|
||
*
|
||
* Interval c(1.0, 3.0), d(2.0, 4.0);
|
||
* int r = c.cmp(d); // r == 0 (重叠，无法确定)
|
||
* @endcode
|
||
*/
|
||
[[nodiscard]] int cmp(const Interval& o) const {
|
||
if (hi < o.lo) return -1; // 确定 this < o
|
||
if (lo > o.hi) return 1; // 确定 this > o
|
||
return 0; // 不确定（区间重叠）
|
||
}
|
||
|
||
// ── 算术运算 ──
|
||
|
||
/**
|
||
* @brief 区间加法
|
||
* @return [this.lo + o.lo, this.hi + o.hi]
|
||
*/
|
||
Interval operator+(const Interval& o) const {
|
||
return {lo + o.lo, hi + o.hi};
|
||
}
|
||
|
||
/**
|
||
* @brief 区间减法
|
||
* @return [this.lo - o.hi, this.hi - o.lo]（保守估计）
|
||
*/
|
||
Interval operator-(const Interval& o) const {
|
||
return {lo - o.hi, hi - o.lo};
|
||
}
|
||
|
||
/**
|
||
* @brief 区间乘法
|
||
*
|
||
* 取四项乘积的最小值和最大值作为结果区间。
|
||
*
|
||
* @return [min(a*c,a*d,b*c,b*d), max(a*c,a*d,b*c,b*d)]
|
||
*/
|
||
Interval operator*(const Interval& o) const {
|
||
double a = lo*o.lo, b = lo*o.hi, c = hi*o.lo, d = hi*o.hi;
|
||
return {std::min({a,b,c,d}), std::max({a,b,c,d})};
|
||
}
|
||
|
||
/**
|
||
* @brief 区间除法
|
||
*
|
||
* 若除区间跨越零，返回极大区间表示不确定性。
|
||
*
|
||
* @return 结果区间或 [-∞, +∞]（除数为零区间）
|
||
*/
|
||
Interval operator/(const Interval& o) const {
|
||
if (o.lo <= 0 && o.hi >= 0) return {-std::numeric_limits<double>::max(),
|
||
std::numeric_limits<double>::max()};
|
||
double a = lo/o.lo, b = lo/o.hi, c = hi/o.lo, d = hi/o.hi;
|
||
return {std::min({a,b,c,d}), std::max({a,b,c,d})};
|
||
}
|
||
|
||
/// 取负
|
||
Interval operator-() const { return {-hi, -lo}; }
|
||
|
||
// ── 数学函数（保守区间）──
|
||
|
||
/**
|
||
* @brief 区间平方根
|
||
*
|
||
* 下界对 0 取 max 以处理负值（保守处理）。
|
||
*
|
||
* @return [sqrt(max(0, lo)), sqrt(max(0, hi))]
|
||
*/
|
||
[[nodiscard]] Interval sqrt() const {
|
||
return {std::sqrt(std::max(0.0, lo)), std::sqrt(std::max(0.0, hi))};
|
||
}
|
||
|
||
/**
|
||
* @brief 从双精度值创建区间（含舍入误差带）
|
||
*
|
||
* @param v 中心值
|
||
* @param error 误差半径（默认 0）
|
||
* @return [v - error, v + error]
|
||
*/
|
||
static Interval from_double(double v, double error = 0) {
|
||
return {v - error, v + error};
|
||
}
|
||
|
||
/**
|
||
* @brief 从坐标值创建区间（含 ULP 误差带）
|
||
*
|
||
* 误差估计为 |v| * 2.22e-16 + 1e-300，即双精度最小单位。
|
||
*
|
||
* @param v 坐标值
|
||
* @return 包含浮点舍入误差的区间
|
||
*/
|
||
static Interval from_coord(double v) {
|
||
double eps = std::abs(v) * 2.22e-16 + 1e-300;
|
||
return {v - eps, v + eps};
|
||
}
|
||
};
|
||
|
||
/**
|
||
* @brief 二维定向行列式的区间验证
|
||
*
|
||
* 用区间算术计算 orient_2d 行列式，结果区间用于判断双精度是否可靠。
|
||
*
|
||
* @param ax,ay 点 a 的坐标
|
||
* @param bx,by 点 b 的坐标
|
||
* @param cx,cy 点 c 的坐标
|
||
* @return 行列式的区间值
|
||
*
|
||
* @see orient_2d_verified
|
||
* @ingroup foundation
|
||
*/
|
||
inline Interval orient_2d_interval(double ax, double ay, double bx, double by, double cx, double cy) {
|
||
Interval acx = Interval::from_coord(ax - cx);
|
||
Interval acy = Interval::from_coord(ay - cy);
|
||
Interval bcx = Interval::from_coord(bx - cx);
|
||
Interval bcy = Interval::from_coord(by - cy);
|
||
return acx * bcy - acy * bcx;
|
||
}
|
||
|
||
/**
|
||
* @brief 使用区间算术验证 orient_2d 结果
|
||
*
|
||
* 若双精度计算结果接近零，调用此函数进行区间验证。
|
||
* 若区间不跨越零，则结果可靠；否则需要升级到精确算术。
|
||
*
|
||
* @param ax,ay 点 a 的坐标
|
||
* @param bx,by 点 b 的坐标
|
||
* @param cx,cy 点 c 的坐标
|
||
* @return -1 确认 Clockwise
|
||
* @return 0 无法确定（需升级到 GMP / 自适应精度）
|
||
* @return 1 确认 CounterClockwise
|
||
*
|
||
* @code{.cpp}
|
||
* int result = orient_2d_verified(0.0, 0.0, 1.0, 1e-16, 2.0, 2e-16);
|
||
* if (result == 0) {
|
||
* // 双精度不可靠，需要精确算术
|
||
* auto ori = GmpExact::orient_2d(a, b, c);
|
||
* }
|
||
* @endcode
|
||
*
|
||
* @see orient_2d_interval, Interval::cmp
|
||
* @ingroup foundation
|
||
*/
|
||
inline int orient_2d_verified(double ax, double ay, double bx, double by, double cx, double cy) {
|
||
Interval det = orient_2d_interval(ax, ay, bx, by, cx, cy);
|
||
return det.cmp(Interval(0.0));
|
||
}
|
||
|
||
} // namespace vde::foundation
|