相关章节: 笛卡尔与齐次坐标系统 | 过滤内核与精确谓词

2.3 代数基础与数类型

理论基础

代数结构层次

CGAL的代数基础建立在抽象的代数结构概念之上,形成清晰的层次结构:

IntegralDomainWithoutDivision  整环(无除法)
    └── IntegralDomain         整环
            └── UniqueFactorizationDomain  唯一分解域
                    └── EuclideanRing      欧几里得环
                            └── Field      域
                                    └── FieldWithSqrt      带平方根的域
                                    └── FieldWithKthRoot   带k次方根的域
                                    └── FieldWithRootOf    带多项式根的域

数学定义

  1. 整环 (Integral Domain)

    • 交换环,含乘法单位元1
    • 无零因子:
    • 示例:整数
  2. 欧几里得环 (Euclidean Ring)

    • 整环,带欧几里得范数
    • 可执行带余除法:,其中
    • 示例:整数 ,多项式环
  3. 域 (Field)

    • 整环,每个非零元都有乘法逆元
    • 支持四则运算(除数非零)
    • 示例:有理数 ,实数 ,复数

Real_embeddable 概念

Real_embeddable 概念描述可以嵌入实数轴的数类型:

要求

  • 可与 double 相互转换
  • 支持比较运算(全序关系)
  • 可计算符号(正/负/零)
  • 可计算绝对值

数学意义:这类数类型可用于几何计算中的坐标表示和谓词评估。

数域扩展

CGAL支持通过代数扩展构造新的数域:

  1. 商域 (Quotient):从整环构造域

    • 示例:
  2. 根式扩展:向域中添加多项式根

架构分析

代数结构特征类

头文件位置CGAL/Algebraic_structure_traits.h

// 代数结构特征主模板
template < typename Type >
struct Algebraic_structure_traits {
    // 默认情况下未定义,需要特化
};
 
// 代数结构标签层次
tag IntegralDomainWithoutDivision_tag;
tag IntegralDomain_tag           : IntegralDomainWithoutDivision_tag;
tag UniqueFactorizationDomain_tag: IntegralDomain_tag;
tag EuclideanRing_tag            : UniqueFactorizationDomain_tag;
tag Field_tag                    : UniqueFactorizationDomain_tag;
tag FieldWithSqrt_tag            : Field_tag;
tag FieldWithKthRoot_tag         : FieldWithSqrt_tag;
tag FieldWithRootOf_tag          : FieldWithKthRoot_tag;
 
// 代数结构特征特化示例:Gmpz(GMP整数)
template <>
struct Algebraic_structure_traits<CGAL::Gmpz> {
    typedef CGAL::Gmpz Type;
    typedef Tag_true Is_exact;
    typedef Tag_true Is_numerical_sensitive;
    typedef EuclideanRing_tag Algebraic_category;
    
    // 整除运算
    struct Integral_division {
        typedef CGAL::Gmpz result_type;
        CGAL::Gmpz operator()(const CGAL::Gmpz& a, const CGAL::Gmpz& b) const {
            CGAL_precondition(b != 0);
            CGAL_precondition(a % b == 0);
            return a / b;
        }
    };
    
    // GCD运算
    struct Gcd {
        typedef CGAL::Gmpz result_type;
        CGAL::Gmpz operator()(const CGAL::Gmpz& a, const CGAL::Gmpz& b) const {
            return ::gcd(a, b);  // 调用GMP的gcd
        }
    };
    
    // 模除运算
    struct Div {
        typedef CGAL::Gmpz result_type;
        CGAL::Gmpz operator()(const CGAL::Gmpz& a, const CGAL::Gmpz& b) const {
            return a / b;
        }
    };
    
    struct Mod {
        typedef CGAL::Gmpz result_type;
        CGAL::Gmpz operator()(const CGAL::Gmpz& a, const CGAL::Gmpz& b) const {
            return a % b;
        }
    };
};
 
// 代数结构特征特化:Gmpq(GMP有理数)
template <>
struct Algebraic_structure_traits<CGAL::Gmpq> {
    typedef CGAL::Gmpq Type;
    typedef Tag_true Is_exact;
    typedef Tag_true Is_numerical_sensitive;
    typedef Field_tag Algebraic_category;
    
