【Codechef】Random Number Generator(多项式除法)
题解
前置技能
1.多项式求逆
求\(f(x)\*g(x) \equiv 1 \pmod {x^{t}}\)
我们在t == 1时,有\(f[0] = frac{1}{g[0]}\)
之后呢,我们倍增一下,假如新的答案是\(g'(x)\)在\(\pmod {x^{2t}}\)意义下,显然有
\(g'(x) - g(x) \equiv 0 \pmod {x^{t}}\)
我们两边平方一下
\(g'^{2}(x) - 2g'(x)g(x) + g^{2}(x) \equiv 0 \pmod {x^{2t}}\)
两边再乘上一个\(f(x)\)
\(g'(x) - 2g(x) + f(x)g^{2}(x) \equiv 0 \pmod {x^{2t}}\)
那么答案就是
\(g'(x) \equiv = g(x)(2 - f(x)g(x))\)
这是我们熟悉的卷积形式,可以直接上FFT
注意指数上界是e * 2因为有三个多项式相乘
2.多项式除法
\(f(x) = g(x)q(x) + r(x)\)
设\(f(x)\)的度数为\(n\),\(g(x)\)的度数为\(m\),\(q(x)\)的度数为\(n - m\)
\(r(x)\)的度数小于等于\(m - 1\)
\(f(\frac{1}{x}) = g(\frac{1}{x})q(\frac{1}{x}) + r(\frac{1}{x})\)
\(x^{n}f(\frac{1}{x}) = x^{n}g(\frac{1}{x})q(\frac{1}{x}) + x^{n}r(\frac{1}{x})\)
设\(R(f)\)为f的每项系数翻转过来
则式子可以变为
\(R(f) = R(g)R(q) + x^{n - m + 1}R(r)\)
再在两边取模\(x^{n - m + 1}\)
则\(R(q) \equiv \frac{R(f)}{R(g)} \pmod {x^{n - m + 1}}\)
3.多项式取模
做完多项式除法后,用减掉商和被除数的乘积
再来看这道题
我们发现,这个式子生成的方式有点难受,我们就把系数翻转过来
我们用指数大小代表是序列的第几项
然后呢,我们发现其实就是如果生成了一个\(x^{k + 1}\)我们要把它分解成那个线性递推式即\(\sum_{i = 1}^{k}a_{i}x^{i}\)
设\(g(x) = x^{k + 1} - \sum_{i = 1}^{k}a_{i}x^{i}\)
答案也就是\(x^{n} \pmod {g(x)}\)的结果
这样只要做多项式快速幂然后取模多项式就好了
代码
#include <iostream>
#include <cstdio>
#include <cstring>
#include <algorithm>
#include <cmath>
#include <ctime>
#include <vector>
#include <set>
//#define ivorysi
#define eps 1e-8
#define mo 974711
#define pb push_back
#define mp make_pair
#define pii pair<int,int>
#define fi first
#define se second
#define MAXN 30005
using namespace std;
typedef long long int64;
typedef unsigned int u32;
typedef double db;
const int MOD = 25 * (1 << 22) + 1,G = 3,L = (1 << 22);
int W[L + 5],F[MAXN],K;
int64 N;
template<class T>
void read(T &res) {
res = 0;char c = getchar();T f = 1;
while(c < '0' || c > '9') {
if(c == '-') f = -1;
c = getchar();
}
while(c >= '0' && c <= '9') {
res = res * 10 - '0' + c;
c = getchar();
}
}
template<class T>
void out(T x) {
if(x < 0) putchar('-');
while(x >= 10) {
out(x / 10);
}
putchar(x % 10 + '0');
}
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;
}
struct poly {
int a[65537],deg;
poly() {
memset(a,0,sizeof(a));
deg = 0;
}
friend void NTT(poly &f,int on,int M) {
for(int i = 1 , j = M / 2 ; i < M - 1 ; ++i) {
if(i < j) swap(f.a[i],f.a[j]);
int k = M / 2;
while(j >= k) {
j -= k;
k >>= 1;
}
j += k;
}
for(int h = 2 ; h <= M ; h <<= 1) {
int wn = W[(L + on * L / h) % L];
for(int k = 0 ; k < M ; k += h) {
int w = 1;
for(int j = k ; j < k + h / 2 ; ++j) {
int u = f.a[j];
int t = 1LL * f.a[j + h / 2] * w % MOD;
f.a[j] = (u + t) % MOD;
f.a[j + h / 2] = (u + MOD - t) % MOD;
w = 1LL * wn * w % MOD;
}
}
}
if(on < 0) {
int Inv = fpow(M,MOD - 2);
for(int i = 0 ; i < M ; ++i) f.a[i] = 1LL * f.a[i] * Inv % MOD;
}
}
friend poly operator + (const poly &f,const poly &g) {
poly h;
int t = max(f.deg,g.deg);
for(int i = 0 ; i <= t ; ++i) {
h.a[i] = (f.a[i] + g.a[i]) % MOD;
}
for(int i = t ; i >= 0 ; --i) {
if(h.a[i] != 0) {h.deg = i;break;}
}
return h;
}
friend poly operator - (const poly &f,const poly &g) {
poly h;
int t = max(f.deg,g.deg);
for(int i = 0 ; i <= t ; ++i) {
h.a[i] = (f.a[i] + MOD - g.a[i]) % MOD;
}
for(int i = t ; i >= 0 ; --i) {
if(h.a[i] != 0) {h.deg = i;break;}
}
return h;
}
friend poly operator * (poly f,poly g) {
int M = 1;while(M <= f.deg + g.deg) M <<= 1;
NTT(f,1,M);NTT(g,1,M);
poly h;
for(int i = 0 ; i < M ; ++i) {
h.a[i] = 1LL * f.a[i] * g.a[i] % MOD;
}
NTT(h,-1,M);
for(int i = M - 1 ; i >= 0 ; --i) {
if(h.a[i] != 0) {
h.deg = i;
break;
}
}
return h;
}
friend void calc_Inverse(const poly &f,poly &g,int e) {
if(e == 1) {
g.a[0] = fpow(f.a[0],MOD - 2);
g.deg = 0;
return;
}
calc_Inverse(f,g,(e + 1) >> 1);
int M = 1;while(M <= g.deg * 2 + e - 1) M <<= 1;
poly t = f;t.deg = e - 1;
for(int i = e ; i < M ; ++i) t.a[i] = 0;
for(int i = g.deg + 1 ; i < M ; ++i) g.a[i] = 0;
NTT(g,1,M);NTT(t,1,M);
for(int i = 0 ; i < M ; ++i) {
g.a[i] = 1LL * g.a[i] * (2 + MOD - 1LL * g.a[i] * t.a[i] % MOD) % MOD;
}
NTT(g,-1,M);
for(int i = e ; i < M ; ++i) g.a[i] = 0;
g.deg = e - 1;
}
friend poly Inverse(const poly &f,int e) {
poly g;
calc_Inverse(f,g,e);
return g;
}
friend poly operator / (const poly &f,const poly &g) {
poly q;
if(f.deg < g.deg) return q;
poly Rf = f,Rg = g;
reverse(Rf.a,Rf.a + Rf.deg + 1);
reverse(Rg.a,Rg.a + Rg.deg + 1);
int t = f.deg - g.deg + 1;
for(int i = t ; i <= Rg.deg ; ++i) Rg.a[i] = 0;
for(int i = t ; i <= Rf.deg ; ++i) Rf.a[i] = 0;
Rf.deg = min(Rf.deg,t - 1);
Rg.deg = min(Rg.deg,t - 1);
poly Rq = Rf * Inverse(Rg,t);
q = Rq;
for(int i = t ; i <= q.deg ; ++i) q.a[i] = 0;
reverse(q.a,q.a + t);
q.deg = t - 1;
return q;
}
friend poly operator % (const poly &f,const poly &g) {
if(f.deg < g.deg) return f;
poly q = f / g;
return f - (g * q);
}
friend poly fpow(const poly &x,int64 c,const poly &Mod) {
poly res,t = x;
res.deg = 0;res.a[0] = 1;
while(c) {
if(c & 1) res = res * t % Mod;
t = t * t % Mod;
c >>= 1;
}
return res;
}
}M,X;
void Init() {
W[0] = 1,W[1] = fpow(G,25);
for(int i = 2 ; i < L ; ++i) W[i] = 1LL * W[i - 1] * W[1] % MOD;
read(K);read(N);
for(int i = 1 ; i <= K ; ++i) read(F[i]);
int a;
M.deg = K + 1;
for(int i = 1 ; i <= K ; ++i) {
read(a);
M.a[K - i + 1]= MOD - a;
}
M.a[K + 1] = 1;
X.deg = 1;X.a[1] = 1;
}
void Solve() {
X = fpow(X,N,M);
int64 ans = 0;
for(int i = 1 ; i <= K ; ++i) {
ans = ans + 1LL * X.a[i] * F[i] % MOD;
}
ans %= MOD;
printf("%lld\n",ans);
}
int main() {
#ifdef ivorysi
freopen("f1.in","r",stdin);
#endif
Init();
Solve();
return 0;
}