Быстрая и точная реализация неполной гамма-функции с двойной точностью

10

Каков современный способ реализации специальных функций двойной точности? Мне нужен следующий интеграл: для и , что можно записать в терминах нижней неполной гамма-функции. Вот моя реализация на Фортране и Си: м=0,1,2,. , , t>0

Fм(T)знак равно01U2ме-TU2dUзнак равноγ(м+12,T)2Tм+12
мзнак равно0,1,2,,,,T>0

https://gist.github.com/3764427

который использует последовательное расширение, суммирует термины до заданной точности, а затем использует рекурсивные соотношения для эффективного получения значений для меньшего . Я хорошо это проверил, и я получаю точность 1e-15 для всех значений параметров, которые мне нужны, подробности см. В комментариях к версии на Fortran.м

Есть ли лучший способ реализовать это? Вот реализация гамма-функции в gfortran:

https://github.com/mirrors/gcc/blob/master/libgfortran/intrinsics/c99_functions.c#L1781

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

Обновление 1 :

Основываясь на комментариях, вот реализация с использованием SLATEC:

https://gist.github.com/3767621

он воспроизводит значения из моей собственной функции, примерно на уровне точности 1e-15. Однако я заметил проблему, заключающуюся в том, что для t = 1e-6 и m = 50 член становится равным 1e-303, а для более высоких «m» он просто начинает давать неправильные ответы. У моей функции нет этой проблемы, потому что я использую отношения расширения / повторения ряда непосредственно для . Вот пример правильного значения: FмTм+12Fм

F100(1e-6)=4.97511945200351715E-003 ,

но я не могу получить это, используя SLATEC, потому что знаменатель взрывается. Как видите, фактическое значение хорошее и маленькое.Fм

Обновление 2 :

Чтобы избежать вышеупомянутой проблемы, можно использовать функцию dgamit(неполная гамма-функция Трикоми), то F(m, t) = dgamit(m+0.5_dp, t) * gamma(m+0.5_dp) / 2есть проблем с больше нет, но, к сожалению, взрывы для . Это , однако , может быть достаточно высокими для моих целей.м 172 мTgamma(m+0.5_dp)м172м

Ондржей Чертик
источник
2
Зачем кодировать свою собственную функцию? GSL, cephes и SLATEC все это реализуют.
Джефф Оксберри
Я обновил вопрос, почему я не использую SLATEC.
Ондржей Чертик
@ OndřejČertík Вы обнаружили ошибку! Проголосовал твой вопрос!
Али
Али --- это не ошибка в SLATEC, но в том факте, что мне действительно нужно разделить на , чтобы получить значение для . Таким образом, численный метод, который работает для может не очень хорошо работать для . t m + 1γ(z,x) Fm(t)γ(z,x)Fm(t)tm+12Fm(t)γ(z,x)Fм(T)
Ондржей Чертик
@ OndřejČertík Хорошо, извините, моя ошибка, я не проверял ваш код, прежде чем сделать свой комментарий.
Али

Ответы:

9

Рассматриваемый интеграл также известен как функция Бойса, в честь британского химика Сэмюэля Фрэнсиса Бойса, который ввел его использование в начале 1950-х годов. Несколько лет назад мне нужно было вычислить эту функцию с двойной точностью, максимально быстро, но точно. Мне удалось добиться относительной ошибки порядка во всем входном домене.10-15

Как правило, выгодно использовать разные приближения для малых и больших аргументов, где оптимальное переключение между «большим» и «малым» лучше всего определяется экспериментально и, как правило, является функцией от . Для моего кода я определил «маленькие» аргументы как те, которые удовлетворяют условию .a m + 1 1мaм+112

Для больших аргументов я вычисляю

Fм(a)знак равно12γ(м+12,a)×п×п,  пзнак равноa-12(м+12)

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

γ(м+12,a)знак равноΓ(м+12)-Γ(м+12,a)

используя табличные значения и вычисляя соответствии с этим ответом , тщательно избегая проблемы вычитания отмены посредством использования слитой операции умножения-сложения. Потенциальная дальнейшая оптимизация состоит в том, чтобы наблюдать, что для достаточно больших , с точностью до заданная точность с плавающей точкой.Γ(m+1Γ(м+12)γ(м+1Γ(м+12,a)aγ(м+12,a)знак равноΓ(м+12)

Для небольших аргументов я начал с разложения в ряд для нижней неполной гамма-функции из

А. Эрдельи, В. Магнус, Ф. Оберхеттингер и Ф. Г. Трикоми, "Высшие трансцендентные функции, том 2". Нью-Йорк, Нью-Йорк: Макгроу-Хилл, 1953

и изменил его, чтобы вычислить функцию Бойса следующим образом (усечение ряда, когда член достаточно мал для заданной точности):Fм(a)

Fм(a)знак равно121м+12ехр(-a)(1+ΣNзнак равно1aN(1+м+12)× ,,, ×(N+м+12))

Существуют также интересные и потенциально важные частные случаи для младших порядков функции Бойса, особенно . Во-первых, у нас есть , где - это функция ошибок, предоставляемая в Fortran 2008 в качестве элементарной функции и в C / C ++ в качестве стандартных функций библиотеки и .мзнак равно0,1,2,3F0(a)знак равноπ4aере(a)ереERFerferff

Для быстрых вычислений, когда , я использую собственные минимаксные полиномиальные приближения для небольших аргументов, скажем , и прямую рекурсию , для больших где проблемы с вычитающей отменой в последнем случае уменьшаются путем использования слитых операций умножения-сложения.мзнак равно1,2,3a<212Fм(a)знак равно12a((2м-1)Fм-1(a)-ехр(-a))

В тех случаях, когда значения функции должны быть вычислены для заданного значения по нескольким порядкам , можно вычислить значение функции непосредственно для самого высокого значения , т. Е. Как обсуждалось выше, а затем использовать численно устойчивую обратную рекурсию для вычисления все остальные значения функций.aммFм-1знак равно12м-1(2a Fм(a)+ехр(-a))

njuffa
источник
Спасибо @njuffa за отличный ответ. Если вы создадите свой код для этого открытого исходного кода, я думаю, что это будет очень полезно для многих людей.
Ондржей Чертик,
1
В настоящее время реализация описанного алгоритма в CUDA доступна для бесплатной загрузки с сайта разработчика NVIDIA (требуется бесплатная регистрация в качестве разработчика CUDA, одобрение обычно в течение одного рабочего дня). Код находится под лицензией BSD, которая должна быть совместима практически с любым проектом.
njuffa
4

Я бы взглянул на книгу Абрамовича и Стегуна, или на новую редакцию, опубликованную NIST пару лет назад, и я думаю, что она доступна онлайн. Они также обсуждают способы стабильной реализации вещей.

Вольфганг Бангерт
источник
Я использовал это: dlmf.nist.gov/8 , при его реализации, но это, вероятно, другой ресурс. Глава 5 «Числовые рецепты» также содержит интересную информацию, но применима только к функциям одной переменной.
Ондржей Чертик
Я не думаю, что вы найдете что-то намного более свежее, чем их справка 2001 года; SLATEC будет старше этого.
Джефф Оксберри
1

Похоже, это не современно, но SLATEC в Netlib предлагает «1400 общих математических и статистических процедур». Неполная гамма доступна по специальным функциям здесь .

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

Али
источник