(in English)

 

БЫСТРЫЕАЛГОРИТМЫИМЕТОДБВЕ

Е.А. Карацуба

 

 

I.═══ БЫСТРЫЕ АЛГОРИТМЫ.

 

 

══════ Область вычислительной математики, которая называется быстрые алгоритмы, появилась в 1960 году.

══════ Под алгоритмом мы будем понимать правило или способ вычисления, не формализуя это понятие. Далее будем считать, что числа записаны в двоичной системе счисления, знаки которой0и1называются битами.

 

Опр.1.Запись знаков══ 0 , 1 , плюс, минус, скобка; сложение, вычитание и умножение двух битов назов╦м одной элементарной или битовой операцией.

 

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

══════ Первые постановки задач о битовой сложности вычислений (1956 г.) принадлежат А.Н. Колмогорову http://mech.math.msu.su/probab/Kolmogorov/kolmogorov.html══ [36], [37],который в то же время подч╦ркивал, что "цикл моих работ по теории информации был создан под большим влиянием публикаций Норберта Винера http://erudite.nm.ru/WienerNorbert.htmи Клода Шеннона http://www.astronet.ru/db/msg/1187395 ."

══════ Прежде чем ввести понятие сложности вычисления, определим, что значит вычислить функцию в заданной точке. Рассмотрим наиболее простой пример вычисления вещественной функцииy = f(x) вещественного переменногоx ,══ axb .

══════ Пусть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" начала развиваться теория быстрых вычислений.

═════ Исследования с целью поиска алгоритма умножения со сложностью, близкой к оптимальной, были

═════ продолжены рядом авторов (среди них были Тоом http://www.de.ufpe.br/~toom/

═════ [48], Кук http://www.cs.toronto.edu/~sacook/

═════ [15], Ш╦нхаге http://www.informatik.uni-bonn.de/~schoe/index.html

═════ [43]), и в 1971 г. Ш╦нхаге

═════ и Штрассен http://www.mathe.uni-konstanz.de/mitarbeiter/strassen.html

═════ [44], [45] построили алгоритм с наилучшей на настоящее время оценкой дляM(n),

 

M(n) = O(n log nlog log n).

 

═════ При построении этого алгоритма кроме "Divide and Conquer" они использовали идею выполнения

═════ арифметических действий по модулю чисел Ферма22 +1, а также быстрое преобразование Фурье.

════════════ На основе "Divide and Conquer" были построены алгоритмы быстрого умножения матриц. Первый

═════ такой алгоритм в 1970 г. наш╦л Штрассен [47], в дальнейшем эти исследования были продолжены

═════ Коперсмитом 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 ,══ 0r<b,

 

═════ сводится к сложению, вычитанию и умножению, и еслиaиbявляютсяn -значными числами,

═════ то сложность деленияO(M(n)).

════════════ Самый простой алгоритм деления заключается в вычислении методомНьютона обратной величины

═════ 1/bс точностью доnзнаков с последующим умножением наaпо алгоритму быстрого умножения.

════════════ Пусть надо вычислитьx = 1/b,где1/2b≤ 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))операций.

 

════════════ В 1976 г. Саламинhttp://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]

═════ предложил первые АГС -алгоритмы для быстрого вычисления простейших трансцендентных функций. В

═════ дальнейшем исследование и использование АГС √алгоритмов было продолжено многими авторами,

═════ среди которых были братья Джонатан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 выполнялось неравенствоRm2⁻ⁿ⁻¹. Это будет,

═════ например, когда══ m4n/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.

════════════ Зигель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¹ <n2k══ и положимN = 2k+¹ (заметим, что

═════ 2n ≤ N < 4n). ПустьxN =α32³ + α424══ + +αN2N,гдеαj = 0или1,══ j = 3, 4, ┘ , N,и

 

| x0 - xN |══ 2N.

 

═════ Тогда

ex0=══ exN + 4θ12N══ ;═══════ | θ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) <2N .

════ Следовательно, (4) можно записать в виде

 

═════ где

 

и постоянныеθν(r),══ 0 <θν(r) < 1 . Легко видеть, что сумма (5) является подходящей для е╦ суммирования с помощью БВЕ. На последнем шаге процесса в (5) имеем

 

 

а затем производим деление целого числаaν ══на целоеbνс точностью доN знаков. Тем самым, получаем такоеην, что

 

ξν = aν /bν = ην═══ +θν 2N,═══ | θν | < 1 .

 

СледовательноexNможет быть представлено в виде

 

| θν |≤ 1 .

Определим целое числоl неравенствами2l¹<k══2l. Предполагая, чтоην=1, θν = 0 для ν > k,мы можем записать последнее выражение в виде

 

| θν |≤ 1 .

Ясно, что══ log kl < log k + 1. Последнее произведение вычисляется за══ lшагов процесса, который аналогичен вычислению══ n! . Перемножим множители(6)последовательно попарно, вычисляя на 1-ом шаге══ 2l¹произведений вида

 

