Алгоритм деления

Эта статья — об алгоритмах деления целых чисел. О делении столбиком см. Деление столбиком; о теореме существования единственного частного и остатка см. Деление с остатком; о делении многочленов см. Деление многочленов.

Алгоритм деления — алгоритм, вычисляющий для двух данных целых чисел ${\displaystyle N}$ и ${\displaystyle D}$ их частное и/или остаток, результат деления с остатком. Некоторые из алгоритмов предназначены для вычислений вручную, другие реализованы в цифровых схемах и программном обеспечении.

Алгоритмы деления разбиваются на две большие категории: медленное и быстрое. Алгоритмы медленного деления дают по одному знаку результата за итерацию. Примерами медленного деления служат алгоритмы деления с восстановлением, без восстановления и SRT. Методы быстрого деления начинаются с аппроксимации конечного частного и дают вдвое больше знаков в конечном результате на каждой итерации. Алгоритмы Ньютона – Рапсона и Гольдшмидта попадают в эту категорию.

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

Обсуждение будет использовать обозначения ${\displaystyle N/D=(Q,R)}$, где

• ${\displaystyle N}$ — числитель (делимое)
• ${\displaystyle D}$ — знаменатель (делитель)

являются входными числами, а

• ${\displaystyle Q}$ — частное
• ${\displaystyle R}$ — остаток

являются выходными данными.

Деление путём вычитаний

Простейший алгоритм, исторически встроенный в алгоритм поиска наибольшего общего делителя и представленный в Началах Евклида, Книга VII Предложение 1, находит остаток деления двух положительных целых чисел с помощью только вычитания и сравнения:

R := N
while R >= D do
  R := R − D
end
return R

Доказательство, что частное и остаток существуют и единственные (описано в статье Деление с остатком) дают полный алгоритм деления, основанный на сложениях, вычитаниях и сравнениях:

function divide(N, D)
  if D = 0 then error(DivisionByZero) end
  if D < 0 then (Q, R) := divide(N, −D); return (−Q, R) end
  if N < 0 then
    (Q,R) := divide(−N, D)
    if R = 0 then return (−Q, 0)
    else return (−Q − 1, D − R) end
  end
  -- Здесь N >= 0 и D >= 0
  return divide_unsigned(N, D)
end

function divide_unsigned(N, D)
  Q := 0; R := N
  while R >= D do
    Q := Q + 1
    R := R − D
  end
  return (Q, R)
end

Эта процедура всегда даёт ${\displaystyle R\geqslant 0}$. Будучи очень простым, алгоритм требует ${\displaystyle \Omega (Q)}$ шагов, а потому экспоненциально медленнее, чем алгоритмы наподобие длинного деления. Алгоритм полезен, если известно, что ${\displaystyle Q}$(шагов) мало (будучи зависимым от объёма вывода).

Деление столбиком

Основная статья: Деление столбиком

Деление столбиком является стандартным алгоритмом многозначных чисел в десятичной системе счисления, используемых для деления с помощью карандаша и бумаги. Он сдвигается постепенно слева направо от делимого, вычитая наибольшее возможное кратное делителя на каждом шаге. Множитель становится цифрой в частном, а конечная разность становится остатком от деления.

Когда алгоритм используется по основанию 2, этот метод образует базис для деления натуральных чисел с остатком. Короткое деление^[англ.] является вариантом деления в столбик, пригодным для деления на одну цифру. Алгоритм уменьшения^[англ.], известный также как метод неполных частных, является менее эффективным видом деления столбиком, но который проще понять. Разрешение вычитать большее кратное число, чем делается на каждом шаге, даёт больше свободы для создания вариантов деления в столбик.

Деление целых чисел (без знака) с остатком

Основная статья: Деление столбиком

Приведённый алгоритм, двоичная версия известного деления в столбик, делит ${\displaystyle N}$ на ${\displaystyle D}$, помещая частное в ${\displaystyle Q}$ и остаток в ${\displaystyle R}$. Все значения трактуются как целые числа без знака.

if D = 0 then error(DivisionByZeroException) end -- Деление на ноль
Q := 0                  -- Начальные значения частного и остатка полагаем равны 0
R := 0                     
for i := n − 1 .. 0 do  -- Здесь n равно числу бит в N
  R := R << 1           -- Сдвиг влево числа R на 1 бит
  R(0) := N(i)          -- Полагаем младший бит R равным биту i делимого
  if R >= D then
    R := R − D
    Q(i) := 1
  end
