问题标签 [double-double-arithmetic]

我正在为仅支持 32 位单精度浮点运算的嵌入式硬件编写程序。然而，我正在实现的算法需要 64 位双精度加法和比较。我正在尝试double使用两个floats 的元组来模拟数据类型。所以 adouble d将被模拟为 astruct包含元组：(float d.hi, float d.low)。

使用字典顺序进行比较应该很简单。然而，添加有点棘手，因为我不确定我应该使用哪个基础。应该是FLT_MAX吗？以及如何检测进位？

如何才能做到这一点？

编辑（清晰）：我需要额外的有效数字而不是额外的范围。

在一段时间内关注双双运算的论文和源代码，我仍然无法找出 dd_real （定义为struct dd_real { double x[2];...}）数字是如何被分成两个双精度数的。假设我用字符串初始化它，anddd_real pi = "3.14159265358979323846264338327950";会是什么？我需要理解它，然后编写一个有希望的小型 Python 函数来完成它。pi.x[0]pi.xi[1]

我不只是想调用 QD 库的原因是我更愿意在 Python 中重新实现正确的拆分，以便将我的 35 位精度常量（以字符串形式给出）double2发送到 CUDA 代码的位置GQD 库将其视为双双实数——似乎是唯一一个处理 CUDA 中扩展精度计算的库。不幸的是，在 Python 方面也排除了 mpmath。

上下文：双双算术

“Double-double”是数字表示为两个双精度数字之和，有效数字没有重叠。这种表示利用现有的双精度硬件实现来进行“近四倍精度”计算。

double-double 实现中的一个典型的低级 C 函数可能需要两个双精度数a，b并|a| ≥ |b|计算(s, e)表示它们和的双精度数：

（改编自这篇文章。）

这些实现通常采用舍入到最近的偶数模式。

在上面的计算中，(s, e)是一个归一化的 double-double 只是因为这个假设。没有它，a == 0x1.0p60在b == 1向上舍入模式下，使用 ,s计算为0x1.0000000000001p60和e略高于-0x0.0000000000001p60。它们的和等于 的数学和，a但b它们的有效数字重叠。

Take和一方面 和 另一方面的a == 0x1.0p120数学和甚至不再重合。abse

问题

有没有办法构建一个具有与典型双双库在舍入到最近偶数（即相对快速且相对准确）中具有相同属性的双双类库，但无论舍入如何都有效模式恰好是？

这样的图书馆是否已经存在？

更一般的上下文：正确舍入的基本函数

双双排序的实现用于正确舍入基本函数库的实现中的中间计算。结果，以这种方式实现的库在 FPU 未处于舍入到最近偶数模式时调用函数时往往会失败。出于性能原因以及在函数执行时到达的信号会使 FPU 处于舍入到最近偶数模式，更改函数内部的舍入模式不是很可口。我认为拥有在任何舍入模式下工作的快速、正确舍入的基本函数的最简单方法是，如果可以以某种方式依赖于在任何舍入模式下工作的双双类型算术。

如果用两个单精度浮点来模拟双精度浮点，性能会怎样，能不能做好？

目前，Nvidia 对启用双精度的 Tesla 卡收取相当高的费用，这使您能够获得三分之一的单精度性能（Titan/Titan Black 例外）。

如果要使用具有 gimped 双精度的 Geforce GPU 并使用 2 个单精度浮点数模拟双精度，性能会如何？

当我第一次得到一个 Haswell 处理器时，我尝试实现 FMA 来确定 Mandelbrot 集。主要算法是这样的：

这确定n像素是否在 Mandelbrot 集中。因此，对于双浮点，它运行超过 4 个像素（floatn = __m256d, intn = __m256i）。这需要 4 个 SIMD 浮点乘法和 4 个 SIMD 浮点加法。

然后我修改了它以像这样使用 FMA

其中 mul_add 调用_mm256_fmad_pd和 mul_sub 调用_mm256_fmsub_pd。此方法使用 4 个 FMA SIMD 操作和两个 SIMD 乘法，这比没有 FMA 的算术操作少了两次。此外，FMA 和乘法可以使用两个端口，而加法只能使用一个。

为了使我的测试不那么有偏见，我放大了一个完全在 Mandelbrot 集中的区域，因此所有值都是maxiter. 在这种情况下，使用 FMA 的方法大约快 27%。这当然是一个进步，但是从 SSE 到 AVX 使我的表现翻了一番，所以我希望 FMA 可能会再增加两倍。

但后来我找到了关于 FMA 的答案，上面写着

fused-multiply-add 指令的重要方面是中间结果的（实际上）无限精度。这有助于提高性能，但不是因为两个操作被编码在一条指令中 - 它有助于提高性能，因为中间结果的几乎无限精度有时很重要，并且在这种级别的普通乘法和加法中恢复非常昂贵精度确实是程序员所追求的。

后面给出了一个 double*double 到double-double乘法的例子

由此，我得出结论，我实施 FMA 的方式不是最优的，因此我决定实施 SIMD 双双。我根据论文Extended-Precision Floating-Point Numbers for GPU Computation实现了 double-double 。该纸是用于双浮动的，所以我将其修改为双双。此外，我没有在 SIMD 寄存器中打包一个 double-double 值，而是将 4 个 double-double 值打包到一个 AVX 高寄存器和一个 AVX 低寄存器中。

对于 Mandelbrot 集，我真正需要的是双倍乘法和加法。在那篇论文中，这些是df64_add和df64_mult函数。下图显示了我df64_mult的软件 FMA（左）和硬件 FMA（右）功能的程序集。这清楚地表明，硬件 FMA 是对双倍乘法的一大改进。

