19.8 混合精度计算

在本节中你将学到…

混合精度计算的设计动机
Lazy_exact_nt的实现原理
过滤机制的工作方式
精度控制策略
在CGAL中实现鲁棒的几何计算

19.8.1 概念解释

什么是混合精度计算？

混合精度计算是一种结合不同精度数值类型的计算策略。在几何计算中，通常结合：

低精度快速计算（如double）：用于大多数情况
高精度精确计算（如Gmpq）：用于关键决策
区间算术：用于判断何时需要精确计算

类比理解：想象你是一位侦探破案：

快速筛选（double）：快速排除明显无关的线索
精确分析（Gmpq）：对关键证据进行详细分析
不确定性评估（interval）：判断哪些线索需要进一步调查

为什么需要混合精度？

几何计算面临一个基本矛盾：

┌─────────────────────────────────────────────────────────────┐
│                     几何计算的困境                           │
├─────────────────────────────────────────────────────────────┤
│                                                             │
│   精确计算                    vs              快速计算       │
│   ─────────                                  ─────────      │
│   • 有理数运算(Gmpq)                         • 浮点运算     │
│   • 无精度误差                              • 硬件加速      │
│   • 结果可靠                                • 内存高效      │
│                                                             │
│   缺点：                                     缺点：          │
│   • 慢100-1000倍                            • 精度误差      │
│   • 内存开销大                              • 可能出错      │
│                                                             │
└─────────────────────────────────────────────────────────────┘

解决方案：智能地结合两者——大部分时候用快速计算，只在必要时使用精确计算。

19.8.2 核心组件

1. 延迟精确数（Lazy_exact_nt）

Lazy_exact_nt是CGAL混合精度系统的核心。它同时维护：

近似值（Interval_nt）：用于快速判断
精确值（ET）：延迟计算，只在需要时求值

┌─────────────────────────────────────────────────────────────┐
│                  Lazy_exact_nt 结构                          │
├─────────────────────────────────────────────────────────────┤
│                                                             │
│   ┌──────────────┐         ┌──────────────┐                │
│   │  Interval_nt │         │   ET (Gmpq)  │                │
│   │  [inf, sup]  │         │   (optional) │                │
│   │              │         │              │                │
│   │  总是存在    │         │  延迟计算    │                │
│   │  快速访问    │         │  缓存结果    │                │
│   └──────────────┘         └──────────────┘                │
│                                                             │
│   决策流程：                                                │
│   1. 使用Interval判断                                       │
│   2. 如果确定（全正/全负）→ 返回结果                        │
│   3. 如果不确定（包含0）→ 计算精确值                        │
│                                                             │
└─────────────────────────────────────────────────────────────┘

2. 过滤机制（Filtering）

过滤是混合精度的关键策略：

// 三级过滤策略
Sign orientation_filtered(const Point& p, const Point& q, const Point& r) {
    // 第1级：静态过滤器（编译期误差界）
    double static_bound = compute_static_bound();
    double approx = compute_approximate(p, q, r);
    if (std::abs(approx) > static_bound) {
        return approx > 0 ? POSITIVE : NEGATIVE;
    }
    
    // 第2级：区间算术
    Interval interval = compute_interval(p, q, r);
    if (interval.is_positive()) return POSITIVE;
    if (interval.is_negative()) return NEGATIVE;
    
    // 第3级：精确计算
    return compute_exact(p, q, r);
}

3. 区间算术（Interval Arithmetic）

区间算术维护数值的可能范围，而非单一值：

计算 a + b：
  a ∈ [1.0, 2.0]
  b ∈ [3.0, 4.0]
  ─────────────
  a+b ∈ [4.0, 6.0]  （向下/向上舍入）

如果结果区间不包含0，则可以确定符号！

19.8.3 代码示例

基础示例：理解Lazy计算

#include <iostream>
#include <memory>
#include <math>
 
// ============================================
// 简化的Lazy_exact_nt实现
// ============================================
 
// 区间类
class Interval {
    double inf_, sup_;
public:
    Interval(double i, double s) : inf_(i), sup_(s) {}
    explicit Interval(double v) : inf_(v), sup_(v) {}
    
