原题链接 https://www.luogu.org/problemnew/show/P4783
一道模板题,更重要的省选难度.....
题目要求的是一个n*n的逆矩阵,还要对大数取膜。
普通高中生:这…………
来一步一步分析:
(1)怎么求逆矩阵?
首先,我们要开一个二维数组,范围是a[401][801];
why?这就和求逆矩阵有关啦。
先输入n*n的矩阵,紧接着在右边罗一个n*n的单位矩阵
然后我们对左半边的矩阵进行高斯消元消成了单位矩阵,则此时右半边的单位矩阵就被消成了左半边原矩阵的逆矩阵
scanf("%lld",&n); for(int i=1;i<=n;i++) { for(int j=1;j<=n;j++) scanf("%lld",&a[i][j]); a[i][i+n]=1; //在输入矩阵的同时在右边罗列一个n*n的单位矩阵 }
接着就要对左半边的矩阵进行高斯消元,怎么高斯消元呢?其实就是把高斯消元的板子套上去啊。P3389
这是模板高斯消元的代码:
#include<iostream> #include<cstdio> #include<cmath> using namespace std; int n,pl; double a[1001][1001]; int main() { cin>>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++) { pl=i; while(a[pl][i]==0&&pl<=n) {pl++;} // 判断第i列首元素非0的最上行,因为第i行第i列元素不能为0 if(pl==n+1) {cout<<"No Solution";return 0;} //一直判到了n+1行,可是一共才只有n行,说明有一列全为0,无解 for(int j=1;j<=n+1;j++) //将第i行第i列元素不为0的那一行与当前行交换 swap(a[i][j],a[pl][j]); double k=a[i][i]; //让第i行每个元素都除以a[i][i]使得a[i][i]为1 for(int j=1;j<=n+1;j++) a[i][j]=a[i][j]/k; //将第i行第i列的元素消成1,注意同行进行同样的操作 for(int j=1;j<=n;j++) { if(i!=j) //将第i列除了第i行的元素全消成0 { //方法是第j行每个元素a[j][m]都减去a[j][1]*a[i][m] double ki=a[j][i]; for(int m=1;m<=n+1;m++) a[j][m]=a[j][m]-ki*a[i][m]; } } } for(int i=1;i<=n;i++) printf("%.2lf ",a[i][n+1]); return 0; }
(2)怎么对有理数取膜?
这也是道模板题…………P2613
#include<iostream> #include<cstdio> using namespace std; const long long mod=19260817; inline long long read() { long long t=0; char ch=getchar(); while(ch<'0'||ch>'9') ch=getchar(); while(ch>='0'&&ch<='9') { t=(t*10+(ch-'0'))%mod; ch=getchar(); } return t; } int exgcd(long long a,long long b,long long &x,long long &y) { if(b==0) { x=1;y=0; return a; } long long r=exgcd(b,a%b,x,y); long long q=x; x=y; y=q-a/b*y; return r; } int main() { long long a,b,x,y; a=read(); b=read(); if(b==0) { cout<<"Angry!"; return 0; } exgcd(b,mod,x,y); x=(x%mod+mod)%mod; printf("%lld",(a%mod*x%mod)%mod); return 0; }
证明过程:
1.费马小定理
如果p是一个质数,而整数a不是p的倍数,则有a^(p-1)≡1(mod p)。
此题已经明确给出mod数19260817,显然它是一个质数,那么我们就可以用费马小定理转化一下,如下:
因为a^(p-1)≡1(mod p)
所以a^(p-2)≡a^(-1) (mod p) (A)
所以c=a/b=a*b^(-1)≡a*b^(p-2) (mod p)
证毕!
所以我们就可以将在膜p意义下的a/b转化成a*b^(p-2)的形式,所以我们只要求出b^(p-2)就大功告成啦,具体做法用快速幂。
2.扩展欧几里德
上面已经证过求在膜p意义下的a/b就是求a*b^(-1),b^(-1)就是b的逆元
下面给出求b的逆元的一种方法:
若存在一个数x,满足bx≡1 (mod p),那么x就是b的逆元
可将bx≡1 (mod p)进一步转化:
bx-1≡0 (mod p)
bx-1=-yp (注:这里说一下为什么是-y,其实这里是不是正负无所谓,写成负的更便于理解)
bx+py=1
二者一结合,就是此题的代码:
#include<iostream> #include<cstring> #include<cstdio> #include<cmath> #include<math.h> using namespace std; const int mod=1000000007; long long n,a[501][1001]; long long x,y; inline int read() { int t=0; char ch=getchar(); while(ch<'0'||ch>'9') ch=getchar(); while(ch>='0'||ch<='9') { t=t*10+(ch-'0'); ch=getchar(); } } int inv(long long a,long long b) { long long ans=1; while(b) { if(b&1) ans=ans*a%mod; a=a*a%mod; b>>=1; } return ans%mod; } int main() { scanf("%lld",&n); for(int i=1;i<=n;i++) { for(int j=1;j<=n;j++) scanf("%lld",&a[i][j]); a[i][i+n]=1; //在输入矩阵的同时在右边罗列一个n*n的单位矩阵 } for(int i=1;i<=n;i++) //进行高斯消元 { for(int j=i;j<=n;j++) { if(a[j][i]) { for(int q=1;q<=2*n;q++) swap(a[j][q],a[i][q]); break; } } if(!a[i][i]) //判无解 { cout<<"No Solution";return 0; } long long k=inv(a[i][i],mod-2)%mod; //利用费马小定理来求逆元 for(int j=1;j<=2*n;j++) a[i][j]=a[i][j]*k%mod; //利用矩阵性质将a[i][i]消成1,注意同样对右半边的单位矩阵操作 for(int j=1;j<=n;j++) //将第i列的其他行消成0 { if(j!=i) { long long k=a[j][i]; for(int m=i;m<=2*n;m++) //注意同时对右半边的单位矩阵进行操作 { a[j][m]=a[j][m]-k*a[i][m]; a[j][m]=(a[j][m]%mod+mod)%mod; } } } } for(int i=1;i<=n;i++) { for(int j=1+n;j<=2*n;j++) //输出右半边的矩阵就是逆矩阵啦 cout<<a[i][j]<<" "; cout<<endl; } return 0; //完结撒花 }
话说真的这道题就是两个模板题的结合qwq~