    // 域不需要整除、GCD等运算
};

Real_embeddable_traits 实现

头文件位置CGAL/number_type_basic.h, CGAL/Real_embeddable_traits.h

// Real_embeddable_traits 主模板
template < typename Type >
struct Real_embeddable_traits {
    // 默认未定义
};
 
// Gmpz 的特化
template <>
struct Real_embeddable_traits<CGAL::Gmpz> {
    typedef CGAL::Gmpz Type;
    typedef Tag_true Is_real_embeddable;
    
    // 转换为double
    struct To_double {
        typedef double result_type;
        double operator()(const CGAL::Gmpz& x) const {
            return x.to_double();
        }
    };
    
    // 转换为区间(用于过滤)
    struct To_interval {
        typedef std::pair<double, double> result_type;
        std::pair<double, double> operator()(const CGAL::Gmpz& x) const {
            double d = x.to_double();
            // 处理转换误差
            if (x.is_finite()) {
                return std::make_pair(d, d);
            }
            // 溢出处理...
        }
    };
    
    // 符号运算
    struct Sgn {
        typedef CGAL::Sign result_type;
        CGAL::Sign operator()(const CGAL::Gmpz& x) const {
            if (x > 0) return CGAL::POSITIVE;
            if (x < 0) return CGAL::NEGATIVE;
            return CGAL::ZERO;
        }
    };
    
    // 比较运算
    struct Compare {
        typedef CGAL::Comparison_result result_type;
        CGAL::Comparison_result operator()(const CGAL::Gmpz& a, 
                                           const CGAL::Gmpz& b) const {
            if (a < b) return CGAL::SMALLER;
            if (a > b) return CGAL::LARGER;
            return CGAL::EQUAL;
        }
    };
};

GMP 集成架构

头文件位置CGAL/Gmpz.h, CGAL/Gmpq.h, CGAL/Gmpfr.h

// Gmpz - GMP整数包装类
class Gmpz {
protected:
    mpz_t mpz;  // GMP整数对象
    
public:
    // 构造函数
    Gmpz() { mpz_init(mpz); }
    Gmpz(int i) { mpz_init_set_si(mpz, i); }
    Gmpz(long l) { mpz_init_set_si(mpz, l); }
    Gmpz(const char* s) { mpz_init_set_str(mpz, s, 10); }
    Gmpz(const Gmpz& other) { mpz_init_set(mpz, other.mpz); }
    
    // 析构函数
    ~Gmpz() { mpz_clear(mpz); }
    
    // 赋值运算符
    Gmpz& operator=(const Gmpz& other) {
        if (this != &other) mpz_set(mpz, other.mpz);
        return *this;
    }
    
    // 算术运算符
    Gmpz operator+(const Gmpz& other) const {
        Gmpz result;
        mpz_add(result.mpz, mpz, other.mpz);
        return result;
    }
    
    Gmpz operator-(const Gmpz& other) const {
        Gmpz result;
        mpz_sub(result.mpz, mpz, other.mpz);
        return result;
    }
    
    Gmpz operator*(const Gmpz& other) const {
        Gmpz result;
        mpz_mul(result.mpz, mpz, other.mpz);
        return result;
    }
    
    Gmpz operator/(const Gmpz& other) const {
        CGAL_precondition(other != 0);
        Gmpz result;
        mpz_div(result.mpz, mpz, other.mpz);
        return result;
    }
    
    Gmpz operator%(const Gmpz& other) const {
        Gmpz result;
        mpz_mod(result.mpz, mpz, other.mpz);
        return result;
    }
    
    // 比较运算符
    bool operator==(const Gmpz& other) const {
        return mpz_cmp(mpz, other.mpz) == 0;
    }
    
    bool operator<(const Gmpz& other) const {
        return mpz_cmp(mpz, other.mpz) < 0;
    }
    
    // 转换函数
    double to_double() const { return mpz_get_d(mpz); }
    long to_long() const { return mpz_get_si(mpz); }
    
