【翻译】用于科学计算的 C++11/14 新特性之六

原文地址 www.numbercrunch.de

【翻译】用于科学计算的 C++11/14 新特性之六

浮点异常

在数值计算过程中，可能会发生各种运行时错误。在 C++ 中，这类错误可以通过浮点异常或通过（全局但线程局部的）变量 errno 来指示。浮点异常与 C++ 异常完全无关。当一个错误的操作引发浮点异常时，异常仅在浮点状态字中被记录，程序会继续执行。该操作产生一个默认值，该值取决于异常（见下表）。C++11 引入了新的头文件 cfenv，其中包含用于检查和操作浮点环境和浮点状态标志的函数。通过这种方式，程序可以找出哪些异常已经被引发。宏 math_errhandling 的值指示了浮点运算符和函数执行的错误处理类型。

IEEE 754 标准定义了五种浮点异常：

• 无效操作，例如，负数的平方根，默认返回静默 NaN。
• 除以零，即有限操作数的运算给出了一个精确的无限结果，例如 $1/0$ 或 $\log(0)$，默认返回 $\pm\infty$。
• 上溢，即结果太大，无法被浮点数正确表示，默认返回 $\pm\infty$（对于最近邻舍入模式）。
• 下溢，即结果非常小（超出正常范围）并且不精确，默认返回一个非规范化值。
• 不精确，即默认返回正确舍入的结果。

以下程序引发各种浮点异常，并展示如何检测它们。首先，必须告诉编译器将通过 #pragma STDC FENV_ACCESS ON 来调查浮点环境，这会强制编译器避免某些优化，这些优化可能会使浮点状态无效。函数 feclearexcept 清除指定的浮点状态标志，而 fetestexcept 确定哪些指定的浮点状态标志被设置。浮点异常由宏常量 FE_INVALID 等表示。程序的预期输出是：

MATH_ERRNO is set
MATH_ERREXCEPT is set

1.0/0.0 = inf
ERROR: division by zero

2^2000 = inf
ERROR: numerical overflow

1/2^2000 = 0
ERROR: numerical underflow

sqrt(-1) = -nan
ERROR: invalid result

#pragma STDC FENV_ACCESS ON
 
#include <cstdlib>
#include <iostream>
#include <cfenv>
#include <cmath>
 
void test_fp_exceptions() {
  bool error=false;
  if (std::fetestexcept(FE_DIVBYZERO)) {
    error=true;
    std::cerr << "ERROR: division by zero\n";
  }
  if (std::fetestexcept(FE_OVERFLOW)) {
    error=true;
    std::cerr << "ERROR: numerical overflow\n";
  }
  if (std::fetestexcept(FE_UNDERFLOW)) {
    error=true; 
    std::cerr << "ERROR: numerical underflow\n";
  }
  if (std::fetestexcept(FE_INVALID)) {
    error=true;
    std::cerr << "ERROR: invalid result\n";
  }
  if (not error)
    std::cerr << "no error\n";
  std::feclearexcept(FE_ALL_EXCEPT);
  std::cerr << '\n';
}
 
int main() {
  double zero=0.0;
 
  std::cout << "MATH_ERRNO is "
	    << (math_errhandling & MATH_ERRNO ? "set" : "not set") << '\n'
	    << "MATH_ERREXCEPT is "
	    << (math_errhandling & MATH_ERREXCEPT ? "set" : "not set") << '\n'
	    << '\n';
 
  std::feclearexcept(FE_ALL_EXCEPT);
 
  std::cout <<  "1.0/0.0 = " << 1.0/zero << '\n';
  test_fp_exceptions();
 
  std::cout << "2^2000 = " << std::pow(2., 2000) << '\n';
  test_fp_exceptions();
 
  std::cout << "1/2^2000 = " << std::pow(2., -2000) << '\n';
  test_fp_exceptions();
 
  std::cout << "sqrt(-1) = " << std::sqrt(-1.) << '\n';
  test_fp_exceptions();
 
  return EXIT_SUCCESS;
}

