【Codechef】Chef and Bike(二维多项式插值)
something wrong with my new blog!
I can't type matrixs
so I come back. qwq
题目:https://www.codechef.com/problems/BIKE
题解
是我naive了,二维和一维其实差不多
首先,n很小,t很大,什么算法?矩阵乘法!没跑了
然后矩阵里填什么?一条边是两个值啊,还要一个%n一个%(n - 1),怎么搞
我们设计一个多项式xayb,x指数(也就是a)代表前轮加上一条边的值后取模的值,y指数(也就是b)代表后轮加上一条边的值后取模的值
最后每一项前面的系数就是答案
这样我们矩阵里填的就是多项式了,例如一条边的值是f=1,r=2那么这个邻接矩阵里填的就是xy2
好了我们做矩阵乘法并且指数分别取模n和(n - 1)就做完啦!(大雾
啥啥啥,为啥这样就做完了啊,我要在矩阵里填多项式,还要资瓷多项式的加和乘,虽然可行?但是复杂度听起来不可靠(项有n2级别啊,我多项式乘暴力搞要n4我FFT还要n2logn,而我矩阵乘法要做n3次,excuse me?)而且很玄学(我不会写),有没有什么简单优美的办法啊。
考虑插值,我们熟悉的插值都有什么呢,例如FFT,例如拉格朗日插值,都是找一些点值代进去然后反解出系数
考虑继续用FFT的思路插值反解系数,因为单位根有一个很棒的性质就是wn=1(先记住这一点,不过可能大家都知道)
啥,你说这模数不是NTT?不是哇你不是N只有22么,你会发现P−1是2..22所有的数的倍数,所以求出原根来就可以求单位根啦,不是2k而是正好N的(所以我们不能FFT了而要暴力DFT,也就是把每个单位根的值代进去求一遍值,显然你n2和大常数nlogn在22没有什么太大的区别,况且你也不能FFT)
那么我们就可以尝试把x的N种取值,y的N−1种取值结合起来
例如xy2就是wNw2N−1
好惹现在窝们(奇怪的口音?)有了N(N−1)个点值了……我们还要画在三维坐标系里不成?
回忆一下学习FFT时候的矩阵
大概是这个样子的
[(ω0n)0(ω0n)1⋯(ω0n)n−1(ω1n)0(ω1n)1⋯(ω1n)n−1⋮⋮⋱⋮(ωn−1n)0(ωn−1n)1⋯(ωn−1n)n−1][a0a1⋮an−1]=[A(ω0)A(ω1)⋮A(ωn−1)]
感觉没有问题,那么回忆一下是怎么反解出来的
假如有这么一个矩阵
[(ω−0n)0(ω−0n)1⋯(ω−0n)n−1(ω−1n)0(ω−1n)1⋯(ω−1n)n−1⋮⋮⋱⋮(ω−(n−1)n)0(ω−(n−1)n)1⋯(ω−(n−1)n)n−1]
我们让两个矩阵相乘
[(ω0n)0(ω0n)1⋯(ω0n)n−1(ω1n)0(ω1n)1⋯(ω1n)n−1⋮⋮⋱⋮(ωn−1n)0(ωn−1n)1⋯(ωn−1n)n−1][(ω−0n)0(ω−0n)1⋯(ω−0n)n−1(ω−1n)0(ω−1n)1⋯(ω−1n)n−1⋮⋮⋱⋮(ω−(n−1)n)0(ω−(n−1)n)1⋯(ω−(n−1)n)n−1]=E
设
D=[(ω0n)0(ω0n)1⋯(ω0n)n−1(ω1n)0(ω1n)1⋯(ω1n)n−1⋮⋮⋱⋮(ωn−1n)0(ωn−1n)1⋯(ωn−1n)n−1]
V=[(ω−0n)0(ω−0n)1⋯(ω−0n)n−1(ω−1n)0(ω−1n)1⋯(ω−1n)n−1⋮⋮⋱⋮(ω−(n−1)n)0(ω−(n−1)n)1⋯(ω−(n−1)n)n−1]
∑n−1i=0∑n−1j=0Eij=DikVkj
当i==j
Eij=n
当i!=j
Eij=∑n−1k=0DikVkj
Eij=∑n−1k=0ωik−kjn
Eij=∑n−1k=0ωk(i−j)n
Eij=1−(ωi−jn)n1−ωi−jn
根据单位根神奇的性质,我们有
(ωi−jn)n=1!
这样的话
E=nI
I是单位矩阵
所以就有了
1n[(ω0n)0(ω0n)1⋯(ω0n)n−1(ω1n)0(ω1n)1⋯(ω1n)n−1⋮⋮⋱⋮(ωn−1n)0(ωn−1n)1⋯(ωn−1n)n−1][A(ω0)A(ω1)⋮A(ωn−1)]=[a0a1⋮an−1]
妙趣横生,我们只需要把原来的FFT的单位根改成反着旋转,最后除一下数组长度即可
然而???你说你会,因为wys在WC2018讲过了。。。(大雾
其实二维的毕克在WC2017也讲了,显然这道题是BI KE出的。。。= =
我们考虑一个n(n−1)的矩阵!
[(ω0n)0(ω0n−1)0(ω0n)0(ω0n−1)1⋯(ω0n)n−1(ω0n−1)n−3(ω0n)n−1(ω0n−1)n−2(ω0n)0(ω1n−1)0(ω0n)0(ω1n−1)1⋯(ω0n)n−1(ω1n−1)n−3(ω0n)n−1(ω1n−1)n−2⋮⋮⋱⋮⋮(ωn−1n)0(ωn−3n−1)0(ωn−1n)0(ωn−3n−1)1⋯⋮⋮(ωn−1n)0(ωn−2n−1)0(ωn−1n)0(ωn−2n−1)1⋯(ωn−1n)n−1(ωn−2n−1)n−3(ωn−1n)n−1(ωn−2n−1)n−2][a0,0a0,1⋮an−1,n−3an−1,n−2]=[A(ω0nω0n−1)A(ω0nω1n−1)⋮A(ωn−1nωn−3n−1)A(ωn−1nωn−2n−1)]
按照一样的方法把号都取反
D=[(ω0n)0(ω0n−1)0(ω0n)0(ω0n−1)1⋯(ω0n)n−1(ω0n−1)n−3(ω0n)n−1(ω0n−1)n−2(ω0n)0(ω1n−1)0(ω0n)0(ω1n−1)1⋯(ω0n)n−1(ω1n−1)n−3(ω0n)n−1(ω1n−1)n−2⋮⋮⋱⋮⋮(ωn−1n)0(ωn−3n−1)0(ωn−1n)0(ωn−3n−1)1⋯⋮⋮(ωn−1n)0(ωn−2n−1)0(ωn−1n)0(ωn−2n−1)1⋯(ωn−1n)n−1(ωn−2n−1)n−3(ωn−1n)n−1(ωn−2n−1)n−2]
V=[(ω−0n)0(ω−0n−1)0(ω−0n)0(ω−0n−1)1⋯(ω−0n)n−1(ω−0n−1)n−3(ω−0n)n−1(ω−0n−1)n−2(ω−0n)0(ω−1n−1)0(ω−0n)0(ω−1n−1)1⋯(ω−0n)n−1(ω−1n−1)n−3(ω−0n)n−1(ω−1n−1)n−2⋮⋮⋱⋮⋮(ω−(n−1)n)0(ω−(n−3)n−1)0(ω−(n−1)n)0(ω−(n−3)n−1)1⋯⋮⋮(ω−(n−1)n)0(ω−(n−2)n−1)0(ω−(n−1)n)0(ω−(n−2)n−1)1⋯(ω−(n−1)n)n−1(ω−(n−2)n−1)n−3(ω−(n−1)n)n−1(ω−(n−2)n−1)n−2]
E=DV
∑n(n−1)−1i=0∑n(n−1)−1j=0Ei,j=Di,kVk,j
考虑把标号拆成一个点,类似一些的矩阵标号方法?
设数值为i拆成的点是ix,iy
当i==j时
Ei,j=(ωixn)kx(ωiyn−1)ky(ω−kxn)jx(ω−kyn−1)jy
Ei,j=n∗(n−1)
而当i!=j时
Ei,j=(ωixn)kx(ωiyn−1)ky(ω−kxn)jx(ω−kyn−1)jy
Ei,j=ωkx(ix−jx)nωky(iy−jy)n−1
嗯?迷一般的。。。眼熟?
可是这样就不是等比数列了?…………等等!
Ei,j=∑n−1kx=0ωkx(ix−jx)n∑n−2ky=0ωky(iy−jy)n−1
Ei,j=1−ωn(ix−jx)n1−ωix−jxn1−ω(n−1)(iy−jy)n−11−ωiy−jyn−1
Ei,j=0
也就是!我们把等比数列的公比q当成了ωix−jxn
后面就不写了,大家都知道了,最后模长除个n(n−1)就可以了!
说到这,我想我已经想到了一个绝妙的n2logn的二维FFT(在界是2^k的情况下)
我猜你们也想到了,大概二进制平摊反转这一个小技巧就不能用了,可以递归
代码
#include <iostream>
#include <cstdio>
#include <cstring>
#include <algorithm>
#define MAXM 100005
//#define ivorysi
using namespace std;
typedef long long int64;
const int MOD = 1163962801,G = 46;
int N,M,T;
int px[105],py[105];
struct node {int s,t,x,y;} E[MAXM];
int v[25][25][25][25],w[25][25];
struct Matrix {
int f[24][24];
void clear() {
memset(f,0,sizeof(f));
}
friend Matrix operator * (const Matrix &a,const Matrix &b) {
Matrix c;c.clear();
for(int i = 1 ; i <= N ; ++i) {
for(int j = 1 ; j <= N ; ++j) {
for(int k = 1 ; k <= N ; ++k) {
c.f[i][j] = (c.f[i][j] + 1LL * a.f[i][k] * b.f[k][j] % MOD) % MOD;
}
}
}
return c;
}
void set_Matrix(int x,int y) {
clear();
for(int i = 1 ; i <= M ; ++i) {
f[E[i].s][E[i].t] = (f[E[i].s][E[i].t] + 1LL * px[E[i].x * x % N] * py[E[i].y * y % (N - 1)] % MOD) % MOD;
}
}
}A,Ans,I;
int fpow(int x,int c) {
int res = 1,t = x;
while(c) {
if(c & 1) res = 1LL * res * t % MOD;
t = 1LL * t * t % MOD;
c >>= 1;
}
return res;
}
void mul(int c) {
Ans = I;Matrix t = A;
while(c) {
if(c & 1) Ans = Ans * t;
t = t * t;
c >>= 1;
}
}
void Process(int x,int y) {
for(int i = 0 ; i < N ; ++i) {
for(int j = 0 ; j < N - 1 ; ++j) {
w[i][j] = 0;
for(int k = 0 ; k < N ; ++k) {
for(int l = 0 ; l < N - 1 ; ++l) {
w[i][j] = (w[i][j] + 1LL * v[x][y][k][l] * px[N - k * i % N] % MOD
* py[N - 1 - l * j % (N - 1)] % MOD) % MOD;
}
}
w[i][j] = 1LL * w[i][j] * fpow(N * (N - 1),MOD - 2) % MOD;
}
}
}
void Init() {
scanf("%d%d%d",&N,&M,&T);
for(int i = 1 ; i <= M ; ++i) {
scanf("%d%d%d%d",&E[i].s,&E[i].t,&E[i].x,&E[i].y);
E[i].x %= N;
E[i].y %= N - 1;
}
px[0] = 1;py[0] = 1;
px[1] = fpow(G,(MOD - 1) / N);
py[1] = fpow(G,(MOD - 1) / (N - 1));
for(int i = 2 ;i <= N ; ++i) {
px[i] = 1LL * px[i - 1] * px[1] % MOD;
py[i] = 1LL * py[i - 1] * py[1] % MOD;
}
}
void Solve() {
Init();
I.clear();
for(int i = 1 ; i <= N ; ++i) I.f[i][i] = 1;
for(int i = 0 ; i < N ; ++i) {
for(int j = 0 ; j < N - 1 ; ++j) {
A.set_Matrix(i,j);
mul(T);
for(int k = 1 ; k <= N ; ++k) {
for(int l = 1 ; l <= N ; ++l) {
v[k][l][i][j] = Ans.f[k][l];
}
}
}
}
for(int i = 1 ; i <= N ; ++i) {
Process(i,i);
for(int k = 0 ; k < N ; ++k) {
for(int l = 0 ; l < N - 1 ; ++l) {
printf("%d%c",w[k][l]," \n"[l == N - 2]);
}
}
}
}
int main() {
#ifdef ivorysi
freopen("f1.in","r",stdin);
#endif
Solve();
}
【推荐】国内首个AI IDE,深度理解中文开发场景,立即下载体验Trae
【推荐】编程新体验,更懂你的AI,立即体验豆包MarsCode编程助手
【推荐】抖音旗下AI助手豆包,你的智能百科全书,全免费不限次数
【推荐】轻量又高性能的 SSH 工具 IShell:AI 加持,快人一步
· 如何编写易于单元测试的代码
· 10年+ .NET Coder 心语,封装的思维:从隐藏、稳定开始理解其本质意义
· .NET Core 中如何实现缓存的预热?
· 从 HTTP 原因短语缺失研究 HTTP/2 和 HTTP/3 的设计差异
· AI与.NET技术实操系列:向量存储与相似性搜索在 .NET 中的实现
· 10年+ .NET Coder 心语 ── 封装的思维:从隐藏、稳定开始理解其本质意义
· 地球OL攻略 —— 某应届生求职总结
· 提示词工程——AI应用必不可少的技术
· Open-Sora 2.0 重磅开源!
· 周边上新:园子的第一款马克杯温暖上架