Опції для вирішення квадратичних порівнянь. Реалізація в MATLAB

Введення

Для вирішення криптографічних завдань необхідно вміти вирішувати квадратичні порівняння по заданому модулю. Алгоритм рішення квадратичного порівняння досить простий і не викликає труднощів у вирішенні при невеликих значеннях модуля і вільного члена, однак у зв'язку з застосуванням досить великих чисел в криптографії, рішення квадратичних порівнянь вручну є досить копітким і тривалим процесом. Звичайно, для вирішення квадратичних порівнянь можна скористатися онлайн-сервісом. Але так як рішення криптографічного завдання не закінчується на вирішенні квадратичного порівняння, то людині, що займається криптографією, буде зручно мати функцію, здатну вирішувати квадратичні порівняння і вільно взаємодіяти з іншими функціями, які використовуються ним. Саме тому було вирішено написати функцію для вирішення квадратичних порівнянь виду x^2 ≡ a ( mod p ) в MATLAB.



Відразу зауважу, так написання функцій для розв'язування квадратичних порівнянь носила повчальний характер, деякі користувальницькі функції, використовувані при обчисленні, дублюють функції, наявні в середовищі MATLAB.

Спочатку пропоную розглянути код двох головних функцій, одна з яких призначена для вирішення квадратичних порівнянь по складеному модулю, а друга — для вирішення квадратичних порівнянь по простому модулю. При цьому ми ознайомимося, в першу чергу з алгоритмом рішення квадратичних порівнянь, а вже в другу чергу з функціями, які необхідні для виконання самих обчислень.

Функція для вирішення порівнянь по складному модулю

Дана функція дозволяє вирішувати квадратичні порівняння, як щодо простого, так і складного модулю. При виклику функції в неї необхідно передати дві змінні р, в результаті виконання свого функція поверне вектор — рядок з двох протилежних за знаком рішень квадратичного порівняння.

function [ result ] = sqcomdif( a, p )

Наступним кроком у вирішенні квадратичного порівняння буде визначення типу модуля р. Для цього буде використовуватися користувацька функція факторизації, призначена для розкладання числа на прості множники. Крім вектора-рядка з простими множниками функція повертає кількість простих множників. По суті функція факторизації дублює стандартну функцію MATLAB factor.

[ mp, sp ] = факторизації( 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 ] = факторизації( 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

А тепер пропоную ознайомитися з допоміжними функціями, які використовуються як при вирішенні квадратичного порівняння, так можуть використовуватися і окремо.



Розкладання числа на множники

При вирішенні квадратичного порівняння у багатьох випадках необхідно вдаватися до розкладання числа на множники, також до цієї операції доведеться вдаватися в тих випадках, коли необхідно перевірити: чи є число простим або складеним.
У функцію факторизації передається число, яке необхідно розкласти на множники. У результаті функція повертає вектор — рядок з множниками і кількість цих множників.

function [ mnojitel, ind ] = факторизації( 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;

Нижче наведу повний код функції факторизації.

function [ mnojitel, ind ] = факторизації( 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
Для обчислення значення символу Лежандра застосовуються наступні властивості:

  1. L ( 1; p ) = 1
  2. L ( -1; p ) = ( -1 ) ^ ( ( p — 1 ) / 2 )
  3. L ( 2; p ) = ( -1 ) ^ ( ( p^2 — 1) / 8 )
  4. Якщо а≡ b*с ( mod p ), L ( b*с; p ) = L ( b; p ) * L (с; p )
  5. Якщо а≡ b ( mod p ), L ( а; p ) = L ( b; p )
  6. Якщо р — прості числа, то L ( а; p ) = (-1) ^( (( p — 1 ) *(a-1)) / 4 ) * L( p, a ). Остання властивість отримало назву закону взаємності Гаусса.


Обчислення символу Лежандра виконується на основі перерахованих вище властивостей. Як тільки одна з умов здійснилось, починаємо перевіряти властивості отриманої пари р, до тих пір, поки не буде знайдене остаточне значення символу Лежандра.
Тепер розглянемо, як можна реалізувати обчислення значення символу Лежандра програмно.
Функція повертає значення символа Лежандра для пари чисел р, які були передані до неї. Це видно із заголовка функції:

function [ sl ] = symvol( a, p )

Наступний крок по своїй суті є застосуванням властивості 5. Якщо число модуля р, то його можна замінити меншим числом, порівнянним з по модулю р.

a=mod( a, p ); % Взяття а по модулю р

Отримане число , намагаємося розкласти на множники. Для розкладання числа на множники була написана користувацька функція факторизації, яка повертає вектор з простими множниками, з яких складається число , а також їх кількість. Детально ця функція була описана вище.

[ mnoj, ind ] = факторизації( 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 ] = факторизації( 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 ] = факторизації( 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 ] = факторизації( 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. Показник ступеня записується в двійковій системі числення.
  2. Якщо він складається з двох або більше біт, під першою 1 записуємо число m=a.
  3. Дивимося на наступний праворуч біт:
    • Якщо це одиниця то зводимо m в квадрат і беремо його по модулю р, а отриманий результат домножаем на і беремо його по модулю р, остаточний результат записується в змінну m.
    • Якщо це нуль, то зводимо число m в квадрат, беремо результат по модулю р, остаточний результат записується в змінну m.
  4. Переходимо до наступного праворуч биту і повторюємо дії описані в пункті 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


Перелік джерел, які допомогли при створенні функцій і написанні статті:

  1. Конспекти та лекції Плотникової Л. В. ( хто в темі — зрозуміє )
  2. http://math.hashcode.ru
  3. http://mathhelpplanet.com
  4. <a href=«www.wolframalpha.com> www.wolframalpha.com



Джерело: Хабрахабр

0 коментарів

Тільки зареєстровані та авторизовані користувачі можуть залишати коментарі.