65.9K
CodeProject 正在变化。 阅读更多。
Home

区间算术的实践实现

Henrik Vestermark
starIconstarIconstarIconstarIconstarIcon

5.00/5 (5投票s)

2024年3月16日

CPOL

25分钟阅读

viewsIcon

5841

在 C++ 中构建区间算术类的步骤

目录

引言

您可以找到C++编程语言中区间算术模板类的示例，但没有一个易于理解和实现。这就是为什么本文名为“区间算术的实际实现”。本文重点介绍如何实现一个支持C++编程语言中`float`和`double`类型的通用区间模板类。

区间算术的主要优点是它能够提供结果的保证边界，这对于处理舍入误差、测量不确定性或不精确输入特别有用。当执行计算时，结果也是包含所有可能结果的区间。这在数值分析、计算机科学和工程等领域尤为重要，因为在这些领域中，了解计算的精度和可靠性至关重要。它允许在复杂计算中进行更 robust 的错误处理和不确定性管理。

构建一个能够处理IEEE754浮点数所有区间算术的区间模板软件类是简单的数学。本文描述了构建区间模板类的实际实现方面，并讨论了实现背后的技巧和数学。该实现涵盖了从+、-、*、/等基本运算符到`sqrt()`、`log()`、`pow()`、`exp()`等初等函数以及`sin()`、`cos()`、`tan()`、`asin()`、`acos()`和`atan()`等三角函数，以及`sinh()`、`cosh()`、`tanh()`、`asinh()`、`acosh()`、`atanh()`等双曲函数。当然，它还支持π、ln2和ln10等区间常数。区间算术支持`float`和`double`等基本的C++类型。

像往常一样，我们关注模板类的实际实现方面，而不是其背后的理论原理。有关区间算术和IEEE754浮点格式使用的一些更深入的见解，请参阅[4]区间算术：从原理到实现和[5]。

本文将仅限于第一阶段的实现。完整的源代码和文档可在我的网站上找到：区间算术 (hvks.com)。

更改日志

  • 2024年3月13日。这是2015年原始版本的完全重写的第二版。

区间算术的历史

区间算术的历史可以追溯到20世纪中叶，并与计算机科学和数值分析的发展交织在一起。以下是其演变的简要概述：

20世纪50年代左右的早期概念。区间算术的根源可以追溯到阿基米德等数学家的工作。然而，其正式发展始于20世纪50年代数字计算机的出现。

区间算术作为一个形式化领域的发展，主要归功于拉蒙·E·摩尔。摩尔在约翰霍普金斯大学应用物理实验室工作，后来在斯坦福大学工作，他开发了区间算术的基本规则和操作。他的工作旨在解决计算中的误差界限问题，这对于可靠的科学计算至关重要。

区间算术在20世纪70年代至80年代的扩展和形式化。在此期间，该领域显著扩展。研究人员开始形式化算术并探索其理论性质。他们还开始将其应用于科学计算、工程和误差分析中的各种问题。

随着计算机科学的兴起，区间算术在20世纪80年代至90年代开始被集成到计算机算法和软件中。这得益于工程和科学领域对可靠数值计算日益增长的需求。

2000年代以后，区间算术持续发展并找到了新的应用。它被应用于从数值分析到机器人和控制系统等领域。开发专门的软件和硬件以高效执行区间算术也一直是重点。

在其历史中，区间算术一直由对精确可靠数值计算的需求所驱动，特别是在面对不确定数据或舍入误差时。它的发展反映了数学和计算机科学中向更 robust 和抗误差方法发展的更广泛趋势。作为一个有趣的例子：

阿基米德是已知最早使用某种形式的区间算术来估计π值的人物之一。他通过在一个圆形内部和外部绘制多边形并计算它们的周长来做到这一点。通过增加多边形的边数，他能够缩小边界并更接近地近似π。当时，阿基米德无法计算π的精确值。他确定π在（大约3.1429）和（大约3.1408）之间。这种近似方法在当时是相当复杂的，是区间算术原理的早期历史例子，其中值被限制在一个上限和下限之间，以提供一个可能的取值范围。

区间算术的应用

区间算术由于其处理不确定性和为计算提供保证边界的能力，在各个领域都有实际应用。以下是一些主要应用，首先是数值分析，它通过考虑舍入误差和其他数值不确定性来确保计算的准确性和可靠性。
在计算机图形渲染和几何计算中，区间算术有助于处理不精确的数据和计算，从而产生更 robust 的算法。对于需要稳定性和可靠性的系统（控制系统），区间算术可以管理测量和模型参数中的不确定性。在机器人规划和控制中，它用于考虑传感器数据的不精确性，并确保安全的移动和决策。在工程设计中，工程师使用区间算术设计具有材料特性公差和变化的系统。对于科学计算中的模拟和建模，如果输入数据可能不确定或在一定范围内变化，区间算术提供了一种方法来确保结果考虑到所有可能的输入值。这还扩展到优化问题、误差分析和传播，这有助于理解输入数据中的误差如何影响计算结果，这在航空航天和医疗设备等敏感应用中至关重要。

通过提供一个处理和传播不确定性的框架，区间算术增强了这些及其他领域计算任务的鲁棒性和可靠性。

区间算术

为了限制浮点算术中的舍入误差，我们可以使用区间算术来跟踪舍入误差。经过一系列使用+、-、*和/等基本运算符的运算后，我们得到的是一个区间而不是结果的近似值。区间宽度表示结果的不确定性，但我们确切地知道正确的结果将在此区间内。

在区间算术中，区间通常表示为 `[a, b]`，其中 '`a`' 是下界，'`b`' 是上界。通常，假定 `a ≤ b`，定义了所谓的正规区间。相反，如果 `a > b`，则区间称为非正规区间。根据 IEEE 1788 区间算术标准，正规和非正规区间都是有效的，并且必须在计算中正确处理，因为我们不能总是假定 `a ≤ b`。这使得必须通过取 `min(a, b)` 来确定区间的下界（下确界），并通过取 `max(a, b)` 来确定上界（上确界）。此要求在实现中引入了额外的步骤：代码必须包含查找任何给定区间的最小值和最大值的函数，而不是假定 `a ≤ b`。

此外，当 `a = b` 时，该区间称为单点区间。

为了简化公式和解释，我们经常将区间表示为正规区间（`a ≤ b`）。然而，在实践中，特别是在编码实现中，根据区间的正规或非正规性质（`a > b`）正确处理区间至关重要。

此外，定义闭区间`[a, b]`、半闭（或半开）区间（例如，`[a, b)` 或 `(a, b]`）和开区间`(a, b)`也很重要。每种类型对于区间中端点值的包含或排除都有不同的含义。

五种不同的区间类型是

区间类型 定义 符号 注释
Close a≤x≤b [a,b]  
左开 a<x≤b (a,b] 与右闭相同
右开 a≤x<b [a,b) 与左闭相同
打开 a<x<b (a,b)  
Empty   ∅ 空区间

除了∅（空区间）之外，其他4种不同类型并非标准严格要求。然而，拥有它们非常方便，否则定义开放或半开放区间将面临挑战。例如，在C++中创建开放区间`(0,1)`将是这样的：

interval<double> i(nextafter(0,+INFINITY), nextafter(1,-INFINITY));

创建这个区间的方法相当不便和繁琐，而我的实现可以这样做：

interval<double> i(0,1,OPEN);

这种方法确保了对区间算术的全面理解和正确实现，遵循IEEE 1788标准。

区间由两个数字表示，分别代表区间的下限和上限

[a,b] 其中 (对于正规区间)

正规区间的4个基本运算符+,-,*,和/定义如下：

