Functions for solving quadratic comparisons. Implementation in MATLAB

Introduction

To solve cryptographic problems, it is necessary to be able to solve quadratic comparisons of a given module. The algorithm for solving quadratic comparisons is quite simple and does not cause difficulties in solving for small values ​​of the modulus and the free term, however, due to the use of sufficiently large numbers in cryptography, solving quadratic comparisons manually is a very laborious and lengthy process. Of course, to solve quadratic comparisons, you can use the online service. But since the solution of a cryptographic problem does not end with solving a quadratic comparison, it will be convenient for a person involved in cryptography to have a function that can solve quadratic comparisons and freely interact with other functions used by it. That is why it was decided to write a function for solving quadratic comparisons of the form x ^ 2 ≡ a (mod p) , where a and p are mutually prime numbers, in MATLAB.



Immediately, I note that the writing of functions for solving quadratic comparisons was educational in nature; some user functions used in the computation duplicate the functions already existing in the MATLAB environment.

First, I propose to consider the code of two main functions, one of which is designed to solve quadratic comparisons in a composite module, and the second to solve quadratic comparisons modulo a simple module. In this case, we will get acquainted, first of all, with the algorithm for solving quadratic comparisons, and secondly, with the functions that are necessary to perform the calculations themselves.
')

Function to solve complex modulus comparisons

This function allows you to solve quadratic comparisons, both in simple and complex modulus. When calling a function, it is necessary to transfer two variables a and p , as a result of its execution, the function will return a vector - a string of two solutions of quadratic comparison opposite in sign.

function [ result ] = sqcomdif( a, p ) 

The next step in solving a quadratic comparison is to determine the type of the module p . For this, the custom factorization function will be used, designed to decompose a number into prime factors. In addition to the string vector with prime factors, the function returns the number of prime factors. In fact, the factorization function duplicates the standard MATLAB factor function.

 [ mp, sp ] = factorization( p ); 

After the module has been factorized using the conditional operator, the number of factors is checked. If the number of factors exceeds 1 , that is, the module p is a composite number, then we need to solve a system of sp quadratic comparisons (in each of the comparisons the module will be one of the factors of the composite module p ). Before tackling the obtained system of quadratic comparisons, it is necessary to make sure that each of the quadratic comparisons of this system has a solution. To do this, we use the for loop, which will make the transition between the elements of the vector with the factors mp . In the body of the loop, a function will be called that calculates the value of the Legendre symbol for each pair of numbers.

 for i=1:1:sp SL( 1, i ) = symvol( a, mp( 1, i ) ); 

If the value of the Legendre symbol is 1 , then the count variable will be incremented by 1 . This is necessary so that after all iterations of the cycle have been completed, we can check whether all the equations of the system have a solution, because it depends on whether we will solve the quadratic comparisons of which the system is composed, or we will display messages that the original quadratic comparison does not making.

  if SL( 1, i ) == 1 count = count + 1; %      count   1 end 

If the number of equations of the system is equal to the number of equations that have solutions, then row vectors are created to store intermediate results of the calculation.

 if count == sp %       answer1 = zeros ( 1, sp ); %     modul = zeros ( 1, sp ); %         answer2 = zeros ( 1, sp ); %          .  

Using a for loop, a transition is made between the multipliers of the p module. The vector vector answer1 records the results of solving a quadratic comparison modulo a simple module, which is obtained using the sqcom function, to which the value of the variable a was passed, as well as the value of the i-th factor of the module p . In the vector-line modul , the quotient of dividing the composite module p by the i-th factor of the module p is written. The answer2 vector will contain the result of solving a linear inequality (p / p (I)) * y ≡ 1 (p (i)) , which will be necessary to calculate the final answer using the formula obtained from the Chinese theorem.

 %        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 

