Автоматика и телемеханика, № 10, 2022

(ekaratsuba@gmail.com)

(ФИЦ ”Информатика и управление” РАН, Москва)

БЫСТРЫЙ АЛГОРИТМ ВЫЧИСЛЕНИЯ ПСИ-ФУНКЦИИ¹

Памяти Константина Владимировича Рудакова

Построен быстрый алгоритм вычисления логарифмической производ-

ной гамма-функции Эйлера, основанный на БВЕ методе. Сложность ал-

горитма близка к оптимальной. Структура алгоритма допускает его рас-

параллеливание.

Ключевые слова: быстрые алгоритмы, пси-функция, гамма-функция

Эйлера, сложность вычисления, метод БВЕ.

DOI: 10.31857/S0005231022100105, EDN: ALDLAI

1. Введение. Быстрые алгоритмы. Метод БВЕ

Быстрые алгоритмы предназначены для оптимизации высокоточных вы-

числений. Основным критерием оценки быстрых алгоритмов является слож-

ность вычисления. Первые задачи по оценке (битовой) сложности вычисления

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

лированы А.Н. Колмогоровым в 1950-х гг. (см. [1, 2]), первый быстрый ал-

горитм (алгоритм умножения n-значных чисел) был найден А.А. Карацубой

в 1960 г. (см. [1-3]), что привело в дальнейшем к созданию серии алгорит-

мов быстрого вычисления различных функций (алгоритм умножен(ия эк)ива-

лентен алгоритму вычисления функции y = x²) со сложностью O

n^1+ε

для

(

)

любого ε > 0 и n > n₁(ε) (вместо O

n^2+ε

лучшей сложности вычисления

функций до эры быстрых алгоритмов).

Далее считаем, что числа записаны в двоичной системе счисления, знаки

которой 0 и 1 называются битами.

Определение 1. Запись знаков 0, 1, плюс, минус, скобка; сложение,

вычитание и умножение двух битов назовем одной элементарной или би-

товой операцией.

Определение 2. Вычислить (вещественную) функцию y = f(x) в точ-

ке x = x₀ с точностью до n знаков, значит найти такое число A, что

|f(x₀) - A| ≤ 2^-n.

Определение 3. Количество битовых операций, достаточное для вы-

числения функции f(x) в точке x₀ с точностью до n знаков посредством

данного алгоритма, называется сложностью вычисления f(x) в точке x₀ и

обозначается s_f(n) = sf,x0(n).

¹ Работа выполнена при финансовой поддержке Российского научного фонда (грант

№ 22-21-00727).

105

Будем далее предполагать, что для сложности умножения двух n-значных

чисел справедлива оценка

(

)

(1)

M (n) = O

n^1+ε

т.е. сложность используемого алгоритма умножения не хуже, чем M(n) =

= O(n log^C n), где C константа (см. алгоритмы из [4, 5]).

Определение 4. Будем называть быстрыми такие алгоритмы вы-

числения функции f, что для этих алгоритмов

(

)

s_f(n) = O

n^1+ε

для любого ε > 0 и n > n₁(ε).

В теории быстрых алгоритмов основным растущим параметром является точ-

ность вычисления n, n → ∞. Поскольку только запись некоторого числа с

точностью до n знаков требует не менее, чем n + 1 операций, из определе-

ния 4 следует, что быстрым вычислительным алгоритмам соответствует пра-

вильный порядок оценки сверху сложности вычисления s_f (n) по n, n → ∞,

n < s_f(n) < n^1+ε.

В 1975 г. были предложены первые алгоритмы быстрого вычисления эле-

ментарных алгебраических функций [6]. Например, самый простой алгоритм

деления числа a на число b заключается в вычислении методом Ньютона об-

ратной величины^1b с точночтью до n знаков с последующим умножением

на a по алгоритму быстрого умножения.

В 1976 г. (см., например, [7]) были предложены первые быстрые алгорит-

мы вычисления константы π, а затем и простейших трансцендентных функ-

ций, основанные на АГС-методе Гаусса (см., например, [8]). В [9], кроме ал-

горитмов вычисления простейших трансцендентных функций, основанных

на АГС-методе Гаусса, содержатся также алгоритмы вычисления некоторых

высших трансцендентных функций (гамма-функции Эйлера, например). Од-

нако сложность этих алгоритмов несколько хуже определенной выше слож-

ности быстрых алгоритмов и оценивается как

(

)

O n^3/2+ε

В 1991 г. автором был построен новый быстрый метод метод Быстрого

Вычисления Е-функций (БВЕ) (см. [10], см. также [11, 12], о E-функциях

см. [13]).

В настоящее время известны два быстрых метода вычисления простейших

трансцендентных функций метод АГС и метод БВЕ. При этом с помощью

БВЕ удается построить быстрые алгоритмы и для вычисления некоторых

высших трансцендентных функций.

Иногда метод БВЕ по ошибке принимают за вариант метода Байнэри

