Synopsis

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

namespace boost{ namespace math{
template <class T1, class T2>
calculated-result-type gamma_p(T1 a, T2 z);
template <class T1, class T2, class Policy>
calculated-result-type gamma_p(T1 a, T2 z, const Policy&);
template <class T1, class T2>
calculated-result-type gamma_q(T1 a, T2 z);
template <class T1, class T2, class Policy>
calculated-result-type gamma_q(T1 a, T2 z, const Policy&);
template <class T1, class T2>
calculated-result-type tgamma_lower(T1 a, T2 z);
template <class T1, class T2, class Policy>
calculated-result-type tgamma_lower(T1 a, T2 z, const Policy&);
template <class T1, class T2>
calculated-result-type tgamma(T1 a, T2 z);
template <class T1, class T2, class Policy>
calculated-result-type tgamma(T1 a, T2 z, const Policy&);
}} // namespaces

Description

Существует четыренеполных гамма-функции: Две являются нормализованными версиями (также известными какупорядоченныенеполные гамма-функции), которые возвращают значения в диапазоне [0, 1], а две являются ненормализованными и возвращают значения в диапазоне [0, Γ(a)]. Пользователи, заинтересованные в статистических приложениях, должны использоватьнормализованные версии (gamma_p и gamma_q)..

Все эти функции требуютa>0иz>= 0, иначе они возвращают результатdomain_error.

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

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

template <class T1, class T2>
calculated-result-type gamma_p(T1 a, T2 z);
template <class T1, class T2, class Policy>
calculated-result-type gamma_p(T1 a, T2 z, const Policy&);

Возвращает нормализованную низшую неполную гамма-функцию a и z:

Эта функция быстро меняется от 0 до 1 вокруг точки z == a:

template <class T1, class T2>
calculated-result-type gamma_q(T1 a, T2 z);
template <class T1, class T2, class Policy>
calculated-result-type gamma_q(T1 a, T2 z, const Policy&);

Возвращает нормализованную низшую неполную гамма-функцию a и z:

Эта функция быстро меняется от 1 до 0 вокруг точки z == a:

template <class T1, class T2>
calculated-result-type tgamma_lower(T1 a, T2 z);
template <class T1, class T2, class Policy>
calculated-result-type tgamma_lower(T1 a, T2 z, const Policy&);

Возвращает полную (ненормализованную) низшую неполную гамма-функцию a и z:

template <class T1, class T2>
calculated-result-type tgamma(T1 a, T2 z);
template <class T1, class T2, class Policy>
calculated-result-type tgamma(T1 a, T2 z, const Policy&);

Возвращает полную (ненормализованную) верхнюю неполную гамма-функцию a и z:

Accuracy

Следующие таблицы дают пиковые и средние относительные ошибки в различных областях a и z, а также сравнения с библиотекамиGSL-1.9иCephes. Обратите внимание, что только результаты для самого широкого типа плавающей точки в системе приведены, поскольку более узкие типы имеютфактически нулевую ошибку.

Обратите внимание, что ошибки растут по мере того, какувеличивается.

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

Все значения находятся в единицах эпсилона.

Testing

Существует два набора тестов: точечные тесты сравнивают значения, взятые из. Онлайн-оценщик Mathworldс этой реализацией выполняет базовую «проверку безумия». Тесты точности используют данные, генерируемые с очень высокой точностью (используя RR-класс NTL, установленный с точностью 1000 бит), используя эту реализацию с очень высокой точностью 60-кратногоприближения Lanczos, и некоторые, но не все, специальные случаи обработки отключены. Это менее чем удовлетворительно: действительно следует использовать независимый метод, но, по-видимому, имеется полное отсутствие таких методов. Мы даже не можем использовать преднамеренно наивную реализацию без специальной обработки корпуса, поскольку непрерывная фракция Legendre нестабильна для малых a и z.

Implementation

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

1)

2)

3)

Нижняя неполная гамма вычисляется из ее последовательного представления:

4)

Или путем вычитания верхнего интеграла из Γ(a) или 1 когдаx - (1(3x)) >a и x >1.1/.

Верхний интеграл вычисляется из непрерывного представления фракции Legendre:

5)

Когда(x >1.1)или путем вычитания нижнего интеграла из Γ(a) или 1 когдаx - (1(3x))

См. документацию дляpowm1иtgamma1pm1для подробной информации об их осуществлении. Обратите внимание, однако, что точностьtgamma1pm1ограничена либо 35 цифрами, либо приближениемLanczos, связанным с типом T - если есть один - в зависимости от того, какой из двух больше. Это накладывает аналогичный предел на точность этой функции в этом регионе.

Дляx< 1.1точка кроссовера, где результат ~0,5 больше не встречается дляx ~ y. Используяx * 0,75< aв качестве критерия кроссовера для0,5< x<= 1.1сохраняет максимальное значение, вычисленное (будь то верхний или нижний интервал) примерно до 0,75. Аналогично дляx<= 0,5, затем используя-0,4 / log(x)< aв качестве критерия кроссовера сохраняет максимальное значение, рассчитанное примерно до 0,7 (будь то верхний или нижний интервал).

Существуют два особых случая, когда a — целое или половина целого числа, и условия кроссовера, перечисленные выше, указывают на то, что мы должны вычислить верхний интеграл Q. Если a — целое число в диапазоне1<= a< 30, то используется следующая конечная сумма:

9)

В то время как для полуцелых чисел в диапазоне0,5<= a< 30используется следующая конечная сумма:

10)

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

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

11)

12)

13)

14)

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

(a > 20) && (a < 200) && fabs(x-a)/a < 0.4

И:

(a > 200) && (fabs(x-a)/a < 4.5/sqrt(a))

Последний диапазон действителен для всех типов до 128-битных длинных двойников и предназначен для обеспечения того, чтобы результат был больше 10^-6, первый диапазон используется только для типов до 80-битных длинных двойников. Эти домены более узкие, чем рекомендованные Temme или Didonato и Morris. Однако использование более широкого диапазона приводит к тому, что большие и неточные (т.е. вычисленные) значения передаются в функции<exp>и<erfc>, что приводит к значительно большей частоте ошибок. Другими словами, здесь существует тонкая грань между эффективностью и ошибкой. Текущие ограничения должны удерживать количество терминов, требуемых (4) и (5), не более ~20 с двойной точностью.

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

15)

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

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

16)

Когдаa-xмал, а a и x большие. Есть еще вычитание и, следовательно, некоторые ошибки отмены - но термины малы, поэтому абсолютная ошибка будет мала - и это абсолютная, а не относительная ошибка, которая учитывается в аргументе к функцииexp. Обратите внимание, что при достаточно больших а и х ошибки все равно достанутся вам в итоге, хотя это и задерживает неизбежное гораздо дольше, чем другие методы. Использованиеlog(1+x)-xвдохновлено Temme (см. ссылки ниже).

References