ν+ 2 θν 2N) (ην +1+ 2 θν+1 2N)══ =ην ην +1+ 2 θν 2Nην +1 + 2 θν+1 2N ην + 4θνθν+1 22N══ = ην ,ν +1+ 8 θν,ν+1 2N,

ν0(2),══ ν = 2, 3,┘ , 2l + 1,

 

где ην ,ν +1══ естьN -значное число, такое что

 

ην ,ν +1= ην ην +1+ θ̃ 2N ,══ | θ̃ |≤ 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,2nN < 4n,то из последнего соотношения находим

 

ex0=η +θ 2n ,

═════ гдеη=η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θ12N ,

 

 

═════ | θ1|≤ 1 .

═════ Запишем

═══════════════════════════

═════ в виде

 

═════ где

 

 

0        <θν(r)<1.Далее мы производим с суммами (8), (9) те же действия, что и с суммой (4). Пусть ην,

fν естьN-значные числа, такие что

 

ξν =aν /bν = ην═══ +θν 2N══ ,══ | θν |≤ 1,

ψν =aν /dν ══= fν ═══+θν 2N,══ | θν |≤ 1,

 

,

═════ Тогда (7) может быть представлено в виде

 

═════ | θν |≤ 1 .

════════════ Вычисляем произведение (10) посредством процесса, описанного в предыдущем параграфе.

═════ В результате со специально выбраннымиl , N , r═══ находим

 

eix0=f══ +i η ══+2θ 2N══ ,══ | θ |≤ 1 .

═════ и отсюда

cos x0=══ f+θ1 2n══ ,══ | θ1|≤ 1 ,

sin x0=══ η+θ2 2n══ ,══ | θ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.

═════ Предполагая, чтоr3n , разложим подинтегральную функцию в (11)y = et ,в ряд Тейлора

═════ по степеням t:

 

 

═════ где══ Rr =θr tr+1/(r+1)!,| θr |══ 1.Отсюда имеем приr+1 = 6n

 

Γ(x0)=nx0 S +θ 2n ,| θ |≤ 1 , (12)

═════ где

 

       и мы будем вычислять═══ Γ(x0)=Γ(a/b)по формулам (12)--(13). Заметим, что для вычисления══ nx0 = na/b,

═════ сточностью до══ n знаков достаточноO(M(n))операций (методом Ньютона). Чтобы вычислитьS

═════ возьмемr+1 =2k,2k1<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)r1n 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! и производим одно

═════ деление с точностью до22nпо формуле

 

S = S1 (k) =prk(k)/qrk(k)/r! ,

 

═════ что да╦т суммуS с точностью2n1 . Следовательно значение══ Γ(x0) =Γ(a/b) вычислено с точностью

═════ 2n+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-1xt1 + ┘ + g1x + g0 ;════════════ g(α) = 0 ;═══ (14)

 

═════ gt ,gt-1, ┘,g0 --- целые числа,t≥ 2.══ Как и в случае рационального аргумента воспользуемся формулами

═════ (12), (13). Отметим, что вычисление══ nα ════с точностью22n══ по формулеnα = eα log n требует

═════ O(M(n) log2n) операций. СуммаSиз (12) принимает вид

 

 

═════ т.е.

S = S1(0) + S2(0) + ┘ + Sr+1(0) ,(16)

 

═════ где

 

════════════ ВычислениеS══ выполняется заkшагов, r+1 =2k,2k1<6n══ 2k,══ k ≥ 1.Особенностью этого

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

═════ сложности вычисления. До шагаi,гдеiопределяется условиями1≤ ik ,2i1 <t2i,производим

════ с суммой (16), (17) действия, аналогичные вышеописанным для случая рационального аргумента. Однако,

═════ группируя слагаемые суммы (16) таким же образом как и в предыдущем параграфе, теперь мы вычисляем

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

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

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

 

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

 

════════════════════════════ ════════════════════════════

════════════════════════════════════════════════════════

 

════════════════════════════════════════════════════════

════════════════════════════════════════════════════════

═════ ν = 0 , 1 , ┘ , (r+1)/2-1.

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

 

═════ где

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

 

═════════════════════════════════════

═════ 0 ≤ m 12i1-1 ,0 ≤ m 22i1.

════════════════

═════ 0 ≤ l 12i1 ,0 ≤ l 22i1.

════

 

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

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

════════════ Передi+1 -ым шагом(1≤ ik ,2i1 <t2i)мы редуцируем многочлены (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, qi-1 , и вычисляемδrk(k),ρrk(k) иS,и следовательноΓ(α) с точностью2n+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., New York (1987).

 

[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, Harvard University (1966).

 

[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γ . University of Helsinki preprint, 21 pp. (1999); J. of Numerical Algorithms, v. 24, pp.83-97, (2000).

 

[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. Princeton University Press, 102 pp., Princeton (1949).

 

[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).

 

 

 

 


© Екатерина Карацуба
Последняя редакция: 28.06.2005.