Сплиттинг (Binary Splitting). В 2010 г. Билл Роско (Andrew William Roscoe)

глава факультета информатики Оксфордского университета (Department of

Computer Science, University of Oxford), ученик Тони Хоара (Charles Antony

106

Richard Hoare) ученика Андрея Николаевича Колмогорова протестировал оба

метода БВЕ и Байнэри Сплиттинг и, отметив, что ошибки от шага к шагу

в этих методах “копятся” по-разному, сделал объективный вывод о том, что

это разные методы. С другой стороны, ни в одной из публикаций до 1991 г.

не содержалось ни одного быстрого алгоритма вычисления какой-либо транс-

цендентной функции или классической константы, основанного на Байнэри

Сплиттинг. Похоже, метод Байнэри Сплиттинг был неудачной попыткой вос-

пользоваться методом А.А. Карацубы (изложен в [1-3]) для получения ал-

горитмов эффективного вычисления трансцендентных функций. При этом

только ошибки в сложности вычисления могли неправильно соориентировать

тех, кто решил, что Байнэри Сплиттинг является быстрым методом вычис-

ления трансцендентных функций.

Метод БВЕ это метод суммирования рядов специального вида. С помо-

щью БВЕ можно вычислить любую элементарную трансцендентную функ-

цию для любого аргумента, классические константы e, π, постоянную Эйле-

ра γ, постоянные Апери и Каталана, такие высшие трансцендентные функ-

ции как гамма-функцию Эйлера, гипергеометрические функции, сфериче-

ские функции, цилиндрические функции и т.д. для алгебраических значений

аргумента и параметров, такие специальные интегралы, как интеграл вероят-

ности, интегралы Френеля, интегральную экспоненциальную функцию, ин-

тегральные синус и косинус и т.д. с оценкой сложности

(

)

(

)

s_f(n) = O

M (n) log² n

n^1+ε

для любого ε > 0 и n > n₁(ε). Структура метода БВЕ дает возможность рас-

параллеливания основанных на БВЕ алгоритмов.

В 2008 г. профессор Эрик Бах (Eric Bach, Univ. of Wisconsin, Madison)

в письме заметил, что никто не знает, как быстро вычислить пси-функцию

(про ψ-функцию см., например, [14]). Алгоритм вычисления пси-функции,

основанный на БВЕ-методе, казался мне достаточно очевидным, но посколь-

ку за эти годы никто его не описал, настоящая статья призвана восполнить

это упущение.

2. Гамма-функция и пси-функция

Одна из самых широко распространенных в анализе и математической

физике высших трансцендентных функций гамма-функция Эйлера опреде-

ляется соотношением (см., например, [14])

+∞

(2)

Γ(z) =

t^z-1e^-t

dt,

Γ(z + 1) = zΓ(z), Re z > 0,

Γ(N + 1) = N!, где N натуральное число. Логарифмическая производная

гамма-функции называется пси-функцией

(3)

ψ(z) =

log Γ(z), ψ(z + 1) =

+ ψ(z).

107

Представление пси-функции в виде ряда

(

)

∑

(4)

ψ(z) = -γ +

z = 0,-1,-2,...,

k+1

z+k

k=0

где γ есть константа Эйлера

(

)

γ = lim

+···+

- log n

n→∞

не дает возможности быстро вычислить эту функцию ввиду медленной сходи-

мости ряда из правой части (4) (о быстром алгоритме вычисления константы

Эйлера см. [10, 15]).

Далее для простоты будем рассматривать лишь случаи вещественного ар-

гумента (для комплексного аргумента все рассуждения нужно проводить от-

дельно для каждой из частей вещественной и мнимой).

Так же, как и асимптотическое выражение для гамма-функции, получае-

мое потенциированием ряда Стирлинга (см., например, [14]), асимптотиче-

ское выражение для пси-функции (см. доказательство в [15])

ψ(x) = log x -

0 ≤ θ = θ(x) ≤ 1,

12x²

120x⁴

126x⁶

не может служить основой для построения быстрого процесса.

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

вавшись уже существующим БВЕ-алгоритмом вычисления гамма-функции и

построив БВЕ-алгоритм вычисления производной (обычной, не логарифми-

ческой) гамма-функции. Из формулы (3) легко видеть, что

Γ^′(x)

(5)

ψ(x) =

Γ(x)

Таким образом, вычислив значения Γ^′(x) и Γ(x) в точке x = x₀ с точностью

(

)

2-(n+1) за O

n^1+ε

операций, а затем разделив одно число на другое с точ-

ностью 2^-n методом Ньютона за

(6)

O (M(n)log n)

операций, получим с учетом (1) зн)чение пси-функции, вычисленное в точке

x = x₀ с точностью 2^-n за O

n^1+ε

операций.

В настоящее время как для гамма-функции, так и для пси-функции быст-

рые алгоритмы на основе метода БВЕ можно построить только для алгеб-