end

Пример

Возьмём ${\displaystyle N=1100_{2}}$ (${\displaystyle 12_{10}}$) и ${\displaystyle D=100_{2}}$ (${\displaystyle 4_{10}}$)

Шаг 1: Set R = 0 and Q = 0
Шаг 2: Take i = 3 (на единицу меньше числа бит в N)
Шаг 3: R = 00 (сдвиг влево на 1)
Шаг 4: R = 01 (полагаем R(0) равным N(i))
Шаг 5: R < D, так что пропускаем команду

Шаг 2: Set i = 2
Шаг 3: R = 010
Шаг 4: R = 011
Шаг 5: R < D, пропускаем команду

Шаг 2: Set i = 1
Шаг 3: R = 0110
Шаг 4: R = 0110
Шаг 5: R >= D, команда выполняется
Шаг 5b: R = 10 (R−D)
Шаг 5c: Q = 10 (полагаем Q(i) равным 1)

Шаг 2: Set i = 0
Шаг 3: R = 100
Шаг 4: R = 100
Шаг 5: R >= D, команда выполняется
Шаг 5b: R = 0 (R−D)
Шаг 5c: Q = 11 (полагаем Q(i) равным 1)

Конец
${\displaystyle Q=11_{2}}$ (${\displaystyle 3_{10}}$) и ${\displaystyle R=0}$.

Методы медленного деления

Все методы медленного деления основываются на стандартном рекуррентном отношении^[1]

${\displaystyle R_{j+1}=B\times R_{j}-q_{n-(j+1)}\times D,}$

где:

• ${\displaystyle R_{j}}$ является j-м частичным остатком деления
• B − основание счисления (равно обычно 2 в компьютерах и калькуляторах)
• ${\displaystyle q_{n-(j+1)}}$ является цифрой частного в позиции ${\displaystyle n-(j+1)}$, где позиции цифр пронумерованы от менее значащих цифр (0) к более значащим (${\displaystyle n-1}$)
• ${\displaystyle n}$ − число знаков в частном
• ${\displaystyle D}$ − делитель

Восстанавливающее деление

Восстанавливающее деление работает с дробными числами числами с плавающей запятой и зависит от предположения ${\displaystyle 0<D<N}$.

Цифры частного ${\displaystyle q}$ формируются из набора ${\displaystyle \{0,1\}}$.

Основной восстанавливающий алгоритм для двоичного (по основанию 2):

R := N
D := D << n             -- R и D должны быть словами вдвое длиннее, чем N и Q
for i := n − 1 .. 0 do  -- Например, 31..0 для 32 бит
  R := 2 * R − D        -- Пробное вычитание из сдвинутого значения (умножение на 2 есть сдвиг в бинарной интерпретации)
  if R >= 0 then
    q(i) := 1          -- Результат - бит 1
  else
    q(i) := 0          -- Результат - бит 0
    R := R + D         -- Новый частичный остаток равен (восстановленному) сдвинутому значению
  end
end

-- Здесь: N = делимое, D = делитель, n = число бит, R = частичный остаток, q(i) = бит № i частного

Вариант алгоритма, не осуществляющий явно ${\displaystyle 2R}$, сохраняет его, а потому ${\displaystyle D}$ не нужно добавлять обратно в случае ${\displaystyle R<0}$.

Невосстанавливающее деление

Невосстанавливающее деление использует набор цифр ${\displaystyle \{-1,1\}}$ для цифр частного вместо ${\displaystyle \{0,1\}}$. Алгоритм более сложен, но имеет преимущество, когда реализуется в микросхемах, что имеется всего одно принятие решения и одно сложение/вычитание на бит частного. Нет шага восстановления после вычитания, что может сократить число операций вдвое, и позволяет выполнить алгоритм быстрее^[2]. Алгоритм невосстанавливающего деления для двоичных (по основанию 2) положительных чисел:

R := N
D := D << n            -- R и D должны быть словами вдвое длиннее, чем N и Q
for i = n − 1 .. 0 do  -- Например, 31..0 для 32 бит
  if R >= 0 then
    q[i] := +1
    R := 2 * R − D
  else
    q[i] := −1
    R := 2 * R + D
  end if
end
 
-- Здесь: N = делимое, D = делитель, n = число бит, R = частичный остаток, q(i) = бит №i частного.

