高斯消元+矩阵快速幂

主元素+若尔当

#include<iostream>
#include<cstdio>
#include<cstring>
#include<algorithm>
#include<cmath>
using namespace std;
const int N=100;
double a[N][N];
double ans[N];
int n;
void swap_r(int q,int p)
{
    for(int i=1;i<=n+1;i++)
    {
        double t=a[p][i];
        a[p][i]=a[q][i];
        a[q][i]=t;
    }
}
void swap_c(int q,int p)
{
    for(int i=1;i<=n+1;i++)
    {
        double t=a[i][p];
        a[i][p]=a[i][q];
        a[i][q]=t;
    }
}
void prt()
{
    for(int i=1;i<=n+1;i++)
    {
        for(int j=1;j<=n+1;j++)
            printf("%2.5f ",a[i][j]);
        printf("\n");
    }
    printf("\n");
}
void gs()
{
    for(int i=1;i<n;i++)
    {
        double m=fabs(a[i][i]);
        int p=i,q=i;
        for(int j=i+1;j<=n;j++)
            for(int k=i;k<=n;k++)
            if(fabs(a[j][k])>m)
            {
                m=fabs(a[j][k]);
                p=j;
                q=k;
            }
            //cout<<p<<endl;
        if(p!=i)
            swap_r(p,i);
        if(q!=i)
            swap_c(q,i);
        //prt();
        for(int j=i+1;j<=n;j++)
        {
            if(a[j][i]==0.0) continue;
            double t=a[j][i]/a[i][i];
            a[j][i]=0.0;
            for(int k=i+1;k<=n+1;k++)
                a[j][k]-=t*a[i][k];
        }
        //prt();
    }
    prt();
}
void red()
{
    for(int i=n;i>1;i--)
    {
        for(int j=i-1;j>=1;j--)
        {
            if(a[j][i]==0) continue;
            double t=a[j][i]/a[i][i];
            a[j][i]=0.0;
            for(int k=i+1;k<=n+1;k++)
                a[j][k]-=t*a[i][k];
        }
    }
    prt();
}
int main()
{
    scanf("%d",&n);
    for(int i=1;i<=n;i++)
        for(int j=1;j<=n+1;j++)
            scanf("%lf",&a[i][j]);
    for(int i=1;i<=n;i++)
        a[n+1][i]=i;
    gs();
    red();
    for(int i=1;i<=n;i++)
        ans[(int)a[n+1][i]]=a[i][n+1]/a[i][i];
    for(int i=1;i<=n;i++)
        printf("x%d=%f\n",i,ans[i]);
    return 0;
}

求出A,B,C,然后套用这个，带入矩阵快速幂即可算的答案

#include<bits/stdc++.h>

typedef long long ll;
typedef long double ld;
typedef unsigned long long ull;
#define MAXN 0x3f3f3f3f3f3f3f3f

using namespace std;

const ll MOD=1e9+7;
const int N=4;

struct node{
    ll a[10][10];
}tmp,ans,t;

//#define debug

node matrix(node x,node y ){
    node q;
    for(int i=1;i<=N;i++){
        for(int j=1;j<=N;j++){
            q.a[i][j]=0;
            for(int k=1;k<=N;k++){
                q.a[i][j]=(q.a[i][j]+x.a[i][k]*y.a[k][j]+MOD)%MOD;
            }
        }
    }
    return q;
}

void quick_ma(ll n){
    for(int i=1;i<=N;i++){
        for(int j=1;j<=N;j++){
            ans.a[i][j]=0;
        }
    }
    for(int i=1;i<=N;i++) ans.a[i][i]=1;
    t=tmp;
    while(n){
        if(n&1) ans=matrix(ans,t);
        n>>=1;
        t=matrix(t,t);
    }
}

int main()
{
    ll n;
    for(int i=1;i<=N;i++){
        for(int j=1;j<=N;j++){
            tmp.a[i][j]=0;
        }
    }
    tmp.a[1][1]=1;
    tmp.a[1][2]=5;
    tmp.a[1][3]=1;
    tmp.a[1][4]=-1;
    tmp.a[2][1]=1;
    tmp.a[3][2]=1;
    tmp.a[4][3]=1;
    while(scanf("%lld",&n)!=EOF){
        if(n==1) printf("1\n");
        else if(n==2) printf("5\n");
        else if(n==3) printf("11\n");
        else if(n==4) printf("36\n");
        else{
            quick_ma(n-4);
            printf("%lld\n",(ans.a[1][1]*36l%MOD+ans.a[1][2]*11l%MOD+ans.a[1][3]*5l%MOD+ans.a[1][4])%MOD);
        }
    }
    return 0;
}