Форум на Исходниках.RU

Арифметика с длинными числами

Небольшое предисловие

Сразу оговорюсь, что я не ставил своей целью работу с числами любого размера, для простоты взял так называемые long числа, т.е. для 16-ти разрядного кода - это 32-х битные числа, а для 32-х разрядного кода - 64-х битные. Все примеры даны для 16-ти битного кода, но при желании легко переводятся под 32-х битные системы. И если понять теорию, то несложно применить данные алгоритмы для работы с числами любой разрядности, по крайней мере, я старался не просто предоставить код, а пояснить саму суть алгоритмов, дать теоретическое обоснование, да еще и на простом и всем понятном русском языке :) Что нам понадобится? Естественно знание ассемблера, хотя бы на уровне синтаксиса, ну и три класса образования :)
В примерах будут использоваться числа, определенные в сегменте данных.
Для удобства и экономии места объявим их заранее:

ideal

P386

model tiny

locals __

dataseg

label X word

dd 0000FFFEh

label Y word

dd 00010003h

label Z word

qword 0

Shift db 5

Хранение и преобразования

Хранить числа будем так же, как и принято в x86 архитектуре - младший байт по младшему адресу. Таким образом, загрузка числа в регистры выглядит следующим образом:

mov ax, [X] ; AX - младшая часть числа

mov dx, [X+2] ; DX - старшая часть числа

; DX:AX = X

Сохранение числа из DX:AX в память:

mov [X], ax ; записываем младшую часть

mov [X+2], dx ; записываем старшую часть

; X = DX:AX

Зачастую может понадобиться получить длинное число из числа меньшей разрядности. Здесь есть два варианта. Для без знаковых чисел достаточно очистить старшую часть:

mov ax, [X] ; AX - исходное число

xor dx, dx ; у беззнакового старшие биты(DX) равны 0

; DX:AX = X

Расширение чисел со знаком осуществляется немного по-другому, поскольку знак числа определяется его старшим битом, которым и нужно заполнить старшую часть числа. Самый быстрый вариант:

mov ax, [X] ; AX - исходное число

mov dx, ax

sar dx, 15 ; if (AX >= 0) DX=0 else DX=sign, т.е. -1

; DX:AX = X

Менее быстрый, но самый короткий вариант:

mov ax, [X] ; AX - исходное число

cwd ; перенос в DX старшего(знакового) бита AX

; DX:AX = X

Сравнение чисел

Для начала рассмотрим сравнение чисел без знака. Но здесь понадобится некоторое знание команд сравнения и условных переходов. Сравнение производится по простому принципу, сначала сравниваем старшую часть, а затем, если это еще необходимо, младшую. Проще всего понять это на примере. Допустим, есть некое число в регистрах DX:AX и если оно больше числа X, то нужно выполнить какие либо действия:

; if DX:AX > [X]:

cmp dx, [X+2] ; сравниваем старшие части

jne __islater

cmp ax, [X] ; сравниваем младшие части

__islater:

jbe __later ; <- переход при DX:AX <= [X]

; сюда попадем при DX:AX > [X]

; соответственно делаем что нам нужно

....

__later:

; продолжаем программу

А теперь по порядку. Почему прыгаем на второй переход, после первого сравнения? Да все очень просто. Если DX меньше [X+2], то будет установлен флаг переноса CF, поэтому второй переход перенаправит программу на метку __later. Если же DX больше [X+2], то флаг CF будет нулевым и переход не сработает, мы попадем на код, который выполняется при DX:AX > [X]. Ну и если старшие части равны, то первый переход не сработает, и все будет зависеть от сравнения младших частей чисел. Далее приведу только примеры, поскольку все они работают по такому же принципу.

; if DX:AX >= [X]:

cmp dx, [X+2] ; сравниваем старшие части

jne __islater

cmp ax, [X] ; сравниваем младшие части

__islater:

jb __later ; <- переход при DX:AX < [X]

; сюда попадем при DX:AX >= [X]

; соответственно делаем что нам нужно

....

__later:

; продолжаем программу

; if DX:AX < [X]:

cmp dx, [X+2] ; сравниваем старшие части

jne __isgreat

cmp ax, [X] ; сравниваем младшие части

__isgreat:

jae __great ; <- переход при DX:AX >= [X]

; сюда попадем при DX:AX < [X]

; соответственно делаем что нам нужно

....

__great:

; продолжаем программу

; if DX:AX <= [X]

cmp dx, [X+2] ; сравниваем старшие части

jne __isgreat

cmp ax, [X] ; сравниваем младшие части

__isgreat:

