## RSA для программиста

версия 1.1

(x) 2000-2001 Z0MBiE
http://z0mbie.host.sk

### Введение

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

### Алгоритм, в самых общих чертах

Алгоритм из тех, что с открытым и секретным ключом.

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

Подписать сообщение может только владелец секретного ключа, а проверить подпись - владелец открытого.

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

    <<...С  резким  скачком  производительности  вычислительной техники сначала
     столкнулся  алгоритм  RSA,  для вскрытия которого необходимо решать задачу
     факторизации.  В  марте  1994 была закончена длившаяся в течение 8 месяцев
     факторизация  числа  из  129 цифр (428 бит). Для этого было задействовано
     600  добровольцев  и  1600 машин, связанных посредством электронной почты.
     Затраченное машинное время было эквивалентно примерно 5000 MIPS-лет.
         Прогресс  в решении проблемы факторизации во многом связан не только с
     ростом вычислительных мощностей, но и с появлением в последнее время новых
     эффективных алгоритмов. (На факторизацию следующего числа из 130 цифр ушло
     всего  500 MIPS-лет). На сегодняшний день в принципе реально факторизовать
     512-битные   числа.   Если   вспомнить,   что   такие  числа  еще  недавно
     использовались  в программе PGP, то можно утверждать, что это самая быстро
     развивающаяся область криптографии и теории чисел...>>

### Шифровка/расшифровка (теория)

Основная фича в том, что алгоритм оперирует "большими" числами, то есть числами, состоящими из сотен бит.

Открытый ключ (public key) -- это числа e и m. (e от слова encrypt, шифровка)

Секретный ключ (secret key, private key) -- числа d и m. (d от слова decrypt, расшифровка)

(иногда пишут {e,m} и {d,m} - на самом деле это ничего не значит)

Числа m и d - большие. Число e - маленькое, часто принимается равным 3, 9, 17, и т.п. Как видите, число m - одно и то же в обоих ключах (известно всем), а число e, зная m и d, обычно может быть легко найдено. Исходя из этого логично было бы считать m открытым ключом, а d секретным, но, так уж принято.

Согласно RSA, большое число x (причем x < m) зашифровывается в y так:
y = (x ^ e) % m,
а расшифровывается обратно так:
x = (y ^ d) % m

Здесь ^ означает возведение в степень, а % -- остаток от деления, иначе называемый модуль. (не путайте этот модуль с абсолютным значением, оно здесь вообще нигде использоваться не будет) Часто модуль записывается так: a = b (mod m). Читается как "а равно b по модулю m". Значок % используется в Си, 'x mod y' в Паскале, а 'x = y (mod m)' во всякой литературе. Например: 10 mod 3 = 1.

Но вернемся к шифровке. Откуда берется x -- понятно и ежу, это же и есть ваше исходное сообщение, только представленное в виде большого числа. А нахождение чисел m,d и e называется "генерация ключей".

### Генерация ключей (теория)

Для генерации ключей можно использовать старое доброе досовское PGP 2.6.x с параметрами '-kg -d'. Сгенеренные числа m/d/e будут выданы на экран. Если вы твердо решили написать свою реализацию RSA, то лучше, оставив генератор ключей на потом, сразу приступайте к шифровке/расшифровке (подпрограмма modexp).

Итак, генерации RSA-ключей (чисел m/d/e) заключается в следующем:

