为何CGAL Lazy数默认除法求和的区间上界异常巨大？

我有一个形如 s = a₁/b₁ + a₂/b₂ + ... 的求和式，其中a_i和b_i都是CGAL Lazy数（类型为CGAL::Lazy_exact_nt<CGAL::Quotient<CGAL::MP_Float>>）。调用s.interval()返回的区间范围极大：下界接近预期值（也与CGAL::as_double(s)的结果接近），但上界异常巨大，甚至是无穷大。

不过用自定义的mydivision(a_i, b_i)替代默认的a_i/b_i后，区间就会变得紧凑，符合预期。自定义除法的实现如下：

typedef CGAL::Quotient<CGAL::MP_Float> Quotient;
typedef CGAL::Lazy_exact_nt<Quotient>  lazyScalar;

lazyScalar mydivision(lazyScalar x1, lazyScalar x2) {
  Quotient q1 = x1.exact();
  Quotient q2 = x2.exact();
  lazyScalar q = lazyScalar(Quotient(
    q1.numerator() * q2.denominator(), q1.denominator() * q2.numerator())
  );
  return q;
}

这个求和式的收敛值是π/2 - 1，对应的公式为：
$$\sum_{k=1}^\infty \frac{k!}{(3)(5)(7)\cdots(2k+1)} = \frac{\pi}{2} - 1$$

完整测试代码如下：

#include <vector>
#include <CGAL/number_utils.h>
#include <CGAL/Lazy_exact_nt.h>
#include <CGAL/MP_Float.h>
#include <CGAL/Quotient.h>
#include <CGAL/Interval_nt.h>

typedef CGAL::Quotient<CGAL::MP_Float> Quotient;
typedef CGAL::Lazy_exact_nt<Quotient>  lazyScalar;
typedef std::vector<lazyScalar>        lazyVector;

// 自定义除法
lazyScalar mydivision(lazyScalar x1, lazyScalar x2) {
  Quotient q1 = x1.exact();
  Quotient q2 = x2.exact();
  lazyScalar q = lazyScalar(Quotient(
    q1.numerator() * q2.denominator(), q1.denominator() * q2.numerator())
  );
  return q;
}

// 对lazy数向量求和
lazyScalar lazySum(lazyVector lv) {
  const size_t n = lv.size();
  lazyScalar sum(0);
  for(size_t i = 0; i < n; i++) {
    sum += lv[i];
  }
  return sum;
}

// 向量元素默认除法
lazyVector lv1_dividedby_lv2(lazyVector lv1, lazyVector lv2) {
  const size_t n = lv1.size();
  lazyVector lv(n);
  for(size_t i = 0; i < n; i++) {
    lv[i] = lv1[i] / lv2[i];
  }
  return lv;
}

// 向量元素自定义除法
lazyVector lv1_mydividedby_lv2(lazyVector lv1, lazyVector lv2) {
  const size_t n = lv1.size();
  lazyVector lv(n);
  for(size_t i = 0; i < n; i++) {
    lv[i] = mydivision(lv1[i], lv2[i]);
  }
  return lv;
}

// 计算lazy数向量的累积乘积
lazyVector lazyCumprod(lazyVector lvin) {
  const size_t n = lvin.size();
  lazyVector lv(n);
  lazyScalar prod(1);
  for(size_t i = 0; i < n; i++) {
    prod *= lvin[i];
    lv[i] = prod;
  }
  return lv;
}

// 使用默认除法计算求和式
lazyScalar Euler(int n) {
  lazyVector lv1(n);
  for(int i = 0; i < n; i++) {
    lv1[i] = lazyScalar(i + 1);
  }
  lazyVector lv2(n);
  for(int i = 0; i < n; i++) {
    lv2[i] = lazyScalar(2*i + 3);
  }
  return lazySum(lv1_dividedby_lv2(lazyCumprod(lv1), lazyCumprod(lv2))); 
}

// 使用自定义除法计算求和式
lazyScalar myEuler(int n) {
  lazyVector lv1(n);
  for(int i = 0; i < n; i++) {
    lv1[i] = lazyScalar(i + 1);
  }
  lazyVector lv2(n);
  for(int i = 0; i < n; i++) {
    lv2[i] = lazyScalar(2*i + 3);
  }
  return lazySum(lv1_mydividedby_lv2(lazyCumprod(lv1), lazyCumprod(lv2))); 
}

// 测试：当n≥170时出现差异
int main() {
  lazyScalar euler = Euler(171);
  CGAL::Interval_nt<false> interval = euler.approx();
  std::cout << "lower bound: " << interval.inf() << "\n"; // 输出约0.57
  std::cout << "upper bound: " << interval.sup() << "\n"; // 输出inf
  lazyScalar myeuler = myEuler(171);
  CGAL::Interval_nt<false> myinterval = myeuler.approx();
  std::cout << "lower bound: " << myinterval.inf() << "\n"; // 输出约0.57
  std::cout << "upper bound: " << myinterval.sup() << "\n"; // 输出约0.57
  return 0;
}

问题：为何使用默认除法时，求和结果的区间上界会异常巨大？ #

1. CGAL Lazy数默认除法的延迟计算特性 #

CGAL::Lazy_exact_nt的核心设计是延迟精确计算，默认的除法运算符并不会立即将a/b转换为精确的商，而是会把这个除法操作作为一个表达式节点存储起来，直到需要精确值或近似值时才会触发计算。

当调用interval()或approx()时，CGAL需要为这个延迟的除法表达式估算区间范围。对于你的场景，随着n增大，累积乘积的数值会变得极大，此时分母的近似值可能因为浮点数精度限制或溢出，导致其区间下界被估算为一个极小的正数——这会直接让1/分母的区间上界爆炸到无穷大，最终整个求和结果的区间上界被这种极端值带偏。

2. 自定义除法的精确计算优势 #

你的mydivision函数直接调用exact()方法，立即将lazy数转换为精确的Quotient对象，然后通过分子分母交叉相乘的方式构造出精确的商，再重新包装为lazyScalar。这相当于提前完成了除法的精确计算，把a/b转换成了一个单一的精确值节点，而非延迟的除法表达式。

当后续计算区间时，这个单一节点的近似值可以直接基于精确值生成紧凑的区间，完全避免了延迟除法带来的区间膨胀问题。

3. n≥170时才出现问题的原因 #

当n较小时，累积乘积的数值还在普通浮点数的精度范围内，默认除法的近似计算还能维持合理的区间。但当n增大到170左右时：

• 累积乘积的数值已经远超浮点数的有效精度范围，默认除法的延迟近似计算中，分母的区间下界会被严重低估，导致1/分母的上界无限放大。
• 自定义除法通过精确计算，完全规避了近似误差的累积和放大，始终基于精确的商来生成区间，因此不会出现异常。

内容的提问来源于Stack Exchange，提问作者Stéphane Laurent