传说中 \(O(n^{2.807})\) 的矩阵乘法模板
实测大数据会比 \(O(n^3)\) 的矩乘快一点，卡常卡不过可以拿来用，并且博主写的比较丑，所以如果优化一下写法大概能跑得更快。
Code
#include<bits/stdc++.h>
#define ri register int
using namespace std;
typedef long long ll;
typedef vector<int> poly;
typedef pair<int,int> pii;
typedef pair<ll,int> pli;
typedef pair<int,ll> pil;
typedef pair<ll,ll> pll;
#define fi first
#define se second
#define pb push_back
#define rez resize
const ll Inf=1e18;
const int rlen=1<<20;
char buf[rlen],*ib=buf,*ob=buf;
#define gc() (((ib==ob)&&(ob=(ib=buf)+fread(buf,1,rlen,stdin)),ib==ob)?-1:*ib++)
inline int read() {
static int ans;
static char ch;
ans=0,ch=gc();
while(!isdigit(ch)) ch=gc();
while(isdigit(ch)) ans=(ans<<3)+(ans<<1)+(ch^48),ch=gc();
return ans;
}
inline ll readl() {
static ll ans;
static char ch;
ans=0,ch=gc();
while(!isdigit(ch)) ch=gc();
while(isdigit(ch)) ans=(ans<<3)+(ans<<1)+(ch^48),ch=gc();
return ans;
}
namespace modular {
const int mod=1e9+7;
inline int add(int a,int b) { return a+b<mod?a+b:a+b-mod; }
inline int dec(int a,int b) { return a<b?a-b+mod:a-b; }
inline int mul(int a,int b) { return (ll)a*b%mod; }
inline void Add(int&a,int b) { a=a+b<mod?a+b:a+b-mod; }
inline void Dec(int&a,int b) { a=a<b?a-b+mod:a-b; }
inline void Mul(int&a,int b) { a=(ll)a*b%mod; }
inline int ksm(int a,int p) { int ret=1;for(;p;p>>=1,Mul(a,a)) (p&1)&&(Mul(ret,a),1);return ret; }
inline int Inv(int a) { return ksm(a,mod-2); }
inline int sqr(int a) { return (ll)a*a%mod; }
inline int cub(int a) { return (ll)a*a%mod*a%mod; }
} using namespace modular;
template<typename T> inline void ckmax(T&a,T b) { a<b?a=b:0; }
template<typename T> inline void ckmin(T&a,T b) { a>b?a=b:0; }
struct mat {
int a[512][512],l;
mat(int _l=0) { l=_l;for(ri i=0;i<l;++i) for(ri j=0;j<l;++j) a[i][j]=0; }
inline int*operator[](const int&k) { return a[k]; }
inline mat operator*(mat b) {
mat c(l);
if(l<=64) {
for(ri i=0;i<l;++i) for(ri k=0;k<l;++k) if(a[i][k])
for(ri j=0;j<l;++j) Add(c[i][j],mul(a[i][k],b[k][j]));
}
else {
l>>=1;
mat aa(l),bb(l),m1(l),m2(l),m3(l),m4(l),m5(l),m6(l),m7(l);
//m1=(a11+a22)*(b11+b22);
for(ri i=0;i<l;++i) for(ri j=0;j<l;++j) {
aa[i][j]=add(a[i][j],a[i+l][j+l]);
bb[i][j]=add(b[i][j],b[i+l][j+l]);
} m1=aa*bb;
//m2=(a21+a22)*b11;
for(ri i=0;i<l;++i) for(ri j=0;j<l;++j) {
aa[i][j]=add(a[i+l][j],a[i+l][j+l]);
bb[i][j]=b[i][j];
} m2=aa*bb;
//m3=a11*(b12-b22);
for(ri i=0;i<l;++i) for(ri j=0;j<l;++j) {
aa[i][j]=a[i][j];
bb[i][j]=dec(b[i][j+l],b[i+l][j+l]);
} m3=aa*bb;
//m4=a22*(b21-b11);
for(ri i=0;i<l;++i) for(ri j=0;j<l;++j) {
aa[i][j]=a[i+l][j+l];
bb[i][j]=dec(b[i+l][j],b[i][j]);
} m4=aa*bb;
//m5=(a11+a12)*b22;
for(ri i=0;i<l;++i) for(ri j=0;j<l;++j) {
aa[i][j]=add(a[i][j],a[i][j+l]);
bb[i][j]=b[i+l][j+l];
} m5=aa*bb;
//m6=(a21-a11)*(b11+b12);
for(ri i=0;i<l;++i) for(ri j=0;j<l;++j) {
aa[i][j]=dec(a[i+l][j],a[i][j]);
bb[i][j]=add(b[i][j],b[i][j+l]);
} m6=aa*bb;
//m7=(a12-a22)*(b21+b22);
for(ri i=0;i<l;++i) for(ri j=0;j<l;++j) {
aa[i][j]=dec(a[i][j+l],a[i+l][j+l]);
bb[i][j]=add(b[i+l][j],b[i+l][j+l]);
} m7=aa*bb;
// c11=m1+m4-m5+m7
// c12=m3+m5
// c21=m2+m4
// c22=m1-m2+m3+m6
for(ri i=0;i<l;++i) for(ri j=0;j<l;++j) {
if((c[i][j]=((ll)m1[i][j]+m4[i][j]-m5[i][j]+m7[i][j])%mod)<0) c[i][j]+=mod;
c[i+l][j]=add(m2[i][j],m4[i][j]);
c[i][j+l]=add(m3[i][j],m5[i][j]);
if((c[i+l][j+l]=((ll)m1[i][j]-m2[i][j]+m3[i][j]+m6[i][j])%mod)<0) c[i+l][j+l]+=mod;
}
l<<=1;
}
return c;
}
inline void print() { for(ri i=0;i<l;++i,puts("")) for(ri j=0;j<l;++j) cout<<a[i][j]<<‘ ‘; }
};
int main() {
mat a(512);
for(ri i=0;i<512;++i) for(ri j=0;j<512;++j) a[i][j]=i*256+j;
a=a*a;
cout<<a[0][0];
return 0;
}
原文地址：https://www.cnblogs.com/ldxcaicai/p/12701602.html
时间： 2025-01-18 05:44:19