    // GMP原生访问
    const mpz_t& mp() const { return mpz; }
    mpz_t& mp() { return mpz; }
};
 
// Gmpq - GMP有理数包装类
class Gmpq {
protected:
    mpq_t mpq;  // GMP有理数对象
    
public:
    // 构造函数
    Gmpq() { mpq_init(mpq); }
    Gmpq(int n) { mpq_init(mpq); mpq_set_si(mpq, n, 1); }
    Gmpq(long n) { mpq_init(mpq); mpq_set_si(mpq, n, 1); }
    Gmpq(const Gmpz& n) { 
        mpq_init(mpq); 
        mpq_set_z(mpq, n.mp()); 
    }
    Gmpq(const Gmpz& n, const Gmpz& d) {
        mpq_init(mpq);
        mpz_set(mpq_numref(mpq), n.mp());
        mpz_set(mpq_denref(mpq), d.mp());
        mpq_canonicalize(mpq);  // 约分
    }
    
    // 析构函数
    ~Gmpq() { mpq_clear(mpq); }
    
    // 算术运算符(自动约分)
    Gmpq operator+(const Gmpq& other) const {
        Gmpq result;
        mpq_add(result.mpq, mpq, other.mpq);
        return result;
    }
    
    Gmpq operator-(const Gmpq& other) const {
        Gmpq result;
        mpq_sub(result.mpq, mpq, other.mpq);
        return result;
    }
    
    Gmpq operator*(const Gmpq& other) const {
        Gmpq result;
        mpq_mul(result.mpq, mpq, other.mpq);
        return result;
    }
    
    Gmpq operator/(const Gmpq& other) const {
        CGAL_precondition(other != 0);
        Gmpq result;
        mpq_div(result.mpq, mpq, other.mpq);
        return result;
    }
    
    // 访问分子分母
    Gmpz numerator() const {
        Gmpz result;
        mpz_set(result.mp(), mpq_numref(mpq));
        return result;
    }
    
    Gmpz denominator() const {
        Gmpz result;
        mpz_set(result.mp(), mpq_denref(mpq));
        return result;
    }
};

实现细节

Quotient 类模板

头文件位置CGAL/Quotient.h

// Quotient 类 - 从整环构造域
template < typename NT >
class Quotient {
public:
    typedef NT Type;
    typedef typename Algebraic_structure_traits<NT>::Algebraic_category 
        Numerator_algebraic_category;
    
private:
    NT num_;  // 分子
    NT den_;  // 分母(始终为正,且与分子互质)
    
    // 规范化:确保分母为正,且分子分母互质
    void normalize() {
        if (den_ < NT(0)) {
            num_ = -num_;
            den_ = -den_;
        }
        // 可选:约分(使用GCD)
        // NT g = gcd(abs(num_), den_);
        // num_ /= g; den_ /= g;
    }
    
public:
    // 构造函数
    Quotient() : num_(0), den_(1) {}
    Quotient(const NT& n) : num_(n), den_(1) {}
    Quotient(const NT& n, const NT& d) : num_(n), den_(d) {
        CGAL_precondition(d != NT(0));
        normalize();
    }
    
    // 从其他数类型构造
    template < typename T >
    Quotient(T t) : num_(NT(t)), den_(NT(1)) {}
    
    // 访问器
    const NT& numerator() const { return num_; }
    const NT& denominator() const { return den_; }
    
    // 算术运算符
    Quotient operator+(const Quotient& other) const {
        // a/b + c/d = (ad + bc) / bd
        return Quotient(num_ * other.den_ + other.num_ * den_, 
                       den_ * other.den_);
    }
    
    Quotient operator-(const Quotient& other) const {
        return Quotient(num_ * other.den_ - other.num_ * den_, 
                       den_ * other.den_);
    }
    
    Quotient operator*(const Quotient& other) const {
        return Quotient(num_ * other.num_, den_ * other.den_);
    }
    
    Quotient operator/(const Quotient& other) const {
        CGAL_precondition(other.num_ != NT(0));
        return Quotient(num_ * other.den_, den_ * other.num_);
    }
    
    // 一元运算符
    Quotient operator-() const { return Quotient(-num_, den_); }
    Quotient operator+() const { return *this; }
    