    double inf() const { return inf_; }
    double sup() const { return sup_; }
    
    bool is_positive() const { return inf_ > 0; }
    bool is_negative() const { return sup_ < 0; }
    bool contains_zero() const { return inf_ <= 0 && sup_ >= 0; }
    
    // 区间运算
    Interval operator+(const Interval& o) const {
        return Interval(inf_ + o.inf_, sup_ + o.sup_);
    }
    Interval operator-(const Interval& o) const {
        return Interval(inf_ - o.sup_, sup_ - o.inf_);
    }
    Interval operator*(const Interval& o) const {
        double p1 = inf_ * o.inf_, p2 = inf_ * o.sup_;
        double p3 = sup_ * o.inf_, p4 = sup_ * o.sup_;
        return Interval(
            std::min({p1, p2, p3, p4}),
            std::max({p1, p2, p3, p4})
        );
    }
};
 
// 精确数（简化版，实际使用Gmpq）
class ExactRational {
    long num_, den_;
public:
    ExactRational(long n = 0, long d = 1) : num_(n), den_(d) {
        if (den_ < 0) { num_ = -num_; den_ = -den_; }
    }
    
    ExactRational operator+(const ExactRational& o) const {
        return ExactRational(num_ * o.den_ + o.num_ * den_, den_ * o.den_);
    }
    
    ExactRational operator-(const ExactRational& o) const {
        return ExactRational(num_ * o.den_ - o.num_ * den_, den_ * o.den_);
    }
    
    ExactRational operator*(const ExactRational& o) const {
        return ExactRational(num_ * o.num_, den_ * o.den_);
    }
    
    bool is_positive() const { return num_ > 0; }
    bool is_negative() const { return num_ < 0; }
    bool is_zero() const { return num_ == 0; }
    
    double to_double() const { return static_cast<double>(num_) / den_; }
    
    friend std::ostream& operator<<(std::ostream& os, const ExactRational& r) {
        os << r.num_;
        if (r.den_ \!= 1) os << "/" << r.den_;
        return os;
    }
};
 
// 延迟精确数
template <typename ET>
class LazyExact {
    mutable std::unique_ptr<ET> exact_;  // 延迟计算的精确值
    Interval interval_;                   // 近似区间
    mutable bool computed_;
    
public:
    // 从double构造
    explicit LazyExact(double d) 
        : exact_(nullptr), interval_(d), computed_(false) {}
    
    // 从精确值构造
    LazyExact(const ET& e)
        : exact_(std::make_unique<ET>(e)), 
          interval_(e.to_double()),
          computed_(true) {}
    
    // 延迟加法
    LazyExact operator+(const LazyExact& o) const {
        LazyExact result;
        result.interval_ = interval_ + o.interval_;
        // 不立即计算精确值
        result.computed_ = false;
        return result;
    }
    
    // 获取区间（快速）
    const Interval& interval() const { return interval_; }
    
    // 获取精确值（可能触发计算）
    const ET& exact() const {
        if (\!computed_) {
            compute_exact();
        }
        return *exact_;
    }
    
    // 判断符号（使用过滤）
    int sign() const {
        // 第1级：区间判断
        if (interval_.is_positive()) return 1;
        if (interval_.is_negative()) return -1;
        if (\!interval_.contains_zero()) {
            return interval_.inf() > 0 ? 1 : -1;
        }
        
        // 第2级：精确计算
        const ET& e = exact();
        if (e.is_positive()) return 1;
        if (e.is_negative()) return -1;
        return 0;
    }
    
private:
    LazyExact() : exact_(nullptr), interval_(0), computed_(false) {}
    
    void compute_exact() const {
        // 这里简化处理，实际会递归计算操作数
        exact_ = std::make_unique<ET>(ET(static_cast<long>(interval_.inf())));
        computed_ = true;
    }
};
 
