中国剩余定理的练习和求逆元，扩展欧几里得算法matlab小结

偶尔敲代码的创作猿

已于 2022-05-26 17:24:04 修改

阅读量1k

点赞数 3

文章标签：算法线性代数

于 2022-05-26 17:11:15 首次发布

本文链接：https://blog.csdn.net/weixin_55339858/article/details/124985040

版权

1.（k > p）k = a*p + r , k = r(mod p )

2.（k < p）p = a*k + r , k = -r(mod p )

中国剩余定理

运行实例：

小结：

求逆元(inv)

首先要介绍求余运算：a = k*p+m,则mod(a , p) = m,写做a = m(mod p)

[注：mod(a,p)是matlab求余函数]

要求a与p互质时a有逆元a^，a*a^= 1(mod p).

求逆元有两种常用的方法，分别是扩展欧几里得算法和公式线性破解。

扩展欧几里得算法求逆元

扩展欧几里得算法

欧几里得算法又称辗转相除法，是用来求两个数最大公因数（gcd）的算法。

核心数学公式：gcd(a , b) = gcd(b , mod(a , b))，通过迭代求解出最后的结果。

扩展欧几里得是通过对裴蜀定理的证明得出的：

说明了对任何整数 a、b和它们的最大公约数d，关于未知数x和y的线性丢番图方程（称为裴蜀等式）：有解当且仅当 m是d的倍数。

即方程：ax + by = m*gcd(a , b)时，方程有解

特例：当a与b互质时,有ax + by =1成立。（也是求解逆元的核心方程）

接下来把辗转相除法加入，就得到扩展欧几里得算法了：

首先有a*x1 + b*y1 = gcd(a , b)和b*x2 + mod(a , b)*y2 = gcd(b , mod(a , b))

因为欧几里得算法的核心：gcd(a , b) = gcd(b , mod(a , b))

所以有：a*x1 + b*y1 = b*x2 + mod(a , b)*y2

解出：x1 = y2 , y1 = x2 - floor(a /b)y2

[注：因为matlab没有地板除，除法取整只能再加函数，这里用floor，向原点方向取整]

然后通过a ,b的辗转相除法,x ,y的迭代，最后得到：

an *xn + 0*yn = gcd(a , b)

根据高数辗转相除法的计算结果可知，当bn为0时,an = gcd(a，b)

即xn = 1, yn = 0(这里取最小，确保是最小解)，然后通过迭代式就可以求出x1和y1的最小值。当然该方程也有解。

扩欧求逆元

a *a^ = 1(mod p)可改写为：a*a^ + m*p = 1。即求a^的最小值。

梳理一下迭代的核心式子：

an = bn-1 , bn = mod(an-1,bn-1)

xn = yn-1 + floor(a /b)xn-1 , yn = xn-1

最后:xn = 1,yn= 0。

那么，有两种实现代码的方式,第一种是正向迭代，迭代到最后，赋值:xn = 1,yn= 0，然后程序返回x和y的值。matlab自定义函数迭代功能好像不支持，代码无法正常运行。所以不采用。

第二钟就是从结果式子an *xn + 0*yn = 1，逆向迭代，直到回到原始式子。由于辗转相除法的前提，ax + by = 1,a必须大于b。所以，如果遇到a<p的情况，需要a加上一定量的p,并且加后的a必须和p互质。代码如下：

%扩展欧几里得算法
function f = getinv1( a ,b)
while(a<b)
    a = a + b;
    while(gcd(a,b)~=1)
        a = a + b;
    end
end
A = [];B = [];
A(end+1)= a;B(end+1)= b;
x = 1;y = 0;i = 1;c = 0;
while(b ~= 0)
    c = b;b = mod(a , b);i = i + 1;
    a =c ;
    A(end+1)= a;B(end+1)= b;
end
while(i ~= 1)
    a = A(i-1);b = B(i-1);
    i = i - 1;
    c = y;
    y = x - floor(a/b)*y;
    x = c;
end
while(x<0)
    x = x + b;
end
f = x;

代码解释

第一部分是要求适用于a<p的情况；

第二部分先正向推演，把an和bn储存下来，因为每一步的x和y迭代需要用到对于的a和b计算，并且计算出辗转相除法了i次;

第三部分就是逆向推导，从结果式出发，求解原x , y，最后输出x时,要确保输出结果为正，加适量p。

公式线性破解

和扩展欧几里得一样，从公式出发。(需要简单的数学运算)

先看看k>p的情况

1.（k > p）k = a*p + r , k = r(mod p )

暴力破解：两边同时乘以k^和r^,得到k*k^*r^ = a*p*k^*r^ + r*r^*k^

接着，逆元的定义式：k*k^ = m*p + 1,r*r^ = n*p + 1

带入，合并同类项：k^ = r^ + (m*r^ - a*k^r^ + n*K^)p

即：k^ = r^(mod p) , k^ = mod(k , p)^(mod p)

则inv(k) = inv(mod(k , p))[注：inv表示逆，inverse element]

这个式子无法迭代，只能把k缩小到小于p的情况

接下来讨论k<p的情况

2.（k < p）p = a*k + r , k = -r(mod p )

步骤同上:k^ = - a*r^ + (k^r^ - a*m*r^ - k^*n)p，即：k^ = -a*r^(mod p)

inv(k) = p -inv(mod(p , k))*floor(p\k)

这是一个迭代式：inv(kn-1) = p - inv(kn)*floor(p\kn-1)

kn = mod(p , kn-1)

因为k与p互质，所以设置极限inv(1) = 1，因为递归不支持，所以还是回溯。

因为迭代中需要有mod(p , k),所以要求mod(p , kn ）不为0 ，即p为素数。

代码如下：

function[e] = getinv2(a , mod) 
%1.ak+r=mp,k=-r(mod p),k<p的情况
%2.k=ap+r,k=r(mod p),k>p的情况
%3.只适用于p为素数的情况
if(a>mod)
    a = rem(a , mod);
end
b = 1;
c = zeros(1,a);
c(1)=1;
while(b < a)
    b = b + 1;
    if(rem(mod,b)==0)
        c(b)=0;
        continue;
    else
        c(b)=(mod-floor(mod/b))*c(rem(mod,b));
        c(b)=rem(c(b),mod);
    end
end
e = c(a);