那么硬件FMA在双双Mandelbrot集计算中表现如何呢？答案是这仅比使用软件 FMA 快 15%。这比我希望的要少得多。双双 Mandelbrot 计算需要 4 次双双加法和四次双双乘法（x*x、y*y、x*y和2*(x*y)）。但是，2*(x*y)对于 double-double，乘法是微不足道的，因此可以在成本中忽略这种乘法。因此，我认为使用硬件 FMA 的改进如此之小的原因是计算以缓慢的双双加法为主（见下面的汇编）。

过去，乘法比加法慢（程序员使用了几种技巧来避免乘法），但在 Haswell 中，情况似乎正好相反。不仅因为 FMA，还因为乘法可以使用两个端口，而加法只能使用一个。

所以我的问题（最后）是：

1. 当加法比乘法慢时如何优化？
2. 有没有一种代数方法可以改变我的算法以使用更多的乘法和更少的加法？我知道有一些方法可以做相反的事情，例如(x+y)*(x+y) - (x*x+y*y) = 2*x*y使用两个加法来减少一个乘法。
3. 有没有办法简单地使用 df64_add 函数（例如使用 FMA）？

如果有人想知道 double-double 方法比 double 慢十倍左右。我认为这还不错，就好像有一个硬件四精度类型一样，它的速度可能至少是双精度类型的两倍，所以我的软件方法比我对硬件的预期慢五倍（如果它存在的话）。

df64_add部件

我在维基百科中看到，实现四精度的方法是使用双精度运算，即使它的位精度不完全相同：https ://en.wikipedia.org/wiki/Quadruple-precision_floating-point_format

在这种情况下，我们使用两个双精度来存储值。因此，我们进行了两次运算来计算结果，每个运算结果的两倍。

在这种情况下，我们可以在每个双精度数上出现舍入错误，或者它们是避免这种情况的机制？

存在双精度浮点不够充分的工作负载，因此需要四精度。这很少在硬件中提供，因此一种解决方法是使用 double-double，其中 128 位数字由一对 64 位数字表示。这不是真正的 IEEE-754 四倍精度 - 一方面，您不会得到任何额外的指数位 - 但在许多方面都足够接近，并且比纯软件实现要快得多。

许多计算机提供向量浮点运算，最好将它们用于双双运算。这可能吗？特别是，在https://github.com/JuliaMath/DoubleDouble.jl/blob/master/src/DoubleDouble.jl查看双双的实现，在我看来，每个算术运算都需要至少一个条件分支中间，我认为这意味着不能使用 SIMD 矢量运算，除非我遗漏了什么？

将精度提高到双精度之外的一种方法（例如，如果我的应用程序正在做一些与空间相关的事情，需要在许多光年的距离上表示准确的位置）是使用double-double，这是一种由两个双精度组成的结构表示两者之和。对于这种结构上的各种算术运算，算法是已知的，例如双双 + 双双、双 × 双双等，例如本文中给出的。

（请注意，这与 IEEE 754-2008 binary128 的格式不同，也称为四精度，并且不保证往返双精度和 binary128 的转换。）

将这种数量表示为字符串的一种明显方法是使用表示双精度的每个单独组件的字符串，例如“1.0+1.0e-200”。我的问题是，是否有一种已知的方法来转换表示值作为单个小数的字符串？即给定字符串“0.3”，然后提供最接近此表示的双双，或者以相反的方向进行。一种天真的方法是使用连续的乘法/除法 10，但这对于双精度数来说是不够的，所以我有点怀疑它们是否能在这里工作。

我正在玩在 qemu 中模拟的 PPC64 虚拟机，试图模仿 POWER8 CPU。

在这里，该long double类型与 x86 中用于长双精度的 80 位浮点数不同，并且据我所知，它也不符合 IEEE754 的 float128，因为根据 C 宏它有一个 106 位的尾数LDBL_MANT_DIG（对比。 IEEE754 为其 float128 规定的 112 位尾数）。

维基百科说 IEEE754 float128 的机器 epsilon 应该在 1.93e-34 左右，比 80 位 x86 浮点数 (1.08e-19) 要好得多。

然而，当我尝试在这个虚拟机中获取机器 epsilon 时，我得到了一个相当令人惊讶的答案：

它输出以下内容：

我从LDBL_EPSILON和 从得到相同的结果std::numeric_limits<long double>::epsilon()。

这将使它比预期的精确度高出大约 10 倍，逻辑告诉我这应该是不可能的。看到尾数正好是 2x53（IEEE754 的 float64 的），我认为它可能使用双双结构，维基百科还说，它对小数字的精度应该低于 IEEE754 float128。

这里发生了什么？

我在 OpenGL 中实现了一个 mandelbrot-explorer。对于更深的缩放，我正在尝试实现双双算术。我在 C++ 中有工作代码，但是 - 如果我没有忽略某些东西 - 完全相同的代码只会在 GLSL 中产生双精度。

乘法的 C++ 代码：

开放 GL 代码：

出于测试目的，我在 C++ 中都这样做了

...和OpenGL

如果我继续对数字进行平方，在 C++ 中，你会看到它变得越来越大（最初它是太多的零）。GLSL 总是产生一个（或零 - 我在输出之前做了 number-1）。

那么 GLSL 中的双打（或可能）有什么奇怪的地方，dvec2还是我忽略了一些东西？另外：当计算带有双打的 mandelbrot-set 时，CPU 真的很慢。GPU 产生了完全正确的图片（直到我放大到大约 10^-14）并且速度非常快......