哪位大佬能帮我看一下，这个代码为什么会出现寻址错误的问题。调用其他函数都没有问题，但是调用Power和LU_Decomposition这两个函数就会寻址错误。

#include <iostream>
#include <cmath>
#include <vector>
#include <algorithm>
#define equal(a,b) (fabs(a-b)<eps)
using namespace std;
static const int N = 505;
const double eps = 1e-8;
struct Matrix {
    int n, m;
    double A[N][N];
    Matrix(int _n = 0, int _m = 0): n(_n), m(_m){
        memset(A,0,sizeof(A));
    }
    Matrix& operator = (const Matrix& mat){
        n = mat.n;
        m = mat.m;
        for(int i = 1; i <= n; i++){
            for(int j = 1; j <= m; j++){
                A[i][j] = mat.A[i][j];
            }
        }
        return *this;
    }
};

Matrix operator + (const Matrix& mat1, const Matrix& mat2){
    int n = mat1.n, m = mat1.m;
    Matrix res(n,m);
    if( n != mat2.n || m != mat2.m ) throw domain_error("The Matrices with different shapes can't be added.");
    for(int i = 1; i <= n; i++){
        for(int j = 1; j <= m; j++){
            res.A[i][j] = mat1.A[i][j] + mat2.A[i][j];
        }
    }
    return res;
}

Matrix operator * (const Matrix& mat1, const Matrix& mat2){
    int n = mat1.n, k1 = mat1.m , k2 = mat2.n , m = mat2.m;
    if( k1 != k2 ) throw domain_error("Two matrix with different shapes can not be added.");
    Matrix res(n,m);
    for(int i = 1; i <= n; i++){
        for(int j = 1; j <= m; j++){
            for(int k = 1; k <= k1; k++){
                res.A[i][j] += mat1.A[i][k]*mat2.A[k][j];
            }
        }
    }
    return res;
}

double Trace( const Matrix& mat){
    int n = mat.n;
    if( n != mat.m ){
        throw domain_error("Only square matrix has trace.");
    }
    double res = 0.0;
    for(int i = 1; i <= n; i++){
        res += mat.A[i][i];
    }
    return res;
}

void Gauss( Matrix& mat){
    if( mat.n != mat.m ) throw domain_error("Gauss Elimination can only be implemented on square matrix.");
    int n = mat.n;
    for(int i = 1, r; i <= n; i++){
        r = i;
        for(int j = i+1; j <= n; j++){
            if( fabs(mat.A[j][i]) > fabs(mat.A[r][i]) ){
                r = j;
            }
        }
        if( i != r ){
            swap(mat.A[i],mat.A[r]);
        }
        if( equal(mat.A[i][i],0.0) ){
            throw domain_error("The matrix only with full rank can be Gauss eliminated.");
        }
        for(int j = i+1; j <= n; j++){
            if( i == j ) continue;
            double tmp = mat.A[j][i]/mat.A[i][i];
            for(int k = 1; k <= n; k++){
                mat.A[j][k] -= tmp*mat.A[i][k];
            }
        }
    }
}

double Determinant( const Matrix& mat){
    if( mat.n != mat.m ) throw domain_error("Only square matrix can have determinant.");
    Matrix tmp = mat;
    int n = mat.n;
    for(int i = 1, r; i <= n; i++){
        r = i;
        for(int j = i+1; j <= n; j++){
            if( fabs(tmp.A[i][j]) > fabs(tmp.A[r][j]) ){
                r = j;
            }
        }
        if( i != r ){
            swap( tmp.A[i], tmp.A[r] );
        }
        if( equal(tmp.A[i][i],0.0) ){
            throw domain_error("Only the matrix with full rank has determinant.");
        }
        for(int j = 1; j <= n; j++){
            if( i == j ) continue;
            double coe = tmp.A[j][i]/tmp.A[i][i];
            for(int k = 1; k <= n; k++){
                tmp.A[j][k] -= coe*tmp.A[i][k];
            }
        }
    }
    double res = 1.0;
    for(int i = 1; i <= n; i++){
        res *= tmp.A[i][i];
    }
    return res;
}

void Inverse(Matrix mat, Matrix& I){
    if( mat.n != mat.m ){
        throw domain_error("The matrix with two different dimension is not invertible.");
    }
    int n = mat.n;
    I.n = n;
    I.m = n;
    for(int i = 1; i <= n; i++){
        I.A[i][i] = 1.0;
    }
    for(int i = 1, r; i <= n; i++){
        r = i;
        for(int j = i+1; j <= n; j++){
            if( fabs(mat.A[j][i]) > fabs(mat.A[r][i]) ){
                r = j;
            }
        }
        if( i != r ){
            swap(mat.A[i],mat.A[r]);
        }
        if( equal(mat.A[i][i],0.0)){
            throw domain_error("This Matirx is not invertible.");
        }
        for(int j = 1; j <= n; j++){
            if( i == j ) continue;
            double tmp = mat.A[j][i]/mat.A[i][i];
            for(int k = 1; k <= n; k++){
                mat.A[j][k] -= tmp*mat.A[i][k];
                I.A[j][k] -= tmp*I.A[i][k];
            }
        }
    }
    for(int i = 1; i <= n; i++){
        for(int j = 1; j <= n; j++){
            I.A[i][j] /= mat.A[i][i];
        }
    }
}

Matrix Power(Matrix mat, int b){
    if( mat.n != mat.m ) throw domain_error("This operation can only be implemented on sqare matrix.");
    int n = mat.n;
    Matrix res(n,n);
    for(int i = 1; i <= n; i++ ) res.A[i][i] = 1;
    while( b ){
        if( b&1 ){
            res = res*mat;
        }
        b >>= 1;
        mat = mat*mat;
    }
    return res;
}

void LU_Decomposition(Matrix mat, Matrix& L_tri, Matrix& U_tri ){
    U_tri = mat;
    Gauss(U_tri);
    Matrix inv_U;
    Inverse(U_tri,inv_U);
    L_tri = mat*inv_U;
}

int main(){
    int n;
    cin >>n;
    Matrix mat(n,n);
    for(int i = 1; i <= n; i++){//输入矩阵
        for(int j = 1; j <= n; j++){
            double x;
            cin >> x;
            mat.A[i][j] = x;
        }
    }
    Matrix mat1(n,n), mat2(n,n);
    Matrix res = Power(mat,2);
    LU_Decomposition(mat,mat1,mat2);
    for(int i = 1; i <= n; i++){
        for(int j = 1; j <= n; j++){
            cout << mat1.A[i][j] <<" ";
        }cout << endl;
    }
    for(int i = 1; i <= n; i++){
        for(int j = 1; j <= n; j++){
            cout << mat2.A[i][j] <<" ";
        }cout << endl;
    }
    for(int i = 1; i <= n; i++){
        for(int j = 1; j <= n; j++){
            cout << res.A[i][j] <<" ";
        }
    }
    return 0;
}