Арифметические операции с длинными целыми положительными числами
- Стандартные средства языков программирования позволяют оперировать с числами ограниченной длинны. В статье рассмотрены алгоритмы вычисления следующих арифметических операций с целыми положительными числами неограниченной длинны (чисел произвольно высокой точности):
- Сложение
- Вычитание
- Умножение
- Целочисленное деление и нахождение остатка от деления
- Возведение в степень
- Выполнение арифметических операций по модулю
Все рассматриваемые алгоритмы реализованы в модуле Arifmetic . pas с использованием Delphi и языка Assembler .
Формат представления чисел
Вначале необходимо рассмотреть способ представления чисел произвольно высокой точности.
Позиционное представление с основанием (или по основанию) b определяется правилом:
(… а 3 а 2 а 1 а 0 ) b =… +a 3 b 3 +a 2 b 2 +a 1 b 1 +a 0
Традиционная десятичная система счисления – это частный случай, когда b=10 и когда значения а выбираются из “десятичных цифр” 0,1,2,3,4,5,6,7,8,9. Простейшее обобщение десятичной системы счисления получается, когда в качестве b берут любое целое число, большее единицы, и числа а – это целые числа из интервала 0?а к < b.
Цифру а к с большим к называют “более значимой”, чем цифру а к с меньшим к. Соответственно самую левую (ведущую или головную) цифру называют наиболее значимой цифрой, а самую правую (хвостовую) – наименее значимой.
Одним из вариантов хранения чисел произвольно высокой точности является размещение их в массивах. При этом придется воспользоваться позиционным представлением чисел по основанию b. Будем хранить каждую цифру в соответствующем элементе массива, индекс которого соответствует позиции разряда цифры в числе. И далее работать будем с этими цифрами.
Для увеличения скорости обработки и уменьшения количества операторов, нужно выбирать систему счисления с как можно большим основанием b. Причем основание b должно быть целочисленными, а цифры этой системы счисления можно обрабатывать стандартными средствами языков программирования. Т.к. большинство современных компьютеров работают с 32-х разрядными микропроцессорами, оптимальным выбором основания b будет число 2^32.
const NewDl = 64 +1; // число 32-х битных элементов массива //64 внутренних разряда, 64*32 = 2048 бит type VnTip = array [0.. NewDl ] of Cardinal ;
1. Сложение
По данным неотрицательным n - разрядным целым числам ( u n -1 u n -2 … u 0 ) и ( v n -1 v n -2 … v 0 ) по основанию b этот алгоритм вырабатывает сумму ( w n w n -1 w n -2 … w 0 ) b . При этом w n – есть “цифра переноса”, всегда равная 0 или 1.
Алгоритм А:
- Начальная установка: присвоить j :=0, k :=0; (Переменная j будет пробегать позиции различных разрядов, а переменная k будет следить за переносами на каждом шаге.)
- Сложить цифры: присвоить w j :=( u j + v j + k ) mod b и k := ( u j + v j + k ) div b . Т.к. u j + v j + k ?( b -1)+( b -1)+1<2 b , то следовательно, k присваивается значение 0 или 1 в зависимости от того, произошел (1) или не произошел (0) перенос.
- Цикл по j : увеличить j на единицу. Если теперь j < n , то вернуться к шагу А2, иначе – присвоить w j := k и завершить выполнение алгоритма.
Реализация алгоритма А.
Сложение реализовано в трех процедурах:
1. procedure slog( a,b,c:VnTip;var r:VnTip); сложение двух чисел r = a + b .
2. procedure slog2( a,b,c:VnTip;var r:VnTip); сложение двух чисел по модулю r = ( a + b ) mod c , при a < c и b < c .
3. procedure slogMod(a,b,c:VnTip;var r:VnTip); сложение двух чисел по модулю r = ( a + b ) mod c .
procedure slog(a,b,c:VnTip;var r:VnTip); assembler; asm PUSH EDI PUSH ESI PUSH EBX mov ESI , eax //сохраняем смещения входных данных mov EBX , ecx CLD // - регистры DI и SI будут увеличиваться mov edi , r // обнуление массива r xor eax,eax mov ecx,dword ptr NewDl rep stosd // r[1..NewDl]:=0 {-------------------------------------------------} lodsd //razmer:= max из a[0] и b[0] cmp EAX,[EDX] jae @Met1 mov EAX,[EDX] @Met1: mov ECX,eax mov edi , r stosd // записываем в r [0] длину результата (большего из слагаемых) clc {-------------------------------------------------} @ met _ c : // сложение двух чисел inc EDX inc EDX inc EDX inc EDX LODSD ADC EAX ,[ EDX ] // только эта операция из цикла изменяет флаг переноса stosd loop @ met _ c {-------------------------------------------------} JNC @ net _ rash // Обработка возможного увеличения размера результата inc dword ptr r mov EAX,1 stosd @ net _ rash : POP EBX POP ESI POP EDI { Выход из п/п с восстановлением регистров } end;
2. Вычитание
По данным неотрицательным n - разрядным целым числам ( u n -1 u n -2 … u 0 ) b ? ( v n -1 v n -2 … v 0 ) b , записанным по основанию b этот алгоритм находит неотрицательную разность ( w n -1 w n -2 … w 0 ) b .
Алгоритм S:
- Начальная установка: присвоить j :=0, k :=0.
- Вычитание разрядов: присвоить w j :=( u j - v j + k ) mod b и k := ( u j - v j + k ) div b . Т.е. k присваивается -1 или 0 в зависимости от того, происходит заимствование ( u j - v j + k <0) или нет.
- Цикл по j : увеличить j на единицу. Если теперь j < n , то вернуться к шагу А2, противном случае завершить выполнение алгоритма. (Когда выполнение алгоритма завершается, должно получиться k =0; условие k =-1 соответствует единственному случаю, когда ( v n -1 v n -2 … v 0 ) b >( u n -1 u n -2 … u 0 ) b , что противоречит сделанному раньше предположению).
Реализация алгоритма S.
Вычитание реализовано в трех процедурах:
1. procedure razn( a,b,c:VnTip;var r:VnTip); разность двух чисел r = a – b , при a > b .
2. procedure razn2( a,b,c:VnTip;var r:VnTip); разность двух чисел по модулю r = ( a – b ) mod c , при a < c и b < c .
3. procedure raznMod(a,b,c:VnTip;var r:VnTip); разность двух чисел по модулю r = ( a – b ) mod c , при a > b .
procedure razn(a,b,c:VnTip;var r:VnTip); assembler; { условие: a должно быть больше b } asm PUSH EDI PUSH ESI PUSH EBX mov ESI , eax //сохраняем смещения входных данных mov EBX , ecx CLD // - регистры DI и SI будут увеличиваться {-------------------------------------------------} mov edi , r // обнуление массива r xor eax,eax mov ecx,dword ptr NewDl inc ecx rep stosd // r[1..NewDl]:=0 {-------------------------------------------------} lodsd // razmer:=a[0] cmp EAX, dword ptr [EDX] jae @Met1 mov EAX, dword ptr [EDX] @Met1: mov ECX,eax //razmer:= max из a[0] и b[0] mov edi , r stosd // записываем в r[0] длину clc {-------------------------------------------------} @ met _ c : // вычитание двух чисел inc EDX inc EDX inc EDX inc EDX LODSD SBB EAX ,[ EDX ] // только эта операция из цикла // изменяет значение флага переноса stosd loop @ met _ c {-------------------------------------------------} // Учет уменьшения размера разности MOV EBX,r @M_Umen_Razm_Raznosti: CMP DWORD PTR [ EBX ],1 // если размер =1, то уменьшать его нельзя JE @M_PROP_UMEN_Raznosti MOV EAX,DWORD PTR [EBX] CMP DWORD PTR [EBX+EAX*4],0 JNE @M_PROP_UMEN_Raznosti DEC DWORD PTR [EBX] JMP @M_Umen_Razm_Raznosti @M_PROP_UMEN_Raznosti: POP EBX POP ESI POP EDI { Выход из п/п с восстановлением регистров } end;
3. Умножение
По данным неотрицательным n - разрядным целым числам ( u m -1 u m -2 … u 0 ) b и ( v n -1 v n -2 … v 0 ) b , записанным по основанию b этот алгоритм находит их произведение ( w m + n -1 … w 1 w 0 ) b . Общепринятый метод умножения при помощи карандаша и бумаги, основан на нахождении сначала частичных произведений ( u m -1 u m -2 … u 0 ) b * v j для 0? j < n , а затем – суммы этих произведений, сдвинутых на соответствующее число масштабирующих разрядов. Но на компьютере предпочтительнее выполнять сложение параллельно умножению, как это и делается в данном алгоритме.
Алгоритм М :
- Начальная установка: присвоить всем параметрам w m + n -1 , …, w 1 , w 0 значения “нуль”. Присвоить j :=0;
- Начальная установка i : присвоить j :=0, k :=0;
- Умножить и сложить: присвоить сначала t := u i * v j + w i + j + k , а затем w i + j := t mod b и k := t div b . (Здесь разряд переноса k всегда будет принадлежать интервалу 0? k < b );
- Цикл по i : увеличить i на единицу. Если после этого i < m , то вернуться к шагу М3; в противном случае присвоить w j + m := k ;
- Цикл по j : увеличить j на единицу. Если после этого j < n , то вернуться к шагу М2; в противном случае завершить выполнение алгоритма.
Реализация алгоритма М.
Умножение реализовано в трех процедурах:
1. procedure umn( a,b,c:VnTip;var r:VnTip); произведение двух чисел r = a * b , при a [0]+ b [0] <= NewDl (сумма числа разрядов а и b меньше или равна количеству 32-х битных элементов массива).
2. procedure umn2(a,b,c:VnTip;var r:VnTip_x2); произведение двух чисел r = a * b . Результат вычислений может превышать размерность типа VnTip .
3. procedure UmnMod(a,b,c:VnTip;var r:VnTip); произведение двух чисел по модулю r = ( a * b ) mod c .
procedure umn(a,b,c:VnTip;var r:VnTip); assembler; var pp:cardinal; AdrA,AdrB,AdrC:LongWord; RazmA,RazmB:LongWord; asm PUSH EDI PUSH ESI PUSH EBX MOV AdrA,EAX //сохраняем смещения входных данных MOV AdrB,EDX MOV AdrC,ECX MOV EAX,DWORD PTR [EAX] MOV RazmA,EAX MOV EAX,DWORD PTR [EDX] MOV RazmB,EAX CLD // - регистры DI и SI будут увеличиваться {-------------------------------------------------} mov edi,r // обнуление массива r xor eax,eax mov ecx,dword ptr NewDl inc ecx rep stosd // r[1..NewDl]:=0 {-------------------------------------------------} MOV EAX,RazmA // считали размер а ADD EAX,RazmB // предварительное задание размера r MOV EBX, R mov DWORD PTR [EBX], EAX // из размера r будем вычитать // тек. индексы циклов {-------------------------------------------------} MOV EBX,0 @M_Cikl_1: INC EBX MOV ECX,0 MOV PP,0 MOV EDI,EBX @M_Cikl_2: // умножение очередных разрядов множителей INC ECX MOV ESI,AdrA MOV EAX, DWORD PTR [ESI+ ECX*4] MOV ESI,AdrB MUL DWORD PTR [ESI+ EBX*4] ADD EAX,PP ADC EDX,0 MOV PP,0 MOV ESI,R ADD DWORD PTR [ESI+ EDI*4],EAX ADC PP,EDX INC EDI MOV EAX,AdrA CMP ECX, DWORD PTR [EAX] JB @M_Cikl_2 MOV EAX,PP ADD DWORD PTR [ESI+ EDI*4],EAX MOV EAX,AdrB CMP EBX, DWORD PTR [EAX] JB @M_Cikl_1 {-------------------------------------------------} MOV ESI,R MOV ECX,DWORD PTR [ESI] // считали размер @M_PROVERKA_RAZMERA: CMP DWORD PTR [ESI+ECX*4],0 LOOPZ @M_PROVERKA_RAZMERA // ПОКА =0 INC ECX // т.к. ecx умен-ся при выходе из цикла MOV DWORD PTR [ESI],ECX // записали новый размер остатка {-------------------------------------------------} POP EBX POP ESI POP EDI { Выход из п/п с восстановлением регистров } end;
4. Деление
Рассмотрим алгоритм деления в столбик ( m + n )-разрядного целого числа на n -разрядное целое число.
По данным неотрицательным целым числам u =( u m + n -1 … u 1 u 0 ) b и v =( v n -1 … v 1 v 0 ) b , записанным по основанию b , в предположении, что v n -1 ?0 и n >1, находим в системе счисления по основанию b частное n div v = ( q m … q 1 q 0 ) b и остаток n mod v = ( r m … r 1 r 0 ) b .
Алгоритм D:
- Нормализация: присвоить d := b div ( v n -1 +1). Затем присвоить ( u m + n u m + n -1 … u 1 u 0 ) b значение ( u m + n -1 … u 1 u 0 ) b , умноженное на d . Также присвоить ( v n -1 … v 1 v 0 ) b значение ( v n -1 … v 1 v 0 ) b , умноженное на d . В результате нормализации вводится новый разряд u m + n слева от u m + n -1 ; если d =1, то все, что необходимо сделать на этом шаге,- присвоить u m + n :=0.
- Начальная установка j : присвоить j := m . Цикл по j , который состоит из шагов D 2- D 6, представляет собой процесс деления ( u j + n … u j +1 u j ) b на v =( v n -1 … v 1 v 0 ) b , дающий один разряд частного q j .
- Вычислить q ': присвоить q ':= ( u j + n b + u j + n -1 ) div v n -1 , и пусть r ' – остаток, ( u j + n b + u j + n -1 ) mod v n -1 . Проверить выполнение неравенства q '= b или q ' v n -2 > br ' + u j + n -2 . Если оно удовлетворяется, то уменьшить q ' на 1, увеличить r ' на v n -1 и повторить попытку при r '< b .
- Умножить и вычесть: заменить ( u j + n … u j +1 u j ) b на ( u j + n … u j +1 u j ) b - q '( v n -1 … v 1 v 0 ) b .
- Сохранить разряд частного: присвоить q j := q '.
- Цикл j : уменьшить j на единицу. Если теперь j ?0, то вернуться к шагу D 3.
- Денормализация: теперь ( q m … q 1 q 0 ) b есть искомое частное, а для получения искомого остатка достаточно ( u m + n -1 … u 1 u 0 ) b разделить на d .
Реализация алгоритма D.
Операции целочисленного деления и нахождения остатка от деления реализованы в процедурах:
1. procedure otn(a,b:VnTip;var c:VnTip;var ost:VnTip); выполняет целочисленное деление c = a div b и нахождение остатка от деления ost = a mod b .
2. procedure modul(var a:VnTip;b:VnTip); находит остаток от деления а = a mod b .
procedure otn(a,b:VnTip;var c:VnTip;var ost:VnTip); assembler; var vm2:array[0..NewDl] of longword; d,K,P:longword; Bn,Bn1:longword; AdrA,AdrB,AdrC: longword; Ra,Rb,Rc,Rost:LongWord; asm PUSHAD // СОХРАНЕНИЕ РЕГИСТРОВ mov ESI,eax //сохраняем смещения входных данных mov EBX,ecx CLD // - регистры DI и SI будут увеличиваться {-------------------------------------------------------} MOV AdrA,EAX MOV AdrB,EDX // СОХРАНИЛИ АДРЕСА МАССИВОВ MOV AdrC,ECX // Ost - соответственно остаток MOV EAX,DWORD PTR [EAX] // РАЗМЕРНОСТЬ МАССИВОВ MOV Ra,EAX MOV EAX,DWORD PTR [EDX] MOV Rb,EAX {--------------------------------------------------------} mov edi,AdrC // обнуление массива c mov eax,1 STOSD // C[0]:=1 = РАЗМЕР mov ecx,dword ptr NewDl // МАХ размер массива XOR EAX,EAX rep stosd // c[1..NewDl]:=0 mov edi,ost // обнуление массива ost mov eax,1 STOSD // Ost[0]:=1 = РАЗМЕР mov ecx,dword ptr NewDl // МАХ размер массива XOR EAX,EAX rep stosd // ost[1..NewDl]:=0 {----------------------------------------------} // копирование А в ost MOV ESI,AdrA MOV EDI,ost MOV ECX,DWORD PTR [ESI] // считывание размера А INC ECX // число итераций REP MOVSD // ТЕПЕРЬ работаем с OST, а не с А !!!!!! MOV EAX,Rb // РАЗМЕР А или ОСТАТКА CMP Ra,EAX // СРАВНЕНИЕ РАЗМЕРА А(ОСТАТКА) и РАЗМЕРА В JB @M_A_MEN_B //ЕСЛИ РАЗМЕР А МЕНЬШЕ В, //ТО ПРОПУСК ВЫПОЛНЕНИЯ ПРОГРАММЫ {O.K.--------------------------------------------------------} // нормализация // находим d (определяем кол-во сдвигов влево) { X = EAX y = EBX c = ECX n = EDX } MOV EAX,[EDX] LEA EAX,[EDX+EAX*4] MOV EAX,[EAX] mov EDX, 32 MOV ECX, 16 // алгоритм нахождения количества нулей слева от // первого значещего бита @M_C_NOT_ZERO: CMP ECX,0 JE @M_C_ZERO MOV EBX,EAX SHR EBX,CL CMP EBX,0 JE @M_Y_ZERO SUB EDX,ECX MOV EAX,EBX @M_Y_ZERO: SHR ECX,1 JMP @M_C_NOT_ZERO @M_C_ZERO: MOV EBX,EAX SHR EBX,31 SUB EDX,EBX SUB EDX,EAX MOV EAX,EDX MOV D,EAX // Сохраняем d { -------------------------------------------------------} // сдвигаем А на ECX бит влево, при этом появляется // новый разряд A[N+M+1] MOV ESI,OST MOV EDI,ESI LODSD INC EAX // УВЕЛИЧИЛИ РАЗМЕР НА 1, появился новый разряд A[N+M+1] STOSD // и сохранили его { если d=0, то необходимо пропустить сдвиг А и В } CMP D,0 JE @M_PROPYSK_SDVIG_NORMALIZACIYA_A_B_VLEVO MOV EBX , D // Сохранение кол-ва сдвигов DEC EAX // т.к. EDI и ESI увеличины на 4 mov ECX,EAX // число итераций( размер А = N+M+1) -1 LEA EDI,[EDI+EAX*4] // Адрес А[Ra+1] DEC EAX LEA ESI,[ESI+EAX*4] // Адрес А[Ra] STD// ESI EDI уменьшаются @M_NORM_SDVIG_A: PUSH ECX //Вложенные циклы LODSD // считали источник mov ECX,EBX // кол-во сдвигов // EAX - источник SHLD dword ptr [edi],EAX,CL // сдвиг и сохранение A[N+M]>>>A[N+M+1] dec EDI dec EDI dec EDI dec EDI POP ECX LOOP @M_NORM_SDVIG_A MOV ECX,EBX // сдвиг последнего (младшего) разряда А[1] SHL DWORD PTR [EDI],CL {--------------------------------------------------------} // сдвигаем B на ECX бит влево, при этом новый разряд // В[N+1] НЕ появляется! MOV ESI,AdrB LODSD // СЧИТАЛИ РАЗМЕР В = N MOV EBX,D // Сохранение кол-ва сдвигов {Т.Е если размер числа В равен 1, то обычный сдвиг} LEA EDI ,[ ESI + EAX *4+4] // Адрес B [ n ] // +4 т.к // ESI уменьшилось при LODSD LEA ESI ,[ ESI + EAX *4] // Адрес B [ n -1] // источник DEC EAX MOV ECX,EAX // ЧИСЛО 32-БИТНЫХ РАЗРЯДОВ ,( размер B = N) -1 // если размер в равен 1, то есх=0, т.е. число повторений цикла =2^32 // это необходимо учитывать и выполнять проверку перед циклом CMP ECX,0 JE @SDVIG_B_RAZMEROM_RAV_1 @M_NORM_SDVIG_B: PUSH ECX //Вложенные циклы LODSD // считали источник mov ECX,EBX // кол-во сдвигов SHLD dword ptr [edi],EAX,CL // EAX - источник dec EDI // СДВИГ И СОХРАНЕНИЕ A[N+M]>>>A[N+M+1] dec EDI dec EDI dec EDI POP ECX LOOP @M_NORM_SDVIG_B @SDVIG_B_RAZMEROM_RAV_1: MOV ECX,EBX // сдвиг последнего (младшего) разряда B[1] SHL DWORD PTR [EDI],CL //@SDVIG_B_RAV_1: {если размер числа В равен 1, то обычный сдвиг} { если d=0, то необходимо пропустить сдвиг А и В } @M_PROPYSK_SDVIG_NORMALIZACIYA_A_B_VLEVO: {--------------------------------------------------------} // СОХРАНЕНИЕ B[n] И B[n-1] MOV EDX,AdrB MOV EAX,DWORD PTR [EDX] LEA EAX,[EDX+EAX*4] // адрес В [n] MOV EBX,DWORD PTR [EAX] // B[n] MOV Bn,EBX MOV EBX,DWORD PTR [EAX-4] // B[n-1] MOV Bn 1, EBX { --------------------------------------------------------} {задание числа итераций основного цикла} MOV ECX,Ra SUB ECX,Rb INC ECX // J := ( Ra +1) – Rb //J равно размер А минус размер В {-----------------------------------------------------------} // ОСНОВНОЙ ЦИКЛ // цикл повторяется М+1 раз @M_OSNOVNOY_CIKL: // нахождение значения q {======== ШАГ 3 ==========} // q^:= ( A[j+n] * MaxZnach + A[j+n-1] ) div B[n] // r^:= ( A[j+n] * MaxZnach + A[j+n-1] ) mod B[n] MOV ESI,Ost MOV EAX,Rb // считывание n ADD EAX,ECX // = n+j MOV EDX,DWORD PTR [ESI+EAX*4] // А [j+n] MOV EAX,DWORD PTR [ESI+EAX*4-4] // А [j+n-1] {----------------------------------------------------------} {ПРОВЕРКА ВОЗМОЖНОГО ПЕРЕПОЛНЕНИЯ q^} // if q^=>MaxZnach или B[N]<=A[N+J] // то при делении возникает переполнение CMP Bn,EDX // СРАВНИВАЕМ B[N] И A[N+J] JA @M_A_MOGNO_DELIT_NA_B // если B[N] > A[N+J] PUSH ECX MOV EBX,EDX // А [j+n] MOV ECX,EAX // А [j+n-1] XOR EAX,EAX // q^:= 0 DEC EAX // q^:= MaxZn -1 MUL Bn // EDX:EAX = q^ * Bn // EDX д.б. = А[j+n] , ПОЭТОМУ (EBX-EDX)=0 И ОСТАТОК ОЦЕНКИ РАВЕН // r^ = (EBX-EDX)*MaxZn + (ECX-EAX) , ГДЕ EBX=А[j+n-1] // ЕСЛИ (EBX-EDX)>0 И r^>MaxZn , ТО ПРОПУСТИТЬ ОЦЕНКУ q^ SUB ECX,EAX SUB EBX,EDX MOV EDX,EBX // = r^ XOR EAX,EAX // q^:= 0 ////MaxZn -1 DEC EAX // q^:= MaxZn -1 POP ECX CMP EBX,0 JNE @M_SHAG_4 // ЕСЛИ (EBX-EDX)>0 И r^>MaxZn, ТО ПРОПУСТИТЬ ОЦЕНКУ q^ JMP @M_NACHALO_OCENKI_q @M_A_MOGNO_DELIT_NA_B: DIV Bn // деление // теперь q^=EAX // r^=EDX {----------------------------------------------------------} { если размер В = 1, то нужно пропустить оценку q^, при этом необходимо помнить, что Ra+1>Rb и следовательно Ra+1=2 (в этом случае нет разряда A[n+j-2]) возможно только когда Rb=1} CMP Rb,1 JE @M_PROPUSK_OCEKI_q_TK_RAZMER_B_RAV_1 {------- оценка значения q^ ---------------} @M_NACHALO_OCENKI_q: PUSH ECX PUSH EAX // СОХРАНЯЕМ q^ MOV EBX,AdrB MOV ESI,Ost MOV EAX,Rb //DWORD PTR [EBX] // = Rb DEC EAX DEC EAX // = Rb - 2 ADD EAX,ECX // = Rb + J -2 MOV ECX,DWORD PTR [ESI+EAX*4] // = A [Rb+J-2] POP EAX PUSH EAX // = q^ MOV EBX,EDX // = r^ MUL Bn1 // EDX:EAX = q^ * B[ Rb-1 ] // СРАВНИВАЕМ ( q^*B[Rb-1] ) И ( r^*MaxZ + A [Rb+J-2] ) // EDX:EAX И EBX:ECX // если ( q^*B[Rb-1] ) меньше или равно ( r^*MaxZ + A [Rb+J-2] ), // то q - верно CMP EDX,EBX JA @q_NE_VERNO // EDX > EBX JB @M_q_Ok // EDX < EBX CMP EAX,ECX // EDX = EBX JA @q_NE_VERNO // EAX > ECX JMP @M_q_Ok // EAX <= ECX @q_NE_VERNO: POP EAX POP ECX DEC EAX // q^:= q^-1 MOV EDX,EBX // =r^ ADD EDX,Bn // r^:= r^+B[Rb] JNC @M_NACHALO_OCENKI_q // ЕСЛИ НЕ БЫЛО ПЕРЕНОСА CF=0 , Т.Е. r^<MaxZ PUSH ECX PUSH EAX @M_q_Ok: POP EAX POP ECX {-----------------------------------------------------------} @M_SHAG_4: @M_PROPUSK_OCEKI_q_TK_RAZMER_B_RAV_1: {========= ШАГ 4 ===============} { Вычитание q^*B из текущего остатка} { A[ J .. J+Rb+1] - q^ * B[ 1 .. Rb] } PUSH ECX MOV ESI,Ost MOV EDI,AdrB MOV EBX,0 // = I =0, но начинается с 1 MOV P,0 // ПЕРЕНОС = 0 @M_CIKL_VICH_q: INC EBX // = I = 1 .. Rb INC ECX // = I+J = j+1 .. j+Rb PUSH EAX // = q^ MUL DWORD PTR [EDI+EBX*4] // EDX:EAX = q^ * B[I] ADD EAX,P // + PERENOS( ЗАЕМ и СТАРШИЙ ADC EDX,0 // РАЗРЯД ПРЕДЫДУЩЕЙ ИТЕРАЦИИ) MOV P,0 // ОБНУЛЯЕМ ПЕРЕНОС SUB DWORD PTR [ESI+ECX*4-4],EAX // A[J+I-1]:= A[J+I-1] // - (q^*B[I] + PERENOS) ADC P,EDX // СОХРАНЕНИЕ НОВОГО ПЕРЕНОСА С УЧЕТОМ ЗАЕМА POP EAX CMP EBX,Rb // I ? Rb JB @M_CIKL_VICH_q //если I < Rb, то продолжить, иначе выход из цикла PUSH EAX MOV EAX,P // Считываем последний перенос MOV P,0 SUB DWORD PTR [ESI+ECX*4],EAX // ВЫЧЕСТЬ ИЗ СТАРШЕГО РАЗРЯДА A[Rb+J] //ПОСЛЕДНИЙ ПЕРЕНОС ADC P,0 // сохраняем последний возможный заем POP EAX // = q^ POP ECX // = J {========== ШАГ 5 =================} // шаг компенсации CMP P ,0 JE @M_PROPUSK_KOMPENSACII // если CF=0 , т.е. вычитание было без заема PUSH EAX // = q^ PUSH ECX // = J LEA ESI,[ESI+ECX*4] // адрес A[J] ИСТОЧНИК MOV EDI,ESI // адрес A[J] ПРИЕМНИК MOV ECX,Rb //[EDI] // = Rb , число итераций MOV EDX,AdrB //EDI // = Adres B clc // сброс флага переноса CLD // увиличение EDI и ESI @M_CIKL_KOMPENSACIYA: inc EDX // сложение двух чисел inc EDX inc EDX inc EDX LODSD // считали ADC EAX,dword ptr [EDX] stosd // записали сумму loop @M_CIKL_KOMPENSACIYA LODSD ADC EAX,0 // перенос в старший разряд А[J+Rb] STOSD POP ECX // = J POP EAX // = q^ // Переносом в A[J+Rb+1] пренебрегаем, т.к. из него заем не вычитали DEC EAX // q^ := q^ - 1 {-------------------------------------------------} {========= ШАГ 6 ============} @M_PROPUSK_KOMPENSACII: MOV EDX,AdrC // сохранение разряда частного ADD DWORD PTR [EDX+ECX*4],EAX // C[J] := C[J] + q^ // MOV DWORD PTR [EDX+ECX*4],EAX // C[J] := q^ {-------------------------------------------------} {========= ШАГ 7 ============} DEC ECX CMP ECX,0 JA @M_OSNOVNOY_CIKL // КОНЕЦ ОСНОВНОГО ЦИКЛА { если d=0, то необходимо пропустить сдвиг А и В } CMP D,0 JE @M_PROPYSK_SDVIG_NORMALIZACIYA_A_B_VPRAVO // сдвигаем остаток на d бит вправо, mov ECX , DWORD PTR [ ESI ] // // ЦИКЛ = РАЗМЕР остатка DEC ECX XOR EAX,EAX CMP ECX,EAX //0 JE @M_RAZMER_A_MEN_1 @M_SDV_VPRAVO: PUSH ECX MOV ECX,D // число сдвигов вправо INC EAX // EAX + 1 MOV EBX,DWORD PTR [ESI+EAX*4+4] SHRD DWORD PTR [ESI+EAX*4],EBX,CL // сдвиг вправо EAX- го разряда POP ECX LOOP @M_SDV_VPRAVO @M_RAZMER_A_MEN_1: MOV ECX,D SHR DWORD PTR [ESI+EAX*4+4],CL // сдвиг вправо старшего разряда MOV ECX,DWORD PTR [ESI] // считали размер остатка @M_RAZMER_A_POSLE_NORMALIZ_OK: CMP DWORD PTR [ESI+ECX*4],0 LOOPZ @M_RAZMER_A_POSLE_NORMALIZ_OK // ПОКА =0 INC ECX // т.к. ecx умен-ся при выходе из цикла MOV DWORD PTR [ESI],ECX // записали новый размер остатка {-----------------------------------------------------------} // сдвигаем делитель на d бит вправо, MOV ESI,AdrB mov ECX,DWORD PTR [ESI] // // ЦИКЛ = РАЗМЕР делителя DEC ECX // уменьшаем на 1 XOR EAX,EAX CMP ECX,EAX //0 JE @M_RAZMER_B_MEN_1 @M_SDV_VPRAVO_DELIT: PUSH ECX MOV ECX,D // число сдвигов вправо INC EAX // EAX + 1 MOV EBX,DWORD PTR [ESI+EAX*4+4] SHRD DWORD PTR [ESI+EAX*4],EBX,CL // сдвиг вправо EAX- го разряда POP ECX LOOP @M_SDV_VPRAVO_DELIT @M_RAZMER_B_MEN_1: MOV ECX,D SHR DWORD PTR [ESI+EAX*4+4],CL // сдвиг вправо старшего разряда { если d=0, то необходимо пропустить сдвиг А и В } @M_PROPYSK_SDVIG_NORMALIZACIYA_A_B_VPRAVO: {-----------------------------------------------------------} {проверка размерности частного} MOV ESI,AdrA MOV ECX,DWORD PTR [ESI] // считали размер делимого MOV ESI,AdrC @M_PROVERKA_RAZMERA_Chastnogo: CMP DWORD PTR [ESI+ECX*4],0 LOOPZ @M_PROVERKA_RAZMERA_Chastnogo // ПОКА =0 INC ECX // т.к. ecx умен-ся при выходе из цикла MOV DWORD PTR [ESI],ECX // записали новый размер остатка {проверка размерности остатка} MOV ESI,Ost MOV ECX,DWORD PTR [ESI] // считали размер остатка @M_PROVERKA_RAZMERA_A: CMP DWORD PTR [ESI+ECX*4],0 LOOPZ @M_PROVERKA_RAZMERA_A // ПОКА =0 INC ECX // т.к. ecx умен-ся при выходе из цикла MOV DWORD PTR [ESI],ECX // записали новый размер остатка {-----------------------------------------------------------} @M_A_MEN_B: // !!! ДЕЛИМОЕ МЕНЬШЕ ДЕЛИТЕЛЯ POPAD CLD end; { конец процедуры Otn}
5. Возведение в степень
Предположим, например, что необходимо вычислить x 16 . Можно просто начать с x и 15 раз умножить его на x . Но тот же ответ можно получить всего за четыре умножения, если несколько раз возвести в квадрат получающийся результат, последовательно вычисляя x 2 , x 4 , x 8 , x 16 .
Рассмотрим алгоритм, который вычисляет значение x n , где х и n – положительные целые числа.
Алгоритм V :
- Инициализация: установить N := n , Y :=1, Z := x .
- Деление N пополам: (в этот момент x n = YZ N ) установить N := N div 2 и одновременно определить, было ли N четно. Если N было четно, перейти к шагу V5.
- Умножение Y на Z : установить Y := Z , умноженное на Y .
- Проверка N =0? Если N =0, то выполнение алгоритма прекращается с Y в качестве результата.
- Возведение Z в квадрат: установить Z := Z * Z и вернуться к шагу V 2.
Реализация алгоритма V.
Операция возведения в степень реализована в двух процедурах:
1. procedure Stepen(a,n,c:VnTip;var r : VnTip);
2. procedure Stepen2(a,n,c:VnTip;var r : VnTip);
Обе процедуры выполняют возведение числа в степень по модулю
r =( a ^ n ) mod c . Но процедура Stepen 2 выполняется на 12-16% быстрее. Это достигается за счет использования оптимизированных для операции возведения в степень процедур нахождения остатка от деления.
procedure Stepen(a,n,c:VnTip;var r : VnTip); {т.к. при возведении в степень резко возрастает размерность числа, то процедура Stepen вычисляет значение по модулю см. пункт Выполнение арифметических операций по модулю.} var i: integer; y,z,prom,nn: VnTip; begin for i:= 1 to NewDl do y[i]:=0; y[0]:=1; y[1]:=1; nn:=n; z:=a; Repeat if (Mod2(nn)=1) then begin UmnMod(y,z,c,prom); y:= prom; end; Div2(nn); if not A_rav_0(nn) then begin UmnMod(z,z,c,prom); z:=prom; end; Until A_rav_0(nn); r := y ; end ;
6. Выполнение арифметических операций по модулю
Арифметические операции по модулю выполняются также как и обычные арифметические операции, но при этом результат вычислений должен принадлежать числовом промежутку [0 , M ). Где M – целое положительное число – модуль, по которому производятся вычисления. Во всех алгоритмах исходные данные должны быть приведены к промежутку [0 , M ).
При сложении, для выполнения этого требования, если результат превышает модуль, достаточно из результата вычесть значение модуля.
При вычитании, результат вычисления может оказаться отрицательным. Чтобы избежать отрицательного результата, перед выполнением операции, нужно проверить неравенство ( u n -1 u n -2 … u 0 ) b ? ( v n -1 v n -2 … v 0 ) b . Если оно не выполняется, то следует число ( u n -1 u n -2 … u 0 ) b увеличить на значение модуля. Или после вычисления разности, проверять результат, и если он отрицательный, то увеличивать его на значение модуля.
При умножении результат, как правило, превышает значение модуля в несколько раз. Поэтому необходимо применять алгоритм нахождения остатка целочисленного деления.
При выполнении целочисленного деления, в качестве модуля выступает делитель, т.е. остаток и целая часть деления будут удовлетворять данному условию.
Для возведения в степень по модулю достаточно вместо обычной процедуры умножения использовать процедуру умножения по модулю.
Модуль Arifmetic.pas и программу демонстрирующую его возможности можно скачать здесь (zip-архив 16 Kb).