C++ 中双双精度扩展精度类的浮点编译器优化问题

Question

我用C++实现了一个“double-double”类。 “双精度”使用两个

doubles

来提高精度。数字被视为

number = hi + lo

，但实际上从未计算出总和，因为

hi + lo == hi

。这是我的代码的缩短版本：

class doubledouble
{
public:
    double hi, lo;

    doubledouble()
    {
        hi = 0.0;
        lo = 0.0;
    }

    doubledouble quickTwoSum(const double a, const double b) const
    {
        int old;
        doubledouble out;
        double s = a + b;
        double e = b - (s - a);
        out.hi = s;
        out.lo = e;
        return (out);
    }

    doubledouble twoSum(const double a, const double b) const
    {
        double s = a + b;
        double v = s - a;
        double e = (a - (s - v)) + (b - v);
        doubledouble out;
        out.hi = s;
        out.lo = e;
        return (out);
    }

    doubledouble operator+(const doubledouble& in) const
    {
        doubledouble s, t;
        s = twoSum(in.hi, this->hi);
        t = twoSum(in.lo, this->lo);
        s.lo += t.hi;
        s = quickTwoSum(s.hi, s.lo);
        s.lo += t.lo;
        s = quickTwoSum(s.hi, s.lo);
        return(s);
    }

    doubledouble split(const double a) const
    {
        const double split = (1 << 12) + 1;
        double t = a * split;
        doubledouble out;
        out.hi = t - (t - a);
        out.lo = a - out.hi;
        return (out);
    }

    doubledouble twoProd(const double a, const double b) const
    {
        double p = a * b;
        doubledouble aS = split(a);
        doubledouble bS = split(b);
        double err = ((aS.hi * bS.hi - p) + aS.hi * bS.lo + aS.lo * bS.hi) + aS.lo * bS.lo;
        aS.hi = p;
        aS.lo = err;
        return (aS);
    }

    doubledouble operator*(const doubledouble& in) const
    {
        doubledouble p;
        p = twoProd(this->hi, in.hi);
        p.lo += this->hi * in.lo + this->lo * in.hi;
        p = quickTwoSum(p.hi, p.lo);
        return(p);
    }
};

我的问题是这段代码在调试模式下工作，但在发布模式下我只能得到正常的

double

精度。我正在使用 MSVC 编译器，如果我在类之前放置一个

#pragma optimize("", off)

并在其后面放置相应的

on

，则代码可以正常工作（当然非常慢）。

我很确定问题出在像

t - (t - a)

这样的东西上，它在数学上简单地等于

a

，所以编译器可能会像这样优化它。为了缩小问题范围，我尝试将

#pragma optimize

放在某些函数周围。但即使我将类中的每个函数都放在编译指示之间，代码也无法正常工作。

无论如何，使用

#pragma optimize

无论如何都不是很有用。该代码速度慢且不可移植。我的第一个问题是：为什么

#pragma optimize

不能在类中的特定函数中工作？这对于缩小问题范围非常有用。

第二个问题：我如何欺骗/说服编译器不优化某些表达式，如

t - (t - a)

？使用临时变量可能不起作用，因为编译器会对其进行优化。

编辑

作为可重现的示例，我使用了以下代码：

double a = 1.23456789012345;
const double split = (1 << 12) + 1;
double t = a * split;
double hi = t - (t - a);
double lo = a - hi;
   
std::cout << "hi: " << hi << ", lo: " << lo << std::endl;

两种情况下的输出（发布模式和调试模式）：

hi: 1.23457, lo: -3.48832e-13

所以

t - (t - a)

似乎不是问题。

至于一个“简单”的例子（如果我找到一个更简单的例子，我会编辑），它会产生错误的行为 - 我计算了缩放比例大于 10^14 的曼德尔布罗特集。在调试模式下我得到了完美的结果，在发布模式下我得到了像素化。我用这个函数来创建图像：

void CalcMandelCPU_doubledouble(std::atomic<unsigned int>* RowCounter)
    {
        std::complex<doubledouble> Z, C;

        int Y = (*RowCounter)++;
        while (Y < ResY)
        {
            for (int X = 0; X < ResX; X++)
            {
                C.real(minX + (X + 0.5) * (maxX - minX) / ResX);
                C.imag(minY + (Y + 0.5) * (maxY - minY) / ResY);

                int i;
                Z = 0;
                for (i = 0; i < MaxIter; i++)
                {
                    Z = Z * Z + C;

                    if (Z.real() * Z.real() + Z.imag() * Z.imag() > 4.0)
                        break;
                }

                if (i == MaxIter)
                {
                    Result[Y * ResX + X] = 0;
                }
                else
                {
                    Result[Y * ResX + X] = i + 1 - log(log(Z.real() * Z.real() + Z.imag() * Z.imag())) / log(2.0);
                }
            }
            Y = (*RowCounter)++;
        }
    }

定义

MaxIter

、

ResX

和

ResY

。

minX

、

maxX

、

minY

和

maxY

将为您提供（缩放）视图。视图变量的示例（如

doubledoubles

）：

minX.hi = -1.9997740601362948

，

minX.lo = 9.2915246622385581e-17

maxX.hi = -1.9997740601362861

，

maxX.lo = 7.5150963188138267e-17

minY.hi = -3.2900430992572130e-09

，

minY.lo = 1.7490722630808694e-25

maxY.hi = -3.2900375437016574e-09

，

maxY.lo = 1.0427104313278710e-25

Answer 1

好的，解决方案如下：

正如多个用户在评论中提到的那样，从

fp:fast

切换到

fp:strict

解决了问题。但是，我不想更改全局优化设置，以免影响我的代码的其余部分。解决方案是在某些变量前加上

volatile

（如果有更好的解决方案，请评论）。事实上只需要两个改变：

功能

twoSum

：

doubledouble twoSum(const double a, const double b) const
    {
        volatile double s = a + b;
        volatile double v = s - a;
        double e = (a - (s - v)) + (b - v);
        doubledouble out;
        out.hi = s;
        out.lo = e;
        return (out);
    }

...和功能

split

:

doubledouble split(const double a) const
    {
        const double split = (1 << 26) + 1;
        volatile double t = a * split;
        doubledouble out;
        out.hi = t - (t - a);
        out.lo = a - out.hi;
        return (out);
    }

谢谢评论者。