显式测试浮点异常可能会变得繁琐。GNU C 库提供了函数 feenableexcept、fedisableexcept 和 fegetexcept，用于启用（和禁用）浮点异常的陷阱。这意味着，浮点异常被转换为信号，可以被信号句柄捕获，如下程序所示，该程序需要 GNU C++ 编译器来编译。

#pragma STDC FENV_ACCESS ON
#include<cstdlib>
#include<iostream>
#include<cfenv>
#include<fenv.h>
#include<signal.h>
 
void floating_point_handler(int signal, siginfo_t *sip, void *uap) {
  std::cerr << "floating point error at " << sip->si_addr << " : ";
  int code=sip->si_code;
  if (code==FPE_FLTDIV)
    std::cerr << "division by zero\n";
  if (code==FPE_FLTUND)
    std::cerr << "underflow\n";
  if (code==FPE_FLTOVF)
    std::cerr << "overflow\n";
  if (code==FPE_FLTINV)
    std::cerr << "invalid result\n";
  std::abort();
}
 
int main() {
  std::feclearexcept(FE_ALL_EXCEPT);
  feenableexcept(FE_DIVBYZERO | FE_UNDERFLOW | FE_OVERFLOW | FE_INVALID);
  struct sigaction act;
  act.sa_sigaction=floating_point_handler;
  act.sa_flags=SA_SIGINFO;
  sigaction(SIGFPE, &act, NULL);
 
  double zero=0.0;
  double one=1.0;
 
  std::cout << "1.0/1.0 = " << one/one << '\n';
  std::cout << "1.0/0.0 = " << one/zero << '\n';
  std::cout << "1.0/1.0 = " << one/one << '\n';
 
  return EXIT_SUCCESS;
}

通常，信号句柄除了发出错误消息、执行一些清理工作以及终止程序外不会做太多事情，特别是当浮点错误是致命的（例如除以零）时候。此外，当程序退出信号句柄时，浮点异常会再次被引发。因此，信号句柄会被再次调用，这会导致信号句柄调用的无限循环。解决这个问题的方法是使用 sigsetjmp 和 siglongjmp 进行非局部跳转，以实现异常机制，如下所示。

#pragma STDC FENV_ACCESS ON
 
#include<cstdlib>
#include<iostream>
#include<cfenv>
#include<fenv.h>
#include<signal.h>
#include<setjmp.h>
 
static sigjmp_buf jmpbuf;
 
void floating_point_handler(int signal, siginfo_t *sip, void *uap) {
  std::cerr << "floating point error at " << sip->si_addr << " : ";
  int code=sip->si_code;
  if (code==FPE_FLTDIV)
    std::cerr << "division by zero\n";
  if (code==FPE_FLTUND)
    std::cerr << "underflow\n";
  if (code==FPE_FLTOVF)
    std::cerr << "overflow\n";
  if (code==FPE_FLTINV)
    std::cerr << "invalid result\n";
  feenableexcept(FE_DIVBYZERO | FE_UNDERFLOW | FE_OVERFLOW | FE_INVALID);
  siglongjmp(jmpbuf, 1);
}
 
int main() {
  std::feclearexcept(FE_ALL_EXCEPT);
  feenableexcept(FE_DIVBYZERO | FE_UNDERFLOW | FE_OVERFLOW | FE_INVALID);
  struct sigaction act;
  act.sa_sigaction=floating_point_handler;
  act.sa_flags=SA_SIGINFO;
  sigaction(SIGFPE, &act, NULL);
 
  double zero=0.0;
  double one=1.0;
 
  if (sigsetjmp(jmpbuf, 1)==0) {  // try
    std::cout << "1.0/1.0 = " << one/one << '\n';
    std::cout << "1.0/0.0 = " << one/zero << '\n';
    std::cout << "1.0/1.0 = " << one/one << '\n';
  } else {  // catch
    std::cerr << "some error occurred\n";
  }
  return EXIT_SUCCESS;
}

使用上述显示的捕捉浮点异常的技术，错误的数值计算变得容易检测，并使得相应地处理它们成为可能。