- Найти большие простые числа p и q, длина произведения (сумма длин) которых (в битах) равна длине искомого RSA-ключа. Простые -- значит не имеющие делителей кроме единицы, например 7, 11, 23. Обычно длины p и q принимаются одинаковыми - например по 128 бит для 256-битного ключа. Хотя есть и модификации алгоритма с существенно разными длинами p/q, для ускорения вычислений и, наверное, затруднения хака.
- Число m принимается равным m = p * q
- Выбирается число e, взаимно простое с произведением (p-1)*(q-1). Взаимно простое -- значит не имеющее общих делителей, кроме единицы. Так как числа p и q - простые, а значит - нечетные (ибо не делятся на 2), то (p-1)*(q-1) -- число четное, и, соответственно, e будет нечетным. При этом e не может быть равно 1, ибо тогда шифровка не имеет смысла, так что минимально допустимое е равно 3.
  Обычно выбирают некоторое начальное e, и пробуют e, e+2, e+4, e+6 и т.п., до тех пор, пока числа e и (p-1)*(q-1) не окажутся взаимно простыми. Записывается это условие так: GCD(e, (p-1)*(q-1)) = 1, что равносильно GCD(e, p-1)=1 и GCD(e, q-1)=1. Здесь GCD (Greatest Common Divisor) есть не что иное как НОД (Наибольший Общий Делитель).
  Должно быть ясно, что если e = 3, то никакого GCD считать нихрена не нужно, а достаточно взять (p-1)*(q-1) и поделить на 3, и если не делится (остаток != 0), то все окей.

  Кроме того, вместо (p-1)*(q-1) можно в разных источниках встретить что-то типа lambda(m) = lcm((p-1)*(q-1)), где LCM (Least Common Multiple) -- НОК (Наименьшее Общее Кратное). Так как множителя два, и друг на друга они, очевидно, не делятся (ибо длины p и q обычно принимаются равными), то никакого lcm считать не надо, достаточно просто умножить p-1 на q-1.

- Число d вычисляется такое, что (e * d) % ((p-1)*(q-1)) = 1.
  Другими словами, d = modinv(e, (p-1)*(q-1)),
  где modinv (modular inverse) -- специальная хитрая функция, кояя пользует так называемый алгоритм Евклида.
  Замечу, что d отнюдь не является уникальным числом; варианты его таковы: d' = d + (p-1)*(q-1) * t, где t=0,1,2,...
  Как видно, t=0 является наилучшим вариантом для достижения наибольшей скорости, однако, если задача, наоборот, максимально замедлить вычисления (бывает и такое, смотри DPGN), то можно выбрать t > 0.

Собственно на этом генерация ключей кончается.

Вникнув в алгоритм, можно видеть, что именно числа p и q и составляют самую секретную фишку. Вы скажете: но ведь m=p*q публикуется? Да. Именно в эффективном алгоритме разложения m на p и q и состоит хак.

### Более быстрая расшифровка

Обычно расшифровку x = (y ^ d) % m заменяют чуть более сложным, но и более быстрым алгоритмом:

     u = modinv(p, q)                -- u,dp,dq -- могут быть вычислены заранее,
     dp = d % (p-1)                  -- при генерации ключей
     dq = d % (q-1)                  --

     a = ((x % p) ^ dp) % p
     b = ((x % q) ^ dq) % q
     if (b < a) b += q
     y = a + p * (((b - a) * u) % q) -- CRT (Chinese Remainder Theorem)

Почему алгоритм более быстрый? -- потому, что оперировать приходится более короткими числами, и в результате производится меньше операций. Хотя, признаюсь, используя этот метод у меня получилось увеличить скорость всего на 10% на C++.

Однако при использовании p и q разных длин, можно получить куда более существенный выигрыш в скорости. Недословно цитирую: "при замене 512-битных p и q на 256 и 768 бит получен выигрыш в скорости в 16 раз".

### шифровка/расшифровка

Прежде всего заметим, что большие числа у нас будут храниться в обычном интелевском формате (младший бит/байт/... сначала).

### modexp()

Теперь рассмотрим вопрос о возведении большого числа в степень:
y = (x ^ e) mod m.
Такая процедура обычно называется modexp.