    // 比较运算符
    bool operator==(const Quotient& other) const {
        return num_ * other.den_ == other.num_ * den_;
    }
    
    bool operator<(const Quotient& other) const {
        return num_ * other.den_ < other.num_ * den_;
    }
    
    // 转换为double
    double to_double() const {
        return num_.to_double() / den_.to_double();
    }
};
 
// Quotient 的代数结构特征
template < typename NT >
struct Algebraic_structure_traits< Quotient<NT> > {
    typedef Quotient<NT> Type;
    typedef typename Algebraic_structure_traits<NT>::Is_exact Is_exact;
    typedef typename Algebraic_structure_traits<NT>::Is_numerical_sensitive 
        Is_numerical_sensitive;
    typedef Field_tag Algebraic_category;
};

MPFR 浮点实现

头文件位置CGAL/Gmpfr.h

// Gmpfr - GMP多精度浮点数(基于MPFR库)
class Gmpfr {
protected:
    mpfr_t mpfr;  // MPFR浮点对象
    
public:
    // 精度控制(全局默认精度)
    static mpfr_prec_t get_default_precision() {
        return mpfr_get_default_prec();
    }
    
    static void set_default_precision(mpfr_prec_t prec) {
        mpfr_set_default_prec(prec);
    }
    
    // 构造函数
    Gmpfr() { mpfr_init(mpfr); }
    
    Gmpfr(double d) { 
        mpfr_init_set_d(mpfr, d, MPFR_RNDN); 
    }
    
    Gmpfr(const Gmpz& z) {
        mpfr_init_set_z(mpfr, z.mp(), MPFR_RNDN);
    }
    
    Gmpfr(const Gmpq& q) {
        mpfr_init_set_q(mpfr, q.mp(), MPFR_RNDN);
    }
    
    // 指定精度构造
    Gmpfr(const Gmpfr& other, mpfr_prec_t prec) {
        mpfr_init2(mpfr, prec);
        mpfr_set(mpfr, other.mpfr, MPFR_RNDN);
    }
    
    // 析构函数
    ~Gmpfr() { mpfr_clear(mpfr); }
    
    // 算术运算(支持舍入模式)
    Gmpfr operator+(const Gmpfr& other) const {
        Gmpfr result;
        mpfr_add(result.mpfr, mpfr, other.mpfr, MPFR_RNDN);
        return result;
    }
    
    Gmpfr operator-(const Gmpfr& other) const {
        Gmpfr result;
        mpfr_sub(result.mpfr, mpfr, other.mpfr, MPFR_RNDN);
        return result;
    }
    
    Gmpfr operator*(const Gmpfr& other) const {
        Gmpfr result;
        mpfr_mul(result.mpfr, mpfr, other.mpfr, MPFR_RNDN);
        return result;
    }
    
    Gmpfr operator/(const Gmpfr& other) const {
        CGAL_precondition(other != 0);
        Gmpfr result;
        mpfr_div(result.mpfr, mpfr, other.mpfr, MPFR_RNDN);
        return result;
    }
    
    // 平方根
    Gmpfr sqrt() const {
        CGAL_precondition(*this >= 0);
        Gmpfr result;
        mpfr_sqrt(result.mpfr, mpfr, MPFR_RNDN);
        return result;
    }
    
    // 获取精度
    mpfr_prec_t get_precision() const {
        return mpfr_get_prec(mpfr);
    }
    
    // 设置精度(可能丢失信息)
    void set_precision(mpfr_prec_t prec) {
        mpfr_prec_round(mpfr, prec, MPFR_RNDN);
    }
};

数类型转换系统

头文件位置CGAL/number_utils.h

// 通用数类型转换
template < typename From, typename To >
struct Number_type_converter {
    To operator()(const From& f) const {
        return To(f);
    }
};
 
// 特化:Gmpz到double
template <>
struct Number_type_converter<CGAL::Gmpz, double> {
    double operator()(const CGAL::Gmpz& z) const {
        return z.to_double();
    }
};
 