int main() {
    // 创建延迟精确数
    LazyExact<ExactRational> a(ExactRational(1, 3));  // 1/3
    LazyExact<ExactRational> b(ExactRational(1, 6));  // 1/6
    
    // 延迟加法（不立即计算）
    auto c = a + b;  // 应该是 1/2
    
    // 使用区间快速判断
    std::cout << "Interval: [" << c.interval().inf() 
              << ", " << c.interval().sup() << "]" << std::endl;
    
    // 获取符号（使用过滤）
    std::cout << "Sign: " << c.sign() << std::endl;
    
    // 获取精确值（触发计算）
    std::cout << "Exact: " << c.exact() << std::endl;
    
    return 0;
}

中级示例：三级过滤系统

#include <iostream>
#include <math>
#include <algorithm>
 
// ============================================
// 三级过滤方向判断
// ============================================
 
enum Sign { NEGATIVE = -1, ZERO = 0, POSITIVE = 1, UNCERTAIN = 2 };
 
struct Point {
    double x, y;
    Point(double x_ = 0, double y_ = 0) : x(x_), y(y_) {}
};
 
// 第1级：静态过滤器（使用预计算的误差界）
class StaticFilter {
    static constexpr double epsilon = 1e-14;  // 误差界
    
public:
    static Sign orientation(const Point& p, const Point& q, const Point& r) {
        double pq_x = q.x - p.x;
        double pq_y = q.y - p.y;
        double pr_x = r.x - p.x;
        double pr_y = r.y - p.y;
        
        double result = pq_x * pr_y - pq_y * pr_x;
        
        // 计算误差界（简化版）
        double abs_sum = std::abs(pq_x * pr_y) + std::abs(pq_y * pr_x);
        double error_bound = epsilon * abs_sum;
        
        if (result > error_bound) return POSITIVE;
        if (result < -error_bound) return NEGATIVE;
        return UNCERTAIN;
    }
};
 
// 第2级：区间过滤器
class IntervalFilter {
public:
    static Sign orientation(const Point& p, const Point& q, const Point& r) {
        // 使用区间算术计算
        // 这里简化，实际使用Interval_nt
        
        double pq_x[2] = { q.x - p.x, q.x - p.x };
        double pq_y[2] = { q.y - p.y, q.y - p.y };
        double pr_x[2] = { r.x - p.x, r.x - p.x };
        double pr_y[2] = { r.y - p.y, r.y - p.y };
        
        // 模拟区间乘法（考虑舍入）
        double prod1[2] = {
            std::nextafter(pq_x[0] * pr_y[0], -INFINITY),
            std::nextafter(pq_x[1] * pr_y[1], INFINITY)
        };
        double prod2[2] = {
            std::nextafter(pq_y[0] * pr_x[0], -INFINITY),
            std::nextafter(pq_y[1] * pr_x[1], INFINITY)
        };
        
        double result_inf = std::nextafter(prod1[0] - prod2[1], -INFINITY);
        double result_sup = std::nextafter(prod1[1] - prod2[0], INFINITY);
        
        if (result_inf > 0) return POSITIVE;
        if (result_sup < 0) return NEGATIVE;
        if (result_inf <= 0 && result_sup >= 0) {
            // 检查是否为精确0
            if (result_inf == 0 && result_sup == 0) return ZERO;
        }
        return UNCERTAIN;
    }
};
 
// 第3级：精确计算（使用扩展精度）
class ExactFilter {
public:
    static Sign orientation(const Point& p, const Point& q, const Point& r) {
        // 使用扩展精度（long double）
        long double pq_x = (long double)q.x - p.x;
        long double pq_y = (long double)q.y - p.y;
        long double pr_x = (long double)r.x - p.x;
        long double pr_y = (long double)r.y - p.y;
        
        long double result = pq_x * pr_y - pq_y * pr_x;
        
        if (result > 0) return POSITIVE;
        if (result < 0) return NEGATIVE;
        return ZERO;
    }
};
 
// 完整的过滤系统
Sign robust_orientation(const Point& p, const Point& q, const Point& r) {
    // 第1级：静态过滤
    Sign s = StaticFilter::orientation(p, q, r);
    if (s \!= UNCERTAIN) {
        std::cout << "[Static filter succeeded] ";
        return s;
    }
    
    // 第2级：区间过滤
    s = IntervalFilter::orientation(p, q, r);
    if (s \!= UNCERTAIN) {
        std::cout << "[Interval filter succeeded] ";
        return s;
    }
    
    // 第3级：精确计算
    std::cout << "[Exact computation required] ";
    return ExactFilter::orientation(p, q, r);
}
 
