Введение
Для решения криптографических задач необходимо уметь решать квадратичные сравнения по заданному модулю. Алгоритм решения квадратичного сравнения достаточно прост и не вызывает сложностей в решении при небольших значениях модуля и свободного члена, однако в связи с применением достаточно больших чисел в криптографии, решение квадратичных сравнений вручную является весьма кропотливым и длительным процессом. Конечно, для решения квадратичных сравнений можно воспользоваться онлайн-сервисом. Но так как решение криптографической задачи не заканчивается на решении квадратичного сравнения, то человеку, занимающемуся криптографией, будет удобно иметь функцию, способную решать квадратичные сравнения и свободно взаимодействовать с другими функциями, которые используются ним. Именно поэтому было решено написать функцию для решения квадратичных сравнений вида x^2 ≡ a ( mod p ), где a и p — взаимно простые числа, в MATLAB.
Сразу замечу, так написание функций для решения квадратичных сравнений носило обучающий характер, некоторые пользовательские функции, используемые при вычислении, дублируют функции, уже имеющиеся в среде MATLAB.
Сначала предлагаю рассмотреть код двух главных функций, одна из которых предназначена для решения квадратичных сравнений по составному модулю, а вторая – для решения квадратичных сравнений по простому модулю. При этом мы ознакомимся, в первую очередь с алгоритмом решения квадратичных сравнений, а уже во вторую очередь с функциями, которые необходимы для выполнения самих вычислений.
Функция для решения сравнений по сложному модулю
Данная функция позволяет решать квадратичные сравнения, как по простому, так и по сложному модулю. При вызове функции в нее необходимо передать две переменные а и р, в результате своего выполнения функция вернет вектор – строку из двух противоположных по знаку решений квадратичного сравнения.
function [ result ] = sqcomdif( a, p )
Следующим шагом в решении квадратичного сравнения будет определение типа модуля р. Для этого будет использоваться пользовательская функция factorization, предназначенная для разложения числа на простые множители. Кроме вектора-строки с простыми множителями функция возвращает количество простых множителей. По сути функция factorization дублирует стандартную функцию MATLAB factor.
[ mp, sp ] = factorization( p );
После того как модуль был разложен на множители с помощью условного оператора осуществляется проверка количества множителей. Если количество множителей превышает 1, то есть модуль р – составное число, то тогда нам необходимо решить систему из sp квадратичных сравнений ( в каждом из сравнений модулем будет служить один из множителей составного модуля р ). Перед тем как браться за решение полученной системы квадратичных сравнений, необходимо убедиться в том, что каждое из квадратичных сравнений этой системы имеет решение. Для этого воспользуемся циклом for, который будет осуществлять переход между элементами вектора с множителями mp. В теле цикла будет вызываться функция, выполняющая вычисление значения символа Лежандра для каждой пары чисел.
for i=1:1:sp
SL( 1, i ) = symvol( a, mp( 1, i ) );
Если значения символа Лежандра будет равно 1, то переменная count будет увеличена на 1. Это нужно для того, чтобы после выполнения всех итераций цикла мы могли проверить, все ли уравнения системы имеют решение, ведь от этого зависит, будем ли мы решать квадратичные сравнения из которых состоит система, или выведем сообщения о том, что исходное квадратичное сравнение не имеет решений.
if SL( 1, i ) == 1
count = count + 1; % Если есть решение квадратичного сравнения count увеличиваем на 1
end
Если количество уравнений системы равняется количеству уравнений, имеющих решения, то создаются векторы-строки для хранения промежуточных результатов вычисления.
if count == sp
% Формируем векторы для хранения результатов рассчета
answer1 = zeros ( 1, sp ); % Хранение решения квадратичного сравнения
modul = zeros ( 1, sp ); % Хранение отношения общего модуля р к его множителю
answer2 = zeros ( 1, sp ); % Хранение решения линейного сравнения для нахождения общего ответа по Кит. теореме
С помощью цикла for, осуществляется переход между множителями модуля р. В вектор-строку answer1 записываются результаты решения квадратичного сравнения по простому модулю, который получен с помощью функции sqcom, в которую были переданы значение переменной а, а также значение i-ого множителя модуля p. В вектор-строку modul записывается частное от деления составного модуля р на i-ый множитель модуля p. В вектор-строку answer2 будет сохраняться результат решения линейного неравенства ( p / p( I ) ) * y ≡ 1 ( p ( i) ), которое будет необходимо для расчета окончательного ответа по формуле полученной из Китайской теоремы.
% Выполнение расчетов для каждого квадратичного сравнения системы
for i=1:1:sp
answer1( 1, i ) = sqcom( a, mp( 1, i ) ) ;
modul( 1, i ) = p / mp( 1, i );
answer2( 1, i ) = lincom ( modul( 1, i ), 1, mp( 1, i ) );
end
После того как цикл закончил свое выполнение, необходимо рассчитать окончательный ответ по формуле: x=( ( p / p( 1)) * b( 1 ) * y( 1 ) + ( ( p / p( 2) ) * b( 2 ) * y( 2 ) + ( ( p / p( i ) ) * b( i ) * y( i ). Для этого воспользуемся поэлементным умножением. В результате будет получена вектор-строка, сумму элементов которой найдем с помощью команды sum. Найдем остаток от деления суммы на составной модуль р – это будет одним из решений квадратичного сравнения по составному модулю. Второй ответ получим, записав первый ответ с противоположным знаком.
% Расчет окончательного ответа
result = zeros ( 1, 2 );
result( 1, 1 ) = mod( sum( ( modul .* answer1 ) .* answer2 ), p );
result( 1, 2 ) = - result( 1, 1 );
В том случае, если изначально модуль р оказался простым числом, то решение квадратичного сравнения будет получено, путем одноразового вызова функции для решения квадратичного сравнения — sqcom. Второе решение будет получено путем взятия первого ответа с противоположным знаком.
else
result( 1, 1 ) = sqcom( a, p );
result( 1, 2 ) = - result( 1, 1 );
end
Ниже привожу код функции sqcomdif полностью.
function [ result ] = sqcomdif( a, p )
% Функция предназначена для решения квадратичных сравнений по заданному модулю
% Даная функция решает квадратичные сравнения как по простому,
% так и по сложному модулю. При решении квадратичного сравнения по сложному модулю,
% общее решение находится по формуле, полученной из Китайской теоремы.
[ mp, sp ] = factorization( p );
if sp > 1 % Если модуль является составным числом
count = 0;
% Выполняется проверка на возможность решения квадратичного сравнения для
% каждого простого множителя модуля
for i=1:1:sp
SL( 1, i ) = symvol( a, mp( 1, i ) );
if SL( 1, i ) == 1
count = count + 1; % Если есть решение квадратичного сравнения count увеличиваем на 1
end
end
% Если все квадратичные сравнения системы имеют решения
if count == sp
% Формируем векторы для хранения результатов рассчета
answer1 = zeros ( 1, sp ); % Хранение решения квадратичного сравнения
modul = zeros ( 1, sp ); % Хранение отношения общего модуля р к его множителю
answer2 = zeros ( 1, sp ); % Хранение решения линейного сравнения для нахождения общего ответа по Кит. теореме
% Выполнение рассчетов для каждого квадратичного сравнения системы
for i=1:1:sp
answer1( 1, i ) = sqcom( a, mp( 1, i ) ) ;
modul( 1, i ) = p / mp( 1, i );
answer2( 1, i ) = lincom ( modul( 1, i ), 1, mp( 1, i ) );
end
% Расчет окончательного ответа
result = zeros ( 1, 2 );
result( 1, 1 ) = mod( sum( ( modul .* answer1 ) .* answer2 ), p );
result( 1, 2 ) = - result( 1, 1 );
else
result = 'net resheniy';
end
else
result( 1, 1 ) = sqcom( a, p );
result( 1, 2 ) = - result( 1, 1 );
end
end
Функция для решения квадратичного сравнения по простому модулю
Эта функция неоднократно вызывалась в функции sqcomdif, и как уже говорилось, функция sqcom применяется для решения квадратичного сравнения по простому модулю и может вызываться независимо от функции sqcomdif, то есть ее можно без проблем вызвать и получить верный ответ в том случае, когда вы уверенны, что модуль – простое число. Так как рассматривается лишь лишение квадратичных сравнений вида x^2 ≡ a ( mod p ), то в функцию необходимо передать числовые значения переменных а и р. В результате работы функции будет получено одно решение квадратичного сравнения.
function [ answer ] = sqcom( a, p )
Так как функция sqcom может использоваться отдельно от функции sqcomdif, то необходимо убедиться в том, что число, записанное в переменную а, является квадратичным вычетом модуля р. Для этого воспользуемся функцией symvol, которая позволяет вычислить значение символа Лежандра для указанной пары чисел.
[ Symvol_Lejandra ] = symvol( a, p );
Если значение переменной Symvol_Lejandra равно 1, то а является квадратичным вычетом по модулю р и будут выполнены дальнейшие действия для нахождения решения квадратичного сравнения. Нам нужно записать число ( р – 1 ) в виде 2^r * q. Начальные значения переменным r, q задаются из расчета, что ( р – 1 ) – нечетное число. Однако если это не так, то они будут изменяться в процессе выполнения цикла ( до тех пор, пока q не станет нечетным числом ).
q = p - 1;
r = 0;
otn = q / 2;
while ( ( q - floor( otn ) * 2 ) == 0 )
q = otn;
r = r + 1;
otn = q / 2;
end
Теперь необходимо найти значение переменной b, которая равняется b = a^q ( mod p ). Так как а и q, обычно выражены достаточно большими числами, то возвести число а в степень q обычным способом не удастся, так как в большинстве случаев возникнет переполнение. Поэтому возведение в степень необходимо выполнять по методу квадрирования. Для того чтобы выполнить это необходимо вызвать функцию kvadrirovanie и передать в нее значение основания, показателя степени, а также модуля по которому будут выполняться вычисления.
b = kvadrirovanie( a, q, p );
Для продолжения вычислений нужно найти наименьшее неотрицательное число f, которое будет квадратичным невычетом по модулю р. Для этого переменной f, присваивается значение 1, и с помощью функции symvol, определяется значение символа Лежандра для пары чисел f и р. Если символ Лежандра для 1 и р равен 1, то переменная f, с помощью цикла while, увеличивается, до тех пор, пока не будет достигнуто значение, являющееся квадратичным невычетом по модулю р.
f = 1;
sym_lej = symvol( f, p );
while sym_lej ~= -1
f = f + 1;
sym_lej = symvol( f, p );
end
Теперь значение f, являющееся квадратичным невычетом по модулю p, необходимо возвести в степень q. Для этого также придется прибегнуть к функции kvadrirovanie, переменной k не обходимо присвоить значение 0.
g = kvadrirovanie( f, q, p);
k = 0;
После вышеперечисленных действий необходимо проверить значение переменной b. Если b сравнимо с 1 по модулю p, то необходимо перейти к расчету ответа, в противном случае найти наименьшее неотрицательное число m, такое что b^( 2^m ) ≡ 1 ( mod p ). Когда такое значение m будет найдено, необходимо пересчитать значения переменных k, g, b, после чего переменной r присвоить значение m. Но это еще не все, необходимо проверить, чтобы новое значение переменной b, также было сравнимо с 1 по модулю р. Если это не так, то необходимо вернуться к подбору числа m. Переменная pok нужна для того, чтобы избежать выполнения определенных математических действий дважды.
if b ~= 1
while b ~= 1
m = 0;
b1 = kvadrirovanie( b, 2^m, p );
while mod( b1, p) ~= 1
m = m + 1;
b1 = kvadrirovanie( b, 2^m, p );
end
pok = 2^(r-m);
g1 = kvadrirovanie( g, pok, p );
b = mod( ( b*g1), p );
k = fix(k + pok);
r = m;
end
end
После того как было найдено m, удовлетворяющее перечисленным выше условиям, можно переходить к непосредственному вычислению ответа. Ответ вычисляется по формуле x = a^( ( q+1) / 2 ) * g^( k / 2 ) ( mod p ). Для вычисления обеих множителей воспользуемся функцией квадрирования, полученный результат возьмем по модулю p.
first = kvadrirovanie( a, ( ( q + 1 ) / 2 ), p );
second = kvadrirovanie( g, ( k / 2 ), p );
answer = mod( ( first * second ), p);
Полученный результат не всегда может быть записан оптимально, по модулю р. Поэтому необходимо выполнить следующую проверку:
delta = p - answer;
if delta < answer
answer = delta;
end
Ниже приведен полный код функции sqcom:
function [ answer ] = sqcom( a, p )
% Функция предназначена для решения квадратичного сравнения по простому
% модулю
% Функция решает квадратичное сравнение вида x^2 = a ( mod p ) , возвращает
% одно решение. Это решение взятое с противоположным знаком будет
% вторым решением.
a=mod(a,p);
% ШАГ 1 проверяем есть ли решение данного квадратичного неравенства
[ Symvol_Lejandra ] = symvol( a, p );
if Symvol_Lejandra == 1
% ШАГ 2 Ищем q, r, вычисляем b
q = p - 1;
r = 0;
otn = q / 2;
while ( ( q - floor( otn ) * 2 ) == 0 )
q = otn;
r = r + 1;
otn = q / 2;
end
b = kvadrirovanie( a, q, p );
% ШАГ 3 Ищем f и квадратичный невычет
f = 1;
sym_lej = symvol( f, p );
while sym_lej ~= -1
f = f + 1;
sym_lej = symvol( f, p );
end
g = kvadrirovanie( f, q, p);
k = 0;
% ШАГ 4 ИТЕРАЦИЯ
if b ~= 1
while b ~= 1
m = 0;
b1 = kvadrirovanie( b, 2^m, p );
while mod( b1, p) ~= 1
m = m + 1;
b1 = kvadrirovanie( b, 2^m, p );
end
pok = 2^(r-m);
g1 = kvadrirovanie( g, pok, p );
b = mod( ( b*g1), p );
k = fix(k + pok);
r = m;
end
end
% ШАГ 5 Рассчет окончательного ответа
first = kvadrirovanie( a, ( ( q + 1 ) / 2 ), p );
second = kvadrirovanie( g, ( k / 2 ), p );
answer = mod( ( first * second ), p);
delta = p - answer;
if delta < answer
answer = delta;
end
else
answer = 'net resheniya';
end
А теперь предлагаю ознакомиться с вспомогательными функциями, которые используются как при решении квадратичного сравнения, так могут использоваться и отдельно.
Разложение числа на множители
При решении квадратичного сравнения во многих случаях необходимо прибегать к разложению числа на множители, также к этой операции придется прибегать в тех случаях, когда будет необходимо проверить: является число простым или составным.
В функцию factorization передается число, которое необходимо разложить на множители. В результате функция возвращает вектор — строку с множителями и количество этих множителей.
function [ mnojitel, ind ] = factorization( delimoe )
Функция состоит из оператора switch, который выполняет различные действия в зависимости от значения входной переменной delimoe. Так если delimoe = 1, то в вектор, хранящий результат разложения на множители, mnojitel записывается 1, в переменную ind, хранящую число множителей, также записывается 1. Аналогичные действия выполняются и в том случае, если delimoe равно -1.
switch delimoe
case { 1 }
mnojitel( 1, 1 )=1;
ind=1;
case { -1 }
mnojitel( 1, 1 )= -1;
ind = 1;
Если эти условия не выполнились, то выполняем проверку знака числа, раскладываемого на множители. Если delimoe меньше 0, то в первый множитель записывается -1, в ind записывается 2, и дальнейшая работа функции будет осуществляться с модулем переменной delimoe.Если же число больше 0, то переменной ind присваивается значение 1.
otherwise
if delimoe < 0
mnojitel( 1, 1 )= -1;
ind = 2;
delimoe = abs ( delimoe );
else
ind = 1;
end
Цикл while работает до тех пор, пока delimoe не равняется delitel, причем изначально переменной delitel присвоено значение 2. При каждой итерации цикла, в переменную ostatok записывается остаток от деления delimoe на delitel. Если остаток равен 0, то есть delitel является множителем delimoe, то это значение записывается в вектор, в котором хранятся множители, а счетчик этого вектора увеличивается на 1. При этом переменной delimoe присваивается частное от выполненного деления. Если же остаток от деления не равен 0, то delitel увеличивается на 1. Когда же происходит выход из цикла, то значение, оставшееся в переменной delimoe, записывается в вектор с множителями, как один из множителей.
while ( delimoe ~= delitel )
ostatok = mod( delimoe, delitel );
if ostatok ~= 0
delitel = delitel + 1;
else
delimoe = delimoe / delitel;
mnojitel( 1, ind ) = delitel;
ind = ind + 1;
end
end
mnojitel( 1, ind ) = delimoe;
Ниже приведу полный код функции factorization.
function [ mnojitel, ind ] = factorization( delimoe )
% Функция для разложения числа на множители
% Начальные значения
delitel = 2;
% Разложение на множители
switch delimoe
case { 1 }
mnojitel( 1, 1 )=1;
ind=1;
case { -1 }
mnojitel( 1, 1 )= -1;
ind = 1;
otherwise
if delimoe < 0
mnojitel( 1, 1 )= -1;
ind = 2;
delimoe = abs ( delimoe );
else
ind = 1;
end
while ( delimoe ~= delitel )
ostatok = mod( delimoe, delitel );
if ostatok ~= 0
delitel = delitel + 1;
else
delimoe = delimoe / delitel;
mnojitel( 1, ind ) = delitel;
ind = ind + 1;
end
end
mnojitel( 1, ind ) = delimoe;
end
end
Вычисление значения символа Лежандра
Для того чтобы проверить является ли число квадратичным вычетом по заданному модулю ( в таком случае квадратичное сравнение имеет решение ) или квадратичным невычетом ( тогда такое квадратичное сравнение не имеет решения ) применяется символ Лежандра, который в отечественной литературе обозначается как L ( a; p ), в зарубежной литературе как L (a / p ).
Символ Лежандра может принимать следующие значения:
L ( a; p ) = 1, в таком случае а принадлежит QR, квадратичное сравнение имеет решение
L ( a; p ) = -1, в таком случае а принадлежит QNR, квадратичное сравнение не имеет решения
Если L ( a; p ) = 0, то а и р не являются взаимно простыми числами, то есть НОД ( а; р ) не равен 1
Для вычисления значения символа Лежандра применяются следующие свойства:
- L ( 1; p ) = 1
- L ( -1; p ) = ( -1 ) ^ ( ( p — 1 ) / 2 )
- L ( 2; p ) = ( -1 ) ^ ( ( p^2 — 1) / 8 )
- Если а≡ b*с ( mod p ), то L ( b*с; p ) = L ( b; p ) * L (с; p )
- Если а≡ b ( mod p ), то L ( а; p ) = L ( b; p )
- Если а и р – простые числа, то L ( а; p ) = (-1) ^( (( p — 1 ) *(a-1)) / 4 ) * L( p; a ). Последнее свойство получило название закона взаимности Гаусса.
Вычисление символа Лежандра выполняется на основе вышеперечисленных свойств. Как только одно из условий выполнилось, начинаем проверять свойства для полученной пары а и р, до тех пор, пока не будет найдено окончательное значение символа Лежандра.
Теперь рассмотрим, как можно реализовать вычисление значения символа Лежандра программно.
Функция будет возвращать значение символа Лежандра для пары чисел а и р, которые были переданы в нее. Это видно из заголовка функции:
function [ sl ] = symvol( a, p )
Следующий шаг по своей сути является применением свойства 5. Если число а больше модуля р, то его можно заменить меньшим числом, сравнимым с а по модулю р.
a=mod( a, p ); % Взятие а по модулю р
Полученное число а, пытаемся разложить на множители. Для разложения числа на множители была написана пользовательская функция factorization, возвращающая вектор с простыми множителями, из которых состоит число а, а также их количество. Подробно эта функция была описана выше.
[ mnoj, ind ] = factorization( a ); % Разложение числа а на множители
Если а не являлось простым числом, то будем действовать по свойству 4. То есть значение символа Лежандра для пары чисел L ( а; p ), найдем как произведение значений символов Лежандра для чисел, являющихся простыми множителями числа а. Для хранения промежуточных результатов создадим вектор-строку, заполненную нулями с размерностью, равной количеству множителей числа а. Если а – простое, то эта строка будет состоять из одного элемента.
sl = zeros( 1, ind ); % Создание матрицы для хранения значений символов Лежандра для каждого множителя
Для вычисления символа Лежандра для каждого множителя поочередно воспользуюсь циклом for, который будет изменять свои значения от 1 до номера последнего множителя. В теле данного цикла происходит непосредственный процесс вычисления значения символа Лежандра с помощью свойств, перечисленных выше.
Код для проверки свойства 1 имеет следующий вид:
% Для символа L(1,p)
if mnoj( 1, i ) == 1
sl( 1, i ) = 1;
end
Так как при проверке значения символа при его виде L(-1,p), по свойству 2, необходимо вычислить значение ( -1 ) ^ ( ( p — 1 ) / 2 ), то необходимо воспользоваться еще одним условным оператором, в котором будет проверяться, является показатель -1 четным или нечетным числом. В зависимости от этого будет варьироваться значение символа Лежандра. Если показатель четный – то символ Лежандра будет равен 1, в противном случае -1. Применение этого условного оператора позволяет избежать возведение -1 в степень ( р — 1)/2, что является весьма затратной операцией.
% Для символа L(-1,p)
if mnoj( 1, i ) == -1
if mod( ( ( p - 1 ) ) / 2, 2 ) == 0
sl( 1, i ) = 1;
else
sl( 1, i ) = -1;
end
end
Аналогичные действия выполняются и в том случае, когда необходимо вычислить значение символа Лежандра при его виде L ( 2; p ). В таком случае оно будет равняться ( -1 ) ^ ( ( p^2 — 1) / 8 ).
% Для символа L(2,p)
if mnoj( 1, i ) == 2
% Если показатель четный то символ Лежанра будет равен 1, иначе -1
if mod( ( ( p^2 - 1 ) / 8 ), 2 ) == 0
sl(1,i)=1;
else
sl(1,i)=-1;
end
end
После проверки этого условия число а, передается в функцию, для того, чтобы проверить, простое оно или нет. В случае если число а – составное ( число его множителей ind1 больше 1 ), то происходит рекурсия и число а передается в эту же функцию для выполнения дальнейших вычислений.
[ mn, ind1 ] = factorization( mnoj( 1, i ) ); % Расскладываем число а на множители
if ind1 > 1 % Если число а - составное
sl( 1, i ) = symvol( mnoj(1,i), p ); % то передаем его в эту же функцию, для дальнейших вычислений
В противном случае число а – простое. Если при этом оно не равно -1, 1, 2, то приходится воспользоваться свойством 6 – законом взаимности Гаусса. Знак перед символом Лежандра определяется путем определения четности множителей его показателя. Он обращается в плюс, если хотя бы один из множителей – четный. После чего происходит рекурсивный вызов функции symvol, а аргументы в нее передаются в другом порядке.
elseif and( mnoj(1,i)~=-1, and( mnoj( 1, i ) ~= 1, mnoj( 1, i ) ~= 2 ) )% если число а - простое, не равно 1, -1, 2
if or( mod( ( ( p - 1 ) / 2 ), 2 ) == 0, mod( ( ( mnoj( 1, i ) - 1 ) / 2 ), 2 ) == 0 ) % если хотя бы один из множителей показателя - четный
sl(1,i)= symvol( p, mnoj( 1, i ) ); % L(p,a)
else
sl(1,i)=-symvol( p, mnoj( 1, i ) ); % -L(p,a)
end
end
В результате проверки вышеперечисленных условий были охвачены все возможные варианты значений числа а.
Если изначально число а было составным, то в переменную sl записывается произведение всех элементов строки — вектора sl, после чего полученное значение возвращается в программу, вызвавшую данную функцию.
if ind ~= 1
sl = prod( sl ); % Вычисление результата L(a,p)
end
Полный код функции symvol:
function [ sl ] = symvol( a, p )
% Функция для вычисления символа Лежандра L(a,p)
% Функция позволяет вычислить символ Лежандра, что дает возможность
% проверить, имеет ли решение квадратичное сравнение по простому модулю
a=mod( a, p ); % Взятие а по модулю р
[ mnoj, ind ] = factorization( a ); % Разложение числа а на множители
sl = zeros( 1, ind ); % Создание матрицы для хранения значений символов Лежандра для каждого множителя
% Цикл перебирающий множители вектора а
for i = 1:ind
% Для символа L(1,p)
if mnoj( 1, i ) == 1
sl( 1, i ) = 1;
end
% Для символа L(-1,p)
if mnoj( 1, i ) == -1
if mod( ( ( p - 1 ) ) / 2, 2 ) == 0
sl( 1, i ) = 1;
else
sl( 1, i ) = -1;
end
end
% Для символа L(2,p)
if mnoj( 1, i ) == 2
% Если показатель четный то символ Лежанра будет равен 1, иначе -1
if mod( ( ( p^2 - 1 ) / 8 ), 2 ) == 0
sl(1,i)=1;
else
sl(1,i)=-1;
end
end
[ mn, ind1 ] = factorization( mnoj( 1, i ) ); % Расскладываем число а на множители,
% чтобы проверить простое число или составное
if ind1 > 1 % Если число а - составное
sl( 1, i ) = symvol( mnoj(1,i), p ); % то передаем его в эту же функцию, для дальнейших вычислений
elseif and( mnoj(1,i)~=-1, and( mnoj( 1, i ) ~= 1, mnoj( 1, i ) ~= 2 ) )% если число а - простое, не равно 1 и не равно 2
if or( mod( ( ( p - 1 ) / 2 ), 2 ) == 0, mod( ( ( mnoj( 1, i ) - 1 ) / 2 ), 2 ) == 0 ) % если хотя бы один из множителей показателя - четный
sl(1,i)= symvol( p, mnoj( 1, i ) ); % L(p,a)
else
sl(1,i)=-symvol( p, mnoj( 1, i ) ); % -L(p,a)
end
end
end
if ind ~= 1
sl = prod( sl ); % Вычисление результата L(a,p)
end
end
Возведение числа в степень по методу квадрирования
При решении квадратичного сравнения часто приходится сталкиваться с необходимостью возведения больших чисел в степень с показателем равным нескольким сотням, а то и тысячам. При обычном возведении числа в такую степень достаточно легко натолкнуться на такое неприятное явление, как переполнение. Для того чтобы избежать его, необходимо возводить такие числа в степень используя метод квадрирования.
Рассмотрим алгоритм данного метода:
- Показатель степени записывается в двоичной системе счисления.
- Если он состоит из двух или более бит, под первой 1 записываем число m=a.
- Смотрим на следующий справа бит:
- Если это единица то возводим m в квадрат и берем его по модулю р, а полученный результат домножаем на а и берем его по модулю р, окончательный результат записывается в переменную m.
- Если это нуль, то возводим число m в квадрат, берем результат по модулю р, окончательный результат записывается в переменную m.
- Переходим к следующему справа биту и повторяем действия описанные в пункте 3, пока не будут выполнены расчеты для последнего бита.
А теперь рассмотрим, как работает функция kvadrirovanie. В функцию передаются значения основания, показателя и модуля, а возвращает функция одно число – конечный результат.
function [ result ] = kvadrirovanie( a, q, p )
Показатель степени q переведем в двоичный вид, и так как полученный результат будет состоять из отдельных символов, то воспользуемся командой size, для определения количества бит, необходимых для записи показателя степени.
q = dec2bin( q );
size_q = size(q);
Если для записи показателя степени необходимо два или более битов, то программа выполнит следующие действия: запишет в переменную m типа uint64 число а. После чего с помощью цикла for, который будет изменять значения переменной i от 2 с шагом 1 до количества бит, необходимых для записи показателя степени q, для того, чтобы иметь возможность осуществлять перебор битов.
if size_q( 1, 2 ) >= 2
m = uint64(a);
for i=2:1:size_q(1,2)
В теле цикла имеется условный оператор, который, в зависимости от значения i-го бита, будет определять действия, которые необходимо произвести с переменной m. Как писалось выше, если бит равен 1, то необходимо возвести m^2 и взять по модулю р, после чего умножить полученное число на а, и результат взять по модулю р.
if q(1,i)=='1'
m = uint64( mod( ( mod( ( m^2 ), p ) * a ), p ));
Иначе, q(1,i)=='0', необходимо возвести m^2 и взять его по модулю р.
else
m = uint64(mod( ( m^2 ), p ));
end
После выполнения всех повторов цикла, окончательное значение переменной m записывается в переменную result.
result =uint64(m);
В показатель степени в двоичной форме равен 1, то в переменную result записывается само же исходное число, которое требовалось возвести в степень.
elseif q(1,1) == '1'
result = uint64( a );
Если же и это условие не выполнилось, то следовательно показатель степени равняется 0. В таком случае в переменную result будет записана 1.
else
result = 1;
end
Полный код функции для возведения в степень по методу квадрирования:
function [ result ] = kvadrirovanie( a, q, p )
% Функция для возведение числа в степень по указанному модулю
% Показатель переводится в двоичный вид, если бит равен 1 то число
% возводится в квадрат и домножается само на себя, если ноль то возводится
% в квадрат. Полученный результат берется по модулю.
q = dec2bin( q );
size_q = size(q);
if size_q( 1, 2 ) >= 2
m = uint64(a);
for i=2:1:size_q(1,2)
if q(1,i)=='1'
m = uint64( mod( ( mod( ( m^2 ), p ) * a ), p ));
else
m = uint64(mod( ( m^2 ), p ));
end
end
result =uint64(m);
elseif q(1,1) == '1'
result = uint64( a);
else
result = 1;
end
end
Решение линейных сравнений
При решении квадратичных сравнений по сложному модулю будет применяться формула, полученная из Китайской теоремы. Для того чтобы рассчитать ответ по этой формуле появляется необходимость решить линейное сравнение k * x ≡ b ( mod p ). Данная функция предназначается для их решения с помощью обратного элемента, то есть такого числа k1 которое при умножении на k, даст результат сравнимый с 1 по заданному модулю р.
В функцию lincom передаются значения коэффициента при неизвестном k, число b, с которым сравнима левая часть, а также модуль p, а возвращает функция соответственно решение линейного сравнения с заданными параметрами.
function [ x ] = lincom ( k, b, p)
Как указывалось выше, неравенство будем решать с помощью обратного элемента. Для того чтобы найти обратный элемент воспользуемся расширенным алгоритмом Эвклида для нахождения НОДа. Напомню, что расширенный алгоритм Эвклида позволяет найти НОД с помощью коэффициентов, которые рассчитываются при выполнении каждой операции деления. Для нахождения НОДа по данному алгоритму необходима рассчитать два ( a, b ) коэффициента, но при нахождении обратного элемента нам потребуется только один – b коэффициент, который как раз и будет обратным элементом для k.
В функции необходимо указать начальные значения коэффициентов b0 и b1, которые будут изменяться по мере выполнения расчетов. Если модуль p, который был записан в переменную pr, делится нацело на коэффициент k, то обратный элемент будет равен 1. Иначе b1 будет рассчитываться с помощью цикла while. В цикле будет пересчитываться значения переменной b0. После того как найдено новое значение b0, с помощью функции swap будет осуществляться обмен значений между переменными b0 и b1. Также при каждой итерации цикла будет происходить присвоение новых значений ряду других переменных.
b0=0;
b1=1;
%Нахождение остатка
ostatok = mod( pr, k );
while ostatok ~= 0
chastnoe = floor( pr / k );
b0 = b0 - b1 * chastnoe;
[ b0, b1 ] = swap( b0, b1 );
pr = k;
k = ostatok;
ostatok = mod( pr, k );
end
После того, как цикл прекратит свое выполнение будет необходимо рассчитать значение решения линейного сравнения. Оно будет равно произведению переменных b1 и b, взятое по модулю p ( именно для этой операции в начале функции значение модуля было продублировано с помощью переменной pr ).
x = mod( b1*b, p );
Ниже приведен код функции целиком:
function [ x ] = lincom ( k, b, p)
% Функция предназначена для решения линейных неравенств вида k*x=b ( mod p )
pr=p;
b0=0;
b1=1;
%Нахождение остатка
ostatok = mod( pr, k );
while ostatok~=0
chastnoe = floor( pr / k );
b0 = b0 - b1 * chastnoe;
[ b0, b1 ] = swap( b0, b1 );
pr = k;
k = ostatok;
ostatok = mod( pr, k );
end
x = mod( b1*b, p );
end
Список источников, которые помогли при создании функций и написании статьи:
- Конспекты и лекции Плотниковой Л.И. ( кто в теме — поймет )
- http://math.hashcode.ru
- http://mathhelpplanet.com
- http://www.wolframalpha.com
Автор: Sany_KENT
