【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
",ans);
}
int main() {
#ifdef ivorysi
    freopen("f1.in","r",stdin);
#endif
    Init();
    Solve();
    return 0;
}