// 特化:double到Gmpq
template <>
struct Number_type_converter<double, CGAL::Gmpq> {
    CGAL::Gmpq operator()(double d) const {
        // 将double转换为有理数
        return CGAL::Gmpq(d);
    }
};
 
// 区间转换
template < typename NT >
std::pair<double, double> to_interval(const NT& x) {
    return Real_embeddable_traits<NT>::To_interval()(x);
}
 
// 符号计算
template < typename NT >
CGAL::Sign sign(const NT& x) {
    return Real_embeddable_traits<NT>::Sgn()(x);
}
 
// 比较
template < typename NT >
CGAL::Comparison_result compare(const NT& a, const NT& b) {
    return Real_embeddable_traits<NT>::Compare()(a, b);
}

使用示例

示例1:基本数类型使用

#include <CGAL/Gmpz.h>
#include <CGAL/Gmpq.h>
#include <CGAL/Gmpfr.h>
#include <iostream>
 
int main() {
    // Gmpz - 任意精度整数
    CGAL::Gmpz a("123456789012345678901234567890");
    CGAL::Gmpz b("987654321098765432109876543210");
    
    std::cout << "Gmpz operations:" << std::endl;
    std::cout << "  a = " << a << std::endl;
    std::cout << "  b = " << b << std::endl;
    std::cout << "  a + b = " << (a + b) << std::endl;
    std::cout << "  a * b = " << (a * b) << std::endl;
    
    // Gmpq - 任意精度有理数
    CGAL::Gmpq p(1, 3);  // 1/3
    CGAL::Gmpq q(2, 6);  // 2/6 = 1/3(自动约分)
    
    std::cout << "\nGmpq operations:" << std::endl;
    std::cout << "  p = " << p << std::endl;
    std::cout << "  q = " << q << std::endl;
    std::cout << "  p == q: " << (p == q) << std::endl;
    std::cout << "  p + 1/6 = " << (p + CGAL::Gmpq(1, 6)) << std::endl;
    
    // Gmpfr - 多精度浮点数
    CGAL::Gmpfr::set_default_precision(256);  // 256位精度
    
    CGAL::Gmpfr x(1.0);
    CGAL::Gmpfr y = x / 3;
    
    std::cout << "\nGmpfr operations (256 bits):" << std::endl;
    std::cout << "  1/3 = " << y << std::endl;
    std::cout << "  sqrt(2) = " << CGAL::Gmpfr(2).sqrt() << std::endl;
    
    return 0;
}

示例2:代数结构特征的使用

#include <CGAL/Gmpz.h>
#include <CGAL/Gmpq.h>
#include <CGAL/Algebraic_structure_traits.h>
#include <iostream>
 
// 检查数类型的代数结构
template < typename NT >
void check_algebraic_structure() {
    typedef CGAL::Algebraic_structure_traits<NT> AST;
    
    std::cout << "Type: " << typeid(NT).name() << std::endl;
    std::cout << "  Is exact: " << AST::Is_exact::value << std::endl;
    std::cout << "  Is numerical sensitive: " 
              << AST::Is_numerical_sensitive::value << std::endl;
    
    // 检查代数类别
    typedef typename AST::Algebraic_category Category;
    
    if (std::is_same<Category, CGAL::IntegralDomainWithoutDivision_tag>::value)
        std::cout << "  Category: IntegralDomainWithoutDivision" << std::endl;
    else if (std::is_same<Category, CGAL::IntegralDomain_tag>::value)
        std::cout << "  Category: IntegralDomain" << std::endl;
    else if (std::is_same<Category, CGAL::EuclideanRing_tag>::value)
        std::cout << "  Category: EuclideanRing" << std::endl;
    else if (std::is_same<Category, CGAL::Field_tag>::value)
        std::cout << "  Category: Field" << std::endl;
    else if (std::is_same<Category, CGAL::FieldWithSqrt_tag>::value)
        std::cout << "  Category: FieldWithSqrt" << std::endl;
}
 