int main() {
    // 测试1：明显的左转（静态过滤器应该成功）
    Point p1(0, 0), q1(1, 0), r1(0, 1);
    std::cout << "Test 1 (clear left turn): ";
    std::cout << robust_orientation(p1, q1, r1) << std::endl;
    
    // 测试2：接近共线（可能需要更高级别的过滤）
    Point p2(0, 0), q2(1, 0), r2(2, 1e-15);
    std::cout << "Test 2 (near collinear): ";
    std::cout << robust_orientation(p2, q2, r2) << std::endl;
    
    // 测试3：共线
    Point p3(0, 0), q3(1, 1), r3(2, 2);
    std::cout << "Test 3 (collinear): ";
    std::cout << robust_orientation(p3, q3, r3) << std::endl;
    
    return 0;
}

高级示例：完整的Lazy_exact_nt模拟

#include <iostream>
#include <memory>
#include <vector>
#include <math>
 
// ============================================
// 完整的Lazy_exact_nt模拟
// ============================================
 
// 精确数类型（模拟Gmpq）
class ExactNumber {
    long num_, den_;
public:
    ExactNumber(long n = 0, long d = 1) : num_(n), den_(d) {
        normalize();
    }
    
    void normalize() {
        if (den_ < 0) { num_ = -num_; den_ = -den_; }
        long g = gcd(std::abs(num_), den_);
        num_ /= g; den_ /= g;
    }
    
    static long gcd(long a, long b) {
        return b == 0 ? a : gcd(b, a % b);
    }
    
    ExactNumber operator+(const ExactNumber& o) const {
        return ExactNumber(num_ * o.den_ + o.num_ * den_, den_ * o.den_);
    }
    
    ExactNumber operator-(const ExactNumber& o) const {
        return ExactNumber(num_ * o.den_ - o.num_ * den_, den_ * o.den_);
    }
    
    ExactNumber operator*(const ExactNumber& o) const {
        return ExactNumber(num_ * o.num_, den_ * o.den_);
    }
    
    ExactNumber operator-() const {
        return ExactNumber(-num_, den_);
    }
    
    double to_double() const {
        return static_cast<double>(num_) / den_;
    }
    
    bool is_zero() const { return num_ == 0; }
    int sign() const {
        if (num_ > 0) return 1;
        if (num_ < 0) return -1;
        return 0;
    }
    
    friend std::ostream& operator<<(std::ostream& os, const ExactNumber& n) {
        os << n.num_;
        if (n.den_ \!= 1) os << "/" << n.den_;
        return os;
    }
};
 
// 区间数
class IntervalNumber {
    double inf_, sup_;
public:
    IntervalNumber(double i, double s) : inf_(i), sup_(s) {}
    explicit IntervalNumber(double v) : inf_(v), sup_(v) {}
    
    static IntervalNumber from_exact(const ExactNumber& e) {
        double d = e.to_double();
        return IntervalNumber(d, d);  // 简化，实际应该考虑舍入
    }
    
    double inf() const { return inf_; }
    double sup() const { return sup_; }
    
    bool is_positive() const { return inf_ > 0; }
    bool is_negative() const { return sup_ < 0; }
    bool contains_zero() const { return inf_ <= 0 && sup_ >= 0; }
    bool is_point() const { return inf_ == sup_; }
    
    IntervalNumber operator+(const IntervalNumber& o) const {
        return IntervalNumber(inf_ + o.inf_, sup_ + o.sup_);
    }
    
    IntervalNumber operator-(const IntervalNumber& o) const {
        return IntervalNumber(inf_ - o.sup_, sup_ - o.inf_);
    }
    
    IntervalNumber operator*(const IntervalNumber& o) const {
        double p[4] = { inf_ * o.inf_, inf_ * o.sup_,
                       sup_ * o.inf_, sup_ * o.sup_ };
        return IntervalNumber(
            *std::min_element(p, p + 4),
            *std::max_element(p, p + 4)
        );
    }
};
 
