3

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

finite_difference :: Fractional a => a -> (a -> a) -> a -> a
finite_difference h f x = ((f $ x + h) - (f x)) / h

您可能知道,有一个微妙的问题:必须确保这一点,(x + h)并且x相差一个完全可表示的数字。否则,由于涉及灾难性取消这一事实,结果会产生巨大的错误(f $ x + h) - (f x)(并且必须仔细选择h,但这不是我的问题)。

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

volatile double temp = x + h;
h = temp - x;

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

我对 Haskell 的了解还不够,不知道如何解决这个问题。恐怕

let temp = x + h
    hh = temp - x 
in ((f $ x + hh) - (f x)) / h

将被 Haskell(或它使用的后端)优化掉。我怎样才能得到volatile这里的等价物(如果可能的话不牺牲懒惰)?我不介意 GHC 的具体答案。

4

3 回答 3

6

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

第一个解决方案:您可以保证不会使用两个辅助函数和 NOINLINE pragma 进行优化:

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

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

normh x h = norm2 x (norm1 x h)

这将起作用,但会产生少量成本。

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

现在提出建议:目前数学尚未优化,因此目前可以正常工作。您担心它会在未来的编译器中中断。我认为这不太可能,但也不是不太可能,我也不想防范它。因此,编写一些涵盖相关案例的单元测试。然后,如果它确实在未来发生故障(出于任何原因),您将确切知道原因。

于 2011-04-12T18:45:24.580 回答
2

一种方法是查看核心。

专注于Doubles(最有可能触发一些优化的情况):

finite_difference :: Double -> (Double -> Double) -> Double -> Double
finite_difference h f x = ((f $ x + hh) - (f x)) / h
   where
        temp = x + h
        hh   = temp - x 

编译为:

A.$wfinite_difference h f x =
    case f (case x of
                  D# x' -> D# (+## x' (-## (+## x' h) x'))
           ) of 
        D# x'' -> case f x of D# y -> /## (-## x'' y) h

同样(甚至更少的重写)对于多态版本。

因此,虽然变量是内联的,但数学并没有被优化掉。除了看核心之外,我想不出一种方法来保证你想要的属性。

于 2011-04-12T18:23:35.897 回答
1

我不这么认为

temp = unsafePerformIO $ return $ x + h

会得到优化。只是一个猜测。

于 2011-04-12T18:20:06.210 回答