目录
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);
中国剩余定理
可以参考搜狗百科的内容,写的也是相当详细。
描述的是一个物不知数的问题:有物不知其数,三三数之剩二,五五数之剩三,七七数之剩二。问物几何?
等价为同余方程组:x = 2(mod 3) , x = 3(mod 5) , x = 2(mod 7)
然后有秦九韶解法:70*2 + 21*3 + 15*2 = 233 =2*(3*5*7)+23。23就是最后的解。
推广到一般:有m1, m2, ... , mn两两互质,有任意的整数a1,a2,a3,...,an,同余方程组:
x = a1 (mod m1)....x = an(mod mn)有解
且:M = m1*m2*m3...*mn ,Mi = M/mi,Mi^*Mi = 1(mod mi)
通解为:x = a1*M1*M1^+....+an*Mn*Mn^ + k*M(k为任意整数)
代码如下:
%中国剩余定理1(模数互质)
clc;clear;
n = input('输入参数个数:');
if (n == 0)
return;
end
A = zeros(1,n);MM = zeros(1,n);inv_B = zeros(1,n);
B = zeros(1,n);M = 1;N = 0;
i = 1;j = 1;k = 1;
while(i<n+1)
fprintf('输入第%d组参数:',i);
A(i) = input('输入a的值:');
B(i) = input('输入m的值:');
if(A(i) == 0||B(i) == 0)
return;
end
i = i + 1;
end
%求M
while(j<n+1)
M = M*B(j);
j = j + 1;
end
fprintf('M的值为:%d\n',M);
while(k<n+1)
MM(k) = floor(M/B(k));
k = k + 1;
end
fprintf('Mi:');
disp(MM);
i = 1;j = 1;k = 1;
while(i<n+1)
if(gcd(MM(i),B(i))~=1)
fprintf('Mi的逆元是0\n');
inv_B(i)= 0;
continue;
end
inv_B(i) = getinv1(MM(i),B(i));
fprintf('M的逆元为%d\n',inv_B(i));
N = N+A(i)*MM(i)*inv_B(i);
i = i + 1;
end
N = mod(N , M);
fprintf('方程的解:%d',N);
运行实例:
输入参数个数:3
输入第1组参数:输入a的值:2
输入m的值:15
输入第2组参数:输入a的值:3
输入m的值:17
输入第3组参数:输入a的值:4
输入m的值:19
M的值为:4845
Mi: 323 285 255
M的逆元为2
M的逆元为4
M的逆元为12
方程的解:2417
小结:
算法逻辑的核心还是数学,找到数学核心表达式,也就很简单了。
各位看官。瑞斯拜。