МАТЕМАТИКА

УДК 511

С. Л. Амбарян

Алгоритм деления
с неограниченной точностью

(Представлено академиком В.С.Закаряном 18/I 2001)

   Точность выполнения арифметических операций на вычислительных системах / персональных компьютерах имеет чрезвычайно важное значение. Можно привести множество примеров, когда данная разрядная сетка не обеспечивает требуемую точность решения поставленной задачи. Кроме того, некоторые исследователи науки и техники (теория чисел, криптография, комбинаторика, теория бесконечных рядов, дифференциальные уравнения, интегральное исчисление, атомная физика, разные направления механики, астрономия и др.) в своих исследованиях, требующих высокой точности, просто лишены возможности непосредственно использовать большие преимущества персональных компьютеров из-за ограничений разрядной сетки. Поэтому в связи с широким распространением и усовершенствованием персональных компьютеров разработка эффективных алгоритмов и программно-аппаратных средств реализации арифметики многократной точности (кратко а.м.т.) имеет хорошую перспективу [1-3].
   Вопросы арифметики однократной, двойной и многократной точности исследованы многими авторами, проблемы, история и библиография, посвященные им, великолепно изложены в [1]. А. Шенхаге и Ф. Штрассен разработали изящный алгоритм, который умножение (деление) двух n-разрядных чисел выполняет за 0 (n log n  log log n) шагов [2].
   Представление цифр-чисел в а.м.т. можно осуществить по-разному. Нами принят следующий подход.
    a. Машинное слово процессора вычислительной системы / персонального компьютера  состоит из  четного  числа  двоичных  разрядов,  равного 2(k + 1) (см. рис.1.), один из которых выделен для знака целого числа с фиксированной точкой (точка фиксирована после младшего разряда).

Рис. 1.

    b. Основанием системы счисления p выбрано целое число из отрезка [2, 2k].
    c. Каждая цифра а.м.т. представлена в одном полуслове как целое положительное число.
   d. Числа а.м.т. представлены как числа с плавающей точкой. Первое полуслово выделено для знака числа, второе полуслово - для порядка числа и, наконец, начиная с третьего полуслова в последовательных адресах оперативной памяти представлены цифры числа (дробная часть) а.м.т. (см. рис. 2.).

Рис. 2.

   Ограничение пункта b равносильно тому, что в одном машинном слове «свободно» помещается двузначное число а. м. т. При соблюдении данного ограничения максимальное число, которое может получиться в результате арифметических операций однократной точности в одном машинном слове, не менее чем 2p2 - 1, т.е. 2p2 - 1 Ј max u Ј 2 · 22k - 1.
   Дадим определение допустимых арифметических операций однократной точности.
   Определение. Арифметическую операцию с однократной точностью будем называть допустимой, если исходные, промежуточные и окончательные результаты не превышают  2 · 22k - 1.
   Пусть даны два n-разрядных положительных целых числа a и b по основанию p.
a = a1a2 ј an, a1 № 0  и b = b1b2 ј bn,  b1 № 0;
Ak = p · Ak-1 + ak, A0 = 0  и Bk = p · Bk-1 + bk, B0 = 0; k = 1,2, ј,n.
(1)

   Рассмотрим операцию деления над числами неограниченной (произвольно высокой) точности в а. м. т.
   Определим число c как

c = a / b = = q0, c1,c2јcn;  q0 = [a / b] = [An / Bn];
где 0 Ј ak, bk, ck < p;  k = 1,2,ј,n.

   При Ak > Bk рассмотрим представление Ak через Bk в форме
Ak = qk · Bk + rk,  0 Ј rk < Bk,  1 Ј qk < p; k = 1,2,ј,n.
 (2)

   При Ak < Bk рассмотрим представление Ak+1 через Bk в форме
Ak+1 = qk · Bk + rk,  0 Ј rk < Bk,  1 Ј qk < p; k і 2.
(2a)

   В представлении (2) rk и rk-1 связаны между собой рекуррентным соотношением
rk = prk-1 + ak - qkbk + p(qk-1 - qk)Bk-1,
где qk Ј qk-1,  k = 1,2,ј,n.
(3)

   При qk № qk-1 (т.е. qk-1 = qk + 1, может быть кроме k = 2) будем иметь
