БЫСТРЫЕ═
АЛГОРИТМЫ═ И═ МЕТОД═
БВЕ
Е.А. Карацуба
═
══════ Область вычислительной математики, которая называется быстрые алгоритмы, появилась в 1960 году.
══════ Под алгоритмом мы будем понимать правило или способ вычисления, не формализуя это понятие. Далее будем считать, что числа записаны в двоичной системе счисления, знаки которой═ 0═ и═ 1═ называются битами.
Опр.1.═ Запись знаков══ 0 ,
1 , плюс, минус, скобка; сложение, вычитание и умножение двух битов
назов╦м одной элементарной или битовой операцией.
Быстрые алгоритмы═ ≈═ это область вычислительной математики, которая изучает алгоритмы вычисления заданной функции с заданной точностью с использованием как можно меньшего числа битовых операций. Тем самым, алгоритмы, которые можно назвать быстрыми, являются реальными алгоритмами. Такие алгоритмы, реализованные на ЭВМ в программном (а иногда и аппаратном) обеспечении позволяют существенно увеличить производительность работы компьютера, а иногда и решить задачи, размер которых не позволял найти решение пут╦м применения обычных методов вычисления. Вопрос о размере задачи, которую можно решить за некоторое время с помощью данного компьютера, приводит нас к понятию сложности вычисления.
══════ Первые постановки задач о битовой сложности
вычислений (
══════ Прежде чем ввести понятие сложности вычисления, определим, что значит вычислить функцию в заданной точке. Рассмотрим наиболее простой пример вычисления вещественной функции═ 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²).
══════ В
══════ В 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" начала развиваться теория быстрых вычислений.
═════ Исследования с целью поиска алгоритма умножения со сложностью, близкой к оптимальной, были
═════ продолжены рядом авторов (среди них были Тоом
http://www.de.ufpe.br/~toom/
═════ [48], Кук http://www.cs.toronto.edu/~sacook/
═════ [15], Ш╦нхаге http://www.informatik.uni-bonn.de/~schoe/index.html
═════ [43]), и в
═════ и Штрассен http://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" были построены алгоритмы быстрого умножения матриц. Первый
═════ такой алгоритм в
═════ Коперсмитом
http://www.research.ibm.com/people/c/coppersmith/,
Виноградом
═════ http://researchsmp2.cc.vt.edu/DB/db/indices/a-tree/w/Winograd:Shmuel.html [16], Паном
═════ http://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], Бендерского http://www.bisguard.com/ [8] и Брента
══════ http://web.comlab.ox.ac.uk/oucl/people/richard.brent.html [11] (см. также [35]).
══════ Используя метод Ньютона можно доказать следующую теорему [11] :
═════ Теорема. Если уравнение f(x) = c имеет простой корень═ ζ
≠═ 0,═ f(x)═
непрерывна по Липшицу
═════ в окрестности═ ζ
, и мы можем вычислить═ f(x)═ с
точностью до═ n═ знаков за═ O(M(n)j(n)) операций ,
═════ где══
j(n)═ ≈═ положительная
, монотонно возрастающая функция, для любого═ x═ из
окрестности ζ ,
═════ то═
ζ можно вычислить с точностью до═
n═ знаков
также за═ O(M(n)j(n))═ операций.
════════════ В
═════ предложили первый алгоритм быстрого вычисления константы π, основанный на АГС -методе Гаусса
═════ (см., например, [13], [14]).═ Используя также восходящие преобразования Ландена ═
═════ http://www-groups.dcs.st-and.ac.uk/~history/Mathematicians/Landen.html═ [38],═ Брент═ [11]
═════ предложил первые АГС -алгоритмы для быстрого вычисления простейших трансцендентных функций. В
═════ дальнейшем исследование и использование АГС √алгоритмов было продолжено многими авторами,
═════ среди которых были братья Джонатан═ http://www.cs.dal.ca/~jborwein/
══ и Питер═ http://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-
═════ чисел.
═════ ┘ и так далее.
═════ Последний,═ 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.
════════════ Зигель═ http://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 ═знаков есть
sf1(n)═ =═ O(M(n) log2n) ,
sf2(n)══ =═ O(M(n) log n) .
4. Вычисления экспоненциальной функции методом БВЕ.
════════════ Пусть═ y = f(x) = ex. Вычислим функцию y = ex в точке═ x = x0═ с точностью═ 2⁻ⁿ, предполагая сначала,
═════ что
0═ < x0 < 1/4 ,══ n > 8.═══ (3)
═════ Возьм╦м целое число═ k,═ удовлетворяющее условиям══ 2k⁻¹ <═ n═ ≤═ 2k══ и положим═ N = 2k+¹ (заметим, что
═════ 2n ≤ N < 4n). Пусть═ xN =═ α32⁻³
+ α42⁻4══ +
┘ +═
αN2⁻N,═ где═ αj = 0═ или═ 1,══ j = 3, 4, ┘ , N,═ и
| x0 - xN |══ ≤═ 2⁻N.
═════ Тогда
ex0═ =══ exN + 4θ1═ 2⁻N══ ;═══════ | θ1═ |═ ≤ 1 .
Будем═ вычислять═ exN.═══ Представим═ xN ═в виде
xN══ =══ β2 /24═
+═ β3 /28 +═
β4 /216 +
┘ +═ βk+1/2N
,
═════ где
β2══ = α3α4 ;═══ β3 =═ α5α6α7α8 ;═ ┘ , ══βk+1 =═ αN √ 2k + 1αN √ 2k + 2══ ┘═ αN √ 1 αN ═.
═════ Здесь═ βν═ ═есть═ 2ν⁻1 -значное целое число, прич╦м═ 2═ ≤═ ν═ ≤═ k+1═ и═ k+1 =═ log N.═ Таким образом, число
═════ exN═ можно представить в виде произведения
═════ Разложим каждый множитель этого произведения в ряд Тейлора:
═════ r = N 2⁻ ν+¹.
═════ Для остаточного члена ряда═ Rν(r)═ справедливы следующее неравенство:
═════ Отсюда═ Rν(r) ═<═ 2⁻N .
═════ Следовательно, (4) можно записать в виде
═════ где
и постоянные═ θν(r),══ 0 <═ θν(r) < 1 . Легко видеть, что сумма (5) является подходящей для е╦ суммирования с помощью БВЕ. На последнем шаге процесса в (5) имеем
а затем производим деление целого числа═ aν ══на целое═ bν═ с точностью до═ N ═знаков. Тем самым, получаем такое═ ην═ , что
ξν = aν /bν = ην═══ +═ θ′ν 2⁻N═ ,═══ | θ′ν | < 1 .
Следовательно═ exN═ может быть представлено в виде
| θν |═ ≤ 1 .
Определим целое число═ l ═неравенствами═ 2l⁻¹═ <═ k═ ≤ ══2l. Предполагая, что═ ην═ =1, θν = 0 для ═ν > k,═ мы можем записать последнее выражение в виде
| θν |═ ≤ 1 .
Ясно, что══ log k═ ≤═ l ═< log k + 1. Последнее произведение вычисляется за══ l═ шагов процесса, который аналогичен вычислению══ n! . Перемножим множители═ (6)═ последовательно попарно, вычисляя на 1-ом шаге══ 2l⁻¹═ произведений вида
(ην═ + 2 θν 2⁻N) (ην +1═ + 2 θν+1 2⁻N)══ =═ ην ην +1═ + 2 θν 2⁻Nην +1 + 2 θν+1 2⁻N ην + 4θνθν+1 2⁻2N══ = ην ,ν +1═ + 8 θν,ν+1 2⁻N,
ν═ ≡═
0(2)═ ,══ ν = 2, 3,═ ┘ , 2l + 1,
где ην ,ν +1══ ═есть═ N -значное число, такое что
ην ,ν +1═ = ην ην +1═ + θ̃ 2⁻N ,══ | θ̃ |═ ≤ 1 , | θν ,ν+1|═ ≤═ 1.
Здесь мы учитываем, что═ | ην | <═ exN ═<══ ex0═ <═ 3/2 .
═════ После 1-го шага такого процесса произведение (6) принимает вид:
══════ ν ≡
0(2)
══════ | θν ,ν+1|═ ≤═ 1.
═════ Далее действуем аналогичным образом.
═════ После 2-го шага такого процесса (6) принимает вид:
═════ ν ≡ 2(4)
═════ | θν,ν+1,ν+2,ν+3|═ ≤═ 1.
═════ ┘ и так далее.
════ Через═ l═ шагов получим
═════ | θ2, ┘, 2l + 1═ |═ ≤═
1.
═════ Поскольку══ l <═ log k +1,══ k═ =═ log N - 1,═ 2n═ ≤═ N < 4n,═ то из последнего соотношения находим
ex0═ =═ η ═+═ θ 2⁻n ,
═════ где═ η═ =═ η2, ┘, 2l + 1═ ,═ | θ |═ ≤ 1.
═════ Сложность вычисления══ ex0══ ═с помощью такого процесса есть
O(M(n) log2n) .
═════ Замечание 1. Мы рассмотрели случай═ 0
<═ x0 ═< 1/4. Если═
x0 ═≥ 1/4, возьм╦м наименьшее целое
═════ число═ m , такое что═ m > 4x0.═ В этом случае══ 0 < x0/m < 1/4 ,══ ex0═ = (ex0/ m)m.══ Тем самым
═════ вычисление══ ex0 ═сводится к вычислению═ ex0/ m ═, 0 < x0 /m < 1/4,═ после чего результат
═════ возводится в степень═ m.═ Это не
ухудшает оценку сложности вычисления, поскольку═
m═ есть
═════ фиксированное число.
5. Вычисление тригонометрических функций методом БВЕ.
════════════
Вычисление══ тригонометрических
функций═ y = f(x) =═ sin x═ и═ y = f(x) = cos x═ в точке═
x =═ x0
═════ может быть сведено к вычислению реальной и мнимой частей══ eix0.══ Число══ xN══ ═представим
═════ в═ том же виде, что и в предыдущем случае, при
этом════ eix0
= eixN ══+═ 4θ1═ 2⁻N ,
═════ | θ1═ |═ ≤ 1 .
═════ Запишем
═══════════════════════════
═════ в виде
═════ где
0
<═ θν(r)═ <═ 1.═
Далее мы производим с суммами (8), (9) те же действия, что и с суммой
(4). Пусть ═ην,
fν ═есть═ N═ -значные числа, такие что
ξν =═ a′ν /bν = ην═══ +═ θ′ν 2⁻N══ ,══ | θ′ν |═ ≤ 1,
ψν =═ a″ν /dν ══= fν ═══+═ θ″ν 2⁻N═ ,══ | θ″ν |═ ≤ 1,
,
═════ Тогда (7) может быть представлено в виде
═════ | θν |═ ≤ 1 .
════════════ Вычисляем произведение (10) посредством процесса, описанного в предыдущем параграфе.
═════ В результате со специально выбранными═ ═l , N , r═══ находим
eix0═
=═ f══ +═ i η ══+═ 2θ 2⁻N══ ,══ | θ |═ ≤ 1 .
═════ и отсюда
cos x0═
=══ f═ +═ θ1 2⁻n══ ,══ | θ1|═
≤ 1 ,
sin x0═
=══ η═ +═ θ2 2⁻n══ ,══ | θ2 |═ ≤ 1 .
═════ Сложность вычисления та же, что и для экспоненциальной функции:
O(M(n) log2n).
═════ В═ [15]═ доказана следующая общая теорема:
═════ Теорема.═ Пусть
y = f(x) ≈ простейшая трансцендентная функция, т.е. экспоненциальная функция или
═════ тригонометрическая функция, или
элементарная алгебраическая функция, или их суперпозиция, или
═════ обратная им функция или суперпозиция
обратных функций. Тогда
sf(n)═
=═ O(M(n) log2n)
.
6. Быстрое вычисление гамма-функции Эйлера при рациональном аргументе.
════════════ Гамма функция Эйлера═ y =═ Γ(s)═ относится к тем высшим трансцендентным функциям, которые не
═════ будучи E- функциями Зигеля, могут быть, тем не менее, быстро вычислены с помощью метода БВЕ. Есть
═ ════два алгоритма для вычисления══ Γ(s): 1) для вычисления═ Γ(s)═ при рациональных значениях аргумента═ s;
═════ 2) для вычисления══ Γ(s)═ при алгебраических значениях═ s. Первый из этих алгоритмов является более
═════ "простым", чем второй. Вычислим функцию═ y =═ Γ(x) ,═ при═ x = x0 =═ a/b ;═ (a,b) = 1,═ предполагая сначала,
═════ что═ 0═ <═ x0═ <═ 1 , т.е. что═ 0═ <═ a═ <═ b . Воспользуемся представлением гамма-функции в виде интеграла
═════ Легко видеть, что при═ p = n,══ 0═ <═ x0═ <═ 1 ,
═,═══ 0═
<═ θ1═
<══ 3/4.
═════ Следовательно,
═════ 0═ <═ θ1═ <══ 3/4.
═════ Предполагая, что═ r═ ≥═ 3n , разложим подинтегральную функцию в (11)═ y = e⁻t ,═ в ряд Тейлора
═════ по степеням ═t:
═════ где══ Rr =═
θr tr+1/(r+1)!,═ | θr
|══ ≤═ 1.═ Отсюда имеем при═ r+1 =
6n
Γ(x0)═ =═ nx0 S +═ θ 2⁻n ,═ | θ |═ ≤ 1 , (12)
═════ где
и мы будем вычислять═══ Γ(x0)═ =═ Γ(a/b)═ ═по формулам (12)--(13). Заметим, что для вычисления══ nx0 = na/b,
═════ с═ точностью до══ n ═знаков достаточно═ O(M(n))═ операций (методом Ньютона). Чтобы вычислить═ S
═════ возьмем═ r+1 =═ 2k,═ 2k⁻1═ <═ 6n══ ≤═ 2k,══ k ≥ 1;═ членов ряда (13). Пусть
═════ Тогда (13) можно записать в виде
S = S1(0) + S2(0) + ┘ + Sr+1(0) .
═════ Вычислим сумму═ S ═с помощью БВЕ за══ k══ шагов, объединяя на каждом шаге слагаемые═ S
═════ последовательно попарно и вынося за скобки "очевидный" общий множитель. Однако, в отличие от
═════ вычисления константы═ e══ методом═ БВЕ, теперь в скобках будут находиться не целые числа, а дроби.
═══ ══Поэтому теперь на каждом шаге будут вычисляться целый числитель и целый знаменатель дроби, стоящей
═════ в скобках (деление не производится до последнего шага). На первом шаге имеем
S = S1(1) + S2(1) + ┘ + Sr1(1) ,═ r1═ ═= (r+1)/2
Sr1-ν(1) =═ Sr+1-2ν(0) + Sr-2ν(0) = (-1)r⁻1n r⁻ 2 ν⁻1/(r-2ν)! (pr1-ν(1)/qr1-ν(1)) .
═════ На 1-ом шаге вычисляем целые числа
pr1-ν(1) =═ - bn(br-2bν
-b+a) + (br-2bν)(br-2bν +a) ;
qr1-ν(1) ═=
(br-2bν +a)(br-2bν -b+a) ;
ν = 0 , 1 , ┘ ,
(r+1)/2-1 ; r+1 = 2k , k ≥ 1.
═════ На═ j -ом шаге (j ≤ k ) имеем
S = S1(j) + S2(j) + ┘ + Srj(j) ,═ rj═ ═= (r+1)/2j
═════ где═ Srj-ν (j); ν = 0 , 1 , ┘ , rj -1 ; определяются равенствами
═════ На═ j -ом шаге (j ≤ k) вычисляем целые числа
══════════
════════════════════════════════ ═════
═════ ν = 0 , 1 , ┘ , rj -1 ;═ rj═ = (r+1)/ 2j.═
═════ И так далее.
═════ На k-ом (последнем) шаге вычисляем целые числа══ prk(k) = p1(k),═ qrk(k) = q1(k), ═r! и производим одно
═════ деление с точностью до═ 2⁻2n═ ═по формуле
S = S1 (k) =═ prk(k)/qrk(k)/r! ,
═════ что да╦т сумму═ S ══с точностью═ 2⁻n⁻1 . Следовательно значение══ Γ(x0) =═ Γ(a/b) ═вычислено с точностью
═════ 2⁻n+1 , n═ ≥ 1 . Сложность такого вычисления есть
s Γ (n)═ =═
O(M(n) log2n) = O(n log3n loglog n).
═════ Замечание 2. Мы вычислили═ Γ(x0),═
x0 =═ a/b ,═ предполагая,
что═ 0══ <═ x0 ═<═ 1 .
Если═
x0═ ≥ 1,═ то пользуясь
═════ соотношением═ Γ(x+1)
= x Γ(x),═ сводим вычисление в нужной точке к вычислению в точке из интервала
═════ (0,1).
7. Быстрое вычисление гамма-функции Эйлера при алгебраическом
аргументе.
════════════ Вычислим═ Γ(x0), ═при условии, что x0 есть алгебраическое число. В соответствии с Замечанием 2,
═════ свед╦м это вычисление к вычислению Γ(α) , 0═ < α < 1 , где═ α═ есть алгебраическое число степени═ t,
═════ t═ ≥ 2. Для простоты, рассмотрим случай, когда═ α ═является вещественным числом. Предполагается, что
═════ известен многочлен═ g(x), наименьшей степени с целыми коэффициентами, корнем которого является число
α, т.е.
g(x) = gt xt + gt-1xt⁻1 + ┘ + g1x + g0 ;════════════ g(α) = 0 ;═══ (14)
═════ gt ,gt-1, ┘,g0 --- целые числа,═ t═ ≥ 2.══ Как и в случае рационального аргумента воспользуемся формулами
═════ (12), (13). Отметим, что вычисление══ nα ════с точностью═ 2⁻2n══ ═по формуле═ nα = eα log n ═требует
═════ O(M(n) log2n) ═операций. Сумма═ S═ из (12) принимает вид
═════ т.е.
S = S1(0) + S2(0) + ┘ + Sr+1(0)
,═ (16)
═════ где
════════════ Вычисление═ S══ выполняется за═ k═ шагов, r+1 =═ 2k,═ 2k⁻1═ <═ 6n══ ≤ ═2k,══ k ≥ 1.═ Особенностью этого
═════ алгоритма является использование основного свойства алгебраических чисел для═ ограничения роста
═════ сложности вычисления. До шага═ i,═ где═ i═ определяется условиями═ 1═ ≤ i═ ≤═ k ,═ 2i⁻1 <═ t═ ≤═ 2i,═ производим
═════ с суммой (16), (17) действия, аналогичные вышеописанным для случая рационального аргумента. Однако,
═════ группируя слагаемые суммы (16) таким же образом как и в предыдущем параграфе, теперь мы вычисляем
═════ на каждом шаге не числитель и знаменатель дробей, находящихся в скобках, а лишь целые коэффициенты
═════ при степенях═ α, многочленов, находящихся в числителе и знаменателе дробей "из скобок".
═════ На 1-ом шаге в скобках находятся дроби
═════ вычисляются числа
════════════════════════════ ════════════════════════════
════════════════════════════════════════════════════════
════════════════════════════════════════════════════════
═
════════════════════════════════════════════════════════
═════ ν = 0 , 1 , ┘ , (r+1)/2-1.
═════ На═ i -ом шаге в скобках находятся дроби
═════ где
═════ вычисляются числа
═════════════════════════════════════
═════ 0 ≤ m 1 ≤═ 2i⁻1-1 ,═ 0 ≤ m 2 ≤═ 2i⁻1.
════════════════
═════ 0 ≤ l 1 ≤═ 2i⁻1 ,═ 0 ≤ l 2 ≤═ 2i⁻1.
════
═════ Заметим. что умножение на степень═ α ═и вычисление═ δ и═ ρ, так же как и деление═ δ на ρ, не
═════ производится до последнего шага.
════════════ Перед═ i+1 -ым шагом═ (1═ ≤ i═ ≤═ k ,═ 2i⁻1 <═ t═ ≤═ 2i)═ мы редуцируем многочлены (19) по модулю
═════ многочлена══ g(x)═ при═ x = α.═ Таким образом, если
═══════════════════════════════════
═════ где═ pt-1, ┘,p0, qt ,┘,q0═ целые, t = 2i, то мы делим═ P(x)═ и═ Q(x)═ на═ g(x)═ с остатками:
P(x) = g(x)P0(x) + P1(x) ,═ Q(x) = g(x)Q0(x) + Q1(x) ,
═════ где═ P1(x)═ и═ Q1(x)═ есть многочлены с рациональными коэффициентами, причем степени═ P1(x)═ и
═════ Q1(x)═ не превышают═ t-1.═ Отсюда и из (14) получаем
P(α) = P1(α), Q(α) = Q1(α),
═════ то есть в (18) имеем
βri-ν(i) = P1(α)/Q1(α) .═══ (20)
═════ Домножая, если нужно, числитель и знаменатель дроби (20) на целый общий множитель, получаем в
═════ числителе и знаменателе дробей "из скобок" многочлены с целыми коэффициентами степени не большей,
═════ чем═ t-1.
════════════ Начиная с шага═ i,═ на каждом шаге═ i, i+1, ┘, k, производится редукция многочленов от═ α,
═════ расположенных в числителе и знаменателе βμ(j) ═по модулю═ g(x).═ Поскольку коэффициенты═ g(x) , так же
═════ как и степень═ g(x)═ являются абсолютными постоянными, то значения этих коэффициентов, участвующие
═════ в вышеописанных редукциях могут возрасти лишь не более чем в═ gt══ раз, где
g =═ max 0≤i ≤t |g i|,
═════ что не влияет на оценку сложности проводимых вычислений. На последнем,═ k -ом, шаге такого процесса
═════ получаем
═════ где
═════ p, q═ ≤ i-1 , и вычисляем═ δrk(k),═ ρrk(k) ═и═ S,═ и следовательно═ Γ(α) ═с точностью═ 2⁻n+1. Сложность
═════ вычисления
s Γ (n)═ =═ O(M(n) log2n)
8. Заключение.
════════════ Представленный метод БВЕ, метод быстрого суммирования специального вида рядов, позволяет
═════ вычислить любую элементарную трансцендентную функцию для любого аргумента, классические
═════ константы═ e,══ π, постоянную Эйлера══ γ, постоянные Апери и Каталана, такие высшие трансцендентные
═════ функции, как гамма-функцию Эйлера, гипергеометрические функции, сферические функции,
═════ цилиндрические функции и т.д. для алгебраических значений аргумента и параметров, дзета-функцию
═════ Римана для целых значений аргумента, дзета-функцию Гурвица для целого аргумента и алгебраических
═════ значений параметра, а также такие специальные интегралы, как интеграл вероятности, интегралы
═════ Френеля, интегральную экспоненциальную функцию, интегральные синус и косинус и т.д. при
═════ Алгебраических значениях аргумента с оценкой сложности вычисления, близкой к оптимальной, а именно
sf(n)═
=═ O(M(n) log2n)
.
══════В настоящее время только метод БВЕ да╦т возможность быстро вычислять значения функций из класса
═════ высших трансцендентных функций, некоторые специальные интегралы математической физики и такие
═════ классические константы, как константы Эйлера, Каталана и Апери. Дополнительным преимуществом
═════ метода БВЕ является возможность распараллеливания основанных на БВЕ алгоритмов.
Список
литературы
[1] В.Б.Алексеев, От метода Карацубы для быстрого умножения чисел к быстрым алгоритмам для
дискретных функций. Труды Математического Института им. В.А.Стеклова, т. 218, с.20-27 (1997).
[2] В.Б.Алексеев, С.А.Ложкин,═ Элементы теории графов, схем и автоматов. Изд. ВМиК МГУ, 58 с.,
Москва (2000).
[3] Г.И.Архипов,
В.А.Садовничий, В.Н.Чубариков, Лекции по математическому анализу. Изд. Высшая
школа, 695 с., Москва (1999).
[4] E.Bach,═ The complexity of number-theoretic constants.
Info.Proc.Letters, N 62, pp.145-152 (1997).
[5] D.H.Bailey,
P.B.Borwein and S.Plouffe,═ On the
rapid computation of various polylogarithmic constants. Math. Comp.,v. 66,
pp.903-913 (1997).
[6] J.C.Bajard, J.M.Muller,═ Calcul et arithmйtique des ordinateurs. Hermes Science, 226 pp. (2004).
[7] Л.А.Бассалыго,═ Замечание о быстром умножении многочленов над полями Галуа. Пробл.Пер.Инф., N 1, pp.101-102 (1978).
[8] Ю.В.Бендерский,═ Быстрые вычисления. Доклады Академии Наук СССР, т.223, N 5, с.1041-1043 (1975).
[9] J.M.Borwein
and P.B.Borwein,═ Pi and the AGM.
Wiley, 414 pp.,
[10] J.M.Borwein,
D.M.Bradley and R.E.Crandall,═
Computational strategies for the Riemann zeta function. J. of Comput.
Appl. Math., v. 121 ,N 1-2 , pp.247-296 (2000).
[11] R.P.Brent,
Fast Multiple-Precision Evaluation of Elementary Functions. ACM, v.23, N 2,
pp.242-251 (1976).
[12] R.P.Brent
and E.M.McMillan,═ Some new
algorithms for high-precision computation of Euler's constant. Math. Comp.,
vol.34, pp.305-312 (1980).
[13] B.C.Carlson,
Algorithms involving arithmetic and geometric means. Amer. Math. Monthly, v.78,
pp.496-505 (1971).
[14] B.C.Carlson,
An algorithm for computing logarithms and arctangents. Math.Comp., v.26,
pp.543-549 (1972).
[15] S.A.Cook,
On the minimum computation time of functions. Thesis,
[16] D.Coppersmith and S.Winograd, On the asymptotic complexity of matrix multiplication. SIAM Journal on Computing, v.11, N 3, pp.472-492 (1982).
[17] С.Б.Гашков,═ О сложности интегрирования рациональных дробей. Труды Математического Института им. В.А.Стеклова, т. 218, с.122-133, (1997).
[18] С.Б.Гашков, В.Н.Чубариков, Арифметика. Алгоритмы. Сложность вычислений. Изд. Высшая школа, 2-е изд., 320 с., Москва (2000).
[19] А.Карацуба и Ю.Офман, Умножение многозначных чисел на автоматах. Доклады Академии Наук СССР т.145, 2, с.293-294(1962).
[20] A.Karacuba, Berechnungen und die Kompliziertheit von Beziehungen. EIK, N 11, s.10-12 (1975).
[21] А.А.Карацуба, Сложность вычислений. Труды Математического института им. Стеклова, т.211, с.169-183 (1995).
[22] Е.А.Карацуба, Быстрое вычисление═ exp(x). Проблемы передачи информации, т.26, N 3, с.109, Хроника-17ая Всесоюзная школа по теории информации и е╦ приложениям (1990).
[23] Е.А.Карацуба, О новом методе быстрого вычисления трансцендентных функций. Успехи Математических Наук, т.46, N 2 (278), с.219-220 (1991).
[24] Е.А.Карацуба,═ О быстром вычислении трансцендентных функций}. Доклады Академии Наук СССР , т.319, N 2, с.278-279 (1991).
[25] Е.А.Карацуба, Быстрое вычисление трансцендентных функций .Проблемы передачи информации, т.27, N 4, с.87-110 (1991).
[26] Е.А.Карацуба, Быстрое вычисление ζ(3). Проблемы передачи информации, т.29, N 1, с.68-73 (1993).
[27] Catherine
A.Karatsuba, Fast evaluation of Bessel functions. Integral Transforms and
Special Functions, v.1, N4, pp. 269-276 (1993).
[28] Е.А.Карацуба,═
Быстрое вычисление дзета-функции Римана ζ(s) для целых значений аргумента═
s. Проблемы передачи информации, т.31, N 4, с.69-80 (1995).
[29] Е.А.Карацуба, О быстром вычислении дзета-функции Римана для целых значений аргумента .Доклады Академии Наук СССР, т.349, N 4, с.463 (1996).
[30] Е.А.Карацуба, Быстрое вычисление дзета-функции Гурвица и L-рядов Дирихле. Проблемы
передачи информации, т.34, N 4, с. 342-353 (1998).
[31] Ekatharine
A. Karatsuba, Fast evaluation of hypergeometric function by FEE.
Computational Methods and Function Theory (CMFT'97), N.Papamichael, St.Ruscheweyh
and E.B.Saff, eds., World Sc.Pub., pp. 303-314 (1999).
[32] E.A.Karatsuba,
On the computation of the Euler constant═
γ .
[33] E.A.Karatsuba,
Fast computation of═ ζ(3) and of some special integrals
,using the polylogarithms, the Ramanujan formula and it's generalization. J. of
Numerical Mathematics BIT, v. 41, N 4, pp.722-730 (2001).
[34] E.A.Karatsuba,
Fast computation of some special integrals of mathematical physics. Scientific
Computing, Validated Numerics, Interval Methods , W.Kr\"amer,J.W.von
Gudenberg,eds. ,pp.29-41, (2001).
[35] Д.E.Kнут, Искусство программирования на ЭВМ. т.2, Изд. Мир, 724 с., Москва (1977).
[36] А.Н.Колмогоров,═ О некоторых асимптотических характеристиках вполне ограниченных метрических пространств. Доклады Академии Наук СССР, т.108, N 3, с.385-388 (1956).
[37] А.Н.Колмогоров, Теория информации и теория алгоритмов. Изд. Наука, 303 с., Москва (1987).
[38] J.Landen,
An investigation of a general theorem for finding the length of any arc of any
conic hyperbola, by means of two elliptic arcs, with some other new and useful
theorems deduced therefrom. Philos.Trans.Royal Soc. London v.65 ,pp.283-289
(1775).
[39] D.W.Lozier
and F.W.J.Olver, Numerical Evaluation of Special Functions. Mathematics of
Computation 1943-1993: A Half -Century of Computational Mathematics,
W.Gautschi,eds., Proc. Sympos. Applied Mathematics, AMS, v.48, pp.79-125
(1994).
[40] В.И.Нечаев, Элементы криптографии. Основы теории защиты информации. Изд. Высшая школа, 110 с., Москва (1999).
[41] V.Ya.Pan,
Strassen's algorithm is not optimal. Proc. ACM Symp. on the Found. of
Comp.Science, pp.166-176 (1978).
[42] E.Salamin,═ Computation of π═ using arithmetic-geometric
mean. Math. Comp., vol.30, N 135, pp.565-570 (1976).
[43] A.Schönhage,
Schnelle Multiplikation grosser Zahlen. Computing, v.1, pp.182-196 (1966).
[44] A.Schönhage
und V.Strassen, Schnelle Multiplikation grosser Zahlen. Computing, v.7, pp.281-292
(1971).
[45] A.Schönhage,
A.F.W.Grotefeld and E.Vetter, Fast Algorithms ≈ A Multitape Turing Machine
Implementation. BI-Wiss.-Verl., 300 pp., Zürich (1994).
[46] C.L.Siegel,
Transcendental numbers.
[47] V.Strassen,
Gaussian elimination is not optimal. J. Numer. Math., N 13, pp.354-356 (1969).
[48] А.Л.Тоом,═ О сложности схемы
из функциональных элементов, реализующей умножение целых чисел. Доклады
Академии Наук СССР, т.150, N 2, с.496-498 (1963).
J
═════ Ответы на часто задаваемые (в основном
программистами) вопросы
═════ О методе БВЕ в статеье:
═════ Е.А.Карацуба ,═ ⌠Быстрое вычисление дзета-функции Римана ζ(s) для целых значений аргумента═ s.■
═════ Проблемы передачи
информации, т.31, N 4, с.69-80 (1995).