// 延迟精确数表示基类
template <typename ET>
class LazyRep {
public:
    virtual ~LazyRep() = default;
    virtual ET exact() const = 0;
    virtual int depth() const = 0;
};
 
// 常量表示
template <typename ET>
class LazyConstant : public LazyRep<ET> {
    ET value_;
public:
    explicit LazyConstant(const ET& v) : value_(v) {}
    ET exact() const override { return value_; }
    int depth() const override { return 0; }
};
 
// 一元操作表示
template <typename ET>
class LazyUnary : public LazyRep<ET> {
public:
    mutable std::shared_ptr<LazyRep<ET>> op;
    mutable ET cached_;
    mutable bool computed_;
    
    explicit LazyUnary(std::shared_ptr<LazyRep<ET>> o) 
        : op(std::move(o)), computed_(false) {}
    
    int depth() const override { return op->depth() + 1; }
    
protected:
    const ET& get_op_exact() const {
        if (\!computed_) {
            cached_ = op->exact();
            computed_ = true;
        }
        return cached_;
    }
};
 
// 取反操作
template <typename ET>
class LazyNegate : public LazyUnary<ET> {
public:
    using LazyUnary<ET>::LazyUnary;
    
    ET exact() const override {
        return -this->get_op_exact();
    }
};
 
// 二元操作表示
template <typename ET>
class LazyBinary : public LazyRep<ET> {
public:
    mutable std::shared_ptr<LazyRep<ET>> left, right;
    mutable ET cached_;
    mutable bool computed_;
    
    LazyBinary(std::shared_ptr<LazyRep<ET>> l, std::shared_ptr<LazyRep<ET>> r)
        : left(std::move(l)), right(std::move(r)), computed_(false) {}
    
    int depth() const override { 
        return std::max(left->depth(), right->depth()) + 1; 
    }
    
protected:
    const ET& get_left_exact() const {
        if (\!computed_) {
            cached_left_ = left->exact();
        }
        return cached_left_;
    }
    
    const ET& get_right_exact() const {
        if (\!computed_) {
            cached_right_ = right->exact();
        }
        return cached_right_;
    }
    
private:
    mutable ET cached_left_, cached_right_;
};
 
// 加法操作
template <typename ET>
class LazyAdd : public LazyBinary<ET> {
public:
    using LazyBinary<ET>::LazyBinary;
    
    ET exact() const override {
        return this->get_left_exact() + this->get_right_exact();
    }
};
 
// 延迟精确数主类
template <typename ET>
class LazyExact {
    std::shared_ptr<LazyRep<ET>> rep_;
    IntervalNumber interval_;
    
public:
    // 从精确值构造
    LazyExact(const ET& e) 
        : rep_(std::make_shared<LazyConstant<ET>>(e)),
          interval_(IntervalNumber::from_exact(e)) {}
    
    // 从double构造
    explicit LazyExact(double d)
        : rep_(std::make_shared<LazyConstant<ET>>(ET(static_cast<long>(d * 1000), 1000))),
          interval_(d) {}
    
    // 从表示构造
    LazyExact(std::shared_ptr<LazyRep<ET>> rep, const IntervalNumber& i)
        : rep_(std::move(rep)), interval_(i) {}
    
    // 获取区间（快速）
    const IntervalNumber& interval() const { return interval_; }
    
    // 获取精确值（可能触发计算）
    ET exact() const { return rep_->exact(); }
    
    // 获取表达式深度
    int depth() const { return rep_ ? rep_->depth() : 0; }
    
    // 判断符号（过滤）
    int sign() const {
        if (interval_.is_positive()) return 1;
        if (interval_.is_negative()) return -1;
        if (interval_.is_point() && interval_.inf() == 0) return 0;
        // 需要精确计算
        return exact().sign();
    }
    
    // 运算符
    LazyExact operator+(const LazyExact& o) const {
        auto new_rep = std::make_shared<LazyAdd<ET>>(rep_, o.rep_);
        return LazyExact(new_rep, interval_ + o.interval_);
    }
    
    LazyExact operator-() const {
        auto new_rep = std::make_shared<LazyNegate<ET>>(rep_);
        return LazyExact(new_rep, IntervalNumber(-interval_.sup(), -interval_.inf()));
    }
    
