时间戳

首页 > 解决方案 > 如何在 OpenGL ES 中使用两个浮点数模拟双精度？

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

问题描述

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

标签： opengl-esfloating-pointglsl

解决方案

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

12345600000 + 8.76542995683848E-4 = 12345600000

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

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

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

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等等...基于此，我设法撰写了以下内容：

//---------------------------------------------------------------------------
//--- High Precision float ver: 1.000 ---------------------------------------
//---------------------------------------------------------------------------
#ifndef _GLSL_HP32
#define _GLSL_HP32
//---------------------------------------------------------------------------
// helper functions (internals)
void hp32_nor(vec3 &c)          // bubble sort c coordinates desc by magnitude
    {
    float x;
    if (fabs(c.x)<fabs(c.y)){ x=c.x; c.x=c.y; c.y=x; }
    if (fabs(c.x)<fabs(c.z)){ x=c.x; c.x=c.z; c.z=x; }
    if (fabs(c.y)<fabs(c.z)){ x=c.y; c.y=c.z; c.z=x; }
    }
void hp32_err(vec3 &c,vec3 &e)  // c+=e; apply rounding error e corection to c
    {
    float q;
    q=c.x; c.x+=e.x; e.x=e.x-(c.x-q);
    q=c.x; c.x+=e.y; e.y=e.y-(c.x-q);
    q=c.x; c.x+=e.z; e.z=e.z-(c.x-q);
    q=c.y; c.y+=e.x; e.x=e.x-(c.y-q);
    q=c.y; c.y+=e.y; e.y=e.y-(c.y-q);
    q=c.y; c.y+=e.z; e.z=e.z-(c.y-q);
    q=c.z; c.z+=e.x; e.x=e.x-(c.z-q);
    q=c.z; c.z+=e.y; e.y=e.y-(c.z-q);
    q=c.z; c.z+=e.z; e.z=e.z-(c.z-q);
    hp32_nor(c);
    }
void hp32_split(vec3 &h,vec3 &l,vec3 a) // (h+l)=a; split mantissas to half
    {
    const float n=8193.0; // 10000000000001 bin uses ~half of mantissa bits
    h=a*n;  // this shifts the a left by half of mantissa (almost no rounding yet)
    l=h-a;  // this will round half of mantissa as h,a have half of mantisa bits exponent difference
    h-=l;   // this will get rid of the `n*` part from number leaving just high half of mantisa from original a
    l=a-h;  // this is just the differenc ebetween original a and h ... so lower half of mantisa beware might change sign
    }
//---------------------------------------------------------------------------
// double api (comment it out if double not present)
vec3 hp32_set(double a)                                             // double -> vec2
    {
    vec3 c;
    c.x=a; a-=c.x;
    c.y=a; a-=c.y;
    c.z=a; hp32_nor(c);
    return c;
    }
double hp32_getl(vec3 a){ double c; c=a.z+a.y; c+=a.x; return c; }  // vec2 -> double
//---------------------------------------------------------------------------
// normal api
vec3 hp32_set(float a){ return vec3(a,0.0,0.0); }                   // float -> vec2
float hp32_get(vec3 a){ float c; c=a.z+a.y; c+=a.x; return c; }     // vec2 -> float
vec3 hp32_add(vec3 a,vec3 b)                                        // = a+b
    {
    // c=a+b; addition
    vec3 c=a+b,e; float q;
    // e=(a+b)-c; rounding error
    c.x=a.x+b.x; e.x=c.x-a.x; e.x-=b.x;
    c.y=a.y+b.y; e.y=c.y-a.y; e.y-=b.y;
    c.z=a.z+b.z; e.z=c.z-a.z; e.z-=b.z;
    e=-e; hp32_err(c,e);
    return c;
    }
vec3 hp32_sub(vec3 a,vec3 b)                                        // = a-b
    {
    // c=a-b; substraction
    vec3 c=a-b,e; float q;
    // e=(a-b)-c; rounding error
    c.x=a.x+b.x; e.x=c.x-a.x; e.x+=b.x;
    c.y=a.y+b.y; e.y=c.y-a.y; e.y+=b.y;
    c.z=a.z+b.z; e.z=c.z-a.z; e.z+=b.z;
    e=-e; hp32_err(c,e);
    return c;
    }
vec3 hp32_mul_half(vec3 a,vec3 b)                                   // = a*b where a,b are just half of mantissas !!! internal call do not use this !!!
    {
    //  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)
    vec3 c,e,f; float q,r;
    // c=a*b; (e,f)=(a*b)-c; multiplication
    c.x=(a.x*b.x);
    r=(a.x*b.y); q=c.x; c.x+=r; e.x=r-(c.x-q);
    r=(a.y*b.x); q=c.x; c.x+=r; e.y=r-(c.x-q);
    c.y=(a.x*b.z);
    r=(a.z*b.x); q=c.y; c.y+=r; e.z=r-(c.y-q);
    r=(a.y*b.y); q=c.y; c.y+=r; f.x=r-(c.y-q);
    c.z=(a.y*b.z);
    r=(a.z*b.y); q=c.z; c.z+=r; f.y=r-(c.z-q);
    r=(a.z*b.z); q=c.z; c.z+=r; f.z=r-(c.z-q);
    e=+hp32_add(e,f); hp32_err(c,e);
    return c;
    }
vec3 hp32_mul(vec3 a,vec3 b)                                        // = a*b
    {
    vec3 ah,al,bh,bl,c;
    // split operands to halves of mantissa
    hp32_split(ah,al,a);
    hp32_split(bh,bl,b);
    //  c = (ah+al)*(bh+bl) = ah*bh + ah*bl + al*bh + al*bl
    c=           hp32_mul_half(ah,bh) ;
    c=hp32_add(c,hp32_mul_half(ah,bl));
    c=hp32_add(c,hp32_mul_half(al,bh));
    c=hp32_add(c,hp32_mul_half(al,bl));
    return c;
    }
//---------------------------------------------------------------------------
#endif
//---------------------------------------------------------------------------

现在我只在 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在浮点数上，尾数的结果位宽作为操作数位宽的总和。所以为了避免四舍五入，我们需要首先将操作数分成两半。可以这样做（从您链接的库中分析）：

// split float into h,l
float a,h,l,n;
n=8193.0;   // n =               10000000000001.00000000000000000000b
a=123.4567; // a =                      1111011.01110100111010101000b
h=a*n;      // h =         11110110111100011000.11000000000000000000b
l=h-a;      // l =         11110110111010011101.01010000000000000000b
h-=l;       // h =                      1111011.01110000000000000000b
l=a-h;      // l =                            0.00000100111010101000b

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

现在减半的乘法c=a*b如下所示：

c = (ah+al)*(bh+bl) = ah*bh + ah*bl + al*bh + al*bl

每个半乘法都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。重要的是按大小对术语进行排序，以尽可能防止舍入错误。沿着术语的总和，我只是像加法一样累积错误。

推荐阅读