раического аргумента. Как правило, задача изучения, а также уменьшения

константы в знаке O в сложности вычисления (что приводит к практическому

улучшению эффективности вычисления в частных случаях) редко рассмат-

ривается в теории быстрых алгоритмов (см. [16]), однако по построению алго-

ритма можно заметить, что алгоритм для алгебраического нерационального

аргумента предполагает гораздо большую константу “в O”, чем алгоритм для

рационального аргумента. Поэтому, хотя рациональные числа являются под-

множеством алгебраических, отдельно строится алгоритм для рационального

аргумента.

108

3. БВЕ алгоритмы вычисления гамма-функции

На основе применения метода БВЕ можно доказать следующие утвержде-

ния относительно сложности вычисления гамма-функции Эйлера (рассмот-

рим для простоты случай вещественного аргумента, для комплексного аргу-

мента соответствующие теоремы доказываются отдельно для вещественной

и мнимой частей аргумента).

Теорема 1. Пусть y = Γ(x₀); x₀ = p/q; p, q целые, (p,q) = 1. Тогда

(

)

(7)

s_Γ(p/q)(n) = O

M (n) log² n

Теорема 2. Пусть y = Γ(x₀); x₀ = α; α алгебраическое число, которое

является корнем многочлена с целыми (заранее известными) коэффициен-

тами. Тогда

(

)

(8)

s_Γ(α)(n) = O

M (n) log² n

Для удобства читателя напомним основные моменты доказательства теоре-

мы 1 из [10]. Вычисляем значение Γ(x₀), предполагая сначала, что 0 < x₀ < 1,

т.е. что 0 < p < q. Воспользовавшись представлением гамма-функции в виде

интеграла из (2) и тем, что в разложении этого интеграла на сумму двух

интегралов

∫

e^-ttx0-1dt = e^-ttx0-1dt +

e^-ttx0-1dt,

для последнего интеграла при a = n, 0 < x₀ < 1, справедлива оценка

∫

e^-ttx0-1dt ≤

e^-tdt = θ₁2^-n,

0 < θ₁ < 1,

a^1-x0

получаем отсюда, что

∫n

(9)

Γ(x₀) = e^-ttx0-1dt + θ₁2^-n = nx0 S + θ₂2^-n,

|θ₁| ≤ 1,

|θ₂

| ≤ 1,

где

∑

n^j

(10)

(-1)^j

r + 1 ≥ 4n, n ≥ 8.

j!(x₀ + j)

j=0

Будем вычислять Γ(x₀) = Γ(p/q) по формулам (9)-(10). Заметим, что для

вычисления nx0 = n^p/q, с точностью до n знаков достаточно O(M(n) log n) =

= O(n log² n log log n) операций (методом Ньютона (см. [6]) при выполнении

умножения методом Шенхаге-Штрассена (см. [4])). Чтобы вычислить S,

109

возьмем в (10) r + 1 = 2^k, 2^k-1 < 4n ≤ 2^k, k ≥ 1. Вычисляем сумму S с по-

мощью БВЕ-процесса за k шагов, имея на j-м шаге (j ≤ k)

r+1

S = S₁(j) + S₂(j) + ··· + Srj-1(j) + Sr_j(j), rj =

где Srj -ν (j); ν = 0, 1, . . . , rj - 1; определяются равенствами

r-(2^jν+2^j-1)

prj-ν(j)

Srj-ν(j) = Sr_j-1-2ν(j - 1) + Sr_j-1-2ν-1(j - 1) =

(r - 2^jν)! qrj -ν (j)

При этом на 1-м шаге (j = 1) вычисляются целые числа

pr1-ν(1) = -qn(qr - 2qν - q + p) + (qr - 2qν)(qr - 2qν + p);

qr1-ν(1) = (qr - 2qν + p)(qr - 2qν - q + p);

r+1

ν = 0,1,...,

- 1; r + 1 = 2^k, k ≥ 1;

а на j-м шаге (1 < j ≤ k) вычисляются целые числа

prj-ν(j) = n2j-1 pr_j-1-2ν(j - 1)qr_j-1-2ν-1(j - 1) +

(r - 2^j ν)!

prj-1-2ν-1(j - 1)qr_j-1-2ν(j - 1);

(r - 2^j-1(2ν + 1))!

qrj-ν(j) = qr_j-1-2ν(j - 1)qr_j-1-2ν-1(j - 1);

ν = 0,1,...,r_j - 1; r_j = (r + 1)/2^j; r + 1 = 2^k, k ≥ 1.

На k-м (последнем) шаге вычисляются целые числа prk (k) = p₁(k), qrk (k) =

= q₁(k), r! и производится одно деление с точностью до 2^-2n знаков по фор-

муле

prk (k) 1

S = S₁(k) =

qrk (k) r!

что дает сумму S с точностью 2^-n-1, а следовательно и значение Γ(x₀) =

