二次比較を解くための関数。 MATLABでの実装

はじめに

暗号の問題を解決するには、特定のモジュールで2次比較を解決できる必要があります。 二次比較を解くためのアルゴリズムは非常に単純であり、モジュールの小さな値と自由項を解くのに困難を引き起こしませんが、暗号化で十分に大きな数を使用するため、手動で二次比較を解くことは非常に骨の折れる長いプロセスです。 もちろん、二次比較を解決するには、オンラインサービスを使用できます。 しかし、暗号問題の解決は二次比較の解決で終わるわけではないので、暗号に携わる人が二次比較を解決し、それによって使用される他の関数と自由にやり取りできる機能を持っていると便利です。 そのため、MATLABでx ^ 2≡a（mod p）の形式の2次比較を解く関数を作成することにしました。ここで、aとpは互いに素な数です。



二次比較を解くための関数を書くことは本質的に教育的であり、計算で使用されるユーザー関数の一部は、MATLAB環境で既に利用可能な関数を複製していることにすぐに注意します。

最初に、2つの主要な関数のコードを検討することを提案します。1つは複合モジュールで2次比較を解決するためのもので、2つ目は単純なモジュールで2次比較を解決するためのものです。 同時に、まず二次比較を解くためのアルゴリズムに精通し、次に、計算自体を実行するために必要な関数に精通します。

複雑なモジュール比較を解決するための関数

この関数を使用すると、単純モジュラスと複素モジュラスの両方で2次比較を解決できます。 関数を呼び出すとき、2つの変数aとpをそれに渡す必要があります;その実行の結果として、関数はベクトル-符号が反対の2次比較の2つの解の文字列を返します。

function [ result ] = sqcomdif( a, p ) 

二次比較を解く次のステップは、モジュールpのタイプを決定することです。 これを行うには、数値を素因数に分解するように設計されたユーザー定義関数因数分解を使用します。 素因数を持つ行ベクトルに加えて、関数は素因数の数を返します。 実際、 因子分解関数は標準のMATLAB 因子関数を複製します。

 [ mp, sp ] = factorization( p ); 

条件演算子を使用してモジュールが因子分解された後、因子の数がチェックされます。 因子の数が1を超える場合、つまりモジュールpが合成数である場合、 sp二次比較のシステムを解く必要があります（各比較では、合成モジュールpの因子の1つがモジュールとして機能します）。 得られた2次比較のシステムの解法を実行する前に、このシステムの2次比較のそれぞれに解があることを確認する必要があります。 これを行うには、乗数mpでベクトルの要素間を遷移するforループを使用します。 ループの本体では、数値の各ペアのルジャンドル記号の値を計算する関数が呼び出されます。

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

Legendreシンボルの値が1に等しい場合、変数カウントは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ループを使用して、モジュールpの因子間で遷移が行われます。 answer1の 2次比較の結果は、 sqcom関数を使用して取得される単純なモジュールになります。このモジュールには、変数aの値とモジュールpの i番目の因子の値がanswer1ベクトル行に書き込まれます。 モジュールpのi番目の係数による複合モジュールpの除算の商は、行ベクトルmodulに書き込まれます。 線形不等式（p / p（I））* y≡1（p（i））を解いた結果は、中国の定理から得られた式に従って最終的な答えを計算するために必要であり、 answer2ベクトル行に保存されます。

 %        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コマンドを使用して見つけることができます。合計を複合モジュールpで除算した残りの部分を見つけます-これは複合モジュールによる2次比較の解決策の1つになります。

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

最初にモジュールpが素数であることが判明した場合、二次比較-sqcomを解く関数の1回の呼び出しで二次比較解が得られます。 2番目の解は、反対の符号を持つ最初の回答を取得することにより得られます。

 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関数は単純なモジュールで2次比較を解決するために使用され、 sqcomdif関数に関係なく呼び出すことができます。つまり、問題なく呼び出すことができ、正しい答えを得ることができますモジュールが素数であること。 x ^ 2≡a（mod p）の形式の2次比較の剥奪のみが考慮されるため、変数aおよびpの数値を関数に転送する必要があります。 関数の結果として、2次比較の1つの解が得られます。

 function [ answer ] = sqcom( a, p ) 

