Я пишу библиотеку 256-битной арифметики для реализации эллиптической криптографии на микроконтроллерах. Вычисления должны проводиться в конечном после, которое, для оптимизации, будет выбрано заранее, и во все вычисления будут заточены под это поле. Таким образом, длинные числа у нас фиксированной длины: допустим, 256 или 128 бит. Для простоты отладки я сейчас реализую 32-битную арифметику. Для тестов мне понадобилось переводить данное число в систему счисления с произвольным основанием. Я для этого написал скрипт на хаскеле: Код (Text): xpand :: Integer -> Integer -> [Integer] expand num base | num < base = [num] | otherwise = num `mod` base : expand (num `div` base) base expand 1000000 64 переведет число в массив цифр по основанию 64. Ближе к реализации. Пусть у нас "длинные числа" имеют размерность 32 бита. Разделим 32 бита на 5 частей: каждый кусок числа содержит 6-битную цифру. Вместе получается 30 бит. Вполне достаточно. Это называется редуцированной записью числа и позволяет не использовать инструкцию adc при суммировании длинных чисел. Дело в том, что для эффективной реализации перенос из разрядов только вредит, т.к. инструкция adc в 6 раз медленее инструкции add. Кроме того, вычисления могут производиться на векторном процессоре, тогда 5 разрядов будут суммироваться за такт, при этом carry теряются, т.к. они бы замедлили вычисления. Поэтому используем редуцированную форму записи (уменьшенный диапазон цифр). При переполнении разрядов переносы нужно посчитать вручную (уже после того, как закончилоссь поразрядное сложение за один такт): Код (C): typedef struct { signed char a[5]; } bigint32; void make_bigint32( bigint32* x, signed char a0, signed char a1, signed char a2, signed char a3, signed char a4) { x->a[0] = a0; x->a[1] = a1; x->a[2] = a2; x->a[3] = a3; x->a[4] = a4; } void bigint32_add(bigint32* r, const bigint32* x, const bigint32* y) { // Addition without carry r->a[0] = x->a[0] + y->a[0]; r->a[1] = x->a[1] + y->a[1]; r->a[2] = x->a[2] + y->a[2]; r->a[3] = x->a[3] + y->a[3]; r->a[4] = x->a[4] + y->a[4]; // Carry propagation r->a[1] += r->a[0] >> 6; r->a[0] &= 0x3F; r->a[2] += r->a[1] >> 6; r->a[1] &= 0x3F; r->a[3] += r->a[2] >> 6; r->a[2] &= 0x3F; r->a[4] += r->a[3] >> 6; r->a[3] &= 0x3F; r->a[4] &= 0x3F; } } Имеем 6-битовые цифры числа. При переполнении carry оказывается в старших битах. Чтобы получить carry, сдвигаем число вправо на 6 бит. Переносим этот carry в старший разряд суммированием. Накладывая маску 0x3F, обнуляем все старшие разряды кроме шести младших. Так работает алгоритм суммирования. У меня не получается правильно реализовать умножение по модулю. Кстати, выберем конечное поле и модуль: пусть [math]p = 2^{12}- 5[/math], [math]F_p[/math] - поле, состоящее из остатков по модулю [math]2^{12} - 5[/math]. Функция умножения: Код (C): #include <stdio.h> #include <stdint.h> typedef struct { signed char a[5]; } big; void make_big(big* x, signed char a0, signed char a1, signed char a2) { x->a[0] = a0; x->a[1] = a1; x->a[2] = a2; } void print(big* x) { for (int i = 0; i < 4; ++i) { printf("%02d", x->a[i]); printf(" "); } printf("\n"); } void int32_mul(big* r, const big* x, const big* y) { const signed char* a = x->a; const signed char* b = y->a; unsigned short t[5] = { a[0]*b[0], a[0]*b[1] + a[1]*b[0], a[0]*b[2] + a[1]*b[1] + a[2]*b[0], a[1]*b[2] + a[2]*b[1], a[2]*b[2] }; t[1] += t[0] >> 6; t[0] &= 0x3F; t[2] += t[1] >> 6; t[1] &= 0x3F; t[3] += t[2] >> 6; t[2] &= 0x3F; t[4] += t[3] >> 6; t[3] &= 0x3F; t[0] += t[2] * 5; // 2^0 t[1] += t[3] * 5; // 2^6 t[2] = 0; t[3] = 0; t[1] += t[0] >> 6; t[0] &= 0x3F; t[0] += t[1] >> 6; t[1] &= 0x3F; r->a[0] = t[0]; r->a[1] = t[1]; r->a[2] = t[2]; r->a[3] = t[3]; } void int32_pow(big* r, const big* base, int exp) { big tmp; make_big(&tmp, 1, 0, 0); while (exp--) { int32_mul(&tmp, &tmp, base); } *r = tmp; } int main() { big r, x, y; make_big(&x, 40, 15, 0); make_big(&y, 16, 31, 0); int32_mul(&r, &x, &y); print(&r); make_big(&x, 39, 15, 0); make_big(&y, 9, 12, 0); int32_mul(&r, &x, &y); print(&r); make_big(&x, 23, 17, 0); make_big(&y, 43, 8, 0); int32_mul(&r, &x, &y); print(&r); make_big(&x, 43, 8, 0); int32_pow(&r, &x, 5); print(&r); make_big(&x, 27, 63, 0); make_big(&y, 43, 8, 0); int32_mul(&r, &x, &y); print(&r); } Объясняю принцип работы по частям: Код (C): unsigned short t[5] = { a[0]*b[0], a[0]*b[1] + a[1]*b[0], a[0]*b[2] + a[1]*b[1] + a[2]*b[0], a[1]*b[2] + a[2]*b[1], a[2]*b[2] }; Это массив коэффициентов произведения без учета переносов и редукции цифр. Пусть длинным числам соответствуют векторы их коэффициентов. Построим многочлены от [math]X[/math] с этими коэффициентами и перемножим. Коэффициента многочлена результата вычисляются по тем формулам (можно расписать на бумажке). Если вместо [math]X[/math] подставить 6, получим результат - число, равное произведению двух длинных чисел. Но в результате цифры получились слишком большими, а как мы помним, под цифры мы отвели всего 6 бит. Для нормализации числа нам помогут carry: Код (C): t[1] += t[0] >> 6; t[0] &= 0x3F; t[2] += t[1] >> 6; t[1] &= 0x3F; t[3] += t[2] >> 6; t[2] &= 0x3F; t[4] += t[3] >> 6; t[3] &= 0x3F; Как и в коде сложения, сдвигом получаем carry (то, что переполнило цифру) и переносим в старший разряд, а в текущем разряде обнуляем все, что за пределами младших 6 бит. Код (C): t[0] += t[2] * 5; // 2^0 t[1] += t[3] * 5; // 2^6 t[2] = 0; t[3] = 0; Понять этот код без модулярной арифметики в конечном поле не получится. Как мы помним, мы зафиксировали конечное поле, состоящее из остатков от деления на [math]2^{12} - 5[/math] (это простое число, поэтому данная алгебраическая структура - поле). Результат умножения можно записать как: [math]a_0 + a_1\cdot 2^6 + a_2\cdot 2^{12} + a_3\cdot 2^{18} + a_4\cdot 2^{24}[/math].Избавимся от больших степеней. Можно заметить, что поскольку [math]2^{12} - 5 = 0\; mod (2^{12} - 5)[/math], то [math]2^{12} = 5\; mod (2^{12} - 5)[/math]. Кроме того, [math]2^{18} = 2^{12}\cdot 2^6 = 5\cdot 2^6[/math], а [math]2^{24} = 2^{12}\cdot 2^{12} = 5 \cdot 5 = 25[/math]. Таким образом, число запишется как [math]a_0 + a_1\cdot 2^6 + 5\cdot a_2 + 5\cdot 2^6\cdot a_3 + 5\cdot 5\cdot a_4 = a_0 + 5\cdot a_2 + (a_1 + 5)\cdot 2^6 + 25\cdot a_4[/math]. В последнем коде я просто сделал группировку слагаемых (в выражении [math]t_0=t_0+ 5\cdot t_2[/math]; [math]t_0[/math] и [math]t_2[/math] - это коэффициенты при шестых степенях, а в выражении [math]t_1=t_1+ 5\cdot t_3[/math]; [math]t_3[/math] это коэффициент при 18-й степени, которая редуцировалась до [math]5\cdot 2^6[/math], поэтому коэффициент [math]t_3[/math] мы умножаем на 5 и прибавляем к коэффициенту при шестой степени, то есть к [math]t_1[/math]. Далее Код (C): t[1] += t[0] >> 6; t[0] &= 0x3F; t[0] += t[1] >> 6; t[1] &= 0x3F; Учитываем возможность переполнения разрядов и нормализуем запись числа при помощи переносов, чтобы все цифры стали 6-битными. Вопрос: почему этот алгоритм работает на первых трех умножениях при возведении в степень, но не работает на четвертом и далее?
Код (Text): void int32_pow(big* r, const big* base, int exp) { big tmp; make_big(&tmp, 1, 0, 0); while (exp--) { int32_mul(&tmp, &tmp, base); } *r = tmp; } Ошибка вычисления возникает и в том случае, когда я умножаю число на его предыдущую степень (вроде, на третью оно вчера умножслось уже с ошибкой).
в форме для ответа есть буква f , которая обрамляет текст тэгами [МАТН] и [/МАТН], степень a^{b}=[math]a^{b}[/math], индекс a_{b}=[math]a_{b}[/math], справка по редактору формул здесь
Что-то не похоже, что умножения вообще работают. Код (C): make_big(&x, 40, 15, 0); make_big(&y, 16, 31, 0); int32_mul(&r, &x, &y); print(&x); print(&y); print(&r); printf("\n"); Код (Text): 40 15 00 00 16 31 00 00 08 56 00 00 Должно было вывести 13 15 31 15, разве нет?
А, да. Не заметил. В таком случае результат должен быть 00 16. А не 08 56, не так ли? Кроме того, 56 - это 7 бит, а не шесть.