= Γ(p/q) с точностью 2^-n+1, n ≥ 1. Сложность такого вычисления есть

(см. [10])

(

)

(11)

s_Γ(p/q)(n) = O

M (n) log² n

Представим доказательство теоремы 2 и алгоритм быстрого вычисления

гамма-функции Эйлера при алгебраическом аргументе α в виде, удобном для

его применении для быстрого вычисления пси-функции.

Доказательство. Как и для случая рационального аргумента, в со-

ответствии с замечанием 1 сведем вычисление Γ(x₀) при алгебраическом x₀

к вычислению Γ(α), 0 < α < 1, где α есть алгебраическое число степени N,

N ≥ 2. Для простоты рассмотрим случай, когда α является вещественным

110

числом. Предполагается, что известен многочлен g(x), наименьшей степени

с целыми коэффициентами, корнем которого является число α, т.е.

(12)

g(x) = g_N x^N + g_N-1x^N-1 + · · · + g₁x + g₀

;

g(α) = 0;

g_N,g_N-1,... ,g₀

целые числа, N ≥ 2. Как и в случае рационального ар-

гумента, воспользуемся формулами (9), (10). Отметим, что вычисление n^α

с точностью 2^-2n по формуле n^α = e^αlogn требует O(M(n) log² n) операций.

Сумма S из (10) принимает вид

∑

n^j

(13)

(-1)^j

j!(j + α)

j=0

(14)

S = S₁(0) + S₂(0) + ··· + S_r+1(0),

где

r-ν

(15)

S_r+1-ν(0) = (-1)^r-ν

;

ν = 0,1,...,r.

(r - ν)!(r - ν + α)

Вычисление S выполняется за k шагов, r + 1 = 2^k; 2^k-1 < 4n ≤ 2^k; k ≥ 1,

БВЕ процесса для “алгебраического случая”, особенностью которого являет-

ся использование основного свойства алгебраических чисел для ограничения

роста сложности вычисления.

До шага i, где i определяется условиями 1 ≤ i ≤ k; 2^i-1 < N ≤ 2ⁱ, произво-

дим с суммой (14), (15) действия, аналогичные описанным выше для случая

рационального аргумента. Однако, группируя слагаемые суммы (14) таким

же образом, теперь вычисляем на каждом шаге не числитель и знаменатель

дробей, находящихся в скобках, а лишь целые коэффициенты при степенях α,

многочленов, находящихся в числителе и знаменателе дробей “из скобок”.

На 1-м шаге в скобках находятся дроби

α(r - 2ν - n) + (r - 2ν)² - n(r - 2ν - 1)

βr1-ν(1) =

α² + α(2r - 4ν - 1) + (r - 2ν)(r - 2ν - 1)

r+1

ν = 0,1,2,...;

- 1;

вычисляются числа

δr1-ν(0,1) = (r - 2ν)² - n(r - 2ν - 1); δr₁-ν(1,1) = r - 2ν - n;

ρr1-ν(0,1) = (r - 2ν)(r - 2ν - 1); ρr₁-ν(1,1) = 2r - 4ν - 1.

На i-м шаге (1 < i ≤ k) в скобках находятся дроби

δri-ν(i)

(16)

βri-ν(i) =

ρri-ν(i)

111

где

∑

(17)

δri-ν(i) =

δri-ν(m,i)α^m; ρr_i-ν(i) =

ρri-ν(l,i)α^l;

m=0

l=0

вычисляются числа

(

∑

δri-ν(m,i) =

n2i-1δri-1-2ν(m1,i - 1)ρr_i-1-2ν-1(m2,i - 1) +

m1+m2=m

0≤m1≤2i-1-1

0≤m2≤2i-1

)

(r - 2ⁱν)!

δri-1-2ν-1(m1,i - 1)ρr_i-1-2ν(m2,i - 1)

;

(r - 2ⁱν - 2^i-1)

∑

ρri-ν(l,i) =

ρri-1-2ν(l1,i - 1)ρr_i-1-2ν-1(l2,i - 1);

l1+l2=l

0≤l1≤2i-1

0≤l2≤2i-1

r+1

m = 0,1,2,... ,2ⁱ - 1, l = 0,1,2,... ,2ⁱ, ν = 0,1,2,... ,

- 1.

Заметим, что умножение на степень α и вычисление δ и ρ, так же как и

деление δ на ρ, не производится до последнего шага. Перед i + 1-м шагом

(1 ≤ i ≤ k, 2^i-1 < N ≤ 2ⁱ) многочлены (17) редуцируются по модулю много-

члена g(x) из (12) при x = α. Таким образом, если

δri-ν(i) = P(x) = pN-1x^N-1 + pN-2x^N-2 + ··· + p1x + p0,

ρri-ν(i) = Q(x) = qNx^N + qN-1x^N-1 + ··· + q1x + q0,

где p_N-1, . . . , p₀; q_N , . . . , q₀

