Synopsis

#include <boost/math/special_functions/gamma.hpp>

namespace boost{ namespace math{
template <class T>
calculated-result-type lgamma(T z);
template <class T, class Policy>
calculated-result-type lgamma(T z, const Policy&);
template <class T>
calculated-result-type lgamma(T z, int* sign);
template <class T, class Policy>
calculated-result-type lgamma(T z, int* sign, const Policy&);
}} // namespaces

Description

Функцияlgammaопределяется:

Вторая форма функции принимает указатель на целое число, которое, если ненулевая, устанавливается на выходе на знак tgamma(z).

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

Эффективно существует две версии этой функции: полностью общая версия, которая является медленной, но достаточно точной, и гораздо более эффективное приближение, которое используется там, где число цифр в значении T соответствует определенному приближениюЛанчоса. На практике любой встроенный тип плавающей точки, с которым вы столкнетесь, имеет подходящее значение.Приближение Ланчосаопределено для него. Также возможно, учитывая достаточное машинное время, генерировать дальнейшее приближениеLanczosс использованием программы libs/math/tools/lanczos_generator.cpp.

Тип возврата этих функций вычисляется с помощьюправил расчета типа результата: Результатом является тип<double>, если T является целым типом, или тип T иначе.

Accuracy

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

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

Testing

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

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

Implementation

Общая версия этой функции реализована с использованием приближения Стерлинга для больших аргументов:

Для небольших аргументов используется логарифм тгаммы.

Для отрицательногоzиспользуется логарифмический вариант формулы отражения:

Для типов известной точности используется аппроксимацияЛанчоса, класс признаков<boost::math::lanczos::lanczos_traits>отображает тип Т в соответствующее приближение. Логарифмическая версия приближения Ланчоса:

Где L_e,g& #160; это сумма Ланчоса, масштабированная на e^g.

Как и раньше, формула отражения используется дляz< 0.

Когда z очень близко к 1 или 2, тогда логарифмическая версияПриближение Ланчосаочень сильно страдает от ошибки отмены: действительно, для значений, достаточно близких к 1 или 2, могут быть получены произвольно большие относительные ошибки (хотя абсолютная ошибка крошечная).

Для типов с точностью до 113 бит (до и включая 128-битные длинные двойники) используются корнесохраняющие рациональные приближения, разработанные JM, через интервалы [1,2] и [2,3]. В интервале [2,3] используется форма приближения:

lgamma(z) = (z-2)(z+1)(Y + R(z-2));

Где Y является константой, а R(z-2) является рациональным приближением: оптимизировано так, что абсолютная ошибка крошечна по сравнению с Y. Кроме того, малые значения z, превышающие 3, могут быть обработаны путем уменьшения аргумента с использованием отношения повторения:

lgamma(z+1) = log(z) + lgamma(z);

В интервале [1,2] необходимо использовать два приближения, одно для малых z использует:

lgamma(z) = (z-1)(z-2)(Y + R(z-1));

Опять же, Y является константой, а R(z-1) оптимизирован для низкой абсолютной ошибки по сравнению с Y. Для z >1.5 вышеприведенная форма не сходится с решением minimax, но эта аналогичная форма:

lgamma(z) = (2-z)(1-z)(Y + R(2-z));

Наконец, для z< 1 отношение повторения может быть использовано для перехода к z >1:

lgamma(z) = lgamma(z+1) - log(z);

Обратите внимание, что, хотя это включает вычитание, оно, по-видимому, не страдает от ошибки отмены: по мере уменьшения z с 1 термин<-log(z)>становится положительным гораздо быстрее, чем термин<lgamma(z+1)>становится отрицательным. Так что в данном конкретном случае значительные цифры сохраняются, а не отменяются.

Для других типов, которые имеют приближениеLanczos, определенное для них, текущее решение выглядит следующим образом: представьте, что мы уравновешиваем два термина в приближенииLanczos, разделяя термин мощности на его значение приz = 1, а затем умножаем коэффициенты Lanczos на то же значение. Теперь каждый термин будет принимать значение 1 вz = 1, и мы можем переставить условия мощности в термины log1p. Аналогично, если вычесть 1 из части суммы Ланчоса (алгебраически, вычитая значение каждого термина приz = 1), мы получим новое суммирование, которое также можно подавать в log1p. Важно отметить, что все термины имеют тенденцию к нулю, какz - >1:

C_k& #160; термины, приведенные выше, такие же, как и в приближении Ланчоса.

Аналогичная перегруппировка может быть выполнена приz = 2: