1

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

乘法的 C++ 代码：

igdd quickTwoSum(double a, double b)
    {
        igdd Temp;
        Temp.hi = a + b;
        Temp.lo = b - (Temp.hi - a);
        return Temp;
    }

igdd Split64(double d)
    {
        const double SPLITTER = (1 << 29) + 1;
        double t = d * SPLITTER;
        igdd Temp;
        Temp.hi = t - (t - d);
        Temp.lo = d - Temp.hi;
        return Temp;
    }

igdd twoProd(double a, double b)
    {
        igdd p;
        p.hi = a * b;
        igdd aS = Split64(a);
        igdd bS = Split64(b);
        p.lo = ((aS.hi * bS.hi - p.hi) + aS.hi * bS.lo + aS.lo * bS.hi + aS.lo * bS.lo);
        return p;
    }

igdd operator*(igdd a)
    {
        igdd Temp;
        Temp = twoProd(a.hi, this->hi);
        Temp.lo += a.hi * this->lo + a.lo * this->hi;
        Temp = quickTwoSum(Temp.hi, Temp.lo);
        return Temp;
    }

开放 GL 代码：

dvec2 Split64(double d)
{
    const double SPLITTER = (1 << 29) + 1;
    double t = d * SPLITTER;
    dvec2 result;
    result.x = t - (t - d);
    result.y = d - result.x;

    return result;
}

dvec2 quickTwoSum(double a, double b)
{
    dvec2 result;
    result.x = a + b;
    result.y = b - (result.x - a);
    return result;
}

dvec2 twoProd(double a, double b)
{
    dvec2 p;
    p.x = a * b;
    dvec2 aS = Split64(a);
    dvec2 bS = Split64(b);
    p.y = (aS.x * bS.x - p.x) + aS.x * bS.y + aS.y * bS.x + aS.y * bS.y;
    return p;
}

dvec2 df128_mul(dvec2 a, dvec2 b)
{
    dvec2 p;

    p = twoProd(a.x, b.x);
    p.y += a.x * b.y + a.y * b.x;
    p = quickTwoSum(p.x, p.y);
    return p;
}

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

igdd Testdd; // The double-double class in C++
Testdd.lo = 0.0000000000000000000000001;
Testdd.hi = 1.0;
std::cout.precision(17);
for (int i = 0; i < 95; i++)
{
    Testdd = Testdd * Testdd;
    std::cout << std::fixed << Testdd.hi << std::endl;
}

...和OpenGL

dvec2 XR;
XR.y = 0.0000000000000000000000001;
XR.x = 1.0;
for (i = 0; i < 95; i++) // I tried different values for 95 here
{
    XR = df128_mul(XR, XR);
}

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

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

4

1 回答 1

2

我终于找到了问题所在。编译器优化了所需的精度。我发现了另一个产生完全相同结果的乘法代码（双精度 - 不是“双双”）。出于某种原因，该代码似乎对这个人有效，但对我无效。

解决问题的方法是使用“精确”关键字。我用它涂抹了所有东西，现在它可以工作了 - 虽然当然要慢得多：

precise dvec2 quickTwoSum(precise double a, precise double b)
{
    precise dvec2 result;
    result.x = a + b;
    result.y = b - (result.x - a);
    return result;
}

precise dvec2 twoProd(precise double a, precise double b)
{
    precise dvec2 p;
    p.x = a * b;
    precise dvec2 aS = Split64(a);
    precise dvec2 bS = Split64(b);
    p.y = (aS.x * bS.x - p.x) + aS.x * bS.y + aS.y * bS.x + aS.y * bS.y;
    return p;
}

precise dvec2 df128_mul(precise dvec2 a, precise dvec2 b)
{
    precise dvec2 p;
    p = twoProd(a.x, b.x);
    p.y += a.x * b.y + a.y * b.x;
    p = quickTwoSum(p.x, p.y);

    return p;
}
于 2022-01-01T18:18:50.987 回答