Совершенно ясно, что показатель степени может быть очень большой, и никто не собирается выполнять соответствующее количество умножений. Вместо этого мы воспользуемся следующей формулой:
x^(a + b) = x^a * x^b
То есть будем рассматривать число e как степени двойки:
пусть e = 32+8+1, тогда x^e = x^32 * x^8 * x^1.
Кроме того, очевидно следующее:
(a*b) % m = (a * (b % m)) % m = ((a % m) * b) % m = ((a % m) * (b % m)) % m
Поэтому (x^(32+8+1)) % m может быть вычислено как
(x^32 * ((x^8 * ((x^1) % m)) % m)) % m
Откуда брать x^i (i-ю степень x) ? - в цикле умножать число x само на себя. В результате алгоритм нашей подпрограммы будет выглядеть так:

     y = 1                         -- результат, (x ^ e) % m
     i = 0                         -- номер текущего бита в показателе степени e
     t = x                         -- временная переменная, t = (x ^ i) % m
     while (1)                     -- повторять всегда, выход по break
     {
         if (e.getbit[i] == 1)     -- если i-й бит числа e установлен в 1
         {
             y = (y * t) % m       -- умножим результат на t = x ^ i
         }
         i = i + 1                 -- следующий бит
         if (i > maxbit(e)) break  -- от 0 и до старшего бита e
         t = (t * x) % m           -- t = t * x, вычисляем следующую степень x^i
     }

Как видно, все, что делает процедура modexp() -- это вызывает процедуру modmult() -- выполняющую c = (a * b) % m -- N раз, где N = e.общее_число_бит + e.число_единичных_бит - 1.

### modmult() - вариант 1

Теперь рассмотрим процедуру modmult():
c = (a * b) % m
По аналогии с modexp()-ом, разложим b на степени двойки, т.е. воспользуемся формулой
a * b = a * (... + b.bit[3]<<3 + b.bit[2]<<2 + b.bit[1]<<1 + b.bit[0]<<0),
здесь << означает shl, двоичный сдвиг влево, т.е. x << y = x * (2^y)
Кроме того, ясно, что:
(a+b) % m = (a + (b % m)) % m = ((a % m) + b) % m = ((a % m) + (b % m)) % m
Таким образом для b=32+8+1, (a * b) % m будет вычислено как
((a<<5)%m + ((a<<3)%m + ((((a<<0)%m) % m)) % m)) % m
Текущее значение a<<i вычисляется в цикле, последовательным сдвигом на 1.

     c = 0                         -- результат, (a * b) % m
     i = 0                         -- номер текущего бита в множителе b
     t = a                         -- временная переменная, t = (x << i) % m
     while (1)
     {
         if (b.getbit[b] == 1)     -- если i-й бит числа b установлен в 1
         {
    //       с = (c + t) % m       -- добавим t = x << i
             с = c + t             -- + t
             if (c >= m) c = c - m -- % m
         }
         i = i + 1                 -- следующий бит
         if (i > maxbit(b)) break  -- от 0 и до старшего бита множителя b
    //   t = (t << 1) % m          -- t=t*2, вычисляем следующую степень a<<i
         t = t << 1                -- << 1
         if (t >= m) t = t - m     -- % m
     }

Вычисление модуля производится последовательным вычитанием m из c, каждый раз когда число с становится больше m. Это становится возможным из-за того, что длина числа c каждый раз увеличивается не больше чем на 1 бит (а значение - в 2 раза), и избавляет от необходимости выполнять нахождение остатка от деления после умножения.

### modmult() - вариант 2

Здесь предлагается другой вариант, умножение столбиком. Такая процедура будет работать намного быстрее за счет встроенной в x86 32-битной команды MUL. Кроме того, что в отличие от предыдущей процедуры потребуется в 32 раза меньше операций (т.к. оперируем двордами, а не битами), насколько я слышал, в x86 процессор заложен еще и так называемый алгоритм быстрого окончания, то есть умножение старших нулевых битов множимого/множителя не происходит. Практика показала, что этот вариант modmultа быстрее в 2 раза.

Вот как выглядит алгоритм умножения столбиком для двух больших чисел:
c = (a * b) % m