Если следовать этому алгоритму, получаем частное в нестандартном формате, состоящем из цифр -1 и +1. Этот вид нужно преобразовывать в двоичную форму. Пример:

Преобразование частного к цифрам ${\displaystyle \{0,1\}}$:
Старт:	${\displaystyle Q=111{\bar {1}}1{\bar {1}}1{\bar {1}}}$
1. Образуем положительный член:	${\displaystyle P=11101010\,}$
2. Образуем отрицательный член:	${\displaystyle M=00010101\,}$
3. Вычитаем: ${\displaystyle P-M}$:	${\displaystyle Q=11010101\,}$
Отрицательный член представлен в бинарном обратном коде, не в дополнительном коде

Если −1 знаки ${\displaystyle Q}$ запомнены как нули (0), как в обычном представлении, то ${\displaystyle P}$ является ${\displaystyle Q}$ и вычисление ${\displaystyle M}$ тривиально: осуществляет побитное дополнение (бит заменяется на дополнение бита) исходного ${\displaystyle Q}$.

Q := Q − bit.bnot(Q)      -- Если цифры −1 в Q представлены как нули.

Наконец, частные, вычисляемые этим алгоритмом, всегда нечётные, а остаток R находится в пределах ${\displaystyle -D\leqslant R\leqslant D}$. Например, ${\displaystyle 5/2=3R-1}$. Для приведения к положительному остатку делаем один шаг восстановления после того, как ${\displaystyle Q}$ приведена из нестандартного вида к стандартному:

if R < 0 then
  Q := Q − 1
  R := R + D  
end if

Истинный остаток равен R >> n. Как и в варианте с восстановлением, младшие биты ${\displaystyle R}$ используются в том же порядке, как образуются биты частного ${\displaystyle Q}$, и имеет смысл использовать единый сдвиг регистра для обоих чисел одновременно.

Деление SRT

Деление названо по первым буквам фамилий создателей (Sweeney, Robertson, Tocher). Деление SRT является популярным методом деления во многих микропроцессорах^[3][4]. Деление подобно делению без восстановления, но использует таблицу поиска на основе делимого и делителя для определения цифры частного.

Наиболее существенная разница в том, что используется избыточное представление для частного. Например, если реализуется деление SRT по основанию 4, каждая цифра частного выбирается из пяти возможностей: ${\displaystyle \{-2,-1,0,+1,+2\}}$. Ввиду этого выбор цифры частного не обязательно должен быть точным. Позднее цифры частного можно откорректировать. Например, пары цифр ${\displaystyle {0,+2}}$ и ${\displaystyle {1,-2}}$ эквивалентны, поскольку ${\displaystyle 0{\times }4+2=1{\times }4-2}$. Это допустимое отклонение позволяет выбирать цифры частного, исходя только из нескольких наиболее значащих бит делимого и делителя, вместо вычитания по полной длине. Это упрощение, в свою очередь, позволяет применять основание для чисел, большее 2.

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

Печально известный баг плавающего деления в процессорах Intel Pentium был вызван неправильно закодированной таблицей поиска. Пять из 1066 ячеек таблицы были ошибочно опущены^[5][6].

Методы быстрого деления

Деление Ньютона – Рапсона

Деление Ньютона – Рапсона использует метод Ньютона для нахождения обратного к ${\displaystyle D}$ числа и умножает это обратное на ${\displaystyle N}$, чтобы найти результирующее частное ${\displaystyle Q}$.

Шаги деления Ньютона – Рапсона:

1. Вычисляем приближение ${\displaystyle X_{0}}$ для обратного числа ${\displaystyle 1/D}$ делителя ${\displaystyle D}$.
2. Последовательно вычисляем более точные приближения ${\displaystyle X_{1},X_{2},\ldots ,X_{S}}$ обратного числа. Это место, где, собственно, и используется метод Ньютона – Рапсона.
3. Вычисляем частное путём умножения делимого на обратное к делителю: ${\displaystyle Q=NX_{S}}$.

Чтобы применить метод Ньютона для нахождения обратного к ${\displaystyle D}$ числа, необходимо найти функцию ${\displaystyle f(X)}$, имеющую нуль в точке ${\displaystyle X=1/D}$. Очевидно, что такой функцией является ${\displaystyle f(X)=DX-1}$, но итерации Ньютона – Рапсона для неё неуспешны, поскольку не могут быть осуществлены без знания обратного к ${\displaystyle D}$ числа (более того, метод пытается вычислить точное обратное за один шаг, а не делать итеративные улучшения). Функция, которая работает – это ${\displaystyle f(X)=(1/X)-D}$, для которой итерация Ньютона – Рапсона даёт