целые, N = 2ⁱ, то мы делим P (x) и Q(x) на g(x)

с остатками:

P (x) = g(x)P₀(x) + P₁(x), Q(x) = g(x)Q₀(x) + Q₁(x),

где P₁(x) и Q₁(x) есть многочлены с рациональными коэффициентами, и

при этом степени многочленов P₁(x) и Q₁(x) не превышают N - 1. Отсюда и

из (12) получаем

P (α) = P₁(α), Q(α) = Q₁(α),

т.е. в (16) имеем

P₁(α)

(18)

βri-ν(i) =

Q₁(α)

Домножая, если нужно, числитель и знаменатель дроби (18) на целый общий

множитель, получаем в числителе и знаменателе дробей “из скобок” много-

члены с целыми коэффициентами степени не большей, чем N - 1.

112

Начиная с шага i, на каждом шаге i, i + 1, . . . , k, производится редукция

многочленов от α, расположенных в числителе и знаменателе β_µ(j), по моду-

лю g(x). Поскольку коэффициенты g(x), так же как и степень g(x) являются

абсолютными постоянными, то значения этих коэффициентов, участвующие

в представленных выше редукциях, могут возрасти лишь не более чем в g^N

раз, где

g = max

|g_i|,

0≤i≤N

что не влияет на оценку сложности проводимых вычислений. На последнем

k-м шаге такого процесса получаем

βrk (k)

δrk (k)

S = S₁(k) =

r!ρrk (k)

где

δrk (k) = δrk (0,k) + δrk (1,k)α + ··· + δrk (j,k)α^j,

ρrk (k) = ρrk (0,k) + ρrk (1,k)α + ··· + ρrk (j,k)α^j;

j,j ≤ N - 1, и вычисляем δ_r

(k), ρrk (k) и S, и следовательно, Γ(α) с точно-

стью 2^-n+1.

Следовательно, сложность вычисления

(

)

s_Γ(α)(n) = O

M (n) log² n

Что и требовалось доказать. Из (1) и оценок (7), (8) следует, что

(

)

sΓ(x0)(n) = O

n^1+ε

для любого ε > 0 и n > n₁(ε), и любого алгебраического аргумента x₀.

Замечание 1. Мы вычислили Γ(x₀), предполагая, что

0 < x₀ < 1.

Если x₀ ≥ 1, то пользуясь соотношением Γ(x + m) = x(x + 1)(x + 2)

...(x + m - 1)Γ(x), сводим вычисление в нужной точке к вычислению в точ-

ке из интервала (0, 1) и вычислению фиксированного (константа) числа пе-

ремножений, что, хотя и увеличивает константу в O из оценок (7), (8), сами

оценки не меняет.

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

При вещественном x, x > 0, дифференцируя (2) по параметру x, получаем

+∞

(19)

Γ^′(x) =

t^x-1e^-t

log tdt.

113

Далее считаем, что 0 < x < 1. Представим интеграл (19) в виде суммы двух

интегралов

∫a

∫

(20)

Γ^′(x) = t^x-1e^-t log tdt +

t^x-1e^-t

log tdt, a ≥ 1.

Оценим последний интеграл в (20). Так как при t ≥ 1, t > log t, то

∫

(21)

t^x-1e^-t log tdt <

t^xe^-tdt ≤

te^-tdt ≤ (a + 1)e^-a.

Рассмотрим теперь первый интеграл правой части формулы (20):

∫a

(22)

I = t^x-1e^-t

log tdt.

Разложим функцию e^-t в знакопеременный ряд Тейлора по степеням t. Для

0 ≤ t ≤ a, m ≥ a,

t³

t^2m

(23)

e^-t = 1 -

+···+

+θ

|θ| ≤ 1.

(2m)!

Подставляя разложение (23) в (22), находим

(24)

I =I₁

+ R(a),

где

∫

∑

(-1)^j

(25)

I₁ =

t^j+x-1

log tdt,

j=0

∫a

t^2m+x-1

(26)

R(a) = θ

log tdt,

|θ| ≤ 1.

(2m)!

Оценим остаточный член R(a). Имеем из (26)

∫a

t^2m+x-1|log t|

(27)

|R(a)| ≤

dt =

(2m)!

∫1

∫

t^2m+x-1

log t

t^2m+x-1 log t

a^2m+x log a

dt +

dt ≤

(2m)!

(2m)!(2m + x)

114

∫_a

Рассмотрим интеграл I₁ из (25). Интегрируя

t^j+x-1 log tdt по частям, по-

лучаем

∑

(-1)^j a^j+x

(28)

I₁ = log a

= ax log aG₁ - a^xG₂,

j! j + x

j! (j + x)²

j=0

где

∑

(-1)^j a^j

(29)

G₁ =

j! j + x

j=0

∑

(-1)^j a^j

(30)

G₂ =

j! (j + x)²

j=0

Из (20), (21), (24), (27)-(30) находим, что

(

)