1. Собственно умножение "столбиком": t = a * b Очевидно, что разрядность числа t (число бит) будет равно сумме бит чисел a и b, то есть при одинаковой длине последних, длина t будет вдвое больше.

     t = 0                         -- произведение
     for i = 0 to maxdword(a)      -- от 0 и до последнего дворда множимого
     for j = 0 to maxdword(b)      -- от 0 и до последнего дворда множителя
     {
    //   ...:t[i+j+2]:t[i+j+1]:t[i+j] += a[i] * b[j]
         EDX:EAX = a[i] * b[j]     -- умножение командой MUL
         k = i + j                 -- позиция с которой добавлять к произведению
         t[k  ] = t[k  ] + EAX
         t[k+1] = t[k+1] + EDX + CF
         while (CF)                -- добавляем, пока перенос
         {
             t[k+2] = t[k+2] + CF
             k = k + 1
         }
     }

2. вычисление остатка от деления

Это, конечно, плохо, но придется. Итак, сейчас мы имеем t = a * b, и надо посчитать c = t % m.

Вот алгоритм:

     c = 0                         -- результат
     for i = maxbit(t) to 0        -- от старшего бита к младшему
     {
         c = c << 1                -- сдвиг влево на 1 бит, shl
         c = c | t.getbit[i]       -- | значит OR, копируем i-й бит t в 0й бит c
         if (c >= m) c = c - m     -- % m
     }

Идея алгоритма в том, что мы последовательно вычитаем из c m<<maxbit(t), ..., m<<i, ..., m<<2, m<<1, m

### сложение, вычитание, сравнение

Эти операции над большими числами релизуются достаточно просто:

     ; input: ESI=big a
     ;        EDI=big b
     ;        ECX=length in DWORDs
     ; output:CF
     ; action:b = b + a

     bn_add:       clc                     ; CF=0
     __cycle:      lodsd
                   adc     eax, [edi]
                   stosd
                   loop    __cycle
                   retn

     ; input: ESI=big a
     ;        EDI=big b
     ;        ECX=length in DWORDs
     ; output:CF
     ; action:b = b - a

     bn_sub:       clc                     ; CF=0
     __cycle:      lodsd
                   sbb     [edi], eax
                   lea     edi, [edi+4]
                   loop    __cycle
                   retn

     ; input: ESI=big a
     ;        EDI=big b
     ;        ECX=length in DWORDs
     ; output:CF==1 (jb) -- a < b
     ;        ZF==1 (jz) -- a = b
     ;        CF==0 (ja) -- a > b

     bn_cmp:
     __cycle:      mov     eax, [esi+ecx*4-4]
                   cmp     eax, [edi+ecx*4-4]
                   loopnz  __cycle
     __exit:       retn

### Генерация ключей: нахождение больших простых чисел p и q

Заполняем число случайными значениями, и проверяем, является ли оно простым. Если не является, то увеличиваем число и проверяем еще раз. И так до тех пор, пока не задымит процессор.

Используется здесь такая оптимизация: начальное число выбирается нечетным, и затем увеличивается не на 1, а на 2. Это делается для того, чтобы заранее выбросить все четные (и соответственно, составные) числа.

Как узнать, является ли текущее число простым?

Есть два способа: быстрый, но ненадежный и медленный, но надежный. Используются оба.

Способ первый. "решето" (sieve). Используется, чтобы отсеять большинство составных чисел (т.е. не являющихся простыми).

Идея в следующем: делим тестируемое число p на сотню первых простых чисел. Если делится (остаток от деления равен нулю) - число составное.

На самом деле, конечно, никто не будет постоянно делить большое число на сотню простых, это очень долго. Поэтому производится следующая оптимизация: составляется таблица остатков remainder[] от деления начального значения тестируемого числа p на нашу сотню простых чисел prime[]. Затем, при поиске простого числа, вместе с самим числом p увеличивается некоторое дельта-значение pdelta. И при поиске остатка от деления p на i-е простое число prime[i] вместо p в качестве делимого берется сумма pdelta+remainder[i]. То есть, используется равенство
p % prime[i] = (pdelta + remainder[i]) % prime[i].
Следует это из уже рассмотренного выше правила:
(a + b) % m = (a + (b % m)) % m
Выигрыш в скорости получается за счет того, что в первом случае мы делили большое число на регистр, а теперь делим регистр на регистр.

