Criska

数论

• 分类: 高斯消元, 扩展欧几里得 • 标签: 扩展欧几里得, 高斯消元

1.高斯消元

高斯消元法主要用于求一个多元一次方程的解。
例如，对于以下的式子：

2x + y - z = 8 (L1)
-3x - y + 2z = -11 (L2)
-2x + y + 2z = -3 (L3

我们可以这样消元：

  • 将L1作为起点，也就是说，我们要通过L1这个式子最终消元后的结果得出x的值，那么就要通过L2，L3将L1中的y，z的系数消为零，同时通过L1将L2，L3中x的系数消为零，这样最后我们就得到了一个关于x的一元一次方程。然后依次类推，L2得到y，L3得到z。
  • 当然，方程组也有可能无解，比如这样的：

x+y+z=1
2x+2y+2z=2
3x+3y+3z=3

  • 如何判断无解呢？其实很简单，假如在消元的过程中我们将一个等式等号左边的所有项的系数都消成零了，那其实就意味着这两个等式是完全相等的。而如果这个方程组中有两个意义一样的方程，这就意味着这个方程组一定是无解或者多解。

具体代码如下：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
#include<iostream>
#include<cstdio>
using namespace std;
int n,mark=0;
double a[200][200],b[120],da[120];
int main()
{
	scanf("%d",&n);
	for(int i=1;i<=n;++i)
	{
		for(int j=1;j<=n+1;++j)
		{
			cin>>a[i][j];
		}
	}
	for(int i=1;i<=n;++i)
	{
		if(mark==1) break;
		for(int j=1;j<=n;++j)
		{
			if(j==i) continue;
			double x=a[j][i]/a[i][i];
			if(x==0) continue;
			for(int k=1;k<=n+1;++k)
			{
				a[i][k]*=x;
				a[j][k]=a[i][k]-a[j][k];
			}
			for(int k=1;k<=n;++k)
			{
				if(a[j][k]!=0) break;
				if(k==n) 
				{
					mark=1;
					break;
				}				
			}
		}
		
	}
	if(mark==1) printf("No Solution");
	else 
	for(int i=1;i<=n;++i) 
	printf("%.2lf\n",a[i][n+1]/a[i][i]);
	return 0;
}

2.扩展欧几里得算法

-同余方程

  • 拿noip的这道题做例子，对于 ax ≡ 1 (mod b)，它完全可以转化为ax=yp+1，也就是 ax+(-y)b=1（y是设出来的一个值）
    显而易见，只要a的系数x=1，只要b的系数是0，这时，我们就会有：a⋅1+b⋅0=gcd(a,b)。这是这个方程的最终结果，而我们需要根据这个结果来推出的初始值。
    假设我们要求出a,b的最大公因数，自然在下一步已经求到了b和a mod b的最大公因数（递归就是这么反人类。。。由未知找已知），并且求出了一组x1和y1使得：b⋅x1+(a%b)⋅y1=gcd
    • 接下来我们可以推出：

ax+by=b*x1+(a%b)*y1
设a÷b=c……d
a=b*c+d
d=a-b*c
c=a/b（由于int是整型，我们正好得出整数部分）
d=a-b*(a/b)
所以
ax+by=b*x1+(a%b)*y1
ax+by=b*x1+(a/b)*b*y1
ax+by=a*y1+b*(x1–[a/b]*y1)

  • x=y1

  • y=x1–[a/b]*y1

    这样我们就得出了递归的方程。
    代码如下：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
#include<iostream>
#include<cstdio>
using namespace std;
int a,b,x,y;
inline int gcd(int a,int b,int &x,int &y)
{
	if(b==0)
	{
		x=1;
		y=0;
	}
	else
	{
		gcd(b,a%b,x,y);
		int t=x;
		x=y;
		y=t-(a/b)*y;		
	}
}
int main()
{
	scanf("%d%d",&a,&b);
	gcd(a,b,x,y);//ax=1(mod b)
	x=x%b;
	if(x<0) x+=b;
	printf("%d",x);
	return 0;
}

-除法取模

  • 我们都知道除法不可以直接取模，但是我们可以对除数的逆元进行取模。

    • 例如 （T/K）%P=(T*K1)%P=(T%P *K1%P)%P
      其中K1就是K的逆元

  • 那么K1该怎么求呢？这里有一个方程：

    • K*K1≡1 (mod P)

    显而易见，这是一个同余方程，我们用上面的方法就可以求出了。