a^2m+x log a

(

)

(31)

Γ^′(x) = a^x log aG₁ - a^xG₂ + O

(a + 1)e^-a

(2m)!(2m + x)

откуда при m = 2n, a = n, n ≥ 8, имеем следующее приближение для произ-

водной гамма-функции в точке x = x₀:

∑

n^j

(32) Γ^′(x₀) = nx0 log n

(-1)^j

j!(j + x₀)

j=0

∑

n^j

-nx0

(-1)^j

+θ₀2^-n,

|θ₀| ≤ 1.

j!(j + x₀)²

j=0

Легко видеть из (13), (32), что сумма (29) совпадает с суммой (13) и поэтому

вычисляется быстро с помощью БВЕ процесса для любого алгебраического

аргумента в соответствии с теоремами 1 и 2. Заметим, что как и ранее, быст-

рое вычисление функций log n и nx0 = ex0 logn требует не более O(M(n) log² n)

операций, что не ухудшает общую оценку сложности вычисления первого сла-

гаемого в (32), которая совпадает с оценкой (8).

Вычислим сумму (30), сначала предполагая, что x₀ = p/q, (p, q) = 1. За-

пишем эту сумму в виде (r = 2m ≥ 4n)

∑

n^j

(33)

(-1)^j

j!(j + p/q)²

j=0

Возьмем

(34)

r + 1 = 2^k,2^k-1 < 4n ≤ 2^k

,k ≥ 1;

членов ряда, из которого выделена сумма S. Пусть

r-ν

S_r+1-ν(0) = (-1)^r-ν

;

ν = 0,1,2,... ,r.

(r - ν)!(r - ν + p/q)²

115

Тогда (33) можно записать в виде

S = S₁(0) + S₂(0) + ··· + S_r+1(0).

Вычислим сумму S с помощью БВЕ-процесса за k шагов, объединяя на каж-

дом шаге слагаемые S последовательно попарно и вынося за скобки “оче-

видный” общий множитель. При этом на каждом шаге будут вычисляться

целый числитель и целый знаменатель дроби, стоящей в скобках (деление не

производится до последнего шага).

На 1-м шаге имеем

r+1

S = S₁(1) + S₂(1) + ··· + Sr1(1); r₁ =

r-2ν-1

pr1-ν(1)

Sr1-ν(1) = Sr+1-2ν(0) + Sr-2ν(0) = (-1)^r-1

(r - 2ν)! qr1-ν (1)

r-2ν-1

pr1-ν(1)

(r - 2ν)! qr1-ν (1)

На 1-м шаге вычисляем целые числа

(

)

(35)

pr1-ν(1) = q²

-n(qr - 2qν - q + p)² + (r - 2ν)(qr - 2qν + p)²

(36)

qr1-ν(1) = (qr - 2qν + p)²(qr - 2qν - q + p)²,

r+1

ν = 0,1,...,

- 1; r + 1 = 2^k, k ≥ 1.

На j-м шаге (j ≤ k) имеем

r+1

S = S₁(j) + S₂(j) + ··· + Srj-1(j) + Sr_j(j); rj =

;

где Srj -ν (j); ν = 0, 1, . . . , rj - 1; определяются равенствами

r-(2^jν+2^j-1)

prj-ν(j)

Srj-ν(j) = Sr_j-1-2ν(j - 1) + Sr_j-1-2ν-1(j - 1) =

(r - 2^jν)! qrj -ν (j)

На j-м шаге (j ≤ k) вычисляем целые числа

(37) prj -ν (j) = n2j-1 pr_j-1-2ν (j - 1)qr_j-1-2ν-1(j - 1) +

(r - 2^j ν)!

prj-1-2ν-1(j - 1)qr_j-1-2ν(j - 1),

(r - 2^j-1(2ν + 1))!

(38)

qrj-ν(j) = qr_j-1-2ν(j - 1)qr_j-1-2ν-1

(j - 1),

ν = 0,1,...,r_j - 1; r_j = (r + 1)/2^j. И так далее.

116

На k-м (последнем) шаге вычисляем целые числа prk (k) = p₁(k), qrk (k) =

= q₁(k), r! и производим одно деление с точностью до 2^-2n знаков по формуле

prk (k) 1

S = S₁(k) =

qrk (k) r!

что дает сумму S с точностью 2^-n-1. Следовательно, значение Γ^′(x₀) =

= Γ^′(p/q) вычислено с точностью 2^-n+1, n ≥ 1.

При алгебраическом x₀ = α первая сумма из (32) совпадает с суммой (13)

и вычисляется так же. Вычисление суммы

∑

n^j

(39)

(-1)^j

j!(j + α)²

j=0

проводится аналогично.

5. Сложность вычисления пси-функции

Оценим сложность вычисления Γ^′(x₀), если x₀ = p/q. Подсчитаем слож-

ность вычисления суммы (33). Сначала, пользуясь формулами (37), (38), под-