    LazyExact operator-(const LazyExact& o) const {
        return *this + (-o);
    }
    
    friend std::ostream& operator<<(std::ostream& os, const LazyExact& l) {
        os << "Lazy[" << l.interval_.inf() << ", " << l.interval_.sup() <> "]";
        return os;
    }
};
 
// 使用示例
int main() {
    typedef LazyExact<ExactNumber> LazyNum;
    
    // 创建延迟精确数
    LazyNum a(ExactNumber(1, 3));  // 1/3
    LazyNum b(ExactNumber(1, 6));  // 1/6
    
    std::cout << "a = " << a << std::endl;
    std::cout << "b = " << b << std::endl;
    
    // 构建表达式树（延迟计算）
    LazyNum c = a + b;  // 不立即计算
    std::cout << "c = a + b = " << c << std::endl;
    std::cout << "Expression depth: " << c.depth() << std::endl;
    
    // 判断符号（使用过滤）
    std::cout << "sign(c) = " << c.sign() << std::endl;
    
    // 获取精确值（触发计算）
    std::cout << "exact(c) = " << c.exact() << std::endl;
    
    // 复杂表达式
    LazyNum d = a + b + a;  // 1/3 + 1/6 + 1/3 = 5/6
    std::cout << "\nd = a + b + a = " << d << std::endl;
    std::cout << "Expression depth: " << d.depth() << std::endl;
    std::cout << "exact(d) = " << d.exact() << std::endl;
    
    return 0;
}

19.8.4 CGAL中的应用

CGAL的Lazy_exact_nt

// 来自CGAL/Lazy_exact_nt.h的核心部分（简化展示）
namespace CGAL {
 
template <class ET>
class Lazy_exact_nt 
  : public Lazy<Interval_nt<false>, ET, To_interval<ET> >
{
public:
    typedef ET ET_type;
    typedef Interval_nt<false> AT_type;
    typedef Lazy<AT_type, ET_type, To_interval<ET_type> > Base;
    
    // 从double构造（近似值）
    Lazy_exact_nt(double d) : Base(new Lazy_exact_Cst<ET>(d)) {}
    
    // 从精确值构造
    Lazy_exact_nt(const ET& e) : Base(new Lazy_exact_Ex_Cst<ET>(e)) {}
    
    // 延迟加法
    Lazy_exact_nt operator+(const Lazy_exact_nt& b) const {
        return Lazy_exact_nt(
            new Lazy_exact_Add<ET>(*this, b),
            this->approx() + b.approx()  // 同时计算近似值
        );
    }
    
    // 获取精确值（触发计算）
    const ET& exact() const {
        if (this->is_lazy()) {
            this->update_exact();  // 计算并缓存
        }
        return *this->ptr();
    }
    
    // 获取近似值（快速）
    const AT_type& approx() const { return this->approx_; }
};
 
} // namespace CGAL

CGAL的Exact_predicates_inexact_constructions_kernel

#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
 
typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
typedef K::Point_2 Point_2;
 
void example() {
    // 构造使用double（快速）
    Point_2 p(0.1, 0.2);  // 使用double存储
    Point_2 q(0.3, 0.4);
    Point_2 r(0.5, 0.6);
    
    // 谓词使用精确计算（鲁棒）
    K::Orientation_2 orientation;
    auto result = orientation(p, q, r);  // 内部使用过滤精确计算
}

CGAL的Exact_predicates_exact_constructions_kernel

#include <CGAL/Exact_predicates_exact_constructions_kernel.h>
 
typedef CGAL::Exact_predicates_exact_constructions_kernel K;
typedef K::Point_2 Point_2;
 
typedef CGAL::Exact_predicates_inexact_constructions_kernel K_fast;
 
void comparison() {
    // 快速内核：构造快，但可能不精确
    K_fast::Point_2 p_fast(0.1, 0.2);
    K_fast::Point_2 q_fast(0.3, 0.6);
    // p_fast + (q_fast - p_fast) * 0.5 可能不精确！
    
    // 精确内核：构造使用Lazy_exact_nt
    K::Point_2 p_exact(0.1, 0.2);
    K::Point_2 q_exact(0.3, 0.6);
    // p_exact + (q_exact - p_exact) * 0.5 是精确的！
}

