zoukankan      html  css  js  c++  java
  • P4783 【模板】矩阵求逆

    原题链接 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~

  • 相关阅读:
    违反了引用完整性约束。Dependent Role 具有多个具有不同值的主体。S级乌龙,自己制造的笑话
    用MVC5+EF6+WebApi 做一个小功能(二) 项目需求整理
    用MVC5+EF6+WebApi 做一个小功能(四) 项目分层功能以及文件夹命名
    用MVC5+EF6+WebApi 做一个小功能(三) 项目搭建
    ASP.NET WebApi总结之自定义权限验证
    用MVC5+EF6+WebApi 做一个小功能(一)开场挖坑,在线答题系统
    Javascript 535种方式!!!实现页面重载
    MVC页面移除HTTP Header中服务器信息
    为什么JavaScript要有null?(翻译)
    可编程渲染管线与着色器语言
  • 原文地址:https://www.cnblogs.com/xcg123/p/10700507.html
Copyright © 2011-2022 走看看