假设$C=AB$, 那么答案就为
$egin{align} otag ans & =sumlimits_{i=0}^{n-1}sumlimits_{j=0}^{n-1}C[i][j]p^{(n-i)n-1-j} \ otag & = sumlimits_{i=0}^{n-1}sumlimits_{j=0}^{n-1}sumlimits_{k=0}^{n-1}A[i][k]B[k][j]p^{(n-i)n-1-j} \ & = sumlimits_{k=0}^{n-1}Big(sumlimits_{i=0}^{n-1}A[i][k]p^{(n-i)n-1}Big)Big(sumlimits_{j=0}^{n-1}B[k][j]p^{-j}Big) otag end{align}$
#include <iostream> #include <sstream> #include <algorithm> #include <cstdio> #include <math.h> #include <set> #include <map> #include <queue> #include <string> #include <string.h> #include <bitset> #define REP(i,a,n) for(int i=a;i<=n;++i) #define PER(i,a,n) for(int i=n;i>=a;--i) #define hr putchar(10) #define pb push_back #define lc (o<<1) #define rc (lc|1) #define mid ((l+r)>>1) #define ls lc,l,mid #define rs rc,mid+1,r #define x first #define y second #define io std::ios::sync_with_stdio(false) #define endl ' ' #define DB(a) ({REP(__i,1,n) cout<<a[__i]<<' ';hr;}) using namespace std; typedef long long ll; typedef pair<int,int> pii; const int P = 1e9+7, P2 = 998244353, INF = 0x3f3f3f3f; ll gcd(ll a,ll b) {return b?gcd(b,a%b):a;} ll qpow(ll a,ll n) {ll r=1%P;for (a%=P;n;a=a*a%P,n>>=1)if(n&1)r=r*a%P;return r;} ll inv(ll x){return x<=1?1:inv(P%x)*(P-P/x)%P;} inline int rd() {int x=0;char p=getchar();while(p<'0'||p>'9')p=getchar();while(p>='0'&&p<='9')x=x*10+p-'0',p=getchar();return x;} //head int n, Aa, Ab, Ac, Ad, Ba, Bb, Bc, Bd, p; uint32_t x, y, z, w; uint32_t xorshift() { uint32_t t = x; t ^= t << 11; t ^= t >> 8; x = y; y = z; z = w; w ^= w >> 19; w ^= t; return w & ((1 << 24) - 1); } void get(uint32_t a, uint32_t b, uint32_t c, uint32_t d) { x = a; y = b; z = c; w = d; } const int N = 7e3+10; int f1[N], f2[N], g1[N], g2[N]; int main() { scanf("%d%d%d%d%d%d%d%d%d%d", &n, &Aa, &Ab, &Ac, &Ad, &Ba, &Bb, &Bc, &Bd, &p); int w1 = inv(qpow(p,n)), w2 = inv(p); int r1 = qpow(p,n*n-1), r2 = 1; REP(i,0,n-1) { f1[i]=r1,r1=(ll)r1*w1%P; f2[i]=r2,r2=(ll)r2*w2%P; } get(Aa,Ab,Ac,Ad); REP(i,0,n-1) REP(j,0,n-1) { g1[j] = (g1[j]+(ll)xorshift()*f1[i])%P; } get(Ba,Bb,Bc,Bd); REP(i,0,n-1) REP(j,0,n-1) { g2[i] = (g2[i]+(ll)xorshift()*f2[j])%P; } int ans = 0; REP(i,0,n-1) ans = (ans+(ll)g1[i]*g2[i])%P; if (ans<0) ans += P; printf("%d ", ans); }