${\displaystyle X_{i+1}=X_{i}-{f(X_{i}) \over f'(X_{i})}=X_{i}-{1/X_{i}-D \over -1/X_{i}^{2}}=X_{i}+X_{i}(1-DX_{i})=X_{i}(2-DX_{i}),}$

что можно вычислить из ${\displaystyle X_{i}}$, используя только умножение и вычитание или два умножения-сложения.

С точки зрения вычислений выражения ${\displaystyle X_{i+1}=X_{i}+X_{i}(1-DX_{i})}$ и ${\displaystyle X_{i+1}=X_{i}(2-DX_{i})}$ не эквивалентны. Чтобы получить точность в 2n бит при использовании второго выражения, нужно вычислить произведение между ${\displaystyle X_{i}}$ и ${\displaystyle (2-DX_{i})}$ с двойной точностью к заданной точности ${\displaystyle X_{i}}$(n бит). В отличие от этого произведение ${\displaystyle X_{i}}$ и ${\displaystyle (1-DX_{i})}$ нужно вычислять только с точностью n бит, поскольку ведущие n бит (после двоичной точки) числа ${\displaystyle (1-DX_{i})}$ нулевые.

Если ошибка определяется как ${\displaystyle \varepsilon _{i}=1-DX_{i}}$, то

${\displaystyle {\begin{aligned}\varepsilon _{i+1}&=1-DX_{i+1}\\&=1-D(X_{i}(2-DX_{i}))\\&=1-2DX_{i}+D^{2}X_{i}^{2}\\&=(1-DX_{i})^{2}\\&={\varepsilon _{i}}^{2}.\\\end{aligned}}}$

Эта ошибка в квадрате на каждом шаге итерации (так называемая квадратичная сходимость метода Ньютона – Рапсона) имеет эффект такой, что число верных знаков в результате, грубо говоря, удваивается для каждой итерации, свойство, которое становится крайне важным, когда встречающиеся числа имеют много знаков. Но это также означает, что начальная сходимость метода может оказаться сравнительно медленной, особенно если значение ${\displaystyle X_{0}}$ выбрано плохо.

Для подзадачи выбора начальной оценки ${\displaystyle X_{0}}$ удобно применить сдвиг делителя ${\displaystyle D}$, чтобы он оказался между ${\displaystyle 0{,}5\leqslant D\leqslant 1}$, с применением того же самого сдвига к делимому ${\displaystyle N}$, чтобы частное не изменилось. Тогда можно использовать линейную Аппроксимацию в виде

${\displaystyle X_{0}=T_{1}+T_{2}D\approx {\frac {1}{D}}\,}$

для инициализации метода Ньютона – Рапсона. Чтобы минимизировать максимальную абсолютную ошибку этого приближения на интервале ${\displaystyle [0{,}5,1]}$, следует использовать

${\displaystyle X_{0}={48 \over 17}-{32 \over 17}D.\,}$

Коэффициенты линейного приближения определяются следующим образом. Абсолютное значение ошибки равно ${\displaystyle |\varepsilon _{0}|=|1-D(T_{1}+T_{2}D)|}$. Минимум максимального абсолютного значения ошибки определяется согласно теореме Чебышева об эквивалентных колебаниях^[англ.], применённой к ${\displaystyle F(D)=1-D(T_{1}+T_{2}D)}$. Локальный минимум функции ${\displaystyle F(D)}$ будет в точке, где ${\displaystyle F'(D)=0}$, что имеет решение ${\displaystyle D=-T_{1}/(2T_{2})}$. Функция в этом минимуме должна иметь значение с противоположным знаком по отношению к крайним точкам интервала, а именно, ${\displaystyle F(1/2)=F(1)=-F(-T_{1}/(2T_{2}))}$. Два равенства от двух неизвестных дают единственное решение ${\displaystyle T_{1}=48/17}$ и ${\displaystyle T_{2}=-32/17}$, а максимальная ошибка равна ${\displaystyle F(1)=1/17}$. При использовании этого приближения абсолютное значение ошибки начального значения меньше, чем

${\displaystyle \vert \varepsilon _{0}\vert \leqslant {1 \over 17}\approx 0{,}059.\,}$

Можно образовать многочлен со степенью большей 1, вычислив коэффициенты согласно алгоритму Ремеза. Тогда начальное приближение требует больших вычислений, что может компенсироваться меньшим числом итераций Ньютона – Рапсона.

Поскольку сходимость этого метода в точности квадратичная, отсюда следует, что достаточно

${\displaystyle S=\left\lceil \log _{2}{\frac {P+1}{\log _{2}17}}\right\rceil \,}$

шагов для вычисления значения с точностью до ${\displaystyle P\,}$ двоичных мест. Это эквивалентно 3 для IEEE чисел одинарной точности и 4 для чисел двойной точности и чисел расширенной двойной точности.

Псевдокод

Следующий псевдокод вычисляет частное от деления N и D с точностью до P двоичных знаков:

Выражаем D как ${\displaystyle M\times 2^{e}}$, где ${\displaystyle 1\leqslant M<2}$ (стандартное представление  с плавающей запятой)
${\displaystyle D^{\prime }{\colon }=D/2^{e+1}}$   // приводим к значению между 0,5 и 1, что можно выполнить битовым сдвигом / вычитанием экспоненты
${\displaystyle N^{\prime }{\colon }=N/2*{e+1}}$
${\displaystyle X{\colon }=48/17-32/17\times D^{\prime }}$   // заранее вычисленные константы с той же точностью, что и D
повторить ${\displaystyle \left\lceil \log _{2}{\frac {P+1}{\log _{2}17}}\right\rceil \,}$ раз   // может быть вычислено заранее для фиксированного P
    ${\displaystyle X{\colon }=X+X\times (1-D^{\prime }\times X)}$
конец
возвратить ${\displaystyle N^{\prime }\times X}$

Например, для деления двойной точности с плавающей запятой этот метод использует 10 умножений, 9 сложений и 2 сдвига.

Варианты деления Ньютона – Рапсона

Метод деления Ньютона – Рапсона можно модифицировать, чтобы сделать чуть быстрее. После сдвига N и D, так что D будет в интервале [0,5, 1,0], инициализируем с

${\displaystyle X:={\frac {140}{33}}+D\cdot \left({\frac {-64}{11}}+D\cdot {\frac {256}{99}}\right)}$.

Это лучшее квадратичное приближение к ${\displaystyle {\frac {1}{D}}}$ и оно даёт значение абсолютной ошибки не больше ${\displaystyle 1/99}$. Параметры выбраны так, чтобы выбрать ошибку равной значению третьего порядка многочлена Чебышёва первого рода. Коэффициенты должны быть вычислены заранее и жёстко закодированы в методе.

Тогда в цикле используем итерацию, которая возводит ошибку в куб.

${\displaystyle E:=1-D\cdot X}$

${\displaystyle Y:=X\cdot E}$

${\displaystyle X:=X+Y+Y\cdot E}$

Если цикл выполняется пока ${\displaystyle X}$ не приблизится к ${\displaystyle 1/D}$ до ведущих ${\displaystyle P}$ бит, то число итераций будет не более

${\displaystyle \left\lceil \log _{3}\left({\frac {P+1}{\log _{2}99}}\right)\right\rceil }$

что равно числу раз возведения 99 в куб, чтобы получить ${\displaystyle 2^{P+1}}$. Тогда

${\displaystyle Q:=N\cdot X}$

равно частному с точностью до ${\displaystyle P}$ бит.

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

Деление Гольдшмидта

Деление Гольдшмидта^[7] (Robert Elliott Goldschmidt) использует итеративный процесс многократного умножения делимого и делителя на один и тот же множитель ${\displaystyle F_{i}}$, выбранный так, что делитель сходится к 1. Это приводит к тому, что делимое сходится к частному ${\displaystyle Q}$:

${\displaystyle Q={\frac {N}{D}}{\frac {F_{1}}{F_{1}}}{\frac {F_{2}}{F_{2}}}{\frac {F_{\ldots }}{F_{\ldots }}}}$

Шаги деления Гольдшмидта:

1. Генерируем приближение для множителя ${\displaystyle F_{i}}$.
2. Умножаем делимое и делитель на ${\displaystyle F_{i}}$.
3. Если делитель достаточно близок к 1, возвращаем делимое, в противном случае возвращаемся к шагу 1.

Предполагаем, что ${\displaystyle N/D}$ масштабировано к значению ${\displaystyle 0<D<1}$, каждый ${\displaystyle F_{i}}$ основан на ${\displaystyle D}$:

${\displaystyle F_{i+1}=2-D_{i}.}$

Умножаем делимое и делитель на множитель и получаем:

${\displaystyle {\frac {N_{i+1}}{D_{i+1}}}={\frac {N_{i}}{D_{i}}}{\frac {F_{i+1}}{F_{i+1}}}.}$

После достаточного числа ${\displaystyle k}$ итераций ${\displaystyle Q=N_{k}}$.

Метод Гольдшмидта используется в процессорах AMD Athlon и более поздних моделях^[8][9]. Он известен также как Anderson Earle Goldschmidt Powers (AEGP) алгоритм и реализован в различных процессорах IBM[10][./Алгоритм_деления#cite_note-_084455512bfde3b8-10 [10]]^[11]. Хотя сходимость метода такая же, как у реализации Ньютона – Рапсона, одно из преимуществ метода Гольдшмидта в том, что умножения в числителе и знаменателе можно делать параллельно^[11].

Биномиальная теорема

Метод Гольдшмидта можно использовать с множителями, которые позволяют упрощение с помощью бинома Ньютона. Предположим, что N/D умножено на степень двойки так, что ${\displaystyle D\in ({\tfrac {1}{2}},1]}$. Мы выбираем ${\displaystyle D=1-x}$ и ${\displaystyle F_{i}=1+x^{2^{i}}}$. Это даёт

${\displaystyle {\frac {N}{1-x}}={\frac {N\cdot (1+x)}{1-x^{2}}}={\frac {N\cdot (1+x)\cdot (1+x^{2})}{1-x^{4}}}=\cdots =Q'={\frac {N'=N\cdot (1+x)\cdot (1+x^{2})\cdot \cdot \cdot (1+x^{2^{(n-1)}})}{D'=1-x^{2^{n}}\approx 1}}}$.

После ${\displaystyle n}$ шагов ${\displaystyle (x\in [0,{\tfrac {1}{2}}))}$ знаменатель ${\displaystyle 1-x^{2^{n}}}$ можно округлить до ${\displaystyle 1}$ с относительной ошибкой

${\displaystyle \varepsilon _{n}={\frac {Q'-N'}{Q'}}=x^{2^{n}}}$,

которая имеет максимальную величину ${\displaystyle 2^{-2^{n}}}$ при ${\displaystyle x={1 \over 2}}$, обеспечивая минимальную точность в ${\displaystyle 2^{n}}$ двоичных цифр.

Методы для больших чисел

Методы, предназначенные для аппаратной реализации, обычно не рассчитаны на целые числа с тысячами или миллионами десятичных цифр, что часто встречается, например, при вычислениях по сравнению по модулю в криптографии. Для этих больших чисел более эффективные алгоритмы преобразуют задачу к использованию малого числа умножений, которые можно сделать с помощью асимптотически эффективных алгоритмов умножения^[англ.], таких как Алгоритм Карацубы, умножение Тума – Кука^[англ.], или алгоритм Шёнхаге — Штрассена. В результате вычислительная сложность деления будет того же порядка (с точностью до умножения на константу) как и умножения. Примеры включают сведение к умножению по методу Ньютона как описанное выше^[12], так и чуть более быстрое деление Бурникеля – Циглера^[13], алгоритмы Барета^[англ.] и Монтгомери^[14]. Метод Ньютона, в частности, эффективен в сценариях, где нужно делить на определённое число несколько раз, поскольку после начального нахождения обратного значения только одно (сокращённое) умножение нужно для каждого умножения.

Деление на константу

Деление на константу ${\displaystyle D}$ эквивалентно умножению на её обратную величину. Поскольку знаменатель постоянен, постоянна и обратная величина ${\displaystyle {\frac {1}{D}}}$. Тогда можно вычислить значение ${\displaystyle 1/D}$ один раз и во время вычислений осуществляем умножение ${\displaystyle N\times {\frac {1}{D}}}$ вместо деления ${\displaystyle N/D}$. В арифметике с плавающей запятой использование ${\displaystyle {\frac {1}{D}}}$ вызывают небольшую проблему, связанную с оптимизацией кода компиляторами^[a], но в арифметике целых чисел остаток всегда будет равен нулю, при условии ${\displaystyle |D|>1}$.

Нет необходимости использовать именно ${\displaystyle 1/D}$, можно использовать любое значение ${\displaystyle X/Y}$ которое сводится к ${\displaystyle 1/D}$. Например, при делении на 3, можно использовать дроби ${\displaystyle 1/3}$, ${\displaystyle 2/6}$, ${\displaystyle 3/9}$ или ${\displaystyle 194/582}$. Следовательно, когда ${\displaystyle Y}$ является степенью двойки, шаг деления можно свести к быстрому правому сдвигу. Как эффект, вычисление ${\displaystyle N/D}$ как ${\displaystyle {\frac {N\times X}{Y}}}$ заменяет деление на умножение и сдвиг. Заметим, что именно такой порядок действий здесь существенен, поскольку ${\displaystyle N\times {\frac {X}{Y}}}$ даст в результате нуль.

Однако, если ${\displaystyle D}$ уже является степенью двойки, нет ${\displaystyle X}$ и ${\displaystyle Y}$, удовлетворяющих условиям выше. К счастью, ${\displaystyle {\frac {N\times X}{Y}}}$ даёт в точности тот же результат в арифметике целых чисел, что и ${\displaystyle N/D}$, даже если ${\displaystyle X/Y}$ не в точности равно ${\displaystyle 1/D}$, но «достаточно близко», так что ошибка, вносимая при аппроксимации находится в битах, которые уйдут после операции сдвига^{[15][16][17][b]}

В качестве конкретного примера арифметики с фиксированной запятой, для 32-битных беззнаковых целых, деление на 3 может быть заменено на умножение на ${\displaystyle {\tfrac {2863311531}{2^{33}}}}$, то есть умножением на 2863311531 (шестнадцатеричное 0xAAAAAAAB) с последующим сдвигом на 33 бита вправо. Значение 2863311531 вычислено как ${\displaystyle {\tfrac {2^{33}}{3}}}$, затем округлено. Аналогично деление на 10 может быть выражено как умножение на 3435973837 (0xCCCCCCCD) с последующим делением на ${\displaystyle 2^{35}}$ (или сдвигом вправо на 35 бит)^[19]. OEIS даёт последовательность констант для умножения как A346495 и для правого сдвига как A346496.

Для общего деления ${\displaystyle x}$-битного беззнакового целого деления, где делитель ${\displaystyle D}$ не является степенью двойки, следующее тождество преобразует деление на два ${\displaystyle x}$-битных сложения/вычитания, одно умножение ${\displaystyle x}$-битного на ${\displaystyle x}$-битное чисел (где только старшая половина результата используется) и несколько сдвигов, предварительного вычисления ${\displaystyle k=x+\lceil \log _{2}{D}\rceil }$ и ${\displaystyle a=\left\lceil {\tfrac {2^{k}}{D}}\right\rceil -2^{x}}$:

${\displaystyle \left\lfloor {\frac {N}{D}}\right\rfloor =\left\lfloor {\frac {\left\lfloor {\frac {N-b}{2}}\right\rfloor +b}{2^{k-x-1}}}\right\rfloor }$ где ${\displaystyle b=\left\lfloor {\frac {Na}{2^{x}}}\right\rfloor }$

В некоторых случаев деление на константу может быть совершено даже за меньшее время путём замены «умножения на константу» в серию сдвигов и сложений или вычитаний^[20]. Особый интерес представляет деление на 10, для которого получается точное частное с остатком, если потребуется^[21].

См. также

• Метод галеры
• Длинная арифметика — Алгоритмы деления
• Ошибка Pentium FDIV

Примечания

1. Несмотря на то, как «мала» проблема, вызываемая оптимизацией, эта оптимизация обычно спрятана за флаг «быстрая математика» в современных компиляторах.
2. Современные Компиляторы обычно осуществляют эту оптимизацию целого умножения-и-сдвига. Для констант, известных только в момент выполнения программа должна реализовать оптимизацию сама^[18]

1. Morris, Iniewski, 2017.
2. Flynn Stanford EE486 (Advanced Computer Arithmetic Division) – Chapter 5 Handout (Division)  (неопр.). Stanford University. Дата обращения: 10 января 2022. Архивировано 18 апреля 2022 года.
3. Harris, Oberman, Horowitz, 1998.
4. McCann, Pippenger, 2005, с. 1279–1301.
5. Statistical Analysis of Floating Point Flaw  (неопр.). Intel Corporation (1994). Дата обращения: 22 октября 2013. Архивировано 23 октября 2013 года.
6. Oberman, Flynn, 1995.
7. Goldschmidt, 1964.
8. Oberman, 1999, с. 106-115.
9. Soderquist, Leeser, 1997, с. 56-66.
10. Anderson, Earle, Goldschmidt, Powers, 1997.
11. Guy, Peter, Ferguson, 2005, с. 118–139.
12. Hasselström, 2003.
13. Burnikel, Ziegler, 1998.
14. Barrett, 1987, с. 311-323.
15. Granlund, Montgomery, 1994, с. 61–72.
16. Möller, Granlund, 2011, с. 165–175.
17. ridiculous_fish. "Labor of Division (Episode III): Faster Unsigned Division by Constants" Архивная копия от 8 января 2022 на Wayback Machine. 2011.
18. ridiculous_fish libdivide, optimized integer division  (неопр.). Дата обращения: 6 июля 2021. Архивировано 23 ноября 2021 года.
19. Warren Jr., 2013, с. 230-234.
20. LaBudde, Robert A.; Golovchenko, Nikolai; Newton, James; Parker, David; Massmind: Binary Division by a Constant Архивная копия от 9 января 2022 на Wayback Machine
21. Vowels, 1992, с. 81–85.

Литература

• R. A. Vowels. Division by 10 // Australian Computer Journal. — 1992. — Т. 24, вып. 3.
• James E. Morris, Krzysztof Iniewski. Nanoelectronic Device Applications Handbook. — CRC Press, 2017. — ISBN 978-1-351-83197-0.
• David L. Harris, Stuart F. Oberman, Mark A. Horowitz. SRT Division: Architectures, Models, and Implementations. — Stanford University, 1998.
• Mark McCann, Nicholas Pippenger. SRT Division Algorithms as Dynamical Systems // SIAM Journal on Computing. — 2005. — Т. 34, вып. 6. — doi:10.1137/S009753970444106X.
• Stuart F. Oberman, Michael J. Flynn. An Analysis of Division Algorithms and Implementations. — Stanford University, 1995.
• Robert E. Goldschmidt. Applications of Division by Convergence. — M.I.T., 1964. — (M.Sc. dissertation).
• Stuart F. Oberman. Floating Point Division and Square Root Algorithms and Implementation in the AMD-K7 Microprocessor // Proceedings of the IEEE Symposium on Computer Arithmetic. — 1999. — С. 106–115. — doi:10.1109/ARITH.1999.762835.
• Peter Soderquist, Miriam Leeser. Division and Square Root: Choosing the Right Implementation // IEEE Micro. — 1997. — July–August (т. 17, вып. 4). — doi:10.1109/40.612224.
• S. F. Anderson, J. G. Earle, R. E. Goldschmidt, D. M. Powers. The IBM 360/370 model 91: floating-point execution unit // IBM Journal of Research and Development. — 1997. — Январь.
• Even Guy, Siedel Peter, Warren Ferguson. A parametric error analysis of Goldschmidt's division algorithm // Journal of Computer and System Sciences. — 2005. — Февраль (т. 70, вып. 1). — doi:10.1016/j.jcss.2004.08.004.
• Karl Hasselström. Fast Division of Large Integers: A Comparison of Algorithms. — Royal Institute of Technology, 2003. — ((thesis) M.Sc. in Computer Science). Архивировано 8 июля 2017 года.
• Christoph Burnikel, Joachim Ziegler. Fast Recursive Division. — Max-Planck-Institut für Informatik, 1998.
• Paul Barrett. Implementing the Rivest Shamir and Adleman public key encryption algorithm on a standard digital signal processor // Proceedings on Advances in cryptology---CRYPTO '86. — London, UK: Springer-Verlag, 1987. — С. 311–323. — ISBN 0-387-18047-8.
• Torbjörn Granlund, Peter L. Montgomery. Division by Invariant Integers using Multiplication // SIGPLAN Notices. — 1994. — Июнь (т. 29, вып. 6). — doi:10.1145/773473.178249.
• Niels Möller, Torbjörn Granlund. Improved Division by Invariant Integers // IEEE Transactions on Computers. — 2011. — Февраль (т. 60, вып. 2). — doi:10.1109/TC.2010.143.
• Henry S. Warren Jr. Hacker's Delight. — 2. — Addison Wesley - Pearson Education, Inc., 2013. — ISBN 978-0-321-84268-8.

Литература для дальнейшего чтения

• John J. G. Savard. Advanced Arithmetic Techniques  (неопр.). quadibloc (2018). Дата обращения: 16 июля 2018. Архивировано 3 июля 2018 года.