After the cycle has completed its execution, it is necessary to calculate the final answer using the formula: x = ((p / p (1)) * b (1) * y (1) + ((p / p (2)) * b (2 ) * y (2) + ((p / p (i)) * b (i) * y (i) . To do this, we use element-wise multiplication. As a result, we obtain a row vector, the sum of the elements of which we will find using the command sum . We will find the remainder of dividing the sum by the composite module p — this will be one of the solutions of the quadratic comparison modulo. The second answer will be obtained by writing the first answer with the opposite sign.

  %    result = zeros ( 1, 2 ); result( 1, 1 ) = mod( sum( ( modul .* answer1 ) .* answer2 ), p ); result( 1, 2 ) = - result( 1, 1 ); 

In the event that initially the module p turned out to be a prime number, then the solution of a quadratic comparison will be obtained by a one-time function call for solving a quadratic comparison - sqcom . The second solution will be obtained by taking the first answer with the opposite sign.

 else result( 1, 1 ) = sqcom( a, p ); result( 1, 2 ) = - result( 1, 1 ); end 

Below is the full sqcomdif function code.

 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 

Function to solve quadratic comparisons modulo simple

This function was repeatedly called in the sqcomdif function, and as already mentioned, the sqcom function is used to solve quadratic comparisons in a simple module and can be called independently of the sqcomdif function, that is, you can call it and get the right answer without any problems when you are sure that the module is a prime number. Since only deprivation of quadratic comparisons of the form x ^ 2 ≡ a (mod p) is considered , it is necessary to transfer the numerical values ​​of the variables a and p to the function. As a result of the function, one solution of the quadratic comparison will be obtained.

 function [ answer ] = sqcom( a, p ) 

Since the sqcom function can be used separately from the sqcomdif function, it is necessary to make sure that the number written in the variable a is a quadratic residue of the module p . To do this, we use the symvol function, which allows us to calculate the value of the Legendre symbol for a specified pair of numbers.

 [ Symvol_Lejandra ] = symvol( a, p ); 

If the value of the Symvol_Lejandra variable is 1 , then a is a quadratic residue modulo p and further steps will be taken to find a solution to the quadratic comparison. We need to write the number (p - 1) in the form 2 ^ r * q . The initial values ​​of the variables r , q are given on the basis that (p - 1) is an odd number. However, if this is not the case, they will change during the execution of the cycle (until q becomes an odd number).

 q = p - 1; r = 0; otn = q / 2; while ( ( q - floor( otn ) * 2 ) == 0 ) q = otn; r = r + 1; otn = q / 2; end 

Now it is necessary to find the value of the variable b , which equals b = a ^ q (mod p) . Since a and q are usually expressed by sufficiently large numbers, it will not be possible to raise the number a to the power of q in the usual way, since in most cases an overflow will occur. Therefore, the exponentiation must be performed according to the quadriration method. In order to accomplish this, it is necessary to call the kvadrirovanie function and transfer to it the value of the base, the exponent, and the module on which the calculations will be performed.

 b = kvadrirovanie( a, q, p ); 

To continue the calculations, you need to find the smallest non-negative number f , which will be a quadratic non-residue modulo p . For this, the variable f is assigned the value 1 , and using the symvol function, the value of the Legendre symbol is determined for a pair of numbers f and p . If the Legendre symbol for 1 and p is 1 , then the variable f , using the while loop , increases until a value that is a quadratic non-residue modulo p is reached.

 f = 1; sym_lej = symvol( f, p ); while sym_lej ~= -1 f = f + 1; sym_lej = symvol( f, p ); end 

Now the value of f , which is a quadratic non-residue modulo p , must be raised to the power of q . For this, it is also necessary to resort to the kvadrirovanie function, the variable k must be assigned the value 0 .

  g = kvadrirovanie( f, q, p); k = 0; 

After the above actions, you need to check the value of the variable b . If b is comparable with 1 modulo p , then it is necessary to proceed to the calculation of the answer, otherwise find the smallest non-negative number m such that b ^ (2 ^ m) 1 (mod p) . When such a value of m is found, it is necessary to recalculate the values ​​of the variables k , g , b , after which the variable r is assigned the value of m . But this is not all, it is necessary to check that the new value of the variable b is also comparable to 1 modulo p . If this is not the case, then it is necessary to return to the selection of the number m . The variable pok is needed to avoid performing certain mathematical actions twice.

 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^(rm); g1 = kvadrirovanie( g, pok, p ); b = mod( ( b*g1), p ); k = fix(k + pok); r = m; end end 

After m has been found that satisfies the above conditions, we can proceed to the direct calculation of the answer. The answer is calculated by the formula x = a ^ ((q + 1) / 2) * g ^ (k / 2) (mod p) . To calculate both factors, we use the quadration function; we take the result modulo p .

  first = kvadrirovanie( a, ( ( q + 1 ) / 2 ), p ); second = kvadrirovanie( g, ( k / 2 ), p ); answer = mod( ( first * second ), p); 

The result obtained cannot always be written optimally, modulo p . Therefore it is necessary to perform the following check:

  delta = p - answer; if delta < answer answer = delta; end 

The following is the full sqcom function code :

 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^(rm); 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 

And now I propose to become familiar with the auxiliary functions that are used both in solving quadratic comparisons, and so can be used separately.

Factor expansion

When solving a quadratic comparison, in many cases it is necessary to resort to factorization of the number, and this operation will also have to be resorted to when it is necessary to check: is the number simple or composite.
The factorization function is passed a number that needs to be factorized. As a result, the function returns a vector - a string with multipliers and the number of these multipliers.

 function [ mnojitel, ind ] = factorization( delimoe ) 

A function consists of a switch statement that performs various actions depending on the value of the input variable delimoe . So if delimoe = 1 , then mnojitel is written 1 to the vector storing the result of the factorization, and the variable ind , which stores the number of factors, is also written 1 . Similar actions are performed in the event that delimoe is -1 .

 switch delimoe case { 1 } mnojitel( 1, 1 )=1; ind=1; case { -1 } mnojitel( 1, 1 )= -1; ind = 1; 

If these conditions are not met, then we check the sign of the number factorized. If delimoe is less than 0 , then -1 is written to the first factor, 2 is written to ind , and the function will continue to work with the module of the delimoe variable . If the number is greater than 0 , then the value of ind is assigned to 1 .

  otherwise if delimoe < 0 mnojitel( 1, 1 )= -1; ind = 2; delimoe = abs ( delimoe ); else ind = 1; end 

The while loop runs as long as delimoe is not equal to delitel , with the initial delitel variable being assigned the value 2 . At each iteration of the loop, the remainder of delimoe divided by delitel is written into the ostatok variable . If the remainder is 0 , that is, the delitel is a delimoe multiplier, then this value is written to the vector in which the multipliers are stored, and the counter of this vector is increased by 1 . In this case, the delimoe variable is assigned the quotient of the division performed. If the remainder of the division is not equal to 0 , then the delitel is increased by 1 . When the exit from the loop occurs, the value remaining in the delimoe variable is written to the vector with factors as one of the factors.

 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; 

Below is the full code of the factorization function.

 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 

Calculating the value of the Legendre symbol

In order to check whether a number is a quadratic residue by a given module (in this case, a quadratic comparison has a solution) or a quadratic non-residue (then such a quadratic comparison has no solution), the Legendre symbol is used, which is denoted in the Russian literature as L (a; p) , in foreign literature as L (a / p) .
The Legendre symbol can take the following values:
L (a; p) = 1 , in this case a belongs to QR, the quadratic comparison has the solution
L (a; p) = -1 , in which case a belongs to QNR, the quadratic comparison has no solution
If L (a; p) = 0 , then a and p are not mutually simple numbers, that is, the gcd (a; p) is not equal to 1
The following properties are used to calculate the value of the Legendre symbol:

1. L (1; p) = 1
2. L (-1; p) = (-1) ^ ((p - 1) / 2)
3. L (2; p) = (-1) ^ ((p ^ 2 - 1) / 8)
4. If ≡ b * (mod p) , then L (b * ; p) = L (b; p) * L (; p)
5. If ≡ b (mod p) , then L (a; p) = L (b; p)
6. If a and p are prime numbers, then L (a; p) = (-1) ^ (((p - 1) * (a-1)) / 4) * L (p; a) . The latter property is called the Gauss law of reciprocity.

The calculation of the Legendre symbol is performed on the basis of the above properties. As soon as one of the conditions is fulfilled, we begin to check the properties for the resulting pair a and p , until the final value of the Legendre symbol is found.
Now consider how you can implement the calculation of the value of the Legendre symbol programmatically.
The function will return the value of the Legendre symbol for a pair of numbers a and p that were passed to it. This is evident from the function header:

 function [ sl ] = symvol( a, p ) 

The next step is essentially the use of property 5 . If the number a is greater than the modulus p , then it can be replaced by a smaller number comparable to a modulo p .

 a=mod( a, p ); %      

The resulting number is, trying to factor out. For factoring a factor, a custom function factorization was written that returns a vector with prime factors that make up the number a , as well as their number. This feature has been described in detail above.

 [ mnoj, ind ] = factorization( a ); %      

If a was not a prime number, then we will act by property 4 . That is, the value of the Legendre symbol for a pair of numbers L (a; p) , we find as the product of the values ​​of the Legendre symbols for numbers that are prime factors of a . To store intermediate results, create a row vector filled with zeros with a dimension equal to the number of factors of a . If a is simple, then this string will consist of one element.

 sl = zeros( 1, ind ); %           

To calculate the Legendre symbol for each multiplier, alternately use the for loop, which will change its values ​​from 1 to the number of the last multiplier. In the body of this cycle, the immediate process of calculating the value of the Legendre symbol takes place using the properties listed above.
The code for checking property 1 is as follows:

  %   L(1,p) if mnoj( 1, i ) == 1 sl( 1, i ) = 1; end 

Since when checking the value of a symbol when it is L (-1, p) , by property 2 , it is necessary to calculate the value (-1) ^ ((p - 1) / 2) , it is necessary to use another conditional operator, in which be checked if the indicator -1 is an even or odd number. Depending on this, the value of the Legendre symbol will vary. If the index is even, then the Legendre symbol will be equal to 1 , otherwise -1 . The use of this conditional operator avoids the erection of -1 to the power of (p - 1) / 2 , which is a very costly operation.

  %   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 

Similar actions are carried out in the case when it is necessary to calculate the value of the Legendre symbol in its form L (2; p) . In this case, it will be (-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 

After checking this condition, the number a is passed to the function in order to check whether it is simple or not. If the number a is a compound (the number of its factors ind1 is greater than 1 ), then recursion occurs and the number a is transferred to the same function to perform further calculations.

  [ mn, ind1 ] = factorization( mnoj( 1, i ) ); %      if ind1 > 1 %    -  sl( 1, i ) = symvol( mnoj(1,i), p ); %       ,    

Otherwise, the number a is prime. If at the same time it is not equal to -1 , 1 , 2 , then one has to use property 6 , the Gauss reciprocity law. The sign in front of the Legendre symbol is determined by determining the parity of the factors of its index. He turns in plus if at least one of the factors is even. After that, a recursive call to the symvol function occurs , and the arguments are passed to it in a different order.

  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 

As a result of the verification of the above conditions, all possible variants of the number a were covered.
If initially the number a was composite, then the variable sl is written the product of all elements of the string - the vector sl , after which the resulting value is returned to the program that caused this function.

 if ind ~= 1 sl = prod( sl ); %   L(a,p) end 

The complete code for the symvol function :

 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 

Raising a number to a power by quadrature method

When solving a quadratic comparison, one often comes across the need to build large numbers to a power with an exponent equal to several hundred, or even thousands. With the usual raising of a number to such a degree, it is rather easy to come across such an unpleasant phenomenon as overflow. In order to avoid it, it is necessary to raise such numbers to a power using the quadration method.
Consider the algorithm of this method:

1. .
2. , 1 m=a .
3. :
   
   • m , , m .
   • , m , , m .
4. Go to the next bit on the right and repeat the steps described in paragraph 3, until the calculations for the last bit are performed.

And now let's look at how the kvadrirovanie function works . The values ​​of the base, indicator and module are transferred to the function, and the function returns a single number - the final result.

 function [ result ] = kvadrirovanie( a, q, p ) 

We translate the exponent q into a binary view, and since the result will consist of individual characters, we use the size command to determine the number of bits needed to record the exponent.

 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); 

In the exponent in binary form is equal to 1 , then the initial number itself is written into the result variable , which was required to be raised to a power.

 elseif q(1,1) == '1' result = uint64( a ); 

If this condition is not fulfilled, then the exponent is 0 . In this case, the result variable will be written 1 .

 else result = 1; end 

Full function code for exponentiation by the quadration method:

 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 

Linear comparison solution

When solving quadratic comparisons in complex moduli, the formula obtained from the Chinese theorem will be applied. In order to calculate the answer using this formula, it becomes necessary to solve the linear comparison k * x ≡ b (mod p) . This function is intended to solve them with the help of an inverse element, that is, such a number k1 which, when multiplied by k , will give a result comparable with 1 given a modulus p .
The lincom function passes the coefficient values ​​for an unknown k , the number b with which the left side is comparable, and also the module p , and the function returns respectively a linear comparison solution with the specified parameters.

 function [ x ] = lincom ( k, b, p) 

, . . , , . ( a, b ) , – b , k .
b0 b1 , . p , which was written to the variable pr , is divided completely into the coefficient k , then the inverse element will be equal to 1 . Otherwise, b1 will be calculated using a while loop . The cycle will recalculate the values ​​of the variable b0 . After a new value of b0 is found , the swap function will exchange the values ​​between the variables b0 and b1 . Also, with each iteration of the cycle, new values ​​will be assigned to a number of other variables.

 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 

After the cycle ends, it will be necessary to calculate the value of the linear comparison solution. It will be equal to the product of the variables b1 and b , taken modulo p (for this operation at the beginning of the function the modulus value was duplicated using the variable pr ).

 x = mod( b1*b, p ); 

The following is the entire function code:

 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 

The list of sources that helped in creating features and writing an article:

1. Summaries and lectures Plotnikova L.I. (who in the subject will understand)
2. http://math.hashcode.ru
3. http://mathhelpplanet.com
4. http://www.wolframalpha.com