19.8.5 常见陷阱

陷阱1：过度使用精确计算

// 错误：所有计算都使用精确内核
typedef CGAL::Exact_predicates_exact_constructions_kernel K;
K::Point_2 p(0, 0);
 
for (int i = 0; i < 1000000; ++i) {
    p = p + K::Vector_2(1, 0);  // 每次迭代都创建延迟计算节点！
}
// 结果：内存爆炸，速度极慢
 
// 正确：只在需要时使用精确计算
typedef CGAL::Exact_predicates_inexact_constructions_kernel K_fast;
K_fast::Point_2 p_fast(0, 0);
 
for (int i = 0; i < 1000000; ++i) {
    p_fast = p_fast + K_fast::Vector_2(1, 0);  // 使用double，快速
}

陷阱2：忽略区间不确定性

// 问题：假设区间总是确定的
Lazy_exact_nt<Gmpq> a = ...;
Lazy_exact_nt<Gmpq> b = ...;
 
if (a.approx().is_positive()) {  // 可能不确定！
    // 错误：区间可能包含0
}
 
// 正确：处理不确定性
if (a.approx().is_positive()) {
    // 确定为正
} else if (a.approx().contains_zero()) {
    // 不确定，需要精确值
    if (a.exact().is_positive()) {
        // 实际上是正的
    }
}

陷阱3：比较延迟计算对象

// 危险：直接比较延迟计算对象
Lazy_exact_nt<Gmpq> a = ...;
Lazy_exact_nt<Gmpq> b = ...;
 
if (a == b) {  // 这可能触发精确计算！
    // ...
}
 
// 更好的做法：先比较近似值
if (a.approx().is_same(b.approx())) {
    // 快速路径：近似值相同
} else {
    // 需要精确比较
    if (a.exact() == b.exact()) {
        // ...
    }
}

19.8.6 最佳实践

实践1：选择合适的内核

// 场景1：只需要谓词（方向判断、点在内部等）
// 使用 EPEC（精确谓词，不精确构造）
typedef CGAL::Exact_predicates_inexact_constructions_kernel K1;
 
// 场景2：需要精确构造（交点计算）
// 使用 EPEC（完全精确）
typedef CGAL::Exact_predicates_exact_constructions_kernel K2;
 
// 场景3：需要平方根（距离计算）
// 使用带平方根的精确内核
typedef CGAL::Exact_predicates_exact_constructions_kernel_with_sqrt K3;

实践2：批量处理策略

// 低效：逐个处理，频繁触发精确计算
for (auto v : mesh.vertices()) {
    auto neighbors = get_neighbors(v);
    if (neighbors.size() > 2) {  // 每次判断都可能触发精确计算
        process(v);
    }
}
 
// 高效：先批量收集，再统一处理
std::vector<Vertex> to_process;
for (auto v : mesh.vertices()) {
    auto neighbors = get_neighbors(v);
    to_process.push_back(v);  // 只收集，不判断
}
 
for (auto v : to_process) {
    process(v);  // 统一处理
}

实践3：监控性能

// 使用CGAL的性能计数器
#include <CGAL/Profile_counter.h>
 
void algorithm() {
    CGAL_PROFILER("algorithm_total");
    
    for (...) {
        CGAL_PROFILER("loop_iteration");
        
        // 谓词调用
        CGAL_PROFILER("orientation_test");
        auto orient = orientation(p, q, r);
    }
}
 
// 输出统计信息
// 设置环境变量 CGAL_PROFILE 为 "y" 查看统计

实践4：混合使用不同精度

// 使用不同精度的内核
typedef CGAL::Exact_predicates_inexact_constructions_kernel Fast_kernel;
typedef CGAL::Exact_predicates_exact_constructions_kernel Exact_kernel;
 
// 快速计算
typedef Fast_kernel::Point_2 Fast_point;
// 精确计算
typedef Exact_kernel::Point_2 Exact_point;
 
// 在需要时转换
Fast_point fast_p(0.1, 0.2);
Exact_point exact_p(fast_p.x(), fast_p.y());  // 转换为精确点