rk = prk-1 + ak - qkbk + pBk-1; k = 2,3,ј,n.
 (3a)

   При qk = qk-1 = q (см. теорему 4) получим
rk = prk-1 + ak - qbk; k = 3,4,ј,n.
(3b)

   Tребуется при помощи допустимых арифметических операций эффективно (быстро) определить значение q0.
   В этих обозначениях справедлива следующая лемма.
   Лемма 1. q - 1 Ј q0 Ј q.
   Доказательство. c = a / b = a1a2,a3јan / b1b2,b3јbn =
= (A2 + x) / (B2 + y) = q + (r + x - qy) / (B2 + y), где q = q2, r = r2 и 0 < x, y < 1.
   Очевидно, что |(r + x - qy) / (B2 + y)| < 1, действительно, это следует из следующих цепочек неравенств:
0 Ј r + x - qy Ј r + x Ј B2 - 1 + x Ј B2 + y,
0 і r + x - qy і - qy і -q > - B2 і -B2 - y.
Таким образом,
q0 = м
н
о
q,
если r + x - qy і 0,
q - 1,
если r + x - qy < 0.

   Следствие 1. Если q = 0, то q0 = 0;
                           если q Ј r, то q0 = q.
   Согласно следствию, если r і q, то q0 = q, и поэтому предположим, что 0 Ј r < q. Обозначим через m то минимальное значение k, для которого в (3b) имеет место одно из следующих неравенств:
rm < 0  либо   rm і q,
тогда справедлива следующая теорема.
   Теорема 1(основная теорема). Если rm < 0, то q0 = q - 1, или
                                                          если rm і q, то q0 = q.
   Доказательство. Пусть 0 Ј ri < q, i = 2,3,ј, m - 1, и для m имеет место условие теоремы, т.е.
rm < 0  либо   rm і q.

   В результате поочередной подстановки значений ri, i = 3,4,ј,m - 1, в рекуррентное соотношение (3b) получим
0 Ј Am-2 - q · Bm-2 + am-1 - qbm-1 < q,  т.е. 0 Ј Am-1 - q · Bm-1 < q.

   Рассмотрим случай, когда rm < 0; тогда
p(Am-1 - q · Bm-1) + am - qbm < 0,
Am - q · Bm Ј -1, (Am + 1) / Bm Ј q.
Из последнего неравенства и леммы 1 следует, что q0 = q - 1.
   Пусть теперь rm і q, тогда
p(Am-1 - q · Bm-1) + am - qbm і q,
Am - q(Bm + 1) і 0, Am / (Bm + 1) і q.

из последнего неравенства и леммы 1 следует, что q0 = q. Теорема доказана.
   Если не существует такое m, т.е. 0 Ј ri < q для всех i = 2,3,ј, n, то справедливо следующее следствие.
   Следствие 2. a / b = q + rn / b,  0 Ј rn < q.
   Теперь подробно рассмотрим случай A2 < B2. Так как имеет место неравенство p2 Ј A3 Ј p3 - p - 1, то необходимо разработать алгоритм, который с помощью допустимых арифметических операций позволяет эффективно определить значения неполного частного q = q2 и остатка r = r2.
   Элементарные преобразования приводят к формуле
