БЫСТРЫЕ АЛГОРИТМЫ И МЕТОД БВЕ
Е.А. Карацуба
I. БЫСТРЫЕ АЛГОРИТМЫ.
Область вычислительной математики, которая называется быстрые алгоритмы, появилась в 1960 году.
Под алгоритмом мы будем понимать правило или способ вычисления, не формализуя это понятие. Далее будем считать, что числа записаны в двоичной системе счисления, знаки которой 0 и 1 называются битами.
Опр.1. Запись знаков 0 , 1 , плюс, минус, скобка; сложение, вычитание и умножение двух битов назовём одной элементарной или битовой операцией.
Быстрые алгоритмы — это область вычислительной математики, которая изучает алгоритмы вычисления заданной функции с заданной точностью с использованием как можно меньшего числа битовых операций. Тем самым, алгоритмы, которые можно назвать быстрыми, являются реальными алгоритмами. Такие алгоритмы, реализованные на ЭВМ в программном (а иногда и аппаратном) обеспечении позволяют существенно увеличить производительность работы компьютера, а иногда и решить задачи, размер которых не позволял найти решение путём применения обычных методов вычисления. Вопрос о размере задачи, которую можно решить за некоторое время с помощью данного компьютера, приводит нас к понятию сложности вычисления.
Первые постановки задач о битовой сложности вычислений (1956 г.) принадлежат А.Н. Колмогорову https://mech.math.msu.su/probab/Kolmogorov/kolmogorov.html [36], [37], который в то же время подчёркивал, что "цикл моих работ по теории информации был создан под большим влиянием публикаций Норберта Винера https://erudite.nm.ru/WienerNorbert.htm и Клода Шеннона https://www.astronet.ru/db/msg/1187395 ."
Прежде чем ввести понятие сложности вычисления, определим, что значит вычислить функцию в заданной точке. Рассмотрим наиболее простой пример вычисления вещественной функции y = f(x) вещественного переменного x , a ≤ x ≤ b .
Пусть f(x) на (a,b) удовлетворяет условию Липшица порядка α , 0 < α < 1 , так что при x1, x2 ∊ (a,b) :
|f(x1) - f(x2)| ≤ | x1 - x2 |α .
Пусть n — натуральное число.
Опр.2. Вычислить функцию y = f(x) в точке x = x0 ∊ (a,b) с точностью до n знаков, значит найти
такое число A, что
| f(x0) – A | ≤ 2⁻ⁿ.
Опр.3. Количество битовых операций, достаточное для вычисления функции f(x) в точке x = x0 с точностью до n знаков посредством данного алгоритма, называется сложностью вычисления f(x) в точке x = x0.
Таким образом, сложность вычисления f(x) в точке x = x0 есть функция n; а также f(x) и x = x0 . Эту функцию обозначают символом
sf(n) = sf,x0 (n).
Ясно, что sf зависит также от алгоритма вычисления и при разных алгоритмах будет разной. Сложность вычисления непосредственно связана со временем, затрачиваемым компьютером на это вычисление и потому иногда в литературе обозначается "временной" функцией T(n) [35].
Вопрос о поведении sf(n) при n → ∞ для класса функций или конкретной функции f, был впервые поставлен А.Н. Колмогоровым около 1956 года. Поскольку при вычислениях в первую очередь используются четыре арифметических действия: сложение, вычитание, умножение и деление, то прежде всего нужно знать количество битовых операций, достаточное для выполнения этих действий. Из определений 2 и 3 следует, что числа x0 и A можно представить в виде целой части и n двоичных знаков после запятой, т.е.
A = [A] + 0,ν1ν2 ν3 …νn,
x0 = [x0] + 0,μ1μ2μ3 …μn,
где νj , μj = 0 или 1 , j = 1, 2, … , n.
Так как целые части [A] , [x0] — фиксированные величины, а n → ∞ , то действия производятся по существу с n -значными числами. Отсюда прежде всего возникает вопрос о сложности вычисления суммы, разности, произведения и частного двух n -значных чисел a и b. Заметим, что деление (с остатком) сводится к сложению, вычитанию и умножению чисел.
Для записи чисел a и b требуется по меньшей мере 2n битовых операций. Следовательно, сложность сложения (вычитания) a и b не меньше 3n битовых операций, и не больше (если делать это обычным способом), чем 4n битовых операций. Таким образом, порядок количества битовых операций, необходимых и достаточных для выполнения сложения (вычитания), один и тот же.
В дальнейшем для краткости, будем называть "битовые операции" просто "операциями".
Следующим вопросом является вопрос о количестве операций, достаточных для вычисления произведения ab, или о сложности умножения. Функция сложности умножения получила специальное обозначение M(n).
Перемножая два n -значных числа обычным школьным способом "в столбик," мы при этом фактически складываем n n -значных чисел. Так что для сложности этого "школьного" или "обычного" метода мы имеем оценку сверху M(n) = O(n²).
В 1956 г. А.Н.Колмогоров высказал гипотезу, что нижняя оценка M(n) при любом методе умножения есть также величина порядка n² (так называемая "гипотеза n² Колмогорова"). На правдоподобность "гипотезы n²" указывал тот факт, что метод умножения "в столбик" известен не менее 4-х тысячелетий (например, этим методом пользовались шумеры), и если бы был более быстрый метод умножения, то он, вероятно, уже был бы найден.
В 1960 А.А. Карацуба http://www.mi.ras.ru/~karatsuba/ [19], [20], [21] (см. также [35]) нашёл новый метод умножения двух n -значных чисел с оценкой сложности
M(n) = O(nlog 23 ), log 2 3 = 1,5849… ,
и тем самым опроверг “гипотезу n².” Этот метод впоследствии был назван
"Divide and Conquer" ("Разделяй и властвуй"); другие, используемые в настоящее
время названия этого метода — "Binary Splitting", "Принцип Дихотомии" и т. п.
С момента появления "Divide and Conquer" начала развиваться теория быстрых вычислений.
Исследования с целью поиска алгоритма умножения со сложностью, близкой к оптимальной, были
продолжены рядом авторов (среди них были Тоом https://www.de.ufpe.br/~toom/
[48], Кук https://www.cs.toronto.edu/~sacook/
[15], Шёнхаге https://www.informatik.uni-bonn.de/~schoe/index.html
[43]), и в 1971 г. Шёнхаге
и Штрассен https://www.mathe.uni-konstanz.de/mitarbeiter/strassen.html
[44], [45] построили алгоритм с наилучшей на настоящее время оценкой для M(n),
M(n) = O(n log n log log n).
При построении этого алгоритма кроме "Divide and Conquer" они использовали идею выполнения
арифметических действий по модулю чисел Ферма 22ⁿ +1, а также быстрое преобразование Фурье.
На основе "Divide and Conquer" были построены алгоритмы быстрого умножения матриц. Первый
такой алгоритм в 1970 г. нашёл Штрассен [47], в дальнейшем эти исследования были продолжены
Коперсмитом https://www.research.ibm.com/people/c/coppersmith/, Виноградом
https://researchsmp2.cc.vt.edu/DB/db/indices/a-tree/w/Winograd:Shmuel.html [16], Паном
https://comet.lehman.cuny.edu/vpan/ [41] и др.
Приблизительно в это же время стали появляться и первые быстрые алгоритмы для вычисления
элементарных алгебраических функций. Как уже упоминалось выше, деление числа a на число b
с остатком, т.е. вычисление чисел q и r в формуле
a = qb + r , 0 ≤ r < b,
сводится к сложению, вычитанию и умножению, и если a и b являются n -значными числами,
то сложность деления O(M(n)).
Самый простой алгоритм деления заключается в вычислении методом Ньютона обратной величины
1/b с точностью до n знаков с последующим умножением на a по алгоритму быстрого умножения.
Пусть надо вычислить x = 1/b, где 1/2 ≤ b ≤ 1 (если это не так, то умножая или деля b на 2N
получим выполнение указанного условия). Возьмём в качестве начального приближения x0 = 3/2 и
вычисляем по формуле
xn+1 = 2xn - bxn².
Этот метод обеспечивает достаточно быструю сходимость к 1/b, поскольку из равенства xn = 1/b - ε,
следует, что xn+1 = 1/b - bε².
Для извлечения корня k -ой степени из числа, т.е. для вычисления a1/k можно также воспользоваться
методом Ньютона. Пусть k ≥ 2 — целое число и надо вычислить x = a1/k, где a ≥ (k+1)k. Вычисляем x
по формуле:
xn+1 = xn(k+1)/k - xnk⁺¹/(ka).
Заметим, что умножением или делением на 2kN можно всегда добиться выполнения для a условия
a ≥ (k+1)k.(иногда это условие оказывается лишним, т.е. указанный процесс быстро сходится и при
других a).
За начальное приближение можно взять, например, одно из двух чисел d или d+1, где
d — целое число, удовлетворяющее условию dk < a ≤ (d+1)k. Тогда одно из чисел |d - z1/k| или
|d + 1 - z1/k| не превосходит 1/2.
Сложность вычисления этим методом будет O(M(n)).
Первые быстрые алгоритмы для вычисления элементарных алгебраических функций, основанные
на методе Ньютона, появились в работах Кука [15], Бендерского https://www.bisguard.com/ [8] и Брента
http://web.comlab.ox.ac.uk/oucl/people/richard.brent.html [11] (см. также [35]).
Используя метод Ньютона можно доказать следующую теорему [11] :
Теорема. Если уравнение f(x) = 0 имеет простой корень ζ ≠ 0, f(x) непрерывна по Липшицу
в окрестности ζ , и мы можем вычислить f(x) с точностью до n знаков за O(M(n)j(n)) операций ,
где j(n) — положительная , монотонно возрастающая функция, для любого x из окрестности ζ ,
то ζ можно вычислить с точностью до n знаков также за O(M(n)j(n)) операций.
В 1976 г. Саламин https://www.ams.org/mathscinet-getitem?mr=53+%237928 [41] и Брент [11]
предложили первый алгоритм быстрого вычисления константы π, основанный на АГС -методе Гаусса
(см., например, [13], [14]). Используя также восходящие преобразования Ландена
http://www-groups.dcs.st-and.ac.uk/~history/Mathematicians/Landen.html [38], Брент [11]
предложил первые АГС -алгоритмы для быстрого вычисления простейших трансцендентных функций. В
дальнейшем исследование и использование АГС –алгоритмов было продолжено многими авторами,
среди которых были братья Джонатан https://www.cs.dal.ca/~jborwein/
и Питер https://www.cecm.sfu.ca/~pborwein/ Борвайны, написавшие книгу "The Pi and the AGM"
("Пи и АГС") [9], где собрано наибольшее число АГС -алгоритмов.
Кроме алгоритмов вычисления простейших трансцендентных функций, в этой книге содержатся также
алгоритмы для вычисления некоторых высших трансцендентных функций (гамма-функции Эйлера,
например). Чтобы вычислить такие функции с точностью до n знаков с помощью описанных в книге
алгоритмов требуется
O(n3/2 log3 n log log n)
операций.
О других результатах, связанных с "Divide and Conquer" и быстрыми вычислениями, см. [1], [2], [3], [6],
[7], [17], [18], [35], [40].
II. МЕТОД БВЕ.
Введение. Определение быстрого алгоритма.
Будем называть "быстрым" такой алгоритм вычисления функции f = f(x), что для этого алгоритма
sf(n) = O(M(n) logc n) ,
где c есть константа. Здесь и далее предполагаем, что для M(n) справедлива оценка Шёнхаге-Штрассена
. Ясно, что для любого алгоритма и любой функции f, выполняется неравентство:
sf(n) > n .
(Заметим, что только запись числа с точностью до n знаков требует не меньше, чем n+1 операций.)
Следовательно, для быстрого алгоритма
n < sf(n) < c1 n logc+1n loglog n < n1+ε
для любого ε > 0 и n > n1(ε). Тем самым, быстрым алгоритмам соответствует правильный порядок
оценки сверху sf(n) по n , n → ∞ .
Метод БВЕ (от 1990) [22]—[34] (см. также [4]—[6], [39]) является вторым после метода АГС
известным на настоящее время методом быстрого вычисления классических констант e и π, и
простейших трансцендентных функций. Но, в отличие от АГС, метод БВЕ является также методом
быстрого вычисления высших трансцендентных функций. В настоящее время БВЕ является единственным
методом, позволяющим вычислить быстро такие высшие трансцендентные функции, как гамма-функцию
Эйлера, гипергеометрические функции, цилиндрические, сферические и т.д. функции при алгебраических
значениях аргумента и параметров, дзета-функцию Римана и дзета- функцию Гурвица при целых
значениях аргумента и алгебраических значениях параметра.
2. Вычисление константы e.
Основная идея метода БВЕ проще всего объясняется на примере вычисления классической
постоянной e. Рассмотрим сначала алгоритм вычисления n! за log n шагов со сложностью
O(n log3 n log log n). Для простоты мы предполагаем, что n = 2k , k ≥ 1.
1-й шаг. Вычисляются n/2 произведений вида:
a1(1) = n(n-1) , a2(1) = (n-2)(n-3) , … , an/2(1) = 2*1 ;
2-й шаг. Вычисляются n/4 произведений вида:
a1(2) = a1(1)a2(1) , a2(2) = a3(1)a4(1), … , an/4(2) = an/2(1)an/2-1(1) ,
… и так далее.
k -й, последний шаг. Вычисляется одно произведение:
a1(k) =a1(k-1)a2(k-1) ,
и это даёт результат: n!
Для вычисления константы e возьмём m = 2k , k ≥ 1 , членов ряда Тейлора для e,
e = 1 + 1/1! + 1/2! + … + 1/(m-1)! + Rm .
При этом выбираем m так, чтобы для остатка Rm выполнялось неравенство Rm ≤ 2⁻ⁿ⁻¹. Это будет,
например, когда m ≥ 4n/log n. Таким образом, возьмем m = 2k таким, что натуральное число k
определяется неравенствами: 2k ≥ 4n/log n > 2k⁻¹.
Будем вычислять сумму
S = 1 + 1/1! + 1/2! + … + 1/(m-1)! =
за k шагов следующего процесса.
1-й шаг. Объединяя слагаемые S последовательно попарно и вынося за скобки "очевидный" общий
множитель, получаем
S = (1/(m-1)! + 1/(m-2)!) + (1/(m-3)! + 1/(m-4)!) + … = (1/(m-1)!)(1+m-1) + (1/(m-3)!)(1+m-3) + …
Будем вычислять только целые значения выражений, стоящих в скобках, т.е. значения m , m-2 , m-4, … .
Таким образом, на первом шаге сумма S преобразуется к виду
S = S(1) ;
m1= m/2 , m = 2k , k ≥ 1. На первом шаге вычисляются m/2 целых чисел вида
αm1-j(1) = m - 2j , j = 0, 1, … , m1-1 .
Далее мы действуем аналогично: объединяя на каждом шаге слагаемые суммы S последовательно
попарно, мы выносим за скобки "очевидный" общий множитель и вычисляем только целые значения
выражений в скобках. Пусть сделано i шагов такого процесса.
i+1 -й (i+1 ≤ k) шаг.
S = S(i+1);
mi+1 = mi /2 = m/2i+1 , вычисляем только m/2i+1 целых чисел вида
j = 0 , 1 , … , mi+1-1, m = 2k , k ≥ i+1 . Здесь ((m-1-2i+1j)!/(m-1-2i-2i+1j)!) есть произведение 2i целых
чисел.
… и так далее.
Последний, k -й шаг. Вычисляем одно целое значение α1(k), вычисляем, пользуясь вышеописанным
быстрым алгоритмом, значение (m-1)!, и производим одно деление целого числа α1(k) на число (m-1)!
с точностью до n знаков. Получившийся результат и есть сумма S, или константа e с точностью до
2⁻ⁿ. Сложность всех вычислений есть
O(M(m) log2 m) = O(M(n) log n) .
3. Метод БВЕ и быстрое суммирование рядов.
Метод получил название "БВЕ" (Быстрое Вычисление E -функций) потому, что даёт возможность
вычислять быстро E -функции Зигеля, и в частности, ex.
Зигель https://www-groups.dcs.st-and.ac.uk/~history/Mathematicians/Siegel.html [46] назвал "E -функциями"
класс функций,"похожих на экспоненциальную". К ним принадлежат такие высшие трансцендентные
функции как гипергеометрические, сферические, цилиндрические функции и т.д.
С помощью БВЕ можно вычислить быстро любой подходящий степенной ряд. Например, для
быстрого вычисления константы π, можно воспользоваться формулой Эйлера
π/4 = arctg 1/2 + arctg 1/3,
а метод БВЕ применить к суммированию рядов Тейлора для
arctg 1/2 = 1/(1*2) – 1/(3*2³) + … + (-1)r-1/((2r-1)22r-1) + R1 ,
arctg 1/3 = 1/(1*3) – 1/(3*3³) + … + (-1)r-1/((2r-1)32r-1) + R2 ,
с остаточными членами R1, R2, для которых справедливы оценки
| R1| ≤ 4/5/(2r+1)/22r+1;
| R2| ≤ 9/10/(2r+1)/32r+1;
и при r = n, 4(|R1|+| R1|) < 2⁻ⁿ. Чтобы вычислить π с помощью БВЕ можно воспользоваться также и
другими приближениями, например из [5].
Чтобы вычислить постоянную Эйлера гамма с точностью до n знаков, нужно просуммировать с
помощью БВЕ два ряда. А именно, при m = 6n, k = n, k ≥ 1,
При этом
sγ(n) = O(M(n) log2n) .
Для быстрого вычисления константы γ можно также применить метод БВЕ к приближению из [12].
С помощью БВЕ можно вычислить быстро следующие два вида рядов:
при условии, что a(j) , b(j) — целые числа, |a(j)|+|b(j)| ≤ (Cj)K; |z| < 1; K и C есть константы, и z есть
алгебраическое число. Сложность вычисления рядов (1), (2) с точностью до n знаков есть
следующая страница>