считаем количество операций, достаточное для вычисления на j + 1-м шаге,

j + 1 ≤ k, чисел prj+1-ν(j + 1), qr_j+1-ν(j + 1); ν = 0,1,2,... ,rj+1 - 1, rj+1 =

= (r + 1)/2^j+1, предполагая при этом, что числа pµ1 (j), qµ2 (j); µ₁ = 1, 2, . . . ,

k - 1, µ₂ = 1,2,...,k - 1; уже вычислены. Для этого, прежде всего, оценим

сверху разрядность чисел, с которыми проводятся вычисления на j + 1-м ша-

ге. Пусть

p(j) = maxpµ1 (j), q(j) = maxqµ2 (j).

µ₁

µ₂

Поскольку

(r - 2^j+1ν)!

≤r2j

(r - 2^j+1ν - 2^j )!

из (37), (38) находим

p(j + 1) ≤ p(j)q(j)(n2j + r2j ) ≤ p(j)q(j)(nr)2j ; q(j + 1) ≤ (q(j))².

Отсюда и из (36), (38) получаем

p (j + 1)

≤ q (j)q (j - 1)...q (1)(nr)2j+2j-1+···+2 ≤

p(1)

≤ (q (1))2j-1+2j-2+···+1 (nr)2j+1-2 .

Следовательно,

(40)

p(j + 1) ≤ p(1)(q(1))2j -1(nr)2j+1-2, q(j + 1) ≤ (q(1))2j .

117

Учитывая, что из (35), (36)

p(1) ≤ q⁴r³p², q(1) ≤ q⁴r⁴p⁴,

из (40) имеем

(41)

p(j + 1) ≤ (nqp)2j+2 r3(2j+1-1) ≤ (qprn)2j+3 , q(j + 1) ≤ (qpr)2j+2 .

Рассмотрим формулу (37). Сложность вычисления значений произведений

(r - 2^j+1ν)!

r+1

; ν = 0,1,2,... ,r_j+1 - 1, r_j+1 =

(r - 2^j+1ν - 2^j )!

2j+1

составляет

(

)

∑

(42)

M (2^τ log r)

τ=1

операций. Из (41), (42) получаем, что для вычисления prj+1-ν (j + 1) доста-

точно O(B(j + 1)) операций, где

∑

(

)

B(j + 1) = M(2^τ log r) + M

2j+3 log qprn

τ=1

Чтобы вычислить все αrj+1-ν (j + 1), которых ровно rj+1 = 2^-(j+1)(r + 1), до-

статочно O(r_j+1B(j + 1)) операций. Предполагая для сложности умноже-

ния M(µ) оценку Шенхаге-Штрассена (см. [4]), получаем, что для вычис-

ления α₁(k) на последнем шаге достаточно





∑

O r_j+1B(j + 1) =

j=1



∑

= O  (r + 1)2^-(j+1)

2τ log r(τ + log log r) log(τ + log log r)+

j=1

τ=1



∑

+ (r + 1)2^-(j+1)2^j+3 log qprn(j + 3 + log log qprn) log(j + 3 + log log qprn) =

j=1





∑

= O rlogqprn(j + loglogqprn)log(j + loglogqprn) =

j=1

(43)

= O(r log²

r log qprn log log qprn)

операций. Учитывая, что q и p это фиксированные константы, а параметр r

удовлетворяет условиям (34), получаем из (43) оценку сложности вычисления

суммы S

(44)

s_S(n) = O(n log³

n log log n).

118

Сложность вычисления суммы (30) при алгебраическом аргументе x = α до

момента редукции по модулю многочлена (12), скажем до шага i, оценивается

так же, как и в изложенном выше случае рационального аргумента. Затем,

как и ранее, на каждом шаге: i, i + 1, i + 2, . . . ; проводится редукция мно-

гочленов от α, стоящих в числителе и знаменателе дробей “в скобках” по

модулю g(x) при x = α. Поскольку коэффициенты в g(x) являются абсолют-

ными постоянными, степень g(x) также абсолютная постоянная, при таких

редукциях значения участвующих в вычислениях коэффициентов могут уве-

личиваться лишь в постоянное число раз, что не меняет оценку сложности

вычисления.

Тем самым, учитывая (8), (10), (32), (33), (44), для сложности вычисления

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

лива оценка

(

)

(45)

s_Γ′(α)(n) = O(n log³ n log log n) = O

n^1+ε

Замечание 2. Мы вычисляли Γ^′(x₀), предполагая, что 0 < x₀ < 1. Если

x₀ ≥ 1, то пользуясь соотношением

Γ^′(x + m) = ((x + 1)(x + 2)... (x + m - 1) + x(x + 2)... (x + m - 1) + ... +

+ x(x + 1)... (x + m - 2)) Γ(x) + x(x + 1)(x + 2)... (x + m - 1)Γ^′(x),

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