[a,b] + [c,d] = [a+c,b+d]
[a,b] – [c,d] = [a-d, b-c]
[a,b] * [c,d] = [min(ac,ad,bc,bd),max(ac,ad,bc,bd)]
[a,b] / [c,d] = [min( and [c,d] does not contain 0.

除法也可以改写为

[a,b] / [c,d] = [a,b] * (1/[c,d]) = [a,b] * [1/d,1/c]

我们不喜欢的是区间乘法需要额外的工作。根据定义，一次区间乘法需要8次乘法。然而，我们可以用少得多的次数来完成。

我们首先定义一个区间类如下：

类 [a,b] 条件 a≤b
Zero a=0 ^ b=0
混合 a<0 ^ b≥0
正数 a>0
负片 b<0

根据下表，在大多数情况下，我们可以用2次乘法完成，在最坏情况下用4次乘法，而不是8次乘法。

类 [a,b] 类 [c,d] 乘法
正数 正数 [a*c,b*d]
正数 混合 [b*c,b*d]
正数 负片 [b*c,a*d]
混合 正数 [a*d,b*d]
混合 混合 [min(a*d,b*c),max(a*c,b*d)]
混合 负片 [b*c,a*c]
负片 正数 [a*d,b*c]
负片 混合 [a*d,a*c]
负片 负片 [b*d,a*c]
Zero - [0,0]
- Zero [0,0]

从实际角度来看，我们不能直接使用上述定义。在使用有限表示（如大多数编程语言和大多数供应商（如Intel、AMS等）的CPU中常见的IEEE754标准用于32位二进制浮点算术和64位二进制浮点算术）执行任何基本操作时，我们总是会引入一些舍入误差。

一般来说，当执行任何基本操作（如+,-,*,/）时，我们都会有一个误差，该误差受½β1-t的限制，使用IEEE754中的默认舍入模式或向正无穷大或负无穷大舍入时。在IEEE754二进制算术中，32位浮点数β=2且t=24（23+1），64位浮点数t=53（52+1），这使得执行任何基本操作时都会产生相对误差。

IEEE754 舍入 朝向+∞或-∞
32位浮点数 ~5.96-8 ~1.19-7
64位浮点数 ~2.22-16 ~1.11-16

在执行区间操作时，我们还需要控制舍入模式，以确保结果确实在区间内。 表示向-∞舍入， 表示向+∞舍入。

根据IEEE1788标准，区间也可以是空集。它用符号∅表示。对于空区间有一些规则。

[a,b]+∅=[a,b], same goes for ∅+[a,b]:=[a,b]
[a,b]-∅:=[a,b], ∅-[a,b]:=[-b,-a]
[a,b]*∅:=∅, ∅*[a,b]:=∅
[a,b]/∅:=∅, ∅/[a,b]:=∅

最后，我们还需要注意区间平方（例如[a,b]2）时会发生什么。如果您取区间`[-1,1]`并使用乘法公式将其与自身相乘，您将得到一个不必要的大区间`[-1,1]`。如果您查看x2在`[-1,1]`范围内的图，您将得到：

正如您所看到的，结果从+1下降到0，然后再次上升到+1。这种操作的正确结果是[0,1]。IEEE 1788标准预见到了这个问题，并要求符合标准的实现具有一个名为`sqr()`的函数用于区间的平方。这个问题对于区间的偶次幂仍然存在，需要在实现中正确处理。

如果你看x3的图，你会得到：

如您所见，操作的值贯穿`[-1,1]`的区间。

控制舍入

IEEE754浮点标准中的默认舍入模式是就近舍入（简称RN）。然而，如前一节所示，我们需要访问向下舍入到-∞（RD）和向上舍入到+∞（RU）。

在区间算术中，舍入是一个关键方面，它可以使用基于硬件的方法或基于软件的仿真来实现。基于硬件的舍入提供速度和准确性，但缺乏灵活性和可移植性。相比之下，软件仿真更具可移植性和灵活性，但可能更慢且准确性可能较低。

基于硬件的舍入

以下是使用基于硬件的舍入的一些优点。
基于硬件的舍入通常更快，因为它利用了处理器的内置功能。由于CPU中的直接实现，它通常提供更精确和一致的舍入。减少了软件仿真的开销，使计算更高效。

然而，速度上的优势是有代价的，因为它依赖于特定的硬件特性，这些特性可能在不同的平台或处理器之间不可用或不一致。

硬件实现比软件更难定制或更新。依赖特定硬件功能的软件可能在不同系统之间面临可移植性问题。

根据IEEE754的要求，向-∞（RD）舍入和向+∞（RU）舍入是可用的，并且在这些模式与默认的就近舍入之间切换非常简单。

控制硬件舍入

为了控制舍入过程，IEEE754提供了至少4种舍入方式。

  1. 四舍五入到最近，这是默认的舍入方法。
  2. 向下舍入至-∞
  3. 向上舍入至+∞
  4. 向零舍入（截断）

根据CPU架构的不同，不同的指令可以启用这些不同的模式。然而，遗憾的是，它在架构之间不具有可移植性。我们只依赖于以下三个函数以某种方式可用，它们可以启用这些模式。

Fp_near();       // Set Floating-point mode to round to the nearest

Fp_Down();       // Set floating-point mode to round down towards -∞

Fp_up();         // Set floating-point mode to round towards +∞

作为如何使用的一个例子。我们可以实现两个区间`[a,b]`和`[c,d]`的区间加法，如下所示：

Fp_up();
f=b+d;
Fp_down();
e=a+c;
Fp_near();
return [e,f];

使用硬件控制的舍入似乎是有利的，然而，在[7]中，他们指出在不同的硬件舍入模式之间来回切换是一个相对昂贵的操作，令人惊讶的是，在许多情况下，使用软件舍入方法提供了更高的性能。尽管这高度依赖于实际的CPU架构。

基于软件的舍入（模拟）

基于软件的舍入无疑具有其优势。基于软件的舍入更具可移植性，因为它不依赖于特定的硬件功能。它更容易定制或适应特定需求，例如实现非标准舍入模式，并且它确保了在不同平台和硬件上的一致行为。

基于软件的舍入的缺点是它通常比硬件实现慢，因为存在仿真开销，并且与硬件实现相比，其实现和维护其准确性和一致性可能更复杂，这取决于软件仿真的质量。

基于软件的舍入利用了就近舍入是IEEE754中的默认舍入模式。并且两种算法是有用的。二和算法和二积算法在区间算术中对于保持精确的边界至关重要。这些操作对于解决浮点算术中的舍入误差问题是必不可少的。参见[7]

  1. 二和算法用于在捕获舍入误差的同时添加两个浮点数。二和操作返回两个数字：精确和（尽可能接近浮点表示所允许的）和误差。
  2. 二积算法与二和算法类似，二积算法将两个浮点数相乘，提供乘积和误差。

为了了解我们如何利用这两种算法，我们考虑算术运算`a+b`的和。结果将四舍五入到最近的结果（RN）。就近舍入（RN）要么与向下舍入（RD）相同，要么与向上舍入（RU）相同。为了找出我们使用二和算法的误差。如果误差为负，那么：

RN(a+b)==RU(a+b) and RD(a+b) is predecessor(RN(a+b))

如果e为正，则

RN(a+b)==RD(a+b) and RU(a+b) is successor(RN(a+b))

`successor` 返回下一个可表示的浮点数，朝向 +∞；`predecessor` 返回下一个可表示的浮点数，朝向 -∞。

乘法也可以使用相同的类比。如果误差为负，则：

RN(a*b)==RU(a*b) and RD(a*b) is predecessor(RN(a*b))

如果e为正，则

RN(a*b)==RD(a*b) and RU(a*b) is successor(RN(a*b))

您可以将其扩展到除法、倒数和平方根，请参见下表。

操作 (Operation) 误差<0 误差>0
RN(a+b) RU=RN, RD=前驱 RD=RN, RU=后继
RN(a-b) RU=RN, RD=前驱 RD=RN, RU=后继
RN(a*b) RU=RN, RD=前驱 RD=RN, RU=后继
RN(a/b) RU=RN, RD=前驱 RD=RN, RU=后继
RN(1/b) RU=RN, RD=前驱 RD=RN, RU=后继
RN(sqrt(a)) RU=RN, RD=前驱 RD=RN, RU=后继

Two-sum算法由O. Møller于1965年发明，有两种变体：two-sum和fast two-sum算法，其实现如下。

原始的C++二和算法。

//  Implement the twosum algorithm. Assuming round to nearest mode
//  (default IEEE754)
//  sum=(a+b)
//  a1=sum-b
//  b1=sum-a
//  da=a-a1
//  db=b-b1
//  err=da+db
//  return sum, err
//  if err is negative then sum=Round_up(a+b),
//  and round_down(a+b)=previous(sum)
//  if err is positive then sum=Round_down(a+b),
//  and Round_up(a+b)=successor(sum)
//  The two_sum algorithm requires 6 floating point operations
std::pair<_IT, _IT> two_sum(_IT a, _IT b)
{
       const _IT sum = (a + b);
       const _IT a1 = sum - b;
       const _IT b1 = sum - a;
       const _IT da = a - a1;
       const _IT db = b - b1;
       const _IT err = da + db;

       return std::make_pair(sum, err);
}

该算法大约需要6次浮点运算。

下面是C++中一个更快版本，出人意料地称为快速二和算法。

//  Implement the fast two-sum algorithm Assuming round to nearest mode (default IEEE754)
//  sum=(a+b)
//  a1=sum-a
//  err=b-a1
//  return sum, err
//  if err is negative then sum=Round_up(a+b),
//  and round_down(a+b)=previous(sum)
//  if err is positive then sum=Round_down(a+b),
//  and Round_up(a+b)=successor(sum)
//  The fast two-sum algorithm requires only 3 floating point instructions
//  and a comparison
std::pair<_IT, _IT> fasttwo_sum(_IT a, _IT b)
{
       const _IT sum = a + b;
       _IT err, tmp;
       if (abs(a) > abs(b))
       {
             tmp = sum - a; err = b - tmp;
       }
       else
       {
             tmp = sum - b; err = a - tmp;
       }
       return std::make_pair(sum, err);
}

这只需要3次浮点运算和比较运算。但是，对于`a+b`的和，它要求`a≥b`。

两种乘积由Rump发明，参见下面的C++实现。

// Split argument into a high and low. (Dekker's method)
// double has 53bits in mantissa(52+1) and shifting is therefore (53+1)/2=27
// float has 24bits in mantissa (23+1) and shifting is therefore (23+1)/2=12
_IT split(_IT a)
{
       const _IT C = a*double((1 << 27) + 1);
       _IT xhigh, xlow;
       xhigh = C - (C - a);
       xlow = a - xhigh;
       return std::make_pair(xhigh, xlow);
}

// The standard two-product algorithm (by RUMP's method)
// (xh,xl)=split(a)
// (yh,yl)=split(b)
// p=a*b
// t1=-p+xh*yh
// t2=t1+xh*yl
// t3=t2+xl*yh
// err=t3+xl*yl
// return (p,err)
//     if err is negative then sum=Round_up(a*b),
//     and round_down(a*b)=previous(sum)
//     if err is positive then sum=Round_down(a*b),
//     and Round_up(a*b)=successor(sum)
// The fast two-product algorithm requires 17 floating-point instructions
// and a comparison
std::pair<_IT, _IT> two_prod(_IT a, _IT b)
{
       const std::pair<_IT, _IT> x= split(a);
       const std::pair<_IT, _IT> y = split(b);
       const _IT p = a * b;
       const _IT t1 = -p + x.first * y.first;
       const _IT t2 = t1 + x.first * y.second;
       const _IT t3 = t2 + x.second * y.first;
       const _IT err = t3 + x.second * y.second;
       if (abs(p) > ldexp(1, -1021))  // avoiding overflow
             err = x.second * y.second - ((((p * 0.5 - (x.first * 0.5) * y.first) * 2)
                   -x.second * y.first) - x.first * y.second);
       else
             err=x.second*y.second-(((p-x.first*y.first)-x.second*y.first)
                 -x.first*y.second);
       return std::make_pair(p,err);
}

然而，正如[7]所指出的，如果您可以访问 fma 操作（Fused-Multiply-Add），您可以将 two-product 代码简化为简单：

使用`fma`的二积

//     The fast two product algorithm
//     p=a*b
//     err=fma(a,b,-p) // fma is required in IEEE754 and either implemented 
//     in hardware or software
//     return (p,err)
//     if err is negative then sum=Round_up(a*b),
//     and round_down(a*b)=previous(sum)
//     if err is positive then sum=Round_down(a*b),
//     and Round_up(a*b)=successor(sum)
// The fast two-product algorithm requires only 2 floating point instructions
// and a comparison
std::pair<_IT, _IT> fasttwo_prod(_IT a, _IT b)
{
       const _IT p = a * b;
       const _IT err = fma(a, b, -p);
       return std::make_pair(p, err);
}

支持区间算术的软件库和系统通常实现这些算法，以确保区间内的所有可能值都得到精确的考虑。一些著名的软件和库包括INTLAB：一个用于区间算术的MATLAB工具箱，提供具有舍入误差边界的可靠计算。Boost区间算术库：Boost C++库的一部分，它支持区间算术操作，并具有各种控制舍入和区间边界处理的策略。

这些工具和库被应用于科学计算、工程和金融等各个领域，只要需要准确的误差跟踪和 robust 的区间计算。它们不仅实现了二和和二积，还实现了一整套适用于区间算术的算术和数学运算。

区间类要求

随着区间算术的定义以及两个关键算法二和算法和二积算法的完成，我们现在可以转向区间类的要求。

在开始C++编码过程之前，准备好需求总是很重要的。我通常将其分为几个阶段，所以我从一个基础实现开始，然后通过定义的阶段进一步扩展。对于每个阶段，您都要对该阶段的功能进行测试。

Phase 1: establishing the base functionality
All the essential operators: +,-,*,/,=,+=,-=,*=,/=.
Constructors for a float, double, or long double to interval type.
Conversion operators from interval to regular float, double, long double, short, int, long.
Support of closed, semi-closed, semi-open, and open intervals.
Support the empty interval set ∅.
cout and cin operator.
Comparison operators ==,!=,<,<=,>,>=.
Essentials methods:
            negate: Negate the interval (same as monadic -).
            toString: return the string representation of an interval number.
            width: return the width of the interval. (standard names this wid).
            center: return the center of the interval (standard names this mid).
            radius: return the radius of the interval (standard names this rad).
            infimum: return the infimum of the interval (standard name it inf).
            supremum: return the supremum of the interval. (standard name it sup).
            leftinterval: return or set the lower range of the interval as is.
            rightinterval: return or set the upper range of the interval as is.
            isEmpty: return the Boolean value of emptiness of the interval.
            isProper: return true if the interval is a proper interval otherwise false.
            isImproper: return true if the interval is an improper interval otherwise false.
            isPoint: return true if the interval is a singleton, otherwise false.
            isClass: return the classification of the interval.
           abs: return the absolute value of the interval.
           sqrt: return the square root of the interval.
           sqr: return the square of the interval (x2).

Phase 2: expanding the base functionality
mixed data type for arithmetic operations.
manifest interval constants: LN2, LN10, e, PI
exponential (exp) and logarithm functions. (log,log10)
power x<sup>y</sup>
Union: Union of two intervals
Intersection: intersections of two intervals
Interior: is the interval an interior to another interval
Precedes: Does the interval precede another interval
Subset: is the interval a subset of another interval
Inclusion: is the interval included in another interval
In: is the float or double number within the interval

Phase 3: Completed all the auxiliary functions
Trigonometric interval functions
Hyperbolic interval functions:

Phase 4:
Expand support for arbitrary precision interval type
other useful functions e.g. gamma, beta, etc.

第一阶段实现

一切都关乎建立区间类的基本功能。

定义区间模板类

通过将实现分为几个阶段，我们现在可以开始我们的初步设计。
我们需要为我们的区间功能创建一个模板类。本质上，类似于C++中的标准复数模板库。此外，该类需要保存区间的左侧和右侧以及区间类型，以便能够处理前面讨论的闭区间、半开区间或全开区间。不出所料，我们得到了这个模板类定义以及私有部分和变量。

// Interval class
// Realistically the class Type can be float, or double.
// Any other type is not supported
//
template<class _IT> class interval {
       _IT left;                  // Left interval
       _IT right;                 // Right interval
       enum interval_type type;   // Interval type CLOSE, OPEN, LEFT_OPEN,
                                    RIGHT_OPEN, (RIGHT_CLOSE is synonym for
                                    LEFT_OPEN and LEFT_CLOSE same as RIGHT_OPEN
       public:

   };

我们还定义了`fasttwo_sum()`和`fasttwo_prod()`这两个函数作为类的内部`private`函数，我们得到如下结果：

// Interval class
// Realistically the class Type can be float, or double.
// Any other type is not supported
//
template<class _IT> class interval {
       _IT left;                  // Left interval
       _IT right;                 // Right interval
       enum interval_type type;   // Interval type CLOSE, OPEN, LEFT_OPEN,
                                  // RIGHT_OPEN, (RIGHT_CLOSE is synonym for
                                  // LEFT_OPEN and LEFT_CLOSE same as RIGHT_OPEN
                                  // and EMPTY

       // Split argument into a right and left. (Dekker's method)
       // double has 53bits (52+1) in mantissa and shifting is
       // therefore (53+1)/2=27
       // float has 24bits (23+1) in mantissa and shifting is
       // therefore (24+1)/2=12
       _IT split(_IT a)
       {
             const _IT C = a * double((1 << 27) + 1);
             _IT xright, xleft;
             xright = C - (C - a);
             xleft = a - xright;
             return std::make_pair(xright, xleft);
       }

       //     Implement the fast two-sum algorithm Assuming round to nearest
       //     mode (default IEEE754)
       //     sum=(a+b)
       //     a1=sum-a
       //     err=b-a1
       //     return sum, err
       //     if err is negative then sum=Round_up(a+b), and
       //     round_down(a+b)=previous(sum)
       //     if err is positive then sum=Round_down(a+b), and
       //     Round_up(a+b)=succesor(sum)
       //     The fast two-sum algorithm requires only 3 floating point
       //     instructions and a comparison
       static std::pair<_IT, _IT> fasttwo_sum(_IT a, _IT b)
       {
             const _IT sum = a + b;
             _IT err, tmp;
             if (abs(a) > abs(b))
             {
                    tmp = sum - a; err = b - tmp;
             }
             else
             {
                    tmp = sum - b; err = a - tmp;
             }
             return std::make_pair(sum, err);
       }

       //     The fast two-product algorithm
       //     p=a*b
       //     err=fma(a,b,-p)     fma is required in IEEE754 and either
       //     implemented in hardware or software
       //     return (p,err)
       //     if err is negative then sum=Round_up(a*b), and
       //     round_down(a*b)=previous(sum)
       //     if err is positive then sum=Round_down(a*b), and
       //     Round_up(a*b)=succesor(sum)
       //     The fast two-product algorithm requires only 2 floating-point
       //     instructions and a comparison
       std::pair<_IT, _IT> fasttwo_prod(_IT a, _IT b)
       {
             const _IT p = a * b;
             const _IT err = fma(a, b, -p);
             return std::make_pair(p, err);
       }

public:

   };

那是我们模板类的`private`部分。现在来看类的构造函数。

区间构造函数

根据我们的要求，我们的区间类可以声明为零个、一个、两个或三个参数。例如：

// Declare an interval variable with an empty class ∅
Interval<double> i;

// Declare an interval variable initialized with a point
Interval<float> fi(2.0)

// Declare an interval variable as the interval [1,2]
Interval<double> di(1.0,2.0);

// Declare an interval variable with the open interval (3,4)
Interval<double> dit(3.0,4.0,OPEN);

// Declare an improper closed interval variable [3,2]
Interval<float> dit2(3,2);
// constructor. zero, one or two arguments for type _IT
interval();                    // Empty interval
interval(const _IT& d);    // Singleton interval
// Regular interval with interval_type (default CLOSE)
interval(const _IT& l, const _IT& h, const enum interval_type t=CLOSE);

区间类型定义如下：

// The five different interval types
// # Close   a<=x<=b        [a,b]
// # Left open a<x<=b      (a,b]  same as Right close
// # Right open a<=x<b     [a,b)  same as Left close
// # Open a<x<b            (a,b)
// # EMPTY interval        ∅
enum interval_type { EMPTY, CLOSE, LEFT_OPEN, RIGHT_OPEN, OPEN,
                     RIGHT_CLOSE=LEFT_OPEN, LEFT_CLOSE=RIGHT_OPEN };

区间坐标函数

我们还需要一些协调函数。例如，获取区间的左值或右值，或更改它们。可能我们需要找到区间的宽度（`wid`）、区间的半径（`rad`）以及将区间转换为标量元素的中心点（`mid`），我们还需要定义所需的下确界（`inf`）和上确界（`sup`）以及IEEE 1788标准要求的“`is`”方法函数。

我们还需要一些直接的方法来 `get` 和 `set` 区间的左、右值，以及获取或设置区间类型（`[]`、`(]`、`[)`、`()`），我使用 `leftinterval()` 和 `rightinterval()` 来表示区间端点。

// Coordinate functions.
_IT rightinterval() const;                   // Return rightinterval bound
_IT leftinterval() const;                    // Return leftinterval bound
_IT rightinterval(const _IT&);               // Set and return rightinterval bound
_IT leftinterval(const _IT&);                // Set and return leftinterval bound
enum interval_type intervaltype() const;     // Return interval type
enum interval_type intervaltype(enum interval_type); // Set and return interval type

// IEEE1788 standard functions
_IT inf() const;    // Return infimum of interval
_IT sup() const;    // Return supremum of interval
_IT mid() const;    // Return midpoint of interval
_IT rad() const;    // Return radius of interval
_IT wid() const;    // Return width of interval
_IT mig() const;    // Return Mignitude. inf(|x|)
_IT mag() const;    // Return Magnitude. sup(|x|)

// Is methods as required per IEEE 1788 standard
bool isEmpty() const;
bool isPoint() const;
bool isImproper() const;
bool isProper() const;

此前，我引入了区间类的概念，作为对给定区间类型进行分类的方法。例如，零、正、负、混合等。这也可以是一个有用的支持函数。以及`toString()`函数，它生成区间的字符串表示。

// Miscellaneous but useful coordinate functions
enum int_class isClass() const;
std::string toString() const;   // Convert interval to a string

枚举类型`int_class`定义如下：

// The eight different interval classification
// # ZERO           a=0 && b=0
// # POSITIVE0      a==0 && b>0
// # POSITIVE1      a>0 && b>0
// # POSITIVE       a>=0 && b>0
// # NEGATIVE0      a<0 && b==0
// # NEGATIVE1      a<0 && b<0
// # NEGATIVE       a<0 && b<=0
// # MIXED          a<0 && b>0
enum int_class { NO_CLASS, ZERO, POSITIVE0, POSITIVE1, POSITIVE, NEGATIVE0,
                NEGATIVE1, NEGATIVE, MIXED };

区间转换运算符

所有这些都是关于定义转换运算符，将区间转换为C++中的任何内置类型。

// Conversion Operators
operator short() const                   // Conversion to short
operator int() const                     // Conversion to int
operator long() const                    // Conversion to long
operator long long() const               // Conversion to long
operator unsigned short() const          // Conversion to unsigned short
operator unsigned int() const            // Conversion to unsigned int
operator unsigned long() const           // Conversion to unsigned long
operator unsigned long long() const      // Conversion to unsigned long
operator double() const;                 // Conversion to double
operator float() const;                  // Conversion to float

区间基本运算符

这些都是C++赋值运算符，见下文。

// Essential operators
interval<_IT>& operator= ( const interval<_IT>& );
interval<_IT>& operator+=( const interval<_IT>& );
interval<_IT>& operator-=( const interval<_IT>& );
interval<_IT>& operator*=( const interval<_IT>& );
interval<_IT>& operator/=( const interval<_IT>& );

区间非基本运算符

这些是+,-,*,/以及+,-的单目版本，加上布尔==,!=,>,<,>=,<=运算符。

// Arithmetic +,-,*,/ Binary and Unary
template<class _IT> interval<_IT> operator+( const interval<_IT>&, const interval<_IT>& );
template<class _IT> interval<_IT> operator+( const interval<_IT>& );   // Unary
template<class _IT> interval<_IT> operator-( const interval<_IT>&, const interval<_IT>& );
template<class _IT> interval<_IT> operator-( const interval<_IT>& );   // Unary
template<class _IT> interval<_IT> operator*(const interval<_IT>&, const interval<_IT>&);
template<class _IT> interval<_IT> operator/( const interval<_IT>&, const interval<_IT>& );

有了类声明和外部声明，我们现在可以着手实现它。

实现区间类

大部分实现都很简单，并遵循类定义。

  1. 构造函数
  2. 坐标方法
  3. 转换运算符
  4. 基本运算符

实现构造函数

构造函数很简单，只是用初始值填充区间类。如果没有给定参数，区间将初始化为空（∅）区间。请注意，对于三个参数的构造函数，如果不存在，则第三个参数默认为闭区间。

// EMPTY interval
template<class _IT> inline  interval<_IT>::interval()
{ left = _IT(0); right = _IT(0); type = EMPTY; }

// Singleton interval
template<class _IT> inline  interval<_IT>::interval(const _IT& d)
{ left = _IT(d); right = _IT(d); type = CLOSE; }

// Regular interval with optional interval type argument
template<class _IT> inline  interval<_IT>::interval
        (const _IT& l, const _IT& h, const enum interval_type t)
{ left = l; right = h; type = t; }

实现坐标函数（方法）

这些都是区间类中定义的所有方法。大多数方法都相对容易实现。对于这些方法，我们有：

  • 右区间()
  • 左区间()
  • 区间类型()

这些是简单的函数，只返回区间的左、右或类型。无需计算，它只返回当前值。还有一个重载版本允许您设置这些`private`变量。还有`isClass()`和`toString()`方法。

下一组是所需的IEEE1788标准函数

  • 下确界()
  • 上确界()
  • 中点()
  • 半径()
  • 宽度()
  • 最小值幅度()
  • 最大值幅度()
  • 是否为空()
  • 是否为点()
  • 是否非正规()
  • 是否正规()

`inf()`和`sup()`这两个方法需要一点解释。

`inf()`函数代表“下确界”或“最大下界（GLB）”。它返回区间中可能存在的最低值，有效地表示区间的下界。此函数对于确保涉及区间的任何操作或比较都考虑到区间元素可能取到的最低值至关重要，从而保持数学计算的准确性和可靠性。这意味着无论区间是正规（左≤右）还是非正规（左＞右），它仍然返回区间的下确界。这在区间是闭区间、开区间还是半开区间的情况下也都有效。

`sup()`函数代表“上确界”。它返回区间中可能存在的最高值，有效地表示区间的上界。同样，`sup()`函数确保在操作或比较期间考虑区间内的最高可能值。这对于涵盖区间内所有可能的取值范围以及准确的区间算术计算至关重要。同样，无论区间是正规（左≤右）还是非正规（左＞右），它仍然返回区间的上确界。这在区间是闭区间、开区间还是半开区间的情况下也都有效。

在区间算术中，这些函数用于定义区间（例如，对于区间`x`，定义为`[inf(x), sup(x)]`）以及执行区间之间的算术运算。当在区间之间执行加法、减法、乘法或除法等操作时，`inf()`和`sup()`函数通过计算基于所执行操作的新边界，并考虑到区间类型和区间的正规性，来帮助确定结果区间。

// Coordinate functions.
// Return the right interval "as is"
template<class _IT> inline _IT interval<_IT>::rightinterval() const
{ return right; }

// Return the left interval "as is"
template<class _IT> inline _IT interval<_IT>::leftinterval() const
{ return left; }

// Set the right interval "as is".
// If interval is empty set the interval type to CLOSE
template<class _IT> inline _IT interval<_IT>::rightinterval(const _IT& r)
{
       if (type == EMPTY)
             type = CLOSE;
       right = r;
       return right;
}

// Set the left interval "as is".
// If interval is empty set the interval type to CLOSE
template<class _IT> inline _IT interval<_IT>::leftinterval(const _IT& l)
{
       if (type == EMPTY)
             type = CLOSE;
       left = l;
       return left;
}

// Return the current intervaltype
template<class _IT> inline enum interval_type interval<_IT>::intervaltype() const
{ return type; }

// Set the interval type
// The below table is for a proper interval
//                  CLOSE []     OPEN ()      LEFT_OPEN (]       RIGHT_OPEN [)
//     CLOSE  []     #            -,+          -,#           #,+
//     OPEN   ()     +,-          #            #,-           +,#
//     LEFT_OPEN (]  +,#          #,+          #             +,+
//     RIGHT_OPEN [) #,-          -,#          _,_          #
//
//     For improper intervals we preserve the improperness in the result.
//
template<class _IT> inline enum interval_type interval<_IT>::intervaltype
                    (const enum interval_type to)
{
    interval<_IT> x = *this;
    const _IT infi(INFINITY);

    if (x.type != to)
    {
        if (this->isImproper())
        {    // If improper swap the interval and switch the half open intervals
            std::swap(x.left, x.right);
            if (x.type == LEFT_OPEN)
                x.type = RIGHT_OPEN;
            else
                if (x.type == RIGHT_OPEN)
                    x.type = LEFT_OPEN;
        }

        switch (x.type)
        {
        case CLOSE:
            // if the interval is already CLOSE then do nothing
            if (to == LEFT_OPEN || to == OPEN)
                x.left = nextafter(x.left, -infi);
            if (to == RIGHT_OPEN || to == OPEN)
                x.right = nextafter(x.right, +infi);
            break;
        case OPEN:
            // if the interval is already Open then do nothing
            if (to == RIGHT_OPEN || to == CLOSE)
                x.left = nextafter(x.left, +infi);
            if (to == LEFT_OPEN || to == CLOSE)
                x.right = nextafter(x.right, -infi);
            break;
        case LEFT_OPEN:
            // If the interval is already RIGHT_CLOSE same as LEFT_OPEN then do nothing
            if (to == RIGHT_OPEN || to == CLOSE)
                x.left = nextafter(x.left, +infi);
            if (to == RIGHT_OPEN || to == OPEN)
                x.right = nextafter(x.right, +infi);
            break;
        case RIGHT_OPEN:
            // If the interval is already LEFT_CLOSE then do nothing
            if (to == LEFT_OPEN || to == OPEN)
                x.left = nextafter(x.left, -infi);
            if (to == LEFT_OPEN || to == CLOSE)
                x.right = nextafter(x.right, -infi);
            break;
        }
        x.type = to;
        if (this->isImproper())
        {
            std::swap(x.left, x.right);
        }
        *this = x;
    }
    return x.type;
}

// compute the infimum(greatest lower bound) of an interval, 
// taking into account the type of interval
// (closed, open, left-open, or right-open) and whether the interval is proper
// (left endpoint is less than or equal to right endpoint) or improper
// (left endpoint is greater than the right endpoint).
template<class _IT> inline _IT interval<_IT>::inf() const
{
    // For a closed interval, directly return the minimum of left and right.
    if (type == CLOSE)
        return min(left, right);

    _IT adjustedLeft = left;
    _IT adjustedRight = right;
    const _IT infi(INFINITY);

    // Adjust left boundary for open intervals
    if (type == LEFT_OPEN || type == OPEN)
        adjustedLeft = nextafter(left, (left <= right) ? +infi : -infi);

    // Adjust right boundary for open intervals, taking into account improper intervals
    if (type == RIGHT_OPEN || type == OPEN)
        adjustedRight = nextafter(right, (left <= right) ? -infi : +infi);

    // Return the minimum of the adjusted boundaries
    return std::min(adjustedLeft, adjustedRight);
}

// Optimizing the function for calculating the supremum(least upper bound) 
// of an interval follows
// a similar approach to optimizing the infimum calculation
template<class _IT> inline _IT interval<_IT>::sup() const
{
    // For a closed interval, directly return the maximum of left and right.
    if (type == CLOSE)
        return max(left, right);

    _IT adjustedLeft = left;
    _IT adjustedRight = right;
    const _IT infi(INFINITY);

    // Adjust left boundary for open intervals
    if (type == LEFT_OPEN || type == OPEN)
        adjustedLeft = nextafter(left, (left <= right) ? +infi : -infi);

    // Adjust right boundary for open intervals, 
    // considering proper and improper intervals
    if (type == RIGHT_OPEN || type == OPEN)
        adjustedRight = nextafter(right, (left <= right) ? -infi : +infi);

    // Return the maximum of the adjusted boundaries
    return std::max(adjustedLeft, adjustedRight);
}

// Return interval midpoint
template<class _IT> inline _IT interval<_IT>::mid() const
{ if (right == left) return left; else  return (right + left) / _IT(2); }

// Return interval radius
template<class _IT> inline _IT interval<_IT>::rad() const
{ _IT r; r = (right - left) / _IT(2); return r; }

// Return interval width
template<class _IT> inline _IT interval<_IT>::wid() const
{ _IT r; r = right - left; if (r < _IT(0)) r = -r; return r; }

// Return mignitude of class
template<class _IT> inline _IT interval<_IT>::mig() const
{
       _IT l = inf();
       _IT r = sup();
       l = abs(l);
       r = abs(r);
       return std::min(l, r);
}

// Return magnitude of interval
template<class _IT> inline _IT interval<_IT>::mag() const
{
       _IT l = inf();
       _IT r = sup();
       l = abs(l);
       r = abs(r);
       return std::max(l, r);
}

// Required is... methods
template<class _IT> inline bool interval<_IT>::isProper() const
{ return left<=right; }
template<class _IT> inline bool interval<_IT>::isImproper() const
{ return right<left; }
template<class _IT> inline bool interval<_IT>::isPoint() const
{ return left==right; }
template<class _IT> inline bool interval<_IT>::isEmpty() const
{ return type==EMPTY; }

// Return interval classification
template<class _IT> inline enum int_class interval<_IT>::isClass() const
{
       if (left == _IT(0) && right == _IT(0)) return ZERO;
       if (left == _IT(0) && right >  _IT(0)) return POSITIVE0;
       if (left >  _IT(0) && right >  _IT(0)) return POSITIVE1;
       if (left >= _IT(0) && right >  _IT(0)) return POSITIVE;
       if (left <  _IT(0) && right == _IT(0)) return NEGATIVE0;
       if (left <  _IT(0) && right <  _IT(0)) return NEGATIVE1;
       if (left <  _IT(0) && right <= _IT(0)) return NEGATIVE;
       if (left <  _IT(0) && right >  _IT(0)) return MIXED;
       return NO_CLASS;
}

// Return the interval as a String representation
template<class _IT> inline std::string interval<_IT>::toString() const
{
       std::string s;
       enum interval_type t = intervaltype();
       std::ostringstream strs;

       strs.precision(25);
       strs << (t == LEFT_OPEN || t == OPEN ? "(" : "[");
       strs << left;
       strs << ",";
       strs << right;
       strs << (t == RIGHT_OPEN || t == OPEN ? ")" : "]");
       return strs.str();
}

实现转换运算符

实现转换运算符是微不足道的，因为它只需要在转换为任何整数类型或`float`和`double`类型时返回区间的中心点。如果区间为空，则转换是未定义的。

// Conversion Operators
template<class _IT> inline interval<_IT>::operator short() const
{      // Conversion to short
       return static_cast<short>(mid());
}

template<class _IT> inline interval<_IT>::operator int() const
{      // Conversion to int
       return static_cast<int>(mid());
}

template<class _IT> inline interval<_IT>::operator long() const
{      // Conversion to long
       return static_cast<long>(mid());
}

template<class _IT> inline interval<_IT>::operator long long() const
{      // Conversion to long long
       return static_cast<long long>(mid());
}

template<class _IT> inline interval<_IT>::operator unsigned short() const
{      // Conversion to unsigned short
       return static_cast<unsigned short>(mid());
}

template<class _IT> inline interval<_IT>::operator unsigned int() const
{      // Conversion to unsigned int
       return static_cast<unsigned int>(mid());
}

template<class _IT> inline interval<_IT>::operator unsigned long() const
{      // Conversion to unsigned long
       return static_cast<unsigned long>(mid());
}

template<class _IT> inline interval<_IT>::operator unsigned long long() const
{      // Conversion to long long
       return static_cast<unsigned long long>(mid());
}

template<class _IT> inline interval<_IT>::operator double() const
{      // Conversion to double
       return static_cast<double>(mid());
}

template<class _IT> inline interval<_IT>::operator float() const
{      // Conversion to float
       return static_cast<float>(mid());
}

实现基本运算符

现在我们可以通过软件仿真来控制舍入，以可移植的方式实现区间算术就变得非常简单。赋值只是将`=`运算符右侧的值复制到左侧。

// Assignment operator. Works for all class types
//
template<class _IT> inline interval<_IT>& interval<_IT>::operator=( const interval<_IT>& a )
{
       left = a.left;
       right = a.right;
       type = a.type;
       return *this;
}

区间加法

区间加法为

使用`fastwo_sum`进行加法的代码。其中`.inf()`返回区间的最小值，`.sup()`返回区间的最大值。

// += operator. Works for all classes.
// Always return a "proper" and closed [] interval
// a:=a+[EMPTY] or b:=[EMPTY]+b or [EMPTY]:=[EMPTY]+[EMPTY]
template<class _IT> inline interval<_IT>& interval<_IT>::operator+=(const interval<_IT>& a)
{
    // Handle EMPTY interval first
    if (a.type == EMPTY)
        return *this;
    if (type == EMPTY)
        return (*this = a);

    const _IT infi(INFINITY);
    // Neither a or b is [EMPTY]
    std::pair<_IT, _IT> xlow = fasttwo_sum(inf(), a.inf());
    std::pair<_IT, _IT> xright = fasttwo_sum(sup(), a.sup());
    left = xlow.first;
    right = xright.first;
    // Any adjustment?
    if (xlow.second < _IT(0))
        left = nextafter(left, -infi);
    if (xright.second > _IT(0))
        right = nextafter(right, +infi);
    type = CLOSE;
    return *this;
}

区间减法

区间减法是

使用`fasttwo_sum`进行减法的代码。

// -= operator. Works for all classes.
// Always return a "proper" and closed [] interval
// a=a-[EMPTY] or -b=[EMPTY]-b or [EMPTY]=[EMPTY]-[EMPTY]
//
template<class _IT> inline interval<_IT>& interval<_IT>::operator-=(const interval<_IT>& a)
{
    // Handle EMPTY interval first
    if (a.type == EMPTY)
        return *this;
    if (type == EMPTY)
        return (*this = -a);

    const _IT infi(INFINITY);
    // Neither a or b is [EMPTY]
    std::pair<_IT, _IT> xlow = fasttwo_sum(inf(), -a.sup());
    std::pair<_IT, _IT> xright = fasttwo_sum(sup(), -a.inf());
    left = xlow.first;
    right = xright.first;
    if (xlow.second < _IT(0))
        left = nextafter(xlow.first, -infi);
    if (xright.second > _IT(0))
        right = nextafter(xright.first, +infi);
    type = CLOSE;
    return *this;
}

区间乘法

区间乘法是

使用`fasttwo_product`算法和一些优化进行乘法的代码。

// Works for all classes.
// [EMPTY]:=a*[EMPTY] or [EMPTY]:=[EMPTY]*b or [EMPTY]:=[EMPTY]*[EMPTY]
// Please note that this is for all integer classes. interval<int>, interval<long>,
// were there is no loss of precision
// Instead of doing the mindless low = MIN(low*a.right, low*a.low,right*a.low,right*a.right) and
// right = MAX(low*a.right, low*a.low,right*a.low,right*a.right) requiring a total of 8 multiplication
//
//   low, right, a.low, a.right    result
//    +     +     -     +        -  +  [ right*a.low, right*a.right ] 2205
//    +     +     -     -        -  -  [ right*a.low, low*a.right ]
//    +     +     +     +        +  +  [ low*a.low, right*a.right ]
//    -     +     +     +        -  +  [ low*a.right, right*a.right ]
//    -     +     -     +        -  +  [ MIN(low*a.right,right*a.low), MAX(low*a.low,right*a.right) ]
//    -     +     -     -        -  -  [ right*a.low, low*a.low ]
//    -     -     +     +        -  -  [ low*a.right, right,a.low ]
//    -     -     -     +        -  -  [ low*a.right, low*a.low ]
//    -     -     -     -        +  +  [ right*a.right, low * a.low ]
//
template<class _IT> inline interval<_IT>& interval<_IT>::operator*=(const interval<_IT>& a)
{
    // Handle EMPTY interval first  ∅
    if (type == EMPTY)
        return *this;
    if (a.type == EMPTY)
        return (*this = a);

    // Neither a or b is ∅
    // Extract intervals through inf() and sup()
    _IT al = inf();
    _IT ah = sup();
    _IT bl = a.inf();
    _IT bh = a.sup();

    auto multiply = [&](const _IT& x, const _IT& y)
    {
        std::pair<_IT, _IT> tmp = interval<_IT>::fasttwo_prod(x, y);
        interval<_IT> res;
        const _IT infi(INFINITY);
        res.left = res.right = tmp.first;
        if (tmp.second < _IT(0))
            res.left = nextafter(tmp.first, -infi);
        if (tmp.second > _IT(0))
            res.right = nextafter(tmp.first, +infi);
        return res;
    };

    interval<_IT> itmp;
    type = CLOSE;
    // Shortcuts
    if (al >= _IT(0) && bl >= _IT(0))
    {// Both intervals positive
        itmp = multiply(al, bl);
        left = itmp.left;
        itmp = multiply(ah, bh);
        right = itmp.right;
        return *this;
    }
    if (ah < _IT(0) && bh < _IT(0))
    {// Both intervals negative
        itmp = multiply(al, bl);
        right = itmp.right;
        itmp = multiply(ah, bh);
        left = itmp.left;
        return *this;
    }
    if (al >= _IT(0) && bh < _IT(0))
    {// [A] interval positive, [B] interval negative
        itmp = multiply(ah, bl);
        left = itmp.left;
        itmp = multiply(al, bh);
        right = itmp.right;
        return *this;
    }
    if (ah < _IT(0) && bl >= _IT(0))
    {// [A] interval negative, [B] interval positive
        itmp = multiply(al, bh);
        left = itmp.left;
        itmp = multiply(ah, bl);
        right = itmp.right;
        return *this;
    }
    // Otherwise, we have a mixed interval. Make all 4 combinations
    itmp = multiply(al, bl);
    left = itmp.left;
    right = itmp.right;
    itmp = multiply(al, bh);
    left = min(left, itmp.left);
    right = max(right, itmp.right);
    itmp = multiply(ah, bl);
    left = min(left, itmp.left);
    right = max(right, itmp.right);
    itmp = multiply(ah, bh);
    left = min(left, itmp.left);
    right = max(right, itmp.right);

    return *this;
}

区间除法

区间除法是

// Works for all classes
// [EMPTY]:=a/[EMPTY] or [EMPTY]:=[EMPTY]/b or [EMPTY]:=[EMPTY]/[EMPTY]
// Please note that this is for all integer classes. interval<int>,interval<long>
// where there is no loss of precision
template<class _IT> inline interval<_IT>& interval<_IT>::operator/=(const interval<_IT>& b)
 {
     const _IT infi(INFINITY);
     // Handle EMPTY interval first
     if (type == EMPTY)
         return *this;
     if (b.type == EMPTY)
         return (*this = b);

     interval<_IT> tmp = *this;  // Save a copy
     // Compute the reverse of y e.g. 1/y
     auto inverse = [&](const _IT& y, const bool up)
     {
         _IT res = _IT(1) / y, r;

         r = -fma(res, y, _IT(-1));
         if (up == false)
         {
             if (r < _IT(0))
                 res = nextafter(res, -infi);
         }
         else
         {
             if (r > _IT(0))
                 res = nextafter(res, +infi);
         }

         return res;
     };

     _IT bl = b.inf();
     _IT bh = b.sup();
     if (bl == _IT(0) && bh == _IT(0))
     {
         left = -infi;
         right = +infi;
         *this *= tmp;
         return *this;
     }
     if (bh == _IT(0))
     {    // b.low is !=0
         right = inverse(bl, true);
         left = -infi;
         *this *= tmp;
         return *this;
     }
     if (bl == _IT(0))
     {    // b.right is !=0
         left = inverse(bh, false);
         right = +infi;
         *this *= tmp;
         return *this;
     }
     // neither b.low or b.right is zero
     left = inverse(bh, false);
     right = inverse(bl, true);
     *this *= tmp;
     return *this;
 }

一元和二元算术运算符的实现

我们只需使用基本运算符来计算二元运算符以简化代码。

// Binary + operator specialization for only interval<_IT> arguments
// Works for all classes
//
template<class _IT> inline interval<_IT> operator+(const interval<_IT>& a, const interval<_IT>& b)
{
       interval<_IT> result(a);
       result += b;
       return result;
}

// Unary + operator
// Works for all classes
//
template<class _IT> inline interval<_IT> operator+(const interval<_IT>& a)
{
       return a;
}

// Binary - operator
// Works for all classes
//
template<class _IT> inline interval<_IT> operator-(const interval<_IT>& a, const interval<_IT>& b)
{
       interval<_IT> result(a);
       result -= b;
       return result;
}

// Unary - operator
// Works for all classes
//
template<class _IT> inline interval<_IT> operator-(const interval<_IT>& a)
{
       interval<_IT> result(0);
       result -= a;
       return result;
}

// Binary * operator
// Works for all classes
//
template<class _IT> inline interval<_IT> operator*(const interval<_IT>& a, const interval<_IT>& b)
{
       interval<_IT> result(a);
       result *= b;
       return result;
}

//Binary / operator
// Works for all classes
//
template<class _IT> inline interval<_IT> operator/(const interval<_IT>& a, const interval<_IT>& b)
{
       interval<_IT> result(a);

       if (result == b)
       {
             enum int_class intclass = b.isClass();
             if (intclass != ZERO && intclass != POSITIVE0 && intclass != NEGATIVE0)
             {
                    result = interval<_IT>(1, 1);
                    return result;
             }
       }
       result /= b;
       return result;
}

比较运算符的实现

在本节中，我们处理所有6个比较运算符：`==`、`!=`、`>=`、`<=`、`>`、`<`。请注意，我们还需要能够处理空区间∅。如果两个区间都是∅，则返回`false`。如果其中一个区间是∅，则返回`true`。请注意，有些实现当两个区间都是∅时返回未定义。

// == operator
// Works for all classes
//
template<class _IT> inline bool operator==(const interval<_IT>& a, const interval<_IT>& b)
{
       if (a.intervaltype() == EMPTY && b.intervaltype() == EMPTY)
             return false; // Both EMPTY=> return false
       if (a.intervaltype() == EMPTY || b.intervaltype() == EMPTY)
             return true;  // One but not both are EMPTY => return true
       // Check for equality. By using inf() and sup() we do not have to worry
       about Improper interval
       // or interval type != CLOSE
       return a.inf() == b.inf() && a.sup() == b.sup();
}

// != operator
// Works for all classes
//
template<class _IT> inline bool operator!=(const interval<_IT>& a, const interval<_IT>& b)
{
       return !(a == b);
}

// >= operator
// Works for all classes
//
template<class _IT> inline bool operator>=(const interval<_IT>& a, const interval<_IT>& b)
{
       if (a.intervaltype() == EMPTY && b.intervaltype() == EMPTY)
             return false;
       if (a.intervaltype() == EMPTY || b.intervaltype() == EMPTY)
             return true;
       if (a.inf() >= b.sup())
             return true;
       return false;
}

// <= operator
// Works for all classes
//
template<class _IT> inline bool operator<=(const interval<_IT>& a, const interval<_IT>& b)
{
       if (a.intervaltype() == EMPTY && b.intervaltype() == EMPTY)
             return false;
       if (a.intervaltype() == EMPTY || b.intervaltype() == EMPTY)
             return true;
       if (a.sup() <= b.inf())
             return true;
       return false;
}

// > operator
// Works for all classes
//
template<class _IT> inline bool operator>(const interval<_IT>& a, const interval<_IT>& b)
{
       if (a.intervaltype() == EMPTY && b.intervaltype() == EMPTY)
             return false;
       if (a.intervaltype() == EMPTY || b.intervaltype() == EMPTY)
             return true;
       if (a.inf() > b.sup() )
             return true;
       return false;
}

// < operator
// Works for all classes
//
template<class _IT> inline bool operator<(const interval<_IT>& a, const interval<_IT>& b)
{
       if (a.intervaltype() == EMPTY && b.intervaltype() == EMPTY)
             return false;
       if (a.intervaltype() == EMPTY || b.intervaltype() == EMPTY)
             return true;
       if (a.sup() < b.inf())
             return true;
       return false;
}

实现sqr(x)和sqrt(x)函数

我们需要为处理x2添加特殊代码。请参阅前面的章节。对于`sqrt(x)`，我们再次可以使用`fma()`库函数来快速确定区间边界。

// sqr(x)=x^2
template<class _IT> inline interval<_IT> sqr(const interval<_IT>& x)
{
    _IT left = x.inf();
    _IT right = x.sup();
    _IT tmpl = left;
    _IT tmpr = right;
    interval<_IT> r(0);

    left *= left;
    right *= right;
    r.rightinterval(max(left, right));
    // Contained zero?
    if ( tmpl > 0 && tmpr > 0)
        r.leftinterval(min(left, right));
    return r;
}

// sqrt(x)
template<class _IT> inline interval<_IT> sqrt(const interval<_IT>& x)
{
    const _IT infi(INFINITY);
    _IT sq, r;
    _IT left, right;

    // Find leftinterval bound
    sq = sqrt(x.inf());
    left = sq;
    r = -fma(sq, sq, -x);
    if (r < _IT(0))
        left = nextafter(left, -infi);

    // Find rightinterval bound
    sq = sqrt(x.sup());
    right = sq;
    r = -fma(sq, sq, -x);
    if (r > _IT(0))
        right = nextafter(right, +infi);
    return interval<_IT>(left, right);
}

实现abs(x)函数

在区间算术中，区间的绝对值（`abs`）定义为包含区间内任何数字所有可能的绝对值。结果始终是非负区间。定义取决于区间相对于零的位置。

如果区间完全非负（即，其下限和上限都大于或等于零），则区间的绝对值就是区间本身，因为其中的所有值都已经是非负的。

如果区间完全为负（即，其下限和上限都小于零），则区间的绝对值是通过取其边界的绝对值并交换它们而获得的（因为下限将是最大的负数，因此具有最高的绝对值）。

如果区间跨越零（即，其下限为负，上限为正），则区间的绝对值是从零到下限和上限绝对值中的最大值。

`abs(x)`的源代码

// abs([a,b])
// if a>=0 in [a,b] then |[a,b]|==[a,b]
// if b<0 in [a,b] then |[a,b]|=[-b,-a]
// if a<0 & b>0 in [a,b] then |[a,b]|=[0,max(-a,b)]
template<class _IT> inline interval<_IT> abs( const interval<_IT>& a )
       {
       if (a.inf() >= _IT(0) ) // entirely positive
             return a;
       else
             if (a.sup() < _IT(0) ) // Entirely negative
                    return -a;

       return interval<_IT>(_IT(0), max(-a.inf(),a.sup()));
       }

cin和cout运算符的实现

标准输出`cout`和输入`cin`都很直接。`cout`以以下格式输出区间：

[ 左, 右 ]

当区间是闭合的

( 左, 右 ] 或 [ 左, 右 ) 用于半开区间

并且

( 左, 右 ) 用于全开区间。

如果区间为空，则输出∅。

而`cin`可以处理上述的一些变化。任何输入语法都是合法的。

            [ left, right ]
            [ interval ] which is equivalent to [ interval, interval ]

您可以使用'('或')'代替方括号'['或']'来表示半开或全开区间。

第一阶段总结

这是第一阶段的目标，定义并实现了下一阶段所需的所有基本功能。

然而，我们确实有一个不足之处。结果并不完全准确。这与构造函数有关。构造函数参数是模板类型`_IT`（可以是`float`、`double`或`long double`）。根据C++标准规则，当实际参数的类型小于`_IT`时，会自动转换为所需的`_IT`类型。但是，转换发生在执行构造函数之前，因此转换中的任何不准确性都不会被区间类捕获。一个例子将阐明这个问题。

Interval<float> x(16777216);
Interval<float> y(16777217);

将为x和y创建相同的区间[16777216.0,16777216.0]。问题在于整数16777217大于可以精确转换为浮点数23位尾数的值，而16777216恰好可以转换。

同样地，将64位双精度浮点数转换为32位浮点数类型也是如此。即使对于`interval`变量，大于9,007,199,254,740,992的64位整数也无法正确转换为`double`（双精度浮点数的尾数为52位）。

在第二阶段，我们将引入混合数据类型，它将神奇地帮助我们解决和修复上述问题。

第一阶段的区间模板库的实际实现到此结束。我们现在已经建立了基本功能，可以继续实施第二阶段和第三阶段。

要阅读包含第二阶段和第三阶段实现的完整论文，完整的源代码和文档可在我的网站上找到：区间算术 (hvks.com)。

参考文献

  1. 英特尔64和IA32架构软件开发手册。第2卷。2014年6月。
  2. Wilkinson, J H, Algebraic Processes中的舍入误差, Prentice-Hall Inc, Englewood Cliffs, NJ 1963
  3. Richard Brent & Paul Zimmermann, 现代计算机算术, 版本 0.5.9 2010年10月17日; http://maths-people.anu.edu.au/~brent/pd/mca-cup-0.5.9.pdf
  4. T. Hickey & M.H. van Emden。区间算术：从原理到实现。
  5. Hend Dawood，区间算术理论，数学基础与应用，Lambert Academic Publishing。2011年。ISBN：978-3-8465-0152-2。
  6. 维基百科，区间算术，https://en.wikipedia.org/wiki/Interval_arithmetic，
  7. Daniel Pfeffer，使用就近舍入的区间算术，https://codeproject.org.cn/Articles/1040839/Interval-arithmetic-using-round-to-nearest-mode-pa
  8. N.T. Hayes，2013年2月，模态区间算术标准。
  9. J.G. Role，区间算术与分析：导论
© . All rights reserved.