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

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

$$F_m(t) = \int_0^1 u^{2m} e^{-tu^2} d u = {\gamma(m+{1\over 2}, t)\over 2 t^{m+{1\over 2}}}$$ $m=0, 1, 2, ...$$t>0$

https://gist.github.com/3764427

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

Есть ли лучший способ реализовать это? Вот реализация гамма-функции в 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^{m+{1\over2}}$$F_m$

$F_{100}$(1e-6)=4.97511945200351715E-003 ,

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

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

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

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

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

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

$$\mathrm{F}_m(a) = {1\over 2}\gamma\left(m + {1\over 2}, a\right) \times p \times p, \space \space p = a^{-{1\over 2}\left(m+ {1\over 2}\right)}$$

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

$$\gamma \left(m + {1\over 2}, a\right) = \Gamma \left(m + {1\over 2}\right) - \Gamma\left(m + {1\over 2}, a\right)$$

используя табличные значения и вычисляя соответствии с этим ответом , тщательно избегая проблемы вычитания отмены посредством использования слитой операции умножения-сложения. Потенциальная дальнейшая оптимизация состоит в том, чтобы наблюдать, что для достаточно больших , с точностью до заданная точность с плавающей точкой.Γ(m+$\Gamma \left(m + {1\over 2}\right)$γ(м+$\Gamma\left(m + {1\over 2}, a\right)$$a$$\gamma \left(m + {1\over 2}, a\right) = \Gamma \left(m + {1\over 2}\right)$

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

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

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

$$\mathrm{F}_{m}(a) = {1\over 2}\frac{1}{m + {1\over 2}}\exp(-a)\left(1+\sum_{n=1}^{\infty}\frac{a^{n}}{(1 + m + {1\over 2}) \times\space ...\space \times (n + m + {1\over 2})}\right)$$

Существуют также интересные и потенциально важные частные случаи для младших порядков функции Бойса, особенно . Во-первых, у нас есть , где - это функция ошибок, предоставляемая в Fortran 2008 в качестве элементарной функции и в C / C ++ в качестве стандартных функций библиотеки и .$m = 0, 1, 2, 3$$\mathrm{F}_{0}(a) = \sqrt{\frac{\pi}{4a}}\mathrm{erf}\left(\sqrt{a}\right)$$\mathrm{erf}$ERFerferff

Для быстрых вычислений, когда , я использую собственные минимаксные полиномиальные приближения для небольших аргументов, скажем , и прямую рекурсию , для больших где проблемы с вычитающей отменой в последнем случае уменьшаются путем использования слитых операций умножения-сложения.$m = 1, 2, 3$$a \lt {2{1\over 2}}$$\mathrm{F}_{m}(a) = \frac{1}{2a}\left(\left(2m-1\right)\mathrm{F}_{m-1}(a) - \exp(-a)\right)$

В тех случаях, когда значения функции должны быть вычислены для заданного значения по нескольким порядкам , можно вычислить значение функции непосредственно для самого высокого значения , т. Е. Как обсуждалось выше, а затем использовать численно устойчивую обратную рекурсию для вычисления все остальные значения функций.$a$$m$$m$$\mathrm{F}_{m-1} = \frac{1}{2m-1} \left(2a \space \mathrm{F}_{m}(a) + \exp\left(-a\right)\right)$

Вы можете взглянуть на Численные методы для специальных функций Ампаро Гила, Хавьера Сегуры и Нико М. Темме.

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

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

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