Конечная разница в Haskell, или как отключить потенциальные оптимизации

Я хотел бы реализовать следующую наивную (первого порядка) функцию конечных разностей:

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.


optimization numerical-methods haskell rounding-error
person Alexandre C.    schedule 12.04.2011    source источник
comment
Почему жертвование ленью является проблемой?   -  person Don Stewart    schedule 12.04.2011
comment
@Don: На самом деле, скорее всего, нет. Я планирую использовать это довольно надуманными способами (например, с процедурами экстраполяции Ричардсона) в контекстах, где h будет происходить из ленивых бесконечных списков. Но что касается именно этой функции, вы правы, я не сразу это понял (и пытаюсь познакомиться с языком).   -  person Alexandre C.    schedule 12.04.2011

Ответы (3)


arrow_upward
6
arrow_downward

У меня есть два решения и предложение:

Первое решение: вы можете гарантировать, что это не будет оптимизировано с помощью двух вспомогательных функций и прагмы NOINLINE:

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

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

normh x h = norm2 x (norm1 x h)

Это будет работать, но потребует небольших затрат.

Второе решение: написать функцию нормализации на C с использованием volatile и вызвать ее через FFI. Штраф за производительность будет минимальным.

Теперь предложение: в настоящее время математика не оптимизирована, поэтому в настоящее время она будет работать правильно. Вы боитесь, что это сломается в будущем компиляторе. Я думаю, что это маловероятно, но не настолько маловероятно, чтобы я также не хотел остерегаться этого. Поэтому напишите несколько модульных тестов, которые охватывают рассматриваемые случаи. Тогда, если он сломается в будущем (по какой-либо причине), вы будете точно знать, почему.

person John L    schedule 12.04.2011
comment
Я скорее боюсь, что бэкэнд LLVM оптимизирует математику (конечно, если мне будет предоставлена ​​какая-то быстрая математическая опция, которую я хотел бы включить и отключить только для этих двух строк). Изменчивые переменные в C выполняют свою работу, потому что есть гарантия, что переменная будет считываться из памяти каждый раз, когда она потребуется. Но я не хочу писать функцию C, так как мне нужен тип Fractional a (а не только Double). Итак, написание функции NOINLINE (достаточно только одной) действительно решает мою проблему. - person Alexandre C.; 12.04.2011
comment
NOINLINE работает на стороне Haskell. И стратегия компиляции GHC такова, что LLVM или GCC вероятно не будут видеть привязки для оптимизации математики каким-либо небезопасным способом. - person Don Stewart; 12.04.2011
comment
Если функция NOINLINE является полиморфной и не экспортируется, я уверен, что LLVM или GCC ее не разгадают. Однако мое предложение по модульному тестированию остается в силе. - person John L; 13.04.2011

arrow_upward
2
arrow_downward

Один из способов — посмотреть на Ядро.

Специализируясь на 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

И аналогично (с еще меньшим переписыванием) для полиморфной версии.

Таким образом, хотя переменные встроены, математика не оптимизирована. Помимо просмотра ядра, я не могу придумать способ гарантировать желаемое свойство.

person Don Stewart    schedule 12.04.2011
comment
Да, но у этого есть тот недостаток, что в будущих версиях компилятора могут быть введены более агрессивные оптимизации. Таким образом, черный ящик для предотвращения оптимизации был бы хорош. - person Robin Green; 12.04.2011
comment
Единственный способ «черного ящика» — добавить {-# GHC_OPTIONS -Onot #-} - person Don Stewart; 12.04.2011
comment
Спасибо, но я использую GHC 7.02 и могу поспорить, что серверная часть LLVM будет оптимизировать это (может быть опция придерживаться математики IEEE, но я не хочу активировать ее для всей библиотеки) - person Alexandre C.; 12.04.2011
comment
@Don: -Онот довольно агрессивен. Я просто хочу отключить оптимизацию для двух конкретных строк. - person Alexandre C.; 12.04.2011
comment
Любое преобразование, выполняемое компилятором, которое не учитывает семантику чисел с плавающей запятой, является ошибкой. Вам вообще не нужно ничего делать, чтобы получить правильный код. - person augustss; 13.04.2011
comment
@augustss: у основных компиляторов C есть переключатель для отключения небезопасной математической оптимизации (т.е. не учитывающей семантику FP). Иногда это то, что вы хотите, а здесь это то, чего вы не хотите. Иногда я просто хочу, чтобы компилятор предположил, что не будет ни NaN, ни бесконечностей, ни денормализованных чисел. - person Alexandre C.; 13.04.2011
comment
@Alexandre C: Да, это разумно. У генератора кода LLVM есть возможность включить небезопасные операции FP (при генерации кода, а не оптимизации, насколько я знаю), но я сомневаюсь, что ghc может передать ему такие вещи, как NOINLINE. Так что может вам не повезло. - person augustss; 13.04.2011

arrow_upward
1
arrow_downward

я так не думаю

temp = unsafePerformIO $ return $ x + h

будет оптимизирован. Просто предположение.

person Robin Green    schedule 12.04.2011