sqcom関数はsqcomdif関数とは別に使用できるため、変数aに書き込まれる数値がモジュールpの 2次剰余であることを確認する必要があります。 これを行うには、 symvol関数を使用します。これにより、指定された数値のペアのルジャンドル記号の値を計算できます。

 [ Symvol_Lejandra ] = symvol( a, p ); 

変数Symvol_Lejandraの値が1の場合、aはpを法とする2次剰余であり、2次比較の解を見つけるためにさらにステップが実行されます。 2 ^ r * qの形式で数字（p-1）を書く必要があります。 変数r 、 qの初期値は、 （p-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 = a ^ q（mod p）に等しい変数bの値を見つける必要があります。 通常、aとqは十分に大きい数で表されるため、ほとんどの場合オーバーフローが発生するため、通常の方法でaをq乗することはできません。 したがって、べき乗は平方の方法で実行する必要があります。 これを実現するには、 kvadrirovanie関数を呼び出して、ベース、指数、および計算を実行するモジュールの値を渡す必要があります。

 b = kvadrirovanie( a, q, p ); 

計算を続行するには、最小の非負数fを見つける必要があります。これは、 pを法とする2次の非剰余数になります。 この変数fには値1が割り当てられ、関数symvolを使用して、数値fとpのペアのルジャンドル記号の値が決定されます。 1とpのルジャンドル記号が1の場合、変数fは、 whileループの助けを借りて、値に達するまで増加します。これは、 pを法とする2次非剰余です。

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

ここで、 pを法とする2次の非剰余であるfの値をqの累乗に上げる必要があります。 これを行うには、関数kvadrirovanieを使用する必要があります。変数kには値0を割り当てる必要があります。

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

上記の手順の後、変数bの値を確認する必要があります。 bが 1モジュロpに匹敵する場合、答えの計算に進む必要があります;さもなければ、 b ^（2 ^ m）≡1（mod p）となる最小の非負数mを見つけます。 このようなmの値が見つかった場合、変数k 、 g 、 bの値を再計算し、値mを変数rに割り当てる必要があります。 しかし、これだけではありません。変数bの新しい値もp modulo 1と比較可能であることを確認する必要があります。 そうでない場合は、番号mの選択に戻る必要があります。 変数pokは 、特定の数学演算を2回実行することを避けるために必要です。

 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 

上記の条件を満たすmが見つかった後、答えの直接計算に進むことができます。 答えは、式x = a ^（（q + 1）/ 2）* g ^（k / 2）（mod p）によって計算されます。 両方の係数を計算するには、2乗関数を使用し、 pを法とする結果を取得します。

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

得られた結果は常に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^(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 

そして今、二次比較を解く際に使用される補助関数に精通することを提案し、それらを別々に使用できるようにします。

数値の因数分解

二次比較を解く場合、多くの場合、数値の因数分解に頼る必要があり、数値が素数であるか複合であるかを確認する必要がある場合にも、この操作を使用する必要があります。
因数分解関数には、 因数分解する必要がある数値が渡されます。 結果として、関数はベクトルを返します-因子とこれらの因子の数を持つ行。

 function [ mnojitel, ind ] = factorization( delimoe ) 

この関数は、入力変数delimoeの値に応じてさまざまなアクションを実行するswitchステートメントで構成されています。 したがって、 delimoe = 1の場合、 mnojitelは因子分解の結果を格納するベクトルの1に書き込まれ、 1は因子の数を格納する変数indに も書き込まれます。 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に等しくなるまで実行され、値deltaは最初に2に 設定されます。 ループの各反復で、 デリモをデリテルで除算した余りがostatok変数に 書き込まれ ます 。 剰余が0の場合、つまり、 デリテルがデリモファクターである場合、この値は、ファクターが格納されているベクトルに書き込まれ、このベクトルのカウンターは1増加します。 この場合、変数delimoeには 、実行された除算の商が割り当てられます。 除算の残りが0に等しくない場合、 デリテルは1 ずつ増加します。 ループが終了すると、 delimoe変数に残っている値が、要因の1つとして、要因とともにベクトルに書き込まれます。

 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 ] = 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 

