Synopsis

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

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

Description

Description Дигамма определяется как логарифмическая производная гамма-функции:

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

Тип возврата этой функции вычисляется с помощью правил расчета Результат : результат по типу ДП, когда T является целым типом, и тип T в противном случае.

Accuracy

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

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

Testing

Существует два набора тестов: спотовые значения вычисляются с помощью онлайн-калькулятора на functions.wolfram.com, в то время как случайные тестовые значения генерируются с использованием высокоточной контрольной реализации (дифференциальная Lanczos approximation см. ниже).

Implementation

Реализация подразделяется на следующие области:

Для отрицательных аргументов формула отражения:

digamma(1-x) = digamma(x) + pi/tan(pi*x);

используется для того, чтобы сделать x положительным.

Для аргументов в диапазоне [0,1] отношение повторения:

digamma(x) = digamma(x+1) - 1/x

используется для переноса оценки на [1,2].

Для аргументов в диапазоне [1,2] используется рациональное приближение , указанное JM (см. ниже).

Для аргументов в диапазоне [2,BIG] отношение рецидива:

digamma(x+1) = digamma(x) + 1/x;

используется для переноса оценки на диапазон [1,2].

Для аргументов > BIG асимптотическое расширение:

можно использовать. Однако это расширение расходится после нескольких терминов: точно, сколько терминов зависит от размера x. Поэтому значение BIG должно быть выбрано таким образом, чтобы серия могла быть усечена в термине, который слишком мал, чтобы оказать какое-либо влияние на результат при оценке по BIG. Выбор BIG=10 для 80-битных реалов, а BIG=20 для 128-битных реалов позволяет усечь серию после достаточно небольшого количества терминов и оценить как многочлен в 1/(x<315>*<317>xx).

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

Рациональное приближение , указанное JM в диапазоне [1,2], происходит следующим образом.

Сначала высокоточное приближение к дигамме было построено с использованием 60-терминной дифференцированной Lanczos approximation, используемой формой:

В тех случаях, когда P(x) и Q(x) являются полиномиями из рациональной формы суммы Lanczos, а P'(x) и Q'(x) являются их первыми производными. Часть Ланзоса этого приближения имеет теоретическую точность ~100 десятичных цифр. Тем не менее, отмена в приведенной выше сумме сократит это примерно до 99-(1/y) цифры, если результат y. Это приближение было использовано для расчета положительного корня дигаммы, и было установлено, что оно согласуется с значением, используемым Коди, до 25 цифр (см. Math. Comp. 27, 123-127 (1973) Коди, Strecok и Thacher) и значением, используемым Моррисом, до 35 цифр (см. TOMS Algorithm 708).

Точно так же несколько спотовых тестов согласовывались с значениями, рассчитанными с использованием функций. wolfram.com to >40 digits. Это достаточно точно, чтобы убедиться, что приближение ниже является точным, чтобы удвоить точность. Достижение 128-битной двойной точности требует, чтобы местоположение корня было известно ~70 цифр, и не ясно, соответствует ли значение, рассчитанное этим методом, этому требованию: сложность заключается в независимой проверке полученного значения.

Рациональное приближение , указанное JM, было оптимизировано для абсолютной ошибки с использованием формы:

digamma(x) = (x - X0)(Y + R(x - 1));

Там, где X0 является положительным корнем дигаммы, Y является константой, а R(x - 1) является рациональным приближением. Обратите внимание, что, поскольку X0 иррационален, нам нужно в два раза больше цифр в X0, чем в x, чтобы избежать ошибки отмены во время вычитания (это предполагает, что x является точным значением, если это не все ставки отключены). Это означает, что даже когда x является значением корня, округленного до ближайшего репрезентативного значения, результат дигаммы (x) не будет нулевым.