// 使用GCD运算(仅对欧几里得环有效)
template < typename NT >
NT compute_gcd(NT a, NT b) {
    typedef CGAL::Algebraic_structure_traits<NT> AST;
    typedef typename AST::Algebraic_category Category;
    
    // 编译时检查是否是欧几里得环
    static_assert(
        std::is_base_of<CGAL::EuclideanRing_tag, Category>::value,
        "GCD only available for Euclidean rings"
    );
    
    typename AST::Gcd gcd;
    return gcd(a, b);
}
 
int main() {
    std::cout << "=== Gmpz (Euclidean Ring) ===" << std::endl;
    check_algebraic_structure<CGAL::Gmpz>();
    
    CGAL::Gmpz a(48), b(18);
    std::cout << "  gcd(48, 18) = " << compute_gcd(a, b) << std::endl;
    
    std::cout << "\n=== Gmpq (Field) ===" << std::endl;
    check_algebraic_structure<CGAL::Gmpq>();
    
    // 以下会导致编译错误(Gmpq是域,没有GCD)
    // CGAL::Gmpq p(1, 2), q(1, 3);
    // compute_gcd(p, q);  // 编译错误!
    
    return 0;
}

示例3:自定义数类型

#include <CGAL/Algebraic_structure_traits.h>
#include <CGAL/Real_embeddable_traits.h>
#include <iostream>
 
// 自定义数类型:模p整数
class ModularInteger {
    int value_;
    static int modulus_;
    
    void normalize() {
        value_ %= modulus_;
        if (value_ < 0) value_ += modulus_;
    }
    
public:
    ModularInteger(int v = 0) : value_(v) { normalize(); }
    
    static void set_modulus(int p) { modulus_ = p; }
    static int modulus() { return modulus_; }
    
    // 算术运算(模p)
    ModularInteger operator+(const ModularInteger& other) const {
        return ModularInteger(value_ + other.value_);
    }
    
    ModularInteger operator-(const ModularInteger& other) const {
        return ModularInteger(value_ - other.value_);
    }
    
    ModularInteger operator*(const ModularInteger& other) const {
        return ModularInteger(value_ * other.value_);
    }
    
    ModularInteger inverse() const {
        // 扩展欧几里得算法求逆元
        int a = value_, m = modulus_;
        int x = 1, y = 0;
        while (a > 1) {
            int q = a / m;
            int t = m;
            m = a % m;
            a = t;
            t = y;
            y = x - q * y;
            x = t;
        }
        return ModularInteger(x);
    }
    
    ModularInteger operator/(const ModularInteger& other) const {
        return *this * other.inverse();
    }
    
    // 比较
    bool operator==(const ModularInteger& other) const {
        return value_ == other.value_;
    }
    
    bool operator!=(const ModularInteger& other) const {
        return value_ != other.value_;
    }
    
    int value() const { return value_; }
};
 
int ModularInteger::modulus_ = 17;  // 默认模数
 
// 特化代数结构特征
namespace CGAL {
template <>
struct Algebraic_structure_traits<ModularInteger> {
    typedef ModularInteger Type;
    typedef Tag_true Is_exact;
    typedef Tag_false Is_numerical_sensitive;
    typedef Field_tag Algebraic_category;  // 模素数时是域
};
 
// ModularInteger 不是 Real_embeddable(无法嵌入实数轴)
// 因此不需要特化 Real_embeddable_traits
}
 
int main() {
    ModularInteger::set_modulus(17);
    
    std::cout << "Modular arithmetic (mod 17):" << std::endl;
    
    ModularInteger a(5);
    ModularInteger b(7);
    
    std::cout << "  a = " << a.value() << std::endl;
    std::cout << "  b = " << b.value() << std::endl;
    std::cout << "  a + b = " << (a + b).value() << std::endl;
    std::cout << "  a * b = " << (a * b).value() << std::endl;
    std::cout << "  a / b = " << (a / b).value() << std::endl;
    
    // 验证:5/7 ≡ 5 * 7^(-1) (mod 17)
    // 7^(-1) ≡ 5 (mod 17) 因为 7*5 = 35 ≡ 1 (mod 17)
    // 所以 5/7 ≡ 5*5 = 25 ≡ 8 (mod 17)
    std::cout << "  Verification: 7 * 5 = " << (b * ModularInteger(5)).value() << " (mod 17)" << std::endl;
    
    return 0;
}

示例4:Quotient 与数域扩展

#include <CGAL/Gmpz.h>
#include <CGAL/Quotient.h>
#include <CGAL/number_utils.h>
#include <iostream>
 
// 使用Quotient从Gmpz构造有理数
typedef CGAL::Quotient<CGAL::Gmpz> Rational;
 
// 计算连分数展开
template < typename Rational >
void continued_fraction_expansion(Rational x, int max_terms = 10) {
    std::cout << "Continued fraction of " << x << ": [";
    
    for (int i = 0; i < max_terms && x != 0; ++i) {
        // 整数部分
        typename Rational::Type a = x.numerator() / x.denominator();
        std::cout << a;
        
        // 分数部分
        x = x - Rational(a);
        
        if (x != 0) {
            std::cout << ", ";
            x = Rational(x.denominator(), x.numerator());  // 取倒数
        }
    }
    
    std::cout << "]" << std::endl;
}
 
// 高精度计算:e的近似
typedef CGAL::Quotient<CGAL::Gmpz> ExactRational;
 
ExactRational approximate_e(int n_terms) {
    ExactRational e(0);
    ExactRational term(1);  // 1/0! = 1
    
    for (int i = 0; i < n_terms; ++i) {
        e = e + term;
        term = term / ExactRational(i + 1);  // 除以 (i+1) 得到下一项
    }
    
    return e;
}
 
int main() {
    // 基本Quotient操作
    Rational a(3, 4);  // 3/4
    Rational b(5, 6);  // 5/6
    
    std::cout << "Rational arithmetic:" << std::endl;
    std::cout << "  a = " << a << std::endl;
    std::cout << "  b = " << b << std::endl;
    std::cout << "  a + b = " << (a + b) << std::endl;
    std::cout << "  a * b = " << (a * b) << std::endl;
    std::cout << "  a / b = " << (a / b) << std::endl;
    
    // 连分数展开
    std::cout << "\nContinued fractions:" << std::endl;
    continued_fraction_expansion(Rational(17, 12));
    continued_fraction_expansion(Rational(355, 113));  // π的近似
    
    // 计算e
    std::cout << "\nApproximation of e:" << std::endl;
    for (int n : {5, 10, 15, 20}) {
        ExactRational e = approximate_e(n);
        std::cout << "  n=" << n << ": " << e << " ≈ " << e.to_double() << std::endl;
    }
    
    return 0;
}

复杂度分析

基本运算复杂度

运算GmpzGmpqGmpfr(p位)double
加减O(n)O(n)O(p)O(1)
乘法O(n log n)O(n log n)O(p log p)O(1)
除法O(n log n)O(n log n)O(p log p)O(1)
GCDO(n²)---
平方根O(n log²n)-O(p log p)O(1)

其中 n 是机器字的数量(对于k位整数,n ≈ k/64)。

内存使用

数类型内存开销说明
double8字节固定精度
GmpzO(n)动态数组存储 limbs
GmpqO(n)两个Gmpz(分子+分母)
GmpfrO(p)根据精度动态分配
QuotientO(n)类似Gmpq,可能未约分

精度与性能权衡

应用场景推荐数类型精度相对速度
快速原型double~15位1x
精确谓词Interval_nt可变2-5x
精确有理Gmpq无限10-100x
延迟精确Lazy_exact_nt按需1-10x
高精度浮点Gmpfr可配置10-1000x

关键文件位置

文件路径说明
CGAL/Algebraic_structure_traits.h代数结构特征定义
CGAL/Real_embeddable_traits.hReal_embeddable特征定义
CGAL/Gmpz.hGMP整数包装类
CGAL/Gmpq.hGMP有理数包装类
CGAL/Gmpfr.hMPFR浮点数包装类
CGAL/Quotient.h商域构造模板
CGAL/MP_Float.h多精度浮点数
CGAL/Interval_nt.h区间数类型
CGAL/Lazy_exact_nt.h延迟精确数类型
CGAL/number_type_basic.h数类型基础定义
CGAL/number_utils.h数类型工具函数
CGAL/double.hdouble类型的特征特化
CGAL/int.hint类型的特征特化