Арифметические операции с длинными целыми положительными числами

Стандартные средства языков программирования позволяют оперировать с числами ограниченной длинны. В статье рассмотрены алгоритмы вычисления следующих арифметических операций с целыми положительными числами неограниченной длинны (чисел произвольно высокой точности):
1. Сложение
2. Вычитание
3. Умножение
4. Целочисленное деление и нахождение остатка от деления
5. Возведение в степень
6. Выполнение арифметических операций по модулю

Все рассматриваемые алгоритмы реализованы в модуле 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.

Алгоритм А:

1. Начальная установка: присвоить j :=0, k :=0; (Переменная j будет пробегать позиции различных разрядов, а переменная k будет следить за переносами на каждом шаге.)
2. Сложить цифры: присвоить 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) перенос.
3. Цикл по 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:

1. Начальная установка: присвоить j :=0, k :=0.
2. Вычитание разрядов: присвоить 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) или нет.
3. Цикл по 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 , а затем – суммы этих произведений, сдвинутых на соответствующее число масштабирующих разрядов. Но на компьютере предпочтительнее выполнять сложение параллельно умножению, как это и делается в данном алгоритме.

Алгоритм М :

1. Начальная установка: присвоить всем параметрам w m + n -1 , …, w 1 , w 0 значения “нуль”. Присвоить j :=0;
2. Начальная установка i : присвоить j :=0, k :=0;
3. Умножить и сложить: присвоить сначала t := u i * v j + w i + j + k , а затем w i + j := t mod b и k := t div b . (Здесь разряд переноса k всегда будет принадлежать интервалу 0? k < b );
4. Цикл по i : увеличить i на единицу. Если после этого i < m , то вернуться к шагу М3; в противном случае присвоить w j + m := k ;
5. Цикл по 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:

1. Нормализация: присвоить 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.
2. Начальная установка 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 .
3. Вычислить 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 .
4. Умножить и вычесть: заменить ( 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 .
5. Сохранить разряд частного: присвоить q j := q '.
6. Цикл j : уменьшить j на единицу. Если теперь j ?0, то вернуться к шагу D 3.
7. Денормализация: теперь ( 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 :

1. Инициализация: установить N := n , Y :=1, Z := x .
2. Деление N пополам: (в этот момент x n = YZ N ) установить N := N div 2 и одновременно определить, было ли N четно. Если N было четно, перейти к шагу V5.
3. Умножение Y на Z : установить Y := Z , умноженное на Y .
4. Проверка N =0? Если N =0, то выполнение алгоритма прекращается с Y в качестве результата.
5. Возведение 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).

Автор: Чикин Роман.