如何在 OpenGL ES 中使用两个浮点数模拟双精度？

Question

我正在致力于创建 Mandelbrot 集的深度缩放，您可能知道 OpenGL ES 不支持该double 数据类型。它可以提供的最高精度是 IEEE 754 的精度float。在谷歌搜索上，经过大量搜索后，我发现了这个博客： https: //blog.cyclemap.link/2011-06-09-glsl-part2-emu/，它完全致力于这个主题。但是，不幸的是我无法理解那里提供的加法、减法和乘法的代码。特别是，我无法理解处理纠错和携带的部分。如果您能向我解释错误检查的深度以及从低位部分到较高部分的位转移，那将非常有帮助。到目前为止，我只了解将双精度数拆分为两个浮点数的基本概念。但是，我不清楚基本操作的实施。如果使用二进制数的上下文来完成解释，将会非常有帮助。

Answer 1

首先是用于处理此问题的基本原理。一旦您添加或减去具有高指数差的数字，结果就会四舍五入：

12345600000 + 8.76542995683848E-4 = 12345600000

现在，正如我在这里向您展示的那样，我们可以将数字存储为更多浮点数的总和，例如，vec2, vec3, vec4它们仍然是浮点数，但在一起可以组合成更大的整体尾数位宽。您问题中的链接没有像我一样使用指数范围，而是使用舍入结果和非舍入结果之间的差异。然而，链接库使用的 justvec2仍然比本机64位精确得多double，因为尾数double有53位并且float只有24位，所以24+24 = 48 < 53这就是我决定使用vec3. 现在的技巧是获得舍入误差。对于与上面相同的示例：

a=12345600000 
b=8.76542995683848E-4
c=a+b=12345600000

是运算a,b的操作float数+，c是四舍五入的结果。e所以可以这样得到差异：

e=c-a; // e= 0
e-=b;  // e=-8.76542995683848E-4
e=-e;  // e=+8.76542995683848E-4

为了得到不四舍五入的结果e应该添加到哪里。c

现在，如果我们在 的每个组件中存储数字的某些部分vec3，那么我们可以尝试将其添加e到所有组件中（始终从 中删除添加的部分e）直到e为零。

因此，如果c.x+e确实是圆形，我们将其添加到c.y等等...基于此，我设法编写了此内容：

12345600000 + 8.76542995683848E-4 = 12345600000

目前我仅在 CPU 端（C++）进行了测试。为了在 GLSL 中使用它，只需注释掉或删除我用来验证准确性的双 api 函数。并将 更改fabs为abs或声明：

float fabs(float x){ return abs(x); }

我再次有一些归一化函数hp32_nor，它按大小对组件进行排序，以便fabs(x)>=fabs(y)>=fabs(z)返回float乘法所需的组件。+,-不需要它。

这hp32_err就像上面描述的正常数和舍入误差差之间的加法（有点可怕，但它尽可能地保留了准确性）。

我还没有对此进行广泛的测试！看起来+,-,*操作比较精确double。

乘法实现有点复杂，因为a*b在浮点数上，尾数的结果位宽作为操作数位宽的总和。因此，为了避免舍入，我们需要首先将操作数分成两半。可以这样完成（从您链接的库中分析）：

a=12345600000 
b=8.76542995683848E-4
c=a+b=12345600000

所以 float 有 24 位尾数，8193(24/2)+1=13位宽。因此，一旦将任何浮点数与其相乘，结果就需要比现在多大约一半的尾数位并四舍五入。这只是返回到原始操作数比例并获取另一半作为新一半和原始浮点值之间的差异的问题。所有这些都是在辅助函数中完成的hp32_split。

现在一半的乘法c=a*b看起来像这样：

e=c-a; // e= 0
e-=b;  // e=-8.76542995683848E-4
e=-e;  // e=+8.76542995683848E-4

每个半乘法a?*b?都是这样的：

c = (a.x+a.y+a.z)*(b.x+b.y+b.z)     // multiplication of 2 expresions
c = (a.x*b.x)+(a.x*b.y)+(a.x*b.z)   // expanded
   +(a.y*b.x)+(a.y*b.y)+(a.y*b.z)
   +(a.z*b.x)+(a.z*b.y)+(a.z*b.z)
c = (a.x*b.x)                       // ordered desc by magnitude (x>=y>=z)
   +(a.x*b.y)+(a.y*b.x)
   +(a.x*b.z)+(a.z*b.x)+(a.y*b.y)
   +(a.y*b.z)+(a.z*b.y)
   +(a.z*b.z)

所以我将 的每个组成部分分为 3 个括号项c。重要的是按大小对项进行排序，以尽可能防止舍入错误。在对项求和的过程中，我只是像加法一样累积误差。