C/C++/FORTRAN 中的四倍精度浮点数计算

Shmuel Safonov

5.00/5 (2投票s)

2017年9月18日

CPOL

6分钟阅读

32444

1027

创建了一个以四倍浮点精度进行计算的数值库，并从 MSVC C/C++ 代码中使用该库

引言

本文描述了如何构建一个执行四倍精度浮点数计算的数值库，以及如何从 MSVC C/C++ 代码中访问该库。该技术通过一个示例进行说明。该示例展示了与标准双精度计算相比，计算精度有显著提高。

背景

大多数现有的商业编译器都支持双精度计算作为可能达到的最高精度。这是由于 CPU 的硬件限制。扩展精度需求会涉及到在现有硬件上模拟精确计算，从而导致计算时间大大增加。

GNU MPFR (http://www.mpfr.org/) 等多精度库可以进行任意精度的计算。但是，使用这些库需要用户在很大程度上重写现有的 C 或 FORTRAN 代码。四倍浮点数方法是双精度计算和多精度计算之间的一个良好折衷，因为它不需要重写现有代码（假设它受到编译器的支持）。

建议在 C/C++ 或 FORTRAN 动态库或静态库中实现四倍精度浮点数计算。通过基于 LAPACK 库使用多项式根计算的示例，展示了精度提高的优势。为此，LAPACK 库使用四倍精度浮点数进行编译。

浮点数格式和变量类型

浮点数按照 IEEE 浮点数运算标准 (IEEE 754) 定义的格式表示，请参见 Wiki (https://en.wikipedia.org/wiki/IEEE_754)。浮点数表示如下：

| sign bit (S) |exponent bits (exp) | mantissa bits (man) |

表示的数字是 (-1)^S * 1.mantissa * 2^(exp)，其中 'mantissa' 由二进制数字组成，'exp' 被视为带符号的数字。

我们将考察以下基本格式：

格式	描述	尾数位数	指数位数	十进制数字	十进制指数范围
binary32	单精度	24	8	7.22	38.23
binary64	双精度	53	11	15.95	307.95
binary128	四倍精度	113	15	34.02	4931.77

下表显示了 C++/FORTRAN 中对应的变量类型。

	binary32	binary64	binary128
C/C++	`float`	`double`	`_Quad` (Intel), `__float128` (GNU quadmath)
FORTRAN (real)	REAL*4	REAL*8	REAL*16
FORTRAN (complex)	COMPLEX*8	COMPLEX*16	COMPLEX*32

支持在 x86 和 x86-64 CPU 上进行四倍精度浮点数仿真的流行编译器是 Intel C/C++ 和 Fortran 以及 GNU GCC 和 GFORTRAN（使用 quadlib 库）。Microsoft VC 编译器不支持此类型，并且在可预见的未来似乎也不会支持！

LAPACK 库的四倍精度浮点数编译

LAPACK 库是一个标准的数值包 (http://www.netlib.org/lapack/)，实现了线性代数算法。本文使用该库计算多项式的根。该库是用 FORTRAN 编写的，但它有 C 移植版本 CLAPACK。该库的标准版本仅提供双精度功能。缺乏更高精度的实现很可能是因为双精度算术对大多数应用程序来说已足够。

提出了两种 LAPACK 四倍精度修改方式：

修改 LAPACK FORTRAN 源码
修改 CLAPACK C 源码

在第一种情况下，变量 `REAL*8, REAL, DOUBLE PRECISION, COMPLEX` 的类型声明应替换为 `REAL*16, COMPLEX*32`。

如果使用 FORTRAN90 标准，类型应分别替换为 `real (kind=selected_real_kind (32))` 和 `complex (kind=selected_real_kind (32))`。

将内建函数的 `REAL*8` 和 `COMPLEX*16` 版本替换为其 `REAL*16` 和 `COMPLEX*32` 版本。此过程似乎需要对代码进行大量修改，但相当直接。

在第二种情况下，应修改 include 文件 f2c.h 以重新定义基本类型。标准函数的原型也应替换为四倍精度版本。

多项式根计算

多项式根使用伴随矩阵法计算。伴随矩阵的特征值就是原多项式的根。伴随矩阵由多项式系数计算得出，矩阵特征值使用 LAPACK 函数 `zgeevx` 计算。

Using the Code

以下是关于如何使用本文进行编码的简要说明。

示例代码由以下部分组成：
test_qroot.exe：Win32 可执行文件，由 MSVC 构建
qroot.dll：DLL，四倍精度实用程序，
由 GNU C++（使用 QuadMath）或 Intel C/C++ 或 Intel Visual Fortran 构建

lapack_qd_w32.dll：DLL，四倍精度 LAPACK 库，
由 GNU C++（使用 QuadMath）或 Intel C/C++ 或 Intel Visual Fortran 构建

qroot.dll 中的四倍精度浮点数类型声明如下：

// File: clapack_types.h
#ifdef __INTEL_COMPILER
typedef _Quad real;
typedef _Quad doublereal;
#else
#include <quadmath.h>
typedef __float128 real;
typedef __float128 doublereal;
#endif

qroot.dll 中的导出函数原型定义如下：

//
// Calculate polynomial roots
//
extern "C"
{
__declspec(dllexport) void qroots(doublecomplex *r, doublecomplex *p, long *degP);
}
...

调用应用程序中的相同类型声明如下：

// File: MyQuad.h
struct CMyQuad
{
long a1;
long a2;
long a3;
long a4;
};
typedef struct CMyQuad MyQuad;
typedef MyQuad* PQuad;

struct CMyCQuad
{
 long ra1;
 long ra2;
 long ra3;
 long ra4;
 long ia5;
 long ia6;
 long ia7;
 long ia8;
};
typedef struct CMyCQuad MyCQuad;
typedef MyCQuad* PCQuad;
...

test_qroot.exe 中相同的导入函数原型定义为：

//
// Calculate polynomial roots
//
extern "C"
{
 void  qroots(PCQuad r, PCQuad p, long *degP);
}
...

此处必须保持相同的类型大小！函数参数应通过指针传递，这也保持了与 FORTRAN 的兼容性。

现在，您必须设置额外的编译选项。

对于 Intel C/C++ 编译器，选项如下。它们启用了 `_Quad` 类型的用法。

-Qoption,cpp,--extended_float_type

对于 Intel Visual Fortran 编译器，选项如下。它们将默认的 real、complex、double precision 和 double complex 变量类型和常量解释为 64 位 real、128 位 complex、128 位 real 和 256 位 complex。请注意，确切的类型应手动替换。还建议增加 LAPACK 编译允许的行长。

/real-size:128 /Qimf-precision:high /extend_source:132

对于 GNU C++ 编译器，您只需命令链接 QuadMath 库。

-lquadmath

对于 GNU Gfortran 编译器，选项是：

-fdefault-real-8 -ffixed-line-length-132

`PCQuad` 类型的中间结果可以传递给其他实用程序函数。最后，实用程序函数可以将结果转换为标准的 `double` 类型。

提示：如果 DLL 是由 GCC 构建的，那么在用 MSVC 链接示例时，请在高级 MSVC 链接器选项中设置 /SAFESEH:NO。

数值结果

使用标准双精度进行的计算在 MATLAB 中通过以下脚本完成：

% Create polynom P(x) = ((x-x00)(x-conj(x00)) )^9
% Calculate the roots and compare with the original root x00.
ar = 0.00001;
ai  = 10;

x0 = ar + ai*j;
x1 = ar - ai*j;

va = [1, -2*ar,ar^2+ai^2];
vb = va;
vb2 = conv(va,va);
vb4 = conv(vb2,vb2);
vb8 = conv(vb4,vb4);
vb9 = conv(vb8,va);
vb9r = real(vb9);
aa = vb9r;

xx =  roots( aa);
 disp('Roots=');
disp(xx);
for cnt=1:length(xx)
    yy(cnt) = polyval(aa,xx(cnt));   
end
disp('P(roots)=');
disp(yy.');

结果是：

% Roots=
%   0.112215172144653 +10.131442757995252i
%   0.112215172144653 -10.131442757995252i
%   0.000010235951691 +10.173357532550744i
%   0.000010235951691 -10.173357532550744i
%  -0.112194813184584 +10.131443062600157i
%  -0.112194813184584 -10.131443062600157i
%  -0.169232753684938 +10.027588249761465i
%  -0.169232753684938 -10.027588249761465i
%   0.169252831553145 +10.027587788028235i
%   0.169252831553145 -10.027587788028235i
%  -0.146212177557002 + 9.913350066605950i
%  -0.146212177557002 - 9.913350066605950i
%   0.146231942318667 + 9.913349664139163i
%   0.146231942318667 - 9.913349664139163i
%  -0.057077609506728 + 9.840940518344002i
%  -0.057077609506728 - 9.840940518344002i
%   0.057097171965100 + 9.840940359975059i
%   0.057097171965100 - 9.840940359975059i

% P(roots)=
%   1.0e+004 *
%
%  -3.584000000000000 + 0.515000000000000i
%  -3.584000000000000 - 0.515000000000000i
%  -7.155200000000000 + 0.000062243652344i
%  -7.155200000000000 - 0.000062243652344i
%  -5.324800000000000 - 0.705800000000000i
%  -5.324800000000000 + 0.705800000000000i
%  -4.300800000000000 - 1.023400000000000i
%  -4.300800000000000 + 1.023400000000000i
%  -4.838400000000000 + 0.984000000000000i
%  -4.838400000000000 - 0.984000000000000i
%  -6.412800000000000 - 1.089000000000000i
%  -6.412800000000000 + 1.089000000000000i
%  -6.464000000000000 + 1.140400000000000i
%  -6.464000000000000 - 1.140400000000000i
%  -3.520000000000000 - 0.276900000000000i
%  -3.520000000000000 + 0.276900000000000i
%  -5.030400000000000 + 0.330800000000000i
%  -5.030400000000000 - 0.330800000000000i
%

本示例中的相同计算得出根为：

/// Roots:
A(0,0)=(1.00761e-005,-10.0016)
A(1,0)=(-0.000986753,-10.0012)
A(2,0)=(0.00100687,-10.0012)
A(3,0)=(0.00153725,-10.0003)
A(4,0)=(-0.00151723,-10.0003)
A(5,0)=(-0.0013331,-9.99922)
A(6,0)=(0.00135303,-9.99922)
A(7,0)=(-0.000520498,-9.99854)
A(8,0)=(0.000540355,-9.99854)
A(9,0)=(1.08464e-005,10.0014)
A(10,0)=(0.000940125,10.0011)
A(11,0)=(-0.000918832,10.0011)
A(12,0)=(0.00143429,10.0003)
A(13,0)=(-0.001414,10.0003)
A(14,0)=(0.00126206,9.99928)
A(15,0)=(-0.0012429,9.99928)
A(16,0)=(0.000503879,9.99864)
A(17,0)=(-0.00048546,9.9986)

以及它们在原始多项式 P(roots) 上的求值：

polyval(0.000010,-10.001551)=-5.75096e-014+ i-1.21989e-017
polyval(-0.000987,-10.001188)=-4.41869e-014+ i1.55312e-017
polyval(0.001007,-10.001188)=-5.32907e-014+ i-4.62412e-017
polyval(0.001537,-10.000269)=-5.00711e-014+ i-5.98751e-017
polyval(-0.001517,-10.000269)=-5.63993e-014+ i5.07407e-017
polyval(-0.001333,-9.999225)=-3.36398e-014+ i1.62088e-017
polyval(0.001353,-9.999225)=-3.90799e-014+ i-4.17824e-017
polyval(-0.000520,-9.998543)=-5.80647e-014+ i8.72783e-018
polyval(0.000540,-9.998543)=-5.973e-014+ i-3.54263e-017
polyval(0.000011,10.001446)=-2.81997e-014+ i7.48016e-017
polyval(0.000940,10.001107)=-3.66374e-014+ i9.91367e-017
polyval(-0.000919,10.001108)=-1.12133e-014+ i6.8237e-017
polyval(0.001434,10.000250)=-4.25215e-014+ i1.1563e-016
polyval(-0.001414,10.000252)=-3.06422e-014+ i4.59973e-017
polyval(0.001262,9.999276)=-4.70735e-014+ i1.17243e-016
polyval(-0.001243,9.999278)=-3.43059e-014+ i4.35578e-017
polyval(0.000504,9.998641)=-4.07452e-014+ i8.93112e-017
polyval(-0.000485,9.998641)=-3.9968e-014+ i5.99157e-017

我们可以看到，绝对根误差从 0.16 降低到 0.002；多项式求值从 20000 降低到 1e-14。

关注点

在我选择上述技术之前，我尝试了多精度 LAPACK 库 (MLAPACK)。我没有成功，因为函数 `zgeevx` 没有实现，只有一些处理对称矩阵的函数能够工作。

我还尝试使用新的面向对象的 FORTRAN 标准，并通过使用 MFPR 库的实现来封装标准类型。由于 GNU GFORTRAN 实现中析构函数的错误功能，我也失败了。

事实证明，Intel 编译器在支持 `_Quad` 变量类型方面比 MSVC 具有优势。此外，Intel Visual Fortran 中的数值编程比 Intel C/C++ 简单得多。它支持对四倍精度变量进行调试、查看和打印，而 C++ 环境则没有这些调试功能。`_Quad` 类型变量的正确代码高亮显示也缺失。

但 FORTRAN 的主要优势在于对复数进行直接操作以及处理矩阵和数组。这似乎是令人怀念的 DEC 编译器的伟大遗产，但现在大多数开发人员都对 FORTRAN 敬而远之了 :）。

而且，如果您不想为商业编译器付费，请使用 GNU C/C++/GFORTRAN。

附注：有趣的是，多个
多项式根的计算精度可以通过使用
ZHONGGANG ZEN 开发的代数技术并在其 MultRoot 包中实现来提高。
在其 MultRoot 包中实现。

[1] ZHONGGANG ZENG, "COMPUTING MULTIPLE ROOTS OF INEXACT POLYNOMIALS",
MATHEMATICS OF COMPUTATION,
Volume 74, Number 250, Pages 869-903
S 0025-5718(04)01692-8
Article electronically published on July 22, 2004

[2] ZHONGGANG ZENG, "MultRoot - A Matlab package computing polynomial roots and multiplicities"
Journal ACM Transactions on Mathematical Software (TOMS), Volume 30 Issue 2, June 2004, Pages 218-236