Мне нужно очень быстро рассчитать приближение Стирлинга

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

const double AL[8] =
{ 0.0, 0.0, 0.6931471806, 1.791759469, 3.178053830, 4.787491743,
    6.579251212, 8.525161361 };
void HGD::mpfr_afc(mpfr_t &ret, const mpfr_t &val){

    if(mpfr_cmp_ui(val, 7) <= 0){
        mpfr_set_d(ret, AL[mpfr_get_ui(val, MPFR_RNDN)], MPFR_RNDN);
    }else{
        mpfr_set(ret, val, MPFR_RNDN);
        mpfr_add_d(ret, ret, 0.5, MPFR_RNDN);
        mpfr_log(LV, val, MPFR_RNDN);
        mpfr_mul(ret, LV, ret, MPFR_RNDN);
        mpfr_sub(ret, ret, val, MPFR_RNDN);
        mpfr_add_d(ret, ret, 0.399089934, MPFR_RNDN);
    }
}

У меня есть пара разных идей:

  • Предварительно вычислите больше, чем первые 8 входных параметров функции.
  • Оптимизируйте математику (используйте более грубое приближение для меньшей точности)
  • Используйте несколько потоков для параллельных вычислений на разных входных данных
  • Переключитесь на собственную арифметику, когда числа могут соответствовать типам машинных данных.

Есть ли другие подходы, которые я мог бы использовать?


c c++ mpfr arbitrary-precision numerical
person pg1989    schedule 12.02.2014    source источник
comment
Я заранее извиняюсь, если это наивно, но если вы разберете свое профилирование дальше, какая часть (части) этой функции займет больше всего времени?   -  person Russ Clarke    schedule 12.02.2014
comment
Вычисление журнала определенно доминирует.   -  person pg1989    schedule 12.02.2014
comment
Я с подозрением отношусь к вашим константам, так как для double я ожидаю 15-18 цифр точности. Но тогда, возможно, конечный результат не требует ответа высокой точности.   -  person chux - Reinstate Monica    schedule 12.02.2014
comment
как часто вы вызываете эту функцию с одним и тем же значением val?   -  person Glenn Teitelbaum    schedule 12.02.2014
comment
Может быть, один раз из десятка обращений к функции в среднем.   -  person pg1989    schedule 12.02.2014
comment
Вызывать mpfr_afc() реже? С какими размерами он работает?   -  person brian beuning    schedule 12.02.2014
comment
От 64 до 256 бит.   -  person pg1989    schedule 12.02.2014
comment
Каков диапазон и распространение val? Как часто val меньше 8?   -  person pat    schedule 12.02.2014
comment
Мне пришлось посмотреть на scipy реализацию формулы Стирлинга, чтобы понять проблему: github.com/scipy/scipy/blob/   -  person iljau    schedule 12.02.2014
comment
val в целом меньше 8 очень редко, за исключением случаев, когда входные данные для функции выборки сами по себе малы.   -  person pg1989    schedule 13.02.2014

Ответы (2)


arrow_upward
2
arrow_downward

Переключитесь на собственную арифметику, когда числа могут соответствовать типам машинных данных.

Это будет моя первая попытка. MPFR, вероятно, станет убийцей производительности.

Мне кажется, вы хотите вычислить логарифм n! которую вы уже аппроксимируете формулой Стирлинга.

Обратите внимание, что n!=Гамма(n+1). Существуют (на первый взгляд) высокооптимизированные функции для вычисления как гамма-функции, так и ее логарифма. Например:

Я бы свернул свое собственное более грубое приближение только в том случае, если все вышеперечисленное не работает с точки зрения производительности.

person Ali    schedule 12.02.2014

arrow_upward
1
arrow_downward

Пара мыслей здесь. Во-первых, мне кажется, что использование MPFR для этого может быть излишним. Любая из библиотек с высокой точностью имеет огромные накладные расходы. Не просто большие накладные расходы, а огромные накладные расходы. Вторая мысль заключается в том, что, возможно, вам не нужно использовать функцию многоточечного журнала. Может быть, вы могли бы обойтись стандартным журналом?

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

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

person cassius    schedule 12.02.2014