Haskell的有限差异,或如何禁用潜在的优化

Ale*_* C. 3 optimization haskell rounding-error numerical-methods

我想实现以下天真(一阶)有限差分函数:

finite_difference :: Fractional a => a -> (a -> a) -> a -> a
finite_difference h f x = ((f $ x + h) - (f x)) / h
Run Code Online (Sandbox Code Playgroud)

正如您可能知道的那样,存在一个微妙的问题:必须确保(x + h)并且可以x用完全可表示的数字来区分.否则,结果有一个巨大的错误,由于(f $ x + h) - (f x)涉及灾难性取消的事实(并且必须仔细选择h,但这不是我的问题).

在C或C++中,问题可以像这样解决:

volatile double temp = x + h;
h = temp - x;
Run Code Online (Sandbox Code Playgroud)

并且volatile修饰符禁用与变量有关的任何优化temp,因此我们确信"聪明"的编译器不会优化这两行.

我还不知道Haskell还不知道如何解决这个问题.恐怕

let temp = x + h
    hh = temp - x 
in ((f $ x + hh) - (f x)) / h
Run Code Online (Sandbox Code Playgroud)

将被Haskell(或它使用的后端)优化掉.我如何得到相当于volatile这里(如果可能的话,不牺牲懒惰)?我不介意GHC的具体答案.

Joh*_*n L 6

我有两个解决方案和建议:

第一个解决方案:您可以保证使用两个辅助函数和NOINLINE编译指示不会优化它:

norm1 x h = x+h
{-# NOINLINE norm1 #-}

norm2 x tmp = tmp-x
{-# NOINLINE norm2 #-}

normh x h = norm2 x (norm1 x h)
Run Code Online (Sandbox Code Playgroud)

这将有效,但会引入一小笔费用.

第二种解决方案:使用volatile在C中写入规范化函数,并通过FFI调用它.性能损失将是最小的.

现在提出建议:目前数学没有优化,所以它目前可以正常工作.你担心它会在未来的编译器中崩溃.我认为这不太可能,但也不是不太可能我也不想防范它.因此,请编写一些涵盖相关案例的单元测试.然后,如果它在将来确实破裂(出于任何原因),你就会确切知道原因.