ja __great ; <- переход при DX:AX > [X]

; сюда попадем при DX:AX <= [X]

; соответственно делаем что нам нужно

....

__great:

; продолжаем программу

; if DX:AX = [X]:

cmp dx, [X+2] ; сравниваем старшие части

jne __isnoteq

cmp ax, [X] ; сравниваем младшие части

__isnoteq:

jne __noteq ; <- переход при DX:AX <> [X]

; сюда попадем при DX:AX = [X]

; соответственно делаем что нам нужно

....

__noteq:

; продолжаем программу

; if DX:AX <> [X]:

cmp dx, [X+2] ; сравниваем старшие части

jne __isequ

cmp ax, [X] ; сравниваем младшие части

__isequ:

je __equ ; <- переход при DX:AX = [X]

; сюда попадем при DX:AX <> [X]

; соответственно делаем что нам нужно

....

__equ:

; продолжаем программу

Для сравнения чисел со знаком применяется та же схема, только используются другие флаги, а именно SF и OF. Тот же пример:

; if DX:AX > [X]:

cmp dx, [X+2] ; сравниваем старшие части

jne __islater

cmp ax, [X] ; сравниваем младшие части

__islater:

jle __later ; <- переход при DX:AX <= [X]

; сюда попадем при DX:AX > [X]

; соответственно делаем что нам нужно

....

__later:

; продолжаем программу

Остальной код приводить не буду, он практически идентичный, только вместо команд условного перехода для чисел без знака используются команды условного перехода для знаковых чисел. Т.е. достаточно выполнить следующую замену:

jb -> jl

jbe -> jle

ja -> jg

jae -> jge

Логические и арифметические сдвиги

Логический сдвиг чисел можно реализовать несколькими способами.
1. совместимый со всеми моделями процессоров.
Сдвиг влево:

mov ax, [X] ; AX - младшая часть

mov dx, [X+2] ; DX - старшая часть

mov cl, [Shift] ; CL - величина сдвига

cmp cl, 16

jae __short

shl dx, cl ; сдвигаем старшую часть числа

mov bx, ax ; сохраняем младшую часть, для ее сдвига в старшую

shl ax, cl ; сдвигаем младшую часть числа

sub cl, 16

neg cl ; CL = 16-CL

shr bx, cl ; выделяем старшие биты младшего слова

or dx, bx ; и заносим их в младшие биты старшего слова

ret ; DX:AX = DX:AX shl CL

__short:

and cl, 0fh ; CL = CL-16, но можно и убрать, процессор сам обрежет

mov dx, ax ; SHL DX:AX, 16

xor ax, ax ; младшая часть будет нулевой

shl dx, cl ; сдвигаем старшую часть на CL-16 бит

ret ; DX:AX = DX:AX shl CL

Сдвиг вправо:

mov ax, [X] ; AX - младшая часть

mov dx, [X+2] ; DX - старшая часть

mov cl, [Shift] ; CL - величина сдвига

cmp cl, 16

jae __short

shr ax, cl ; сдвигаем младшую часть числа

mov bx, dx ; сохраняем старшую часть, для ее сдвига в младшую

shr dx, cl ; сдвигаем старшую часть числа

sub cl, 16

neg cl ; CL = 16-CL

shl bx, cl ; выделяем младшие биты старшего слова

or ax, bx ; и заносим их в старшие биты младшего слова

__short:

and cl, 0fh ; CL = CL-16, но можно и убрать, процессор сам обрежет

mov ax, dx ; SHR DX:AX, 16

xor dx, dx ; старшее слово будет нулевым

shr ax, cl ; сдвигаем младшее слово на CL-16 бит

ret ; DX:AX = DX:AX shr CL

2. быстрый вариант, для 386+
Сдвиг влево:

mov ax, [X] ; AX - младшая часть

mov dx, [X+2] ; DX - старшая часть

mov cl, [Shift] ; CL - величина сдвига

cmp cl, 16

jae __short

shld dx, ax, cl ; сдвигаем старшее слово с выдвижением в него старших бит младшего

shl ax, cl ; теперь сдвигаем младшее слово

ret ; DX:AX = DX:AX shl CL

__short:

and cl, 0fh ; CL = CL-16, но можно и убрать, процессор сам обрежет

mov dx, ax ; SHL DX:AX, 16

xor ax, ax ; младшая часть будет нулевой

shl dx, cl ; сдвигаем старшую часть на CL-16 бит

ret ; DX:AX = DX:AX shl CL

Сдвиг вправо:

mov ax, [X] ; AX - младшая часть

mov dx, [X+2] ; DX - старшая часть