лению в точке из интервала (0, 1) и вычислению фиксированного (констан-

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

функции на сумму произведений фиксированных чисел, что, хотя и увели-

чивает константу в O из оценки (45), саму оценку не меняет.

Таким образом, доказаны следующие утверждения

Теорема 3. Пусть y = Γ^′(x₀); x₀ = p/q; p, q - целые, (p,q) = 1. Тогда

(

)

s_Γ′(p/q)(n) = O

M (n) log² n

Теорема 4. Пусть y = Γ^′(x₀); x₀ = α; α алгебраическое число, кото-

рое является корнем многочлена с целыми (заранее известными) коэффици-

ентами. Тогда

(

)

M (n) log² n

s_Γ′(α)(n) = O

На основе формулы (5) и теорем 1, 2, 3, 4, а также оценки (8) получаем

доказательство следующего утверждения

Теорема 5. Пусть y = ψ(x₀); x₀ = α; α алгебраическое число. Тогда

(

)

(46)

s_ψ(α)(n) = O

M (n) log² n

119

6. Заключение

Отметим, что аналогичные быстрые алгоритмы можно построить для ло-

гарифмических производных гамма функции любого порядка. Например, для

производной пси-функции находим

d Γ^′(x)

Γ^′′(x)

(Γ^′(x))2

(47)

dx Γ(x)

Γ(x)

Легко видеть, что для второй производной гамма функции

+∞

Γ^′′(x) =

t^x-1e^-t log² tdt,

можно построить БВЕ-процесс, аналогичный представленному выше. Другие

функции и другое слагаемое из (47) у нас уже вычислены со сложностью (46).

Таким образом, последовательно вычислив все производные до порядка K и

переходя к порядку K + 1, можно вычислять логарифмические производные

любого порядка со сложностью (46). Поскольку K здесь константа, общая

сложность вычисления не изменится, обеспечивая тем самым справедливость

утверждения

Теорема 6. Пусть F(z₀) =dK

log Γ(z₀); z₀ = α; α

алгебраическое

dz^K

число. Тогда

(

)

s_F(α)(n) = O

M (n) log² n

СПИСОК ЛИТЕРАТУРЫ

1. Dynkin E.B., Kolmogorov A.N., Kostrikin A.I., et. al. Six Lectures Delivered at the

International Congress of Mathematicians in Stockholm, 1962 (American Mathemat-

ical Society Translations-series 2). 31. Publ. by AMS, 1963.

2. Карацуба А.А. Сложность вычислений // Тр. МИАН. 1995. № 211. С. 186-202.

3. Карацуба А., Офман Ю. Умножение многозначных чисел на автоматах // Докл.

Академии наук СССР. 1962. Т. 145. № 2. С. 293-294.

4. Schönhage A., Strassen V. Schnelle Multiplikation grosser Zahlen // Computing.

1971. No. 7. P. 281-292.

5. Fürer M. Faster Integer Multiplication // SIAM J. Comput. 2009. V. 39. No. 3.

P. 979-1005.

6. Бендерский Ю.В. Быстрые вычисления // Докл. Академии наук СССР. 1975. Т.

223. № 5. C. 1041-1043.

7. Salamin Е. Computation of π using arithmetic-geometric mean // Math. Comp.

1976. V. 30. No. 135. P. 565-570.

8. Carlson B.C. Algorithms involving arithmetic and geometric means // Amer. Math.

Monthly. 1971. No. 78. P. 496-505.

9. Borwein J.M., Borwein P.B. Pi and the AGM-A Study in Analytic Number Theory

and Computational Complexity. New York: Wiley, 1987.

120

10. Карацуба Е.А. Быстрые вычисления трансцендентных функций // Пробл. пе-

редачи информ. 1991. Т. 27. № 4. С. 76-99.

11. Карацуба Е.А. О вычислении функции Бесселя путем суммирования рядов //

Сиб. журн. вычисл. мат. 2019. Т. 27. № 4. С. 76-99.

12. Lozier D.W., Olver F.W.J., Numerical Evaluation of Special Functions // Mathe-

matics of Computation 1943-1993: A Half -Century of Computational Mathemat-

ics, Gautschi W., eds., Proc. Sympos. Applied Mathematics. AMS. 1994. No. 48.

P. 79-125.

13. Siegel C.L. Transcendental numbers. Princeton: Princeton University Press, 1949.

14. Temme N.M. Special Functions. An Introduction to the Classical Functions of Math-

ematical Physics. New York: J. Wiley and Sons, 1996.

15. Karatsuba Е.А. On the computation of the Euler constant γ // Numerical Algo-

rithms. 2000. No. 24. P. 83-97.

16. Bach Е. The complexity of number-theoretic constants // Info. Proc. Letters. 1997.

No. 62. P. 145-152.

Статья представлена к публикации членом редколлегии А.А. Лазаревым.

Поступила в редакцию 01.02.2022

После доработки 28.04.2022

Принята к публикации 29.06.2022

121