Способ второй. Сюда подаются на проверку только числа, прошедшие предыдущую проверку. Используется так называемая теорема Ферма: для любого a, если (a ^ (p-1)) % p != 1, то x - число составное. В качестве a обычно пробуют числа 2,3,5,7. Как видим, здесь просто вызывается описанная раньше процедура modexp().

### Вычисление GCD

Итак, вот мы сгенерили простые числа p и q.

Теперь перебираем e = 3,5,7,9,... пока GCD(e, (p-1)*(q-1)) не окажется равным единице. (можно начинать с любого другого значения числа e, обычно тройка используется из соображений максимальной скорости)

Как посчитать GCD(x,y)? - используя так называемый алгоритм Евклида

     while (1)
     {
       if (y == 0) return x
       x = x % y
       if (x == 0) return y
       y = y % x
     }

Поскольку число x - большое, а y - обычный дворд, оптимизируем алгоритм так: первый раз считаем остаток от деления большого числа на дворд, а затем уже оперируем только двордами. Вычисление остатка от деления приведено в описании процедуры modmult() - 2.

### Вычисление modinv() (modular inverse)

Процедура modinv(a,m) возвращает такое i (0<i<m), что (i*a) mod m = 1. Число a должно быть 0<a<m.

Вот ее алгоритм из исходников PGP 2.6.x: (genprime.c)

     g[0] = m
     g[1] = a
     v[0] = 0
     v[1] = 1
     i = 1
     while (g[i] != 0)
     {
         g[i+1] = g[i-1] % g[i]
         y      = g[i-1] / g[i]
         t      = y * v(i)
         v[i+1] = v[i-1]
         v[i-1] = v[i-1] - t
         i      = i + 1
     }
     x = v[i-1]
     if (x < 0) x = x + m
     return x

А вот алгоритм из каких-то других исходников:

     b = m
     c = a
     i = 0
     j = 1
     while (c != 0)
     {
       x = b / c
       y = b % с
       b = c
       c = y
       y = j
       j = i - j * x
       i = y
     }
     if (i < 0) i += m;
     return i

Если присмотреться, то можно понять, что обе процедуры по сути одинаковые.

Что все это значит, я пока не представляю. Какая-то хуйня. Хотя на ассемблер в конце концов переписалось и заработало ;-)

### Деление больших чисел

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

     d = 0                         -- частное, d = x / y
     c = 0                         -- остаток от деления, c = x % y
     for i = maxbit(x) to 0        -- от старшего бита к младшему
     {
         d = d << 1                -- сдвиг влево на 1 бит, shl
         c = c << 1                -- сдвиг влево на 1 бит, shl
         c = c | x.bit[i]          -- | значит OR, копируем i-й бит x в 0й бит c
         if (c >= y)
         {
             c = c - y             -- % y
             d = d | 1             -- вычли y из c, добавим 1 к d
         }
     }

### Еще что-нибудь про ассемблер

Удобно оперировать битами в больших числах командами BT/BTS/BTR/BTC. Выглядит это так:

     ebx = указатель на большое число в памяти
     ecx = номер бита

     bt  [ebx], ecx -- копировать бит в CF
     bts [ebx], ecx -- копировать бит в CF, затем установить бит (-->1)
     btr [ebx], ecx -- копировать бит в CF, затем очистить бит (-->0)
     btc [ebx], ecx -- копировать бит в CF, затем инвертировать бит

То же самое, но из обычных команд:

     mov edx, ecx
     shr edx, 3
     and cl, 111b
     mov al, 1
     shl al, cl
     test [ebx+edx], al -- jz bit0, jnz bit1
     xor  [ebx+edx], al -- инвертировать бит
     or   [ebx+edx], al -- установить бит (-->1)
     not  al
     and  [ebx+edx], al -- очистить (-->0)

### сдвиг большого числа влево (x shl 1, x * 2)

     ; input: EDI=bignumber
     ;        ECX=length in DWORDs
     ; output:CF
     ; action:x = x shl 1

     bn_shl1:
     bn_sub:       clc                     ; CF=0
     __cycle:      rcl     dword ptr [edi], 1
                   lea     edi, [edi+4]
                   loop    __cycle
                   retn

