整数行列式计算高斯方法无精度损失

最新推荐文章于 2022-12-02 09:34:39 发布

原创最新推荐文章于 2022-12-02 09:34:39 发布 · 696 阅读

0 ·

CC 4.0 BY-SA版权

文章标签：

#struct

学习笔记专栏收录该内容

14 篇文章

订阅专栏

本文介绍了一种使用高斯消元法计算整数行列式的方法，通过结构体表示矩阵，并使用长整型避免精度损失。程序中实现了矩阵的消元过程，并以分数形式存储中间结果，确保计算的准确性。

摘要生成于 C知道，由 DeepSeek-R1 满血版支持，前往体验 >

没用double，用的分数表示。不能算太大的。

#include <iostream>
#include <cstring>
#include <cstdio>
#include <cmath>
#define MAXN 50
#define MOD 1000000007
using namespace std;
struct mat{
       int size;
       long long mat[MAXN][MAXN];
       };

long long eu(long long a, long long b){
         if(b==0)
         return a;
         else
         return eu(b,a%b);
         }

struct xy{
       long long x,y;
       xy(long long a,long long b){
             x=a;
             y=b;
             }
       xy(void){};
       void init(){
            long long t=eu(x,y);
            x/=t;
            y/=t;
            x%=(MOD*y);
            }
       bool zero(){
            return x==0;
            }
       xy mul(xy other){
            xy ss;
            ss.x=x*other.x;
            ss.y=y*other.y;
            ss.init();
            return ss;
            }
       void del(xy other){
            x=x*other.y;
            y=y*other.x;
            init();
            }
       void sub(xy other){
            long long tx=x*other.y-other.x*y;
            long long ty=other.y*y;
            x=tx;
            y=ty;
            init();
            }
       };

//矩阵高斯消元法
long long det(const mat a){
     xy b[MAXN][MAXN];
     int i,j,k,sign=0;
     xy res(1,1);

     for( i=0; i<a.size; ++i)
     for( j=0; j<a.size; ++j)
     { b[i][j].x=a.mat[i][j]; b[i][j].y=1;}

     for( i=0; i<a.size; ++i){
             if(b[i][i].zero()){///////////////////////////////////////////////////////过会儿改进
                              for(j=i+1; j<a.size; ++j)
                              if(!b[j][i].zero()) break;
                              if(j==a.size)
                              return 0;

                              for(k=0; k<a.size; ++k)
                              {b[i][k].x^=b[j][k].x^=b[i][k].x^=b[j][k].x; b[i][k].y^=b[j][k].y^=b[i][k].y^=b[j][k].y;}
                              ++sign;
             }

             res=res.mul(b[i][i]);

             for(k=i+1; k<a.size; ++k)
             b[i][k].del(b[i][i]);

             for(j=i+1; j<a.size; ++j)
             for(k=i+1; k<a.size; ++k)
             b[j][k].sub(b[j][i].mul(b[i][k]));
             }
             if(sign&1)
             return -res.x/res.y;
             else
             return res.x/res.y;
     }
//
int mem[10000];
int main(){
    int n,lim,i,j;
    mat a,b;
    while(scanf("%d",&n)!=EOF){
    a.size=b.size=n;
    for(i=0; i<n; ++i)
    for(j=0; j<n; ++j)
    scanf("%lld",&a.mat[i][j]);
//    printf("%lld/n",a.mat[0][0]);
//    for(i=0; i<n; ++i)
//    for(j=0; j<n; ++j)
//    scanf("%d",&b.mat[i][j]);
    printf("%lld/n",det(a));

    }
    }