计算科学 - 在损坏的编译器上计算 double 的保守圆角平方 - 吾爱随笔录 - 问答

在损坏的编译器上计算 double 的保守圆角平方

计算科学区间算术

2021-12-22 00:39:39

我正在尝试计算 $x^2$ 在哪里 $x$ 是一个双精度浮点数。我需要保守的四舍五入，这意味着我需要 $x \cdot x$ 四舍五入和 $x \cdot x$ 四舍五入。如果舍入模式设置为 FE_UPWARDS，则代码为

double lo = -(x * -x);
double hi = x * x;

但是，这在 clang 3.4 上对我不起作用：编译器注意到“公共子表达式”并且只做一个乘法。有谁知道如何说服 clang 不要这样做，理想情况下不会损害其他优化？特别是，我想避免使用非内联函数。

3个回答

我相信这可以通过将区间边界打包到 SSE 寄存器中并使用 SSE 执行所有区间操作来解决，如

http://hal.inria.fr/docs/00/28/84/56/PDF/intervals-sse2-long-paper.pdf

无论如何，这应该比我当前的代码快，并且 clang 不应该在同一指令中应用 SSE。

如果您的编译器有nextafter，那么您可以执行以下操作，无论设置哪种舍入模式都可以：

double x2 = x*x;
double lo = nextafter(x2,-1);
double hi = nextafter(x2,x2+1);

但我不知道是否nextafter是内联的（可能不是）。

这是一个悲伤的答案，所以我不会接受它，希望其他人有更好的解决方案：

double hi = x * x;
double lo = (2*epsilon-1)*hi;

这是安全的，因为如果我们四舍五入，

\begin{aligned} - ((2 ϵ - 1) ⊙ (x ⊙ x)) & \leq (1 + ϵ)^{2} (1 - 2 ϵ) x^{2} \\ \leq (1 + 2 ϵ + ϵ^{2}) (1 - 2 ϵ) x^{2} \\ \leq (1 - ϵ^{2} - ϵ^{3}) x^{2} \\ \leq x^{2} \end{aligned}

$\begin{aligned} -\left((2\epsilon-1) \odot (x \odot x)\right) &\le (1 + \epsilon)^2(1-2\epsilon) x^2 \\ &\le (1+2\epsilon+\epsilon^2)(1-2\epsilon)x^2 \\ &\le (1-\epsilon^2-\epsilon^3)x^2 \\ &\le x^2 \end{aligned}$ 所以我们的代码正确地计算了下限。

其它你可能感兴趣的问题

上一篇MATLAB 的开源自动微分？下一篇IBM LSF 是否促进并行计算？