### Как происходит реальная шифровка

В случае, когда длина сообщения подлежащего шифровке больше длины RSA-ключа (длины числа m, в битах), сообщение разбивается на соответствующие блоки: для каждого блока (числа x) должно соблюдаться условие x < m. Можно например посчитать длину числа m в байтах и длину блоков принимать на байт меньше, и т.п.

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

Можно, например, при помощи RSA шифровать некое начальное значение хэш-буфера, а затем, последовательно изменяя этот буфер, шифровать им собственно сообщение.

Перед шифровкой алгоритмом RSA данные рекомендуется зашифровывать с использованием случайного элемента, по причине low-exponent attack.

### low-exponent attack

Пусть некоторое сообщение x посылается e получателям, так, что открытые ключи каждого из них {e,n1}, {e,n2}, {e,n3}, ...
Хак заключается в восстановлении исходного сообщения, зная e открытых ключей и e зашифрованных сообщений.

Пусть e = 3, x -- исходное сообщение, y1,y2,y3 -- зашифрованные сообщения, n1,n2,n3 -- модули открытых ключей. Тогда x^3 может быть найдено следующим образом:

    #define CRT(c1,c2,n1,n2)  (c1+n1*((modinv(n1,n2)*(c2-c1))%n2))

      x3 = CRT(CRT(y1,y2,n1,n2),y3,n1*n2,n3);

Вычисление кубического корня из большого числа осуществляется простейшим бинарным поиском:

      x = 1;
      for (t = x3; t > 1; t = t / 2)
      {
        if (x*x*x < x3) x += t; else
        if (x*x*x > x3) x -= t; else break;
      }

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

### Еще один способ хакнуть RSA

Вот некий способ похакать RSA, предложен толи каким-то Pinoy'ем, то-ли неким Leo de Velez'ом. Способ насколько теоретически прост, настолько же практически не реализуем.

Есть формула: 2^X = 1 (mod m), где m - модуль, присутствующий в обоих ключах, а X необходимо найти. Ясно, что X = e * d - 1, откуда, найдя X и зная открытое e, можно вычислить секретное d по формуле d = (X + 1) / e.

Идея тут в том, чтобы умножить m на некое Y, так, чтобы произведение m * Y состояло из одних двоичных единиц, то бишь и было равно 2^X-1. Отсюда несложно найти X, равное длине полученного произведения в битах.

Однако, способ нахождения Y, коий предложил de Velez, неверен. То есть эффективного способа попросту нет.

А даже если бы он был, то все равно вычислить ни m * Y, ни Y, ни X (для реальных ключей) было бы невозможно, так как _длина_ в битах числа X того же порядка, что и _значение_ e*d, а поскольку число d - большое, то длина X была бы например 2^1024, а значение 2^(2^1024), что есть полнейший бред.

### В помощь изучающим исходники на английском

Рекомендуются исходники от PGP 2.6.3ia -- там все понятно и с комментариями.

dividend | -- делимое
divisor | -- делитель
quotient | -- частное
remainder | -- остаток
multiplicand | -- множимое
multiplier | -- множитель
product | -- произведение
to raise to a power | -- возводить в степень
base | -- основание (степени)
exponent | -- показатель (степени)
module, modulus | -- модуль. здесь - то же что и remainder
prime, prime number | -- простое число
sieve | -- решето (мат. метод нахождения простых чисел)
Fermat's theorem | -- теорема Ферма: (a^(p-1))%p=1, если p-простое
Euclid's algorithm | -- алгоритм Евклида
factoring | -- факторизация, разложение на множители. в частности, разложение m на p и q.
factoring problem | -- вопрос об эффективном алгоритме факторизации больших чисел, не решен до сих пор. кто решит -- поимеет лимон баксов.
LCM, Least(Lowest) Common Multiple | -- НОК, Наименьшее Общее Кратное
GCD, Greatest Common Divisor