mov cl, [Shift] ; CL - величина сдвига

cmp cl, 16

jae __short

shrd ax, dx, cl ; сдвигаем младшее слово с выдвижением в него младших бит старшего

shr dx, cl ; теперь сдвигаем старшее слово

ret ; DX:AX = DX:AX shr CL

__short:

and cl, 0fh ; CL = CL-16, но можно и убрать, процессор сам обрежет

mov ax, dx ; SHR DX:AX, 16

xor dx, dx ; старшее слово будет нулевым

shr ax, cl ; сдвигаем младшее слово на CL-16 бит

ret ; DX:AX = DX:AX shr CL

Но при логическом сдвиге вправо мы получаем только положительные числа, поскольку старшие биты обнуляются. Поэтому для сдвига знаковых чисел необходимо учитывать старший бит. Это делается с помощью команды SAR (Shift Arithmetic operand Right), которая при сдвиге заполняет освобождаемые биты старшим битом операнда.
Арифметический сдвиг вправо:

mov ax, [X] ; AX - младшая часть

mov dx, [X+2] ; DX - старшая часть

mov cl, [Shift] ; CL - величина сдвига

cmp cl, 16

jae __short

shr ax, cl ; сдвигаем младшую часть числа

mov bx, dx ; сохраняем старшую часть, для ее сдвига в младшую

sar dx, cl ; сдвигаем старшую часть числа

sub cl, 16

neg cl ; CL = 16-CL

shl bx, cl ; выделяем младшие биты старшего слова

or ax, bx ; и заносим их в младшее слово

ret ; DX:AX = DX:AX sar CL

__short:

and cl, 0fh ; CL = CL-16, но можно и убрать, процессор сам обрежет

mov ax, dx ; SHR DX:AX, 16

sar dx, 15 ; if (DX >= 0) DX=0 else DX=sign, т.е. -1

sar ax, cl ; сдвигаем младшее слово на CL-16 бит

ret ; DX:AX = DX:AX sar CL

Или покороче, для 386+:

mov ax, [X] ; AX - младшая часть

mov dx, [X+2] ; DX - старшая часть

mov cl, [Shift] ; CL - величина сдвига

cmp cl, 16

jae __short

shrd ax, dx, cl ; сдвигаем младшее слово с выдвижением в него младших бит старшего

sar dx, cl ; теперь сдвигаем старшее слово

ret ; DX:AX = DX:AX sar CL

__short:

and cl, 0fh ; CL = CL-16, но можно и убрать, процессор сам обрежет

mov ax, dx ; SHR DX:AX, 16

sar dx, 15 ; if (DX >= 0) DX=0 else DX=sign, т.е. -1

sar ax, cl ; сдвигаем младшее слово на CL-16 бит

ret ; DX:AX = DX:AX sar CL

Думаю здесь все ясно, и объяснять ничего не надо, сдвиг он и в Африке сдвиг :)

Сложение и вычитание

Сложение реализуется по тому же принципу, по которому производится сложение столбцом в школе. Для тех, кто забыл, как это делается, напомню:

---

165

т.е. сложение можно выразить как два независимых сложения:

8 90

+ +

7 60

--- ---

15 150

150

---

165

Таким образом, сложение можно реализовать достаточно просто - сначала складываем младшие части чисел (единицы), а затем старшие (десятки), но с учетом возможного переноса предыдущего сложения. Допустим, что есть два числа:

DX:AX - 0000FFFEh

CX:BX - 00010003h

Складываем младшие части:

0FFFEh ; AX

00003h ; BX

------

10001h

Результат получился 17-ти битный, но поскольку у нас части чисел не могут превышать 16-ти бит, то в качестве результата получим младшие 16-ть бит, а 17-й отправится в регистр флагов, в бит CF - Carry Flags (флаг переноса). Теперь складываем старшие части, но уже с учетом 17-го бита предыдущего сложения, для этого удобнее всего использовать команду ADC - ADd with Carry, которая аналогична простой ADD, но добавляет к сумме значение флага CF:

00000h ; DX

00001h ; CX

1h ; флаг CF - перенос от предыдущего сложения

-----

00002h

Таким образом, получили правильный результат - 20001h. В теории выглядит немного страшно, но реализуется элементарно:

mov ax, [X] ; AX - младшая часть

mov dx, [X+2] ; DX - старшая часть

add ax, [Y] ; складываем младшие части

adc dx, [Y+2] ; складываем старшие части, с учетом флага переноса CF.

; DX:AX = X+Y

Вычитание чисел производится похожим образом, но с учетом заема, для чего служит команда SBB - SuBtract with Borrow:

mov ax, [X] ; AX - младшая часть

mov dx, [X+2] ; DX - старшая часть

sub ax, [Y] ; вычитаем младшие части

sbb dx, [Y+2] ; вычитаем старшие части, с учетом

; возможного заема предыдущего вычитания.

; DX:AX = X-Y

Изменение знака

Для изменения знака числа на противоположное существует команда NEG, но она не работает с большими числами, поэтому придется реализовать свой вариант NEG, реализованный по тому же принципу что и оригинальная команда. Сама команда работает следующим образом:

not [X] ; инвертируем все биты числа

add [X], 1 ; прибавляем еденицу

Отсюда и получаем несложную реализацию:

mov ax, [X] ; AX - младшая часть

mov dx, [X+2] ; DX - старшая часть

not ax ; инвертируем число

not dx

add ax, 1 ; добавляем еденицу

adc dx, 0 ; учитываем возможный перенос из сложения младшей части

; DX:AX = -DX:AX

А немного подумав, можно написать и покороче:

mov dx, [X+2] ; DX - старшая часть

mov ax, [X] ; AX - младшая часть

not dx

neg ax

sbb dx, -1 ; DX:AX = -DX:AX

Умножение и деление

Вот и подошли к самой интересной части. Прежде чем рассматривать умножение вспомним, как мы это делали в школе, столбцом:

---

686

5880

+ ----

6566

что можно записать как:

98 98

* *

7 60

--- ----

686 5880

686

5880

----

6566

То есть принцип простой, разбиваем 67 на два числа - 60 и 7, после чего перемножаем их на первое и результаты складываем. А теперь представим, что каждая цифра может содержать значение не от 0 до 9, а от 0 до 65535, т.е. целый 16-ти битный регистр, после чего умножение больших чисел сводится к традиционному умножению столбиком:

DX:AX - 00020001h

CX:BX - 00040003h

00020001h

00040003h

---------

00060003h ; 00020001h*0003h (DX:AX * BX)

00080004h ; 00020001h*0004h (DX:AX * CX)

+ -------------

0008000A0003h

Разумеется, складывать нужно с учетом флага переноса. Ну и, учитывая вышесказанное, реализация будет выглядеть следующим образом:

xor dx, dx

mov [Z+4], dx

mov [Z+6], dx ; обнуляем старшие части результата

mov ax, [X] ; сначала перемножаем младшую часть X

mul [Y] ; с обоими частями Y

mov [Z], ax

mov [Z+2], dx

mov ax, [X] ; перемножаем младшую часть X

mul [Y+2] ; со старшей частью Y

add [Z+2], ax ; добавляем к результату

adc [Z+4], dx ; с учетом возможного переноса из предыдущего сложения

adc [Z+6], 0

mov ax, [X+2] ; теперь перемножаем старшую часть X

mul [Y] ; с обоими частями Y

add [Z+2], ax ; с учетом возможного переноса из предыдущего сложения

adc [Z+4], dx

adc [Z+6], 0

mov ax, [X+2] ; перемножаем старшую часть X

mul [Y+2] ; со старшей частью Y

add [Z+4], ax

adc [Z+6], dx

; Z = X*Y

Можно немного упростить умножение, если у одного из множителей старшая часть равна нулю:

xor dx, dx

mov [Z+4], dx ; обнуляем старшую части результата

mov ax, [X] ; сначала перемножаем младшую часть X

mul [Y]

mov [Z], ax

mov [Z+2], dx

mov ax, [X+2] ; теперь перемножаем старшую часть X

mul [Y]

add [Z+2], ax

adc [Z+4], dx ; с учетом возможного переноса из предыдущего сложения

; Z = X*Y

Осталось рассмотреть самую сложную часть - деление. В отличие от остальных операций здесь все немного сложнее. Первым, что приходит на ум - деление вычитанием:

mov ax, [X] ; AX - младшая часть делимого

mov dx, [X+2] ; DX - старшая часть делимого

mov bx, [Y] ; BX - младшая часть делителя

mov cx, [Y+2] ; CX - старшая часть делителя

xor si, si

xor di, di ; DI:SI - результат (частное)

__loop:

cmp dx, cx ; делимое больше делителя?

jne __hinotequ

cmp ax, bx

__hinotequ:

jb __done ; флаг переноса установится если делимое меньше делителя

sub ax, bx ; вычитаем из делимого делитель

sbb dx, cx

add si, 1 ; увеличиваем частное

adc di, 0

jmp __loop ; повторяем, пока делимое не станет меньше делителя

__done:

; DI:SI - частное (результат)

; DX:AX - остаток

Здесь есть один неприятный момент - при большом делимом и маленьком делителе скорость вычисления будет крайне мала, из-за тысячей циклов. Но есть и другой способ, пришедший еще со времен 8-ми битных процессоров, не имевших аппаратной поддержки деления. Он основан на серии сдвигов и вычитаний, неком аналоге деления столбиком, только с двоичной спецификой. Как мы делим столбиком? Допустим:

2565|4

+---

Сначала подбираем наименьшее делимое. Мы видим, что 2 меньше 4, соответственно берем уже 25, а в частное записываем 0 (но поскольку он первый, то его можно проигнорировать). Но 25 намного больше 4, поэтому мы умножаем 4 на такое число, произведение которого будет как можно ближе к делимому ,но не превышая его, в нашем случае подходит 6:

2565|4

24 +---

-- |6

Теперь вычитаем это произведение из делимого и если оставшееся число все еще больше делителя, то повторяем. В результате получаем следующее:

2565|4

24 +---

-- |641 - частное

---

1 - остаток

Но это для десятичной системы счисления. Для двоичной все будет точно так же, но нам не придется подбирать множитель, поскольку основание системы позволяет использовать только ноль или единицу. Опять проще всего посмотреть на примере:

101000000101|100

100 +----

--- |1010000001 - частное

100

---

00000101

100

---

1 - остаток

Легко заметить, что множитель всегда равен единице, что значительно упрощает нашу задачу, поскольку нам достаточно сравнить очередную часть делимого с делителем и если делимое меньше, то в частное добавляем ноль, иначе единицу и вычитаем делитель. Все это легко реализуется с помощью сдвига и вычитания. Собственно реализация:

mov ax, [X]

mov dx, [X+2] ; dx:ax - делимое

mov bx, [Y]

mov cx, [Y+2] ; cx:bx - делитель

or cx, cx ; делитель меньше (2^16-1) ?

jnz big_div ; <- нет, переходим к варианту деления больших чисел

cmp dx, bx ; если старшая часть делимого меньше делителя,

jb one_div ; то можно обойтись одним делением

; придется делить старшую и младшую части отдельно, поскольку

; в противном случае словим переполнение (результат превысит размер слова)

mov cx, ax ; сохраняем младшую часть делимого

mov ax, dx

xor dx, dx

div bx ; делим старшую часть

xchg ax, cx ; ax - младшая часть делимого

one_div:

div bx ; делим младшую часть

mov si, dx ; si - младшая часть остатка

mov dx, cx ; dx - старшая часть результата

xor di, di ; di - старшая часть остатка

ret

; DX:AX - частное

; DI:SI - остаток

big_div: ; сюда попадаем только при делителе большем 0FFFFh

mov bp, 32 ; 32 цикла, по одному на каждый разряд числа

xor si, si

xor di, si ; di:si - остаток деления

__loop:

shl ax, 1 ; сдвигаем делимое влево, на один разряд

rcl dx, 1

rcl si, 1 ; выдвигаемый с делимого бит заносится в остаток

rcl di, 1

cmp di, cx ; остаток больше делителя?

jne __hinotequ

cmp si, bx

__hinotequ:

jc __next ; флаг переполнения установится если остаток меньше делителя

sub si, bx ; вычитаем из накопителя делитель

sbb di, cx ;

inc ax ; и прибавляем еденицу к делимому

__next:

dec bp ; повторяем, пока не пройдемся по всему числу

jne __loop

; DX:AX - частное

; DI:SI - остаток

В самом начале функции мы проверяем размер делителя и если он умещается в одно слово, то используем аппаратное деление, что намного быстрее программного. Тут есть одна тонкость, связанная с реализацией деления на x86. Допустим, нужно разделить 101111h на 1. Естественно получим 101111h. Но это число не влезет в один 16-ти битный регистр и мы получим липовый результат, да еще и процессор сгенерирует исключение. Поэтому мы отлавливаем такие ситуации, деля число в два приема - сначала старшую часть (в примере - 10h), а затем младшую (1111h), что позволяет получить верный 32-х битный результат.
Остается только добавить, что рассмотренные реализации умножения и деления расчитаны на числа без знака. Можно конечно реализовать и для знаковых чисел, но проще поменять знак на положительный, а после умножения/деления вернуть его обратно.

Заключение

Существуют более оптимальные алгоритмы деления больших чисел, но поскольку они специализированны на числах определенной длины, то здесь не рассматриваются. На этом закончу. Надеюсь, я достаточно ясно все разъяснил :)]]>