q2 = q1 - (q1b2 + r2 - (pr1 + a3)) / B2,
q2 = q1 -]q1b2 - (pr1 + a3)[/B2,
(4)
DA = q1b2 - (pr1 + a3)  через B2  в форме
DA = Dq + Dr / B2,  0 Ј Dr < B2,  -1 Ј Dq Ј p-1.
(5)
Следует заметить, что q1b2 Ј (2p - 2)(p - 1) < 2p2 - 1.
   Тогда справедлива следующая теорема.
   Теорема 2 (первая теорема о частном и об остатке). При A2 < B2 значения неполного частного и остатка в (2a) определяются формулами:
q = q1 - Dq  и r = 0,  при Dr = 0;  или
 (6)
q = q1 - (Dq + 1)  и r = B2 - Dr,  при Dr № 0.
 (7)

   Правильность формул (6) и (7) следует из (4) и (5).
   Аналог следствия 1 здесь звучит так: если q2 = 1, то q0 = 1 и если q2 = p, то q0 = p - 1.
   Доказанная теорема занимает ключевое место при разработке алгоритмов операции деления в а.м.т.
   В а.м.т. в алгоритме операции деления [1] после определения неполного частного q0 с точностью единицы (q - 1 Ј q0 Ј q) вычисляется разность
at = at - qb,

где at(at) - текущее (предыдущее) значение делимого, и при at < 0 производится компенсирующее сложение в а.м.т.

at = at + b, 0 < at < b.

   Ниже следует лемма, согласно которой компенсирующее сложение в алгоритме операций деления может отсутствовать.
   Лемма 2 (о компенсирующем сложении). Если неполное частное определено с точностью единицы, то i-ый разряд частного c выражается формулой:
ci = q,  если at > 0  и at і 0,
ci = q - 1,  если   at > 0  и   at < 0,
ci = p + q,  если at < 0  и at і 0,
ci = p + q - 1,  если   at < 0  и   at < 0; i = 1, 2, ј, n.

   Доказательство.  Первые два случая очевидны, докажем третий случай, откуда следует и четвертый.
   Пусть at = at - qb < 0, тогда определение следующей цифры частного достигается при помощи представления
p at / b = q + r / b,  0 Ј r < b, -p < q Ј -1,
так как p(at + b) / b = p + pat / b = p + q + r / b, откуда и следует утверждение леммы.
   Замечание. Отметим, что здесь можно обойтись без использования отрицательных неполных частных (q < 0 при a < 0), положив a = -a после операции a = a - qb.
   Пусть, например, w первоначально принимает значение, равное нулю (w = 0). Далее w принимает значение, равное единице (w = 1), если результат операции a = a - qb окажется отрицательным (a < 0), то w сохраняет это значение до тех пор, пока a заново не станет отрицательным, тогда w принимает значение, равное нулю (w = 0) и т.д. Тогда
ci = q,  если a і 0  и w = 0,
ci = q - 1,  если   a < 0  и   w = 0,
ci = p - q - 1,  если a і 0  и w = 1,
ci = p - q,  если   a < 0  и  w = 1; i = 1, 2, ј, n.

   Согласно основной теореме и первой теореме о неполном частном и об остатке опишем алгоритм операции деления в а.м.т.
   Описание алгоритма операции деления. Работа алгоритма операции деления начинается с проверки на нуль делимого a (a1 = 0 или a1 № 0) и делителя b (b1 = 0 или b1 № 0). Первоначально положим показатель результата e(c) = e(a) - e(b). В дальнейшем показатель делителя b e(b) = 0. При A2 і B2 установить показатель делимого e(a) = 0 и уточнить значение e(c) = e(c) + 1, а при A2 < B2 установить показатель делимого e(a) = 1.
    Алгоритм операции деления неотрицательных целых чисел а.м.т. состоит из следующих шагов.
          Шаг 1. Пр. «Опр. e(c), e(a) = 0/1, e(b) = 0»;
                   2. i = 1;
                   3. (q,r) = A2+e(a) / B2;
                   4. Если r і q, то идти к 17;
                   5. Пр. «Уточнение значения q»;
                   6. Если r і 0, то идти к 17;
                   7. q = q - 1;
                   8. m = k - (3 + e(a));
                   9. Пр. «ci = q, ci+j = p - 1; j = 1,2,ј,m»;
                   10. Если m = 0, то идти к 15;
                   11. q = q + 1;
                   12. a = a - qb;
                   13. e(a) = e(a) + m;
                   14. a = a + b; идти к 16;
                   15. a = a - qb;
                   16. e(a) = e(a) + 1; идти к 3;
                   17. ci = q;
                   18. a = a - qb;
                   19. Пр. «Опр. e(a) и ci+j = 0, j = 1,2,јe(a)».
                   20. i = i + 1, если i Ј n, то идти к 3.
                   21. End.
   Помимо того, что описанный алгоритм исключает компенсирующее сложение (при a - qb < 0) [1], он имеет и другие преимущества. Перейдем к строгому изложению некоторых из них.
   Теорема 3. Если для некоторого i = 1,2,ј,n и некоторого k = m + 3 + e(a) і 3 в (3b) имеет место одно из следующих неравенств: rm < 0 либо rm і q, тогда:
           Случай 1. Если rk < 0, то ci = q - 1 и за этой цифрой следует не менее m цифр ci+j = p - 1; j = 1,2,ј,m.
           Случай 2. Если rk і q, то ci = q и за этой цифрой следует не менее m цифр ci+j = 0; j = 1,2,ј,m.
   Доказательство.
           Случай 1. Пусть rk < 0, тогда 0 Ј rk-1 < qk-1 = q и qk = q - 1;
c = a / b > Ak / (Bk + 1) = (qkBk + rk) / (Bk + 1) = qk + (rk - qk) / (Bk + 1) =
= qk + (prk-1 + ak - qkbk + pBk-1 - qk) / (Bk + 1) і
і qk + (pBk-1 - qk(bk + 1)) / (Bk + 1) = qk + 1 - (qk + 1)(bk + 1) / (Bk + 1) і
і qk + 1 - (qk + 1) / (Bk-1 + 1) = q - q / (Bk-1 + 1) і q - q / (pm+1+e(a) + 1) і
і q - q/(pm+1 + 1);   c > (q - 1) + 1 - q / (pm+1+1),
откуда следует утверждение теоремы по первой части.
   При m = 0 из последнего неравенства следует, что
c > (q - 1) + 1 - (p - 1) / (p + 1) = (q - 1) + 2 / (p + 1),
т.е. ci = q - 1 и за ним следует цифра ci+1 і 1.
           Случай 2. Пусть теперь rk і 0, тогда 0 Ј rk-1 < qk-1 = q и qk = q;
c = a / b < (Ak + 1) / Bk = (qkBk + rk + 1) / Bk = qk + (rk + 1) / Bk =
= qk + (prk-1 + ak - qkbk + 1) / Bk Ј qk + (prk-1 - qkbk + 1) / Bk Ј
Ј qk + (prk-1 + ak + 1) / Bk Ј qk + pqk-1 / Bk Ј qk + qk-1 / Bk =
= q + q / Bk Ј q + q / pm+1+e(a) Ј q + q / pm+1 Ј q + (p - 1)/pm+1;
c < q + (p - 1) / pm+1;
откуда следует утверждение теоремы по второй части.
   При m = 0 из последнего неравенства следует, что
c < q + (p - 1) / p,
т.е. ci = q и за ним следует цифра ci+1 Ј p - 2.
   Теорема полностью доказана.
   Теорема 4 (вторая теорема о частном и об остатке). Если для некоторого k і 3 (в (3b)) имеет место неравенство 0 Ј rk < qk, то 0 Ј rk-1 < qk-1 и qk-1 = qk.
   Доказательство. Прежде всего докажем, что qk-1 = qk = q. Если qk-1 № qk, то из (3a) следует, что
ak = qkbk + rk - p(rk-1 + Bk-1) Ј qkbk + qk - p(rk-1 + Bk-1) Ј
Ј qk(bk + 1) - p(rk-1 + Bk-1) Ј p(qk - rk-1 - Bk-1) Ј 0, т.е. ak < 0,
что противоречит определению ak.
   Предположим, что утверждение теоремы неверно и rk-1 < 0, тогда
ak = qkbk + rk - prk-1,  так как
qbk + rk і 0,  то ak і -prk-1 і p,
что противоречит определению ak.
   Пусть теперь rk-1 і q, тогда
ak = qkbk + rk - prk-1 = -(p(rk-1 - q) + (p - bk)q - rk),  так как
p(rk-1 - q) і 0  и (p - bk)q - rk і 0,  то ak < 0,
что противоречит определению ak. Теорема доказана.

   Ереванский научно-исследовательский
   институт математических машин

Литература

     1. Кнут Д. Искусство программирования для ЭВМ. Получисленные алгоритмы. М.: Мир, 1977.
     2. Schonhage A., Strassen V. Computing, Archiv für elektronisches Rechnen, 7 (1971), Fasc. 3-4.
     3. Амбарян С. Л. Эффективные алгоритмы выполнения арифметических операций в арифметике с многократной точностью. Пакет программ а.м.т. (ПП-АМТ) в среде Visual C++. Ереванский научно-исследовательский институт математических машин (ЕрНИИММ), Ереван, 1998.