Как улучшить производительность этого числового вычисления в Haskell?

Я посреди портирования исходной C реализации David Blei Скрытого распределения Дирихле Haskell, и я пытаюсь решить, оставить ли часть материала низкого уровня в C. Следующая функция является одним примером — это - приближение второй производной lgamma:

double trigamma(double x)
{
    double p;
    int i;

    x=x+6;
    p=1/(x*x);
    p=(((((0.075757575757576*p-0.033333333333333)*p+0.0238095238095238)
         *p-0.033333333333333)*p+0.166666666666667)*p+1)/x+0.5*p;
    for (i=0; i<6 ;i++)
    {
        x=x-1;
        p=1/(x*x)+p;
    }
    return(p);
}

Я перевел это в более или менее идиоматического Haskell следующим образом:

trigamma :: Double -> Double
trigamma x = snd $ last $ take 7 $ iterate next (x' - 1, p')
  where
    x' = x + 6
    p  = 1 / x' ^ 2
    p' = p / 2 + c / x'
    c  = foldr1 (\a b -> (a + b * p)) [1, 1/6, -1/30, 1/42, -1/30, 5/66]
    next (x, p) = (x - 1, 1 / x ^ 2 + p)

Проблема состоит в том, что, когда я выполняю обоих через Критерий, моя версия Haskell в шесть или семь раз медленнее (я компилирую с -O2 на GHC 6.12.1). Некоторые подобные функции еще хуже.

Я практически ничего не знаю о производительности Haskell, и я ужасно не интересуюсь рытьем через Ядро или что-либо как этот, так как я могу всегда просто называть небольшое количество интенсивных математикой функций C через FFI.

Но мне любопытно на предмет того, существует ли низко висящий плод, который я пропускаю — некоторое расширение или библиотека или аннотация, которую я мог использовать для ускорения этого числового материала, не делая это слишком ужасным.

ОБНОВЛЕНИЕ: Вот два лучших решения благодаря Don Stewart и Yitz. Я изменил ответ Yitz немного для использования Data.Vector.

invSq x = 1 / (x * x)
computeP x = (((((5/66*p-1/30)*p+1/42)*p-1/30)*p+1/6)*p+1)/x+0.5*p
  where p = invSq x

trigamma_d :: Double -> Double
trigamma_d x = go 0 (x + 5) $ computeP $ x + 6
  where
    go :: Int -> Double -> Double -> Double
    go !i !x !p
        | i >= 6    = p
        | otherwise = go (i+1) (x-1) (1 / (x*x) + p)

trigamma_y :: Double -> Double
trigamma_y x = V.foldl' (+) (computeP $ x + 6) $ V.map invSq $ V.enumFromN x 6

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

Как camccann сказал в Reddit, мораль истории "Для лучших результатов, используйте Don Stewart в качестве своего генератора кода бэкенда GHC". Запрещая то решение, самая безопасная ставка, кажется, только для перевода управляющих структур C непосредственно в Haskell, хотя сплав цикла может дать подобную производительность в более идиоматическом стиле.

Я, вероятно, закончу тем, что использовал Data.Vector подход в моем коде.