本节要点

混合精度平衡性能与鲁棒性：大部分计算使用快速近似，关键决策使用精确计算。
Lazy_exact_nt延迟计算：构建表达式树，只在需要时计算精确值。
三级过滤策略：静态过滤 → 区间过滤 → 精确计算，逐级提升精度。
区间算术检测不确定性：通过维护数值范围，判断何时需要精确计算。
内核选择影响整体性能：根据应用需求选择合适的内核类型。
避免过度使用精确计算：精确计算有显著开销，应谨慎使用。

进一步阅读

CGAL文档：Number Types模块的Lazy_exact_nt章节
《Robust Geometric Computation》：几何计算鲁棒性的经典教材
CGAL论文：“The Design and Implementation of CGAL”
区间算术资源：IEEE 1788标准

Chunibyo

Explorer

19.8 混合精度计算

19.8 混合精度计算

相关章节

在本节中你将学到…

19.8.1 概念解释

什么是混合精度计算？

为什么需要混合精度？

19.8.2 核心组件

1. 延迟精确数（Lazy_exact_nt）

2. 过滤机制（Filtering）

3. 区间算术（Interval Arithmetic）

19.8.3 代码示例

基础示例：理解Lazy计算

中级示例：三级过滤系统

高级示例：完整的Lazy_exact_nt模拟

19.8.4 CGAL中的应用

CGAL的Lazy_exact_nt

CGAL的Exact_predicates_inexact_constructions_kernel

CGAL的Exact_predicates_exact_constructions_kernel

19.8.5 常见陷阱

陷阱1：过度使用精确计算

陷阱2：忽略区间不确定性

陷阱3：比较延迟计算对象

19.8.6 最佳实践

实践1：选择合适的内核

实践2：批量处理策略

实践3：监控性能

实践4：混合使用不同精度

本节要点

进一步阅读

Graph View

Table of Contents

Backlinks

Chunibyo

Explorer

19.8 混合精度计算

19.8 混合精度计算 §

相关章节 §

在本节中你将学到… §

19.8.1 概念解释 §

什么是混合精度计算？ §

为什么需要混合精度？ §

19.8.2 核心组件 §

1. 延迟精确数（Lazy_exact_nt） §

2. 过滤机制（Filtering） §

3. 区间算术（Interval Arithmetic） §

19.8.3 代码示例 §

基础示例：理解Lazy计算 §

中级示例：三级过滤系统 §

高级示例：完整的Lazy_exact_nt模拟 §

19.8.4 CGAL中的应用 §

CGAL的Lazy_exact_nt §

CGAL的Exact_predicates_inexact_constructions_kernel §

CGAL的Exact_predicates_exact_constructions_kernel §

19.8.5 常见陷阱 §

陷阱1：过度使用精确计算 §

陷阱2：忽略区间不确定性 §

陷阱3：比较延迟计算对象 §

19.8.6 最佳实践 §

实践1：选择合适的内核 §

实践2：批量处理策略 §

实践3：监控性能 §

实践4：混合使用不同精度 §

本节要点 §

进一步阅读 §

Graph View

Table of Contents

Backlinks

19.8 混合精度计算

相关章节

在本节中你将学到…

19.8.1 概念解释

什么是混合精度计算？

为什么需要混合精度？

19.8.2 核心组件

1. 延迟精确数（Lazy_exact_nt）

2. 过滤机制（Filtering）

3. 区间算术（Interval Arithmetic）

19.8.3 代码示例

基础示例：理解Lazy计算

中级示例：三级过滤系统

高级示例：完整的Lazy_exact_nt模拟

19.8.4 CGAL中的应用

CGAL的Lazy_exact_nt

CGAL的Exact_predicates_inexact_constructions_kernel

CGAL的Exact_predicates_exact_constructions_kernel

19.8.5 常见陷阱

陷阱1：过度使用精确计算

陷阱2：忽略区间不确定性

陷阱3：比较延迟计算对象

19.8.6 最佳实践

实践1：选择合适的内核

实践2：批量处理策略

实践3：监控性能

实践4：混合使用不同精度

本节要点

进一步阅读