Motivation

Зачем основывать гамма- и гамма-подобные функции на приближении Ланчоса?

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

Именно сочетание этих двух свойств делает приближение привлекательным: Приближение Стирлинга является очень точным для большого z и имеет некоторые из тех же аналитических свойств, что и приближение Ланчоса, но не может быть легко использовано во всем диапазоне z.

В качестве простейшего примера рассмотрим соотношение двух гамма-функций: можно вычислить результат с помощью lgamma:

exp(lgamma(a) - lgamma(b));

Тем не менее, даже если lgamma одинаково точна для 0,5 ulp, в худшем случае относительная ошибка может быть легко показана:

Erel > a * log(a)/2 + b * log(b)/2

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

Напротив, аналитически комбинируя подобные термины мощности в соотношении приближений Ланчоса, эти ошибки могут быть практически устранены для малыхaиbи держаться под контролем для очень больших (или очень маленьких в этом отношении)aиb. Конечно, вычисление больших мощностей само по себе является печально известной трудной проблемой, но даже в этом случае аналитические комбинации приближений Ланчоса могут сделать разницу между получением действительного результата или просто мусором. См. примечания к реализации функциибетадля примера этого метода на практике. Неполные функцииgamma_p gammaиbetaиспользуют аналогичные аналитические комбинации терминов мощности, чтобы объединить гамма- и бета-функции, разделенные большими полномочиями на единичные (более простые) выражения.

The Approximation

Приближение Ланчоса к гамма-функции определяется:

Если S_g(z) является бесконечной суммой, то есть конвергентной для всех z >0, иgявляется произвольным параметром, который управляет «формой» терминов в сумме, которая дается:

С индивидуальными коэффициентами, определяемыми в закрытом виде:

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

Поэтому приближение Ланчоса часто записывается в виде частичных дробей с ведущими константами, поглощенными коэффициентами в сумме:

где:

Параметрgявляется произвольно выбранной константой, аNявляется произвольно выбранным числом терминов для оценки в части «Сумма Ланчоса».

	Note
	Некоторые авторы определяют сумму от k=1 до N и, следовательно, получают коэффициенты N+1. Это приводит к путанице как следующего обсуждения, так и кода (поскольку C++ имеет дело с полуоткрытыми диапазонами массивов, а не с закрытым диапазоном суммы). Это соглашение согласуется сГодфри, но не сПью, поэтому будьте осторожны, когда ссылаетесь на литературу в этой области.

Computing the Coefficients

Коэффициенты C0..CN-1 необходимо вычислить изNиgс высокой точностью, а затем сохранить как часть программы. Расчет коэффициентов производится методомГодфри; пусть константы содержатся в векторе столбца Р, тогда:

P = D B C F

где B - матрица NxN:

D - матрица NxN:

C - матрица NxN:

и F - вектор столбца элемента N:

Обратите внимание, что матрицы B, D и C содержат все целые термины и зависят только отN, этот продукт следует вычислить сначала, а затем умножить наFв качестве последнего шага.

Choosing the Right Parameters

Хитрость состоит в том, чтобы выбратьNиg, чтобы дать желаемый уровень точности: выбор небольшого значения дляgприводит к строго конвергентному ряду, но тот, который сходится только медленно. Выбор большего значениягприводит к тому, что термины серии являются большими и/или расходящиеся примерно для первыхг-1терминов, а затем внезапно сходятся с «сжатием».

Пьюопределил оптимальное значениегдляNв диапазоне1<= N<= 60: к сожалению, на практике выбор этих значений приводит к ошибкам отмены в сумме Ланчоса, так как наибольший срок в (переменный) ряд примерно в 1000 раз больше результата. Эти оптимальные значения, по-видимому, не будут полезны на практике, если оценка не может быть выполнена с помощью ряда защитных цифри, коэффициенты хранятся с более высокой точностью, чем это требуется в результате. Эти значения лучше всего зарезервировать для, скажем, вычислений с плавающей точностью с двойной точностью арифметики.

Альтернативой, описаннойГодфри, является выполнение исчерпывающего поискаNиgпространства параметров для определения оптимальной комбинации для данногоpцифрового типа с плавающей точкой. Повторение этой работы нашло хорошее приближение для двухточной арифметики (близко к найденномуГодфри), но не нашло действительно хороших приближений для 80 или 128-битных длинных двойников. Кроме того, было отмечено, что полученные приближения имеют тенденцию к оптимизации для малых значений z (1< z< 200), используемых для проверки реализации против факториалов. Было отмечено, что вычислительные соотношения гамма-функций с большими аргументами страдают от ошибок, возникающих в результате усечения ряда Ланкозо.

Пьюопределил все места, где теоретическая ошибка приближения была минимальной, но, к сожалению, опубликовал только самое большое из этих минимумов. Тем не менее, он делает замечание, что минимумы близко совпадают с местом, где первый забытый термин (a_N) в серии Ланчоса S_g(z) меняет знак. Эти места довольно легко найти, хотя и со значительным компьютерным временем. Эти «сладкие пятна» нужно вычислить только один раз, свести в таблицу, а затем искать, когда это необходимо для приближения, которое обеспечивает требуемую точность для некоторого фиксированного типа точности.

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

На этом этапе отметим, что сумма Ланчоса может быть преобразована в рациональную форму (соотношение двух полиномов, полученное из формы частичного дробления с использованием полиномиальной арифметики), и при этом изменяет коэффициенты так, чтовсе они являются положительными. Это означает, что сумма в рациональной форме может быть оценена без ошибки отмены, хотя и с удвоением числа коэффициентов для данного N. Повторяя поиск «сладких пятен», на этот раз оценивая сумму Ланчоса в рациональной форме, и испытывая только те «сладкие пятна», теоретическая ошибка которых меньше, чем у машинного эпсилона для тестируемого типа, давала хорошие приближения для всех тестируемых типов. Оптимальные найденные значения были довольно близки к лучшим случаям, сообщаемымPugh(просто немного большеNи немного меньшеgдля заданной точности, чемPughотчетов), и хотя преобразование в рациональную форму удваивает количество хранимых коэффициентов, следует отметить, что половина из них являются целыми числами (и, следовательно, требуют меньше места для хранения), и приближения требуют меньшегоN, чем в противном случае требовалось бы, поэтому в целом может потребоваться меньше операций с плавающей точкой.

В следующей таблице показаны оптимальные значения дляNиgпри вычислении с фиксированной точностью. Они должны приниматься за работу в процессе: нет значений для 106-битных знаковых и машин (Darwin long doubles & NTL quad_float), и возможна дальнейшая оптимизация значенийг. Ошибки, приведенные в таблице, являются оценками ошибки, обусловленной усечением бесконечного ряда Ланчоса доNтерминов. Они вычисляются из суммы первых пяти забытых терминов - и, как известно, являются довольно пессимистическими оценками - хотя заметно, что лучшие комбинацииNиgпроизошли, когда предполагаемая ошибка усечения почти точно соответствует эпсилону машины для рассматриваемого типа.

Наконец, обратите внимание, что приближение Ланчоса можно записать следующим образом, удалив коэффициент exp(g) из знаменателя, а затем разделив все коэффициенты на exp(g):

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

References