ルジャンドル記号の値の計算

数値が2次剰余モジュロ（この場合、2次比較には解がある）か2次剰余（そのような2次比較には解がない）かどうかを確認するために、ロシア文学ではL（a; p）外国文学ではL（a / p）として 。
ルジャンドル記号は、次の意味をとることができます。
L（a; p）= 1 、この場合aはQRに属し、二次比較には解があります
L（a; p）= -1 、この場合aはQNRに属し、2次比較には解がありません
L（a; p）= 0の場合、aとpは互いに素ではありません。つまり、 GCD（a; p）は 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. a≡b（mod p）の場合、 L（a; p）= L（b; p）
6. aとpが素数の場合、 L（a; p）=（-1）^（（（p-1）*（a-1））/ 4）* L（p; a） 。 最後の特性は、ガウス相反則と呼ばれます。

ルジャンドル記号の計算は、上記のプロパティに基づいています。 条件の1つが満たされるとすぐに、ルジャンドル記号の最終値が見つかるまで、結果のペアaとpのプロパティのチェックを開始します。
次に、ルジャンドル記号の値の計算をプログラムで実装する方法を検討します。
この関数は、転送された数字のペアaとpのルジャンドル記号の値を返します。 これは関数ヘッダーから見ることができます：

 function [ sl ] = symvol( a, p ) 

次のステップは、本質的にプロパティ5を適用することです。 数aがモジュールpより大きい場合、モジュロpに匹敵するより小さい数で置き換えることができます。

 a=mod( a, p ); %      

結果の数a 、因数分解しようとしています。 数値を因数分解するために、数値aとその数値を構成する単純な因子を含むベクトルを返すカスタム因数分解関数が作成されました。 この機能については、上で詳しく説明しました。

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

aが素数でない場合、プロパティ4に進みます。 つまり、数値のペアL（a; p）のルジャンドルsymbol記号の値は、 aの単純な因子である数値のルジャンドルsymbols記号の値の積として求められます。 中間結果を保存するために、aの因子の数に等しい次元を持つゼロで満たされた行ベクトルを作成します。 aが素数の場合、この行は1つの要素で構成されます。

 sl = zeros( 1, ind ); %           

各因子のルジャンドル記号を計算するには、代わりにforループを使用します。これにより、値が1から最後の因子の数に変更されます。 このサイクルの本体には、上記のプロパティを使用してルジャンドル記号の値を計算する直接プロセスがあります。
プロパティ1をチェックするコードは次のとおりです。

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

プロパティ2に従ってL（-1、p）の形式でシンボルの値をチェックするとき、値（-1）^（（p-1）/ 2）を計算する必要があるため、もう1つの条件演算子を使用する必要があります。オンの場合、インジケータ-1は偶数または奇数です。 これに応じて、ルジャンドル記号の意味は異なります。 指数が偶数の場合、ルジャンドル記号は1に等しくなり、そうでない場合は-1になります。 この条件演算子を使用すると、 -1の（p-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 

この条件を確認した後、単純かどうかを確認するために、数値aが関数に渡されます。 数値aが複合の場合（そのind1因子の数が1 より大きい ）、再帰が発生し、数値aが同じ関数に転送されて、さらに計算が実行されます。

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

それ以外の場合、数値aは素数です。 同時に-1、1、2に等しくない場合、プロパティ6-ガウス相反則を使用する必要があります。 ルジャンドル記号の前の記号は、その指標の要素のパリティを決定することにより決定されます。 要因の少なくとも1つが偶数であれば、プラスに変わります。 その後、 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 

上記の条件をチェックした結果、aの値のすべての可能なバリアントがカバーされました。
, 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. .
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であり、指数まで上げる必要があった元の数値自体が結果変数に書き込まれます。

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

この条件が満たされない場合、指数は0です。この場合、結果変数に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. http://www.wolframalpha.com