算法数学笔记-三、多项式
目录
三、多项式
拉格朗日线型插值
x[i] = i; (i = 1 ~ n)
long long ask(long long k){
int n = nn;
pre[0] = k;
fac[0] = suf[n + 1] = 1;
for(int j = 1; j <= n; j ++)
pre[j] = pre[j - 1] * (k - j) % mod, fac[j] = fac[j - 1] * j % mod;
for(int j = n; j >= 1; j --)
suf[j] = suf[j + 1] * (k - j) % mod;
long long re = y[0] * suf[1] % mod * inv(fac[n]) % mod * fuyi(n) % mod;
for(int i = 1; i <= n; i ++)
re = (re + y[i] * pre[i - 1] % mod * suf[i + 1] % mod * fuyi(n - i) % mod * inv(fac[i] * fac[n - i])) % mod;
return (re % mod + mod) % mod;
}
FFT（快速傅里叶变换）及FFT乘法
struct Complex{
double real, fals;
Complex(double a = 0, double b = 0){real = a; fals = b;}
Complex operator + (const Complex & a){ return Complex(real + a.real, fals + a.fals);}
Complex operator - (const Complex & a){ return Complex(real - a.real, fals - a.fals);}
Complex operator * (const Complex & a){ return Complex(real * a.real - fals * a.fals, real * a.fals + fals * a.real);}
Complex operator = (const double & a){ return Complex(a, 0); }
};
void FFT(vector<Complex> &A, int type = 1){
int limit = A.size();
for(int i = 0; i <= limit - 1;i ++)
if(i > rev[i]) swap(A[i], A[rev[i]]);
for(int half = 1; half <= limit - 1; half <<= 1){ //ling : half = 待合并区间的长度的一半
Complex wn(cos(Pi / half) , type * sin(Pi / half));
for(int L = 0; L <= limit - 1; L += (half << 1)){
Complex wn_k(1, 0);
for(int k = 0; k <= half - 1; k ++, wn_k = wn_k * wn){
Complex x = A[L + k], y = wn_k * A[L + half + k];
A[L + k] = x + y;
A[L + half + k] = x - y;
}
}
}
if(type == -1)// ling : type == -1 : 逆变换
for(int i = 0; i <= limit - 1; i ++)
A[i].real /= limit;
}
vector<double > MultiplyFFT(const vector<double> &a, const vector<double> &b){
vector<Complex> A(0), B(0);
for(int i = 0; i <= a.size() - 1; i ++) A.push_back(Complex(a[i], 0));
for(int i = 0; i <= b.size() - 1; i ++) B.push_back(Complex(b[i], 0));
int n = A.size() + B.size() - 1, limit = 1;
while(limit < (A.size() + B.size())) limit <<= 1;
for(int i = 0;i <= limit - 1;i ++)
rev[i] = (rev[i >> 1] >> 1) | ((i & 1) * limit >> 1);
A.resize(limit), FFT(A);
B.resize(limit), FFT(B);
for(int i = 0; i <= limit - 1; i ++)
A[i] = A[i] * B[i];
FFT(A, -1);
vector<double >re(n);
for(int i = 0; i <= n - 1; i ++)
re[i] = A[i].real; // ling : 如果要return整数的话，不要忘了[ + 0.5]
return re;
}
任意模数多项式乘法（MTT）
namespace MTT{
const int N = 5e5+10;
struct complex{
double x,y;
complex operator+(const complex &b)const{return (complex){x+b.x,y+b.y};}
complex operator-(const complex &b)const{return (complex){x-b.x,y-b.y};}
complex operator*(const complex &b)const{return (complex){x*b.x-y*b.y,x*b.y+y*b.x};}
complex operator*(const double &b)const{return (complex){x*b,y*b};}
complex operator~()const{return (complex){x,-y};}
}I,a [N], b[N], c[N], d[N], tmp1[N], tmp2[N], resw[N];
int f[N], g[N], rev[N];
int s;
inline void FFT(complex a[], int n, int ty){
for(int i = 0; i < n; i ++)
if(i < rev[i])
swap(a[i], a[rev[i]]);
resw[0] = (complex){1, 0};
for(int l = 1; l < n; l <<= 1){
complex wn = (complex){cos(pi / l), ty * sin(pi / l)};
for(int j = l - 2; j >= 0; j -= 2)
resw[j] = resw[j >> 1], resw[j | 1] = resw[j] * wn;
for(int i = 0; i < n; i += (l << 1))
for(int j = 0; j < l; j ++){
complex x = a[i + j], y = resw[j] * a[i + j + l];
a[i + j] = x + y; a[i + j + l] = x - y;
}
}
}
inline void work(int f[], int n, int s, int si, complex a[], complex b[]){
for(int i = 0; i < s; i ++)
tmp1[i].x = f[i] / si, tmp1[i].y = f[i] % si;
FFT(tmp1, s, 1);
for(int i = 0; i < s; i ++)
a[i] = (tmp1[i] + (~tmp1[i == 0 ? 0 : s - i])) * 0.5,
b[i] = ((~tmp1[i == 0 ? 0 : s - i]) - tmp1[i]) * 0.5 * I;
}
#define calc(x) (long long)(x / s + 0.5)
inline int get(int si, int s, int i, int p){
return (1ll * si * si * (calc(tmp1[i].x) % p)
+ 1ll * si * ((calc(tmp1[i].y) + calc(tmp2[i].x)) % p)
+ calc(tmp2[i].y)) % p;
}
inline vector<int> MTTMultiply(vector<int> aa, vector<int> bb, int mod){
I = (complex){0, 1};
int n = aa.size() - 1, m = bb.size() - 1;
for(int i = 0; i <= n; i ++) f[i] = aa[i] % mod;
for(int i = 0; i <= m; i ++) g[i] = bb[i] % mod;
s = 1;
int res = 0;
while(s <= n + m)
s <<= 1, res ++;
for(int i = 0; i < s; i ++)
rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << (res - 1));
int si = sqrt(mod);
work(f, n, s, si, a, b);
work(g, m, s, si, c, d);
for(int i = 0; i < s; i ++)
tmp1[i] = a[i] * c[i] + I * b[i] * c[i],
tmp2[i] = a[i] * d[i] + I * b[i] * d[i];
FFT(tmp1, s, -1);
FFT(tmp2, s, -1);
vector<int> ans;
for(int i = 0; i <= n + m; i ++)
ans.push_back(get(si, s, i, mod));
return ans;
}
}
集合幂级数
inline int quickpow(int x, int b, int mod){
int j = 1;
for(; b; b >>= 1, x = 1ll * x * x % mod)
if(b & 1)
j = 1ll * j * x % mod;
return j % mod;
}
namespace Subset{
typedef vector <int> vec;
// ling : !!!!! vector的size（）是ull，所以for的时候size（）不能 为0 , (ull)0 - 1 = inf
const int mod = 998244353;
vec tmpInv, Inv;
void out(vec x){
for(int i = 0; i <= x.size() - 1; i ++)
cout << x[i] << " ";
cout << endl;
}
void FWTor(vec &A, int type = 1){
int limit = A.size();
for(int len = 2, half = 1; half <= limit - 1; len <<= 1, half <<= 1)
for(int now = 0; now <= limit - 1; now += len)
for(int i = 0; i <= half - 1; i ++)
A[now | i | half] = (A[now | i | half] + 1ll * type * A[now | i] + mod) % mod;
}
void FWTand(vec &A, int type = 1){
int limit = A.size();
for(int len = 2, half = 1; half <= limit - 1; len <<= 1, half <<= 1)
for(int now = 0; now <= limit - 1; now += len)
for(int i = 0; i <= half - 1; i ++)
A[now | i] = (A[now | i] + 1ll * type * A[now | i | half] + mod) % mod;
}
void FWTxor(vec &A, int type = 1){
int limit = A.size();
for(int len = 2, half = 1; half <= limit - 1; len <<= 1, half <<= 1)
for(int now = 0; now <= limit - 1; now += len)
for(int i = 0; i <= half - 1; i ++){
int x = A[now | i], y = A[now | i | half];
A[now | i] = 1ll * type * (x + y) % mod;
A[now | i | half] = 1ll * type * (x - y + mod) % mod;
}
}
vec FwtMultiply(vec A, vec B, int type){ // ling : type = 1 or; 2 and; 3 xor
int limit = 1;
while(limit < max(A.size(), B.size())) limit <<= 1;
A.resize(limit), B.resize(limit);
if(type == 1) FWTor(A), FWTor(B);
else if(type == 2) FWTand(A), FWTand(B);
else if(type == 3) FWTxor(A), FWTxor(B);
for(int i = 0; i <= limit - 1; i ++)
A[i] = 1ll * A[i] * B[i] % mod;
if(type == 1) FWTor(A, -1);
else if(type == 2) FWTand(A, -1);
else if(type == 3) FWTxor(A, inv2);
A.resize(limit);
return A;
}
vec SubsetMultiply(int n, vec A, vec B){
int len = 1 << n;
A.resize(len);
B.resize(len);
vector<vec> aa(n + 1), bb(n + 1), cc(n + 1);
vec fra(n + 1), ena(n + 1), frb(n + 1), enb(n + 1);
for(int i = 0; i <= n; i ++)
aa[i].resize(len), bb[i].resize(len), cc[i].resize(len);
for(int i = 0; i < len; i ++)
aa[__builtin_popcount(i)][i] = A[i];
for(int i = 0; i < len; i ++)
bb[__builtin_popcount(i)][i] = B[i];
for(int i = 0; i <= n; i ++)
FWTor(aa[i]), FWTor(bb[i]);
for(int i = 0; i <= n; i ++){
fra[i] = frb[i] = (1 << n);
ena[i] = enb[i] = -1;
for(int S = 0; S < (1 << n); S ++){
if(aa[i][S]){
if(fra[i] == (1 << n))
fra[i] = S;
ena[i] = S;
}
if(bb[i][S]){
if(frb[i] == (1 << n))
frb[i] = S;
enb[i] = S;
}
}
}
for(int x = 0; x <= n; x ++){
for(int i = 0; i <= x; i ++)
for(int j = max(fra[i], frb[x - i]); j <= min(ena[i], enb[x - i]); j ++)
cc[x][j] = (cc[x][j] + 1ll * aa[i][j] * bb[x - i][j]) % mod;
FWTor(cc[x], -1);
}
vec ret(len);
for(int i = 0; i < len; i ++)
ret[i] = cc[__builtin_popcount(i)][i];
return ret;
}
inline vec times(int b, vec a){
for(int i = 0; i < a.size(); i ++)
a[i] = 1ll * a[i] * b % mod;
return a;
}
inline vec dec(vec a, vec b){
int n = a.size() - 1, m = b.size() - 1;
a.resize(max(n, m) + 1);
for(int i = 0; i <= m; i ++)
a[i] = (a[i] - b[i] + mod) % mod;
return a;
}
inline vec cut(vec a, int l){
if(a.size() <= l) return a;
a.resize(l);
return a;
}
inline vec mul(vec a, vec b, int l){
a = SubsetMultiply(__lg(l), a, b);
if(l != -1) a.resize(l);
return a;
}
vec GetInv(int m, vec a){
vec ans(1, quickpow(a[0], mod - 2, mod));
int si = 2;
while((si >> 1) < m){
ans = dec(times(2, ans), mul(mul(ans, ans, si), cut(a, si), si));
si <<= 1;
}
ans.resize(m);
return ans;
}
vec differ(vec A){ // ling: 求导
int n = A.size() - 1;
for(int i = 0; i <= n - 1; i ++)
A[i] = 1LL * A[i + 1] % mod * (i + 1) % mod;
A[n] = 0;
return A;
}
vec integ(vec A){ //ling : 多项式积分
int n = A.size() - 1;
A.push_back(0);
Inv.resize(n + 2);
Inv[1] = 1;
for(int i = 2; i <= n + 1; i ++)
Inv[i] = (long long)(mod - mod / i) * Inv[mod % i] % mod;
for(int i = n; i >= 0; i --)
A[i + 1] = 1LL * A[i] % mod * Inv[i + 1] % mod;
A[0] = 0;
return A;
}
vec GetLin(int n, vec F){
int len = 1 << n;
vector<vec> G(n + 1);
vec a(n + 2), b(n + 2), ret(1 << n), Inv(len + 2);
Inv[1] = 1;
for(int i = 2; i <= len + 1; i ++)
Inv[i] = (long long)(mod - mod / i) * Inv[mod % i] % mod;
for(int i = 0; i <= n; i ++)
G[i].resize(len);
for(int S = 1; S < (1 << n); S ++)
G[__builtin_popcount(S)][S] = F[S];
for(int i = 0; i <= n; i ++) FWTor(G[i]);
for(int S = 1; S < (1 << n); S ++){
for(int i = 0; i <= n; i ++)
a[i] = G[i][S];
for(int i = 0; i < n; i ++){
int res = 0;
for(int j = 0; j < i; j ++)
res = (res + 1ll * b[j] * a[i - j]) % mod;
b[i] = (1ll * a[i + 1] * (i + 1) + mod - res) % mod;
}
for(int i = 1; i <= n; i ++)
G[i][S] = 1ll * b[i - 1] * Inv[i] % mod;
}
for(int i = 0; i <= n; i ++)
FWTor(G[i], -1);
for(int S = 0; S < (1 << n); S ++)
ret[S] = G[__builtin_popcount(S)][S];
return ret;
}
}
\[\ \ FWTor[A]_i = \sum_{j | i = i}A_j \ \ \ \ \ i的子集和
\\
FWTand[A]_i = \sum_{j \& i = i}A_j \ \ \ \ \ i的超集和
\]
多项式全家桶
namespace poly{
#define LL long long
typedef vector <int> vec;
// ling : !!!!! vector的size（）是ull，所以for的时候size（）不能 为0 , (ull)0 - 1 = inf
const int N = 5e5 + 10, mod = 998244353, g = 3, invg = 332748118;
vec xi, yi, f, v, a, InterP, InterPxi;
vec Q[N], P[N];// ling : 多项式快速多点求值和快速插值
long long p = mod, I, out1, out2, inv2 = 499122177, realK;
vec aa, bb, tmpInv, tmpSqrt, F, G;
int rev[N], Inv[N];
int kp;
int gpow[33], invgpow[333];
bool NTTflag;
int GetInv(int x){
return ::quickpow(x, mod - 2, mod);
}
vec GetInvSum(vec a){
int n = a.size() - 1;
vec s(n + 1), sinv(n + 1);
s[0] = a[0];
for(int i = 1; i <= n; i ++)
s[i] = 1ll * s[i - 1] * a[i] % mod;
sinv[n] = GetInv(s[n]);
for(int i = n - 1; i >= 0; i --)
sinv[i] = 1ll * sinv[i + 1] * a[i + 1] % mod;
for(int i = 1; i <= n; i ++)
sinv[i] = 1ll * sinv[i] * s[i - 1] % mod;
return sinv;
}
vec PolyAdd(const vec &A, const vec &B){
vec ret = A;
ret.resize(max(A.size(), B.size()));
for(int i = 0; i <= B.size() - 1; i ++)
ret[i] = (ret[i] + B[i]) % mod;
return ret;
}
void NTT(vec &A, int type = 1){
if(NTTflag == false){
for(int i = 1; i <= 22; i ++){ // ling : 预处理：gpow[i] = p^((mod - 1) / (1 << i))
gpow[i] = quickpow(g, (mod - 1) / (1 << i), mod);
invgpow[i] = quickpow(invg, (mod - 1) / (1 << i), mod);
}
NTTflag = true;
}
int limit = A.size();
for(int i = 0;i <= limit - 1;i ++)
rev[i] = (rev[i >> 1] >> 1) | ((i & 1) * limit >> 1);
for(int i = 0; i <= limit - 1;i ++)
if(i > rev[i]) swap(A[i], A[rev[i]]);
for(int half = 1, num = 1; half <= limit - 1; half <<= 1, num ++){ //ling: half = 待合并区间的长度的一半
int gn = (type == 1) ? gpow[num] : invgpow[num];
for(int L = 0; L <= limit - 1; L += (half << 1)){
long long gn_k = 1;
for(int k = 0; k <= half - 1; k ++, gn_k = gn_k * gn % mod){
int x = A[L + k], y = gn_k * A[L + half + k] % mod;
A[L + k] = (x + y) % mod;
A[L + half + k] = ((x - y) % mod + mod) % mod;
}
}
}
if(type == -1){// ling : type == -1 : 逆变换
int inv = quickpow(limit, mod - 2, mod);
for(int i = 0; i <= limit - 1; i ++)
A[i] = (long long)A[i] * inv % mod;
}
}
void out(vec x){
for(int i = 0; i <= x.size() - 1; i ++)
cout << x[i] << " ";
cout << endl;
}
vec Multiply(vec A, vec B){
int n = A.size() + B.size() - 1, limit = 1;
while(limit < (A.size() + B.size())) limit <<= 1;
A.resize(limit), NTT(A);
B.resize(limit), NTT(B);
for(int i = 0; i <= limit - 1; i ++)
A[i] = 1ll * A[i] * B[i] % mod;
NTT(A, -1), A.resize(n);
return A;
}
vec GetInv(int Deg, vec A){
A.resize(Deg << 1);
vec ret(1, quickpow(A[0], mod - 2, mod));
for(int limit = 4; limit < (Deg << 2); limit <<= 1){
int deg = limit >> 1;
tmpInv.resize(limit), ret.resize(limit);
for(int i = 0; i <= deg - 1; i ++) //cout << limit << endl;
tmpInv[i] = A[i];
for(int i = deg; i <= limit - 1; i ++)
tmpInv[i] = 0;
NTT(tmpInv), NTT(ret);
for(int i = 0; i <= limit - 1; i ++)
ret[i] = (2 - 1LL * tmpInv[i] * ret[i] % mod) * ret[i] % mod;
NTT(ret, -1), ret.resize(deg);
}
ret.resize(Deg);
for(int i = 0; i <= Deg - 1; i ++)
if(ret[i] < 0) ret[i] += mod;
return ret;
}
vec differ(vec A){ // ling: 求导
int n = A.size() - 1;
for(int i = 0; i <= n - 1; i ++)
A[i] = 1LL * A[i + 1] % mod * (i + 1) % mod;
A[n] = 0;
return A;
}
vec integ(vec A){ //ling : 多项式积分
int n = A.size() - 1;
A.push_back(0);
Inv[1] = 1;
for(int i = 2; i <= n + 1; i ++)
Inv[i] = (long long)(mod - mod / i) * Inv[mod % i] % mod;
for(int i = n; i >= 0; i --)
A[i + 1] = 1LL * A[i] % mod * Inv[i + 1] % mod;
A[0] = 0;
return A;
}
vec GetLin(int deg, vec F){ //ling : 多项式取对数 （%n）
F.resize(deg);
F = Multiply(differ(F), GetInv(deg, F)), F.resize(deg - 1);
return integ(F);
}
vec GetSqrt(int Deg, vec A){ //ling : 多项式开方
A.resize(Deg << 1);
vec ret(1, Shengyu2::GetShengyu2(A[0], mod).first); // 得到A[0] 的二次剩余
for(int limit = 4; limit < (Deg << 2); limit <<= 1){
int deg = limit >> 1;
tmpSqrt.resize(deg);ret.resize(deg);
for(int i = 0; i <= deg - 1; i ++)
tmpSqrt[i] = A[i];//ling : tmp is A now
tmpSqrt = Multiply(tmpSqrt, GetInv(deg, ret)); // ling : tmp is A/ret now
for(int i = 0; i <= deg - 1; i ++)
ret[i] = (long long)(tmpSqrt[i] + ret[i]) * inv2 % mod;
ret.resize(deg);
}
ret.resize(Deg);
return ret;
}
vec GetExp(int Deg, vec A){ //ling : 多项式求指数
A.resize(Deg << 1);
vec ret = vec(1, 1);
for(int limit = 4; limit < (Deg << 2); limit <<= 1){
int deg = limit >> 1;
vec linret = GetLin(deg, ret);
for(int i = 0; i <= deg - 1; i ++)
linret[i] = (A[i] - linret[i] + mod) % mod;
linret[0] = (linret[0] + 1) % mod;
ret = Multiply(ret, linret);
ret.resize(deg);
}
ret.resize(Deg);
return ret;
}
vec opR(vec A){
for(int i = 0; i <= (A.size() - 1) >> 1; i ++)
swap(A[i], A[A.size() - 1 - i]);
return A;
}
pair<vec, vec> GetDiv(vec F, vec G){ //ling : 多项式整除和求余
pair<vec, vec> ret;
int n = F.size() - 1, m = G.size() - 1;
ret.first = Multiply(opR(F), GetInv(n - m + 1, opR(G)));
ret.first.resize(n - m + 1);
ret.first = opR(ret.first);
vec GQ = Multiply(G, ret.first);
GQ.resize(m), ret.second.resize(m);
for(int i = 0; i <= m - 1; i ++)
ret.second[i] = (F[i] - GQ[i] + mod) % mod;
return ret; // first: 商, second:余
}
struct big{
vector < int > num;
big(){
num.clear(); num.push_back(0);
}
void scan(){
num.clear();
char c = getchar();
while('0' > c || c > '9')
c = getchar();
while('0' <= c && c <= '9')
num.push_back(c - '0'), c = getchar();
for(int i = 0; i <= (num.size() - 1) >> 1; i ++)
swap(num[i], num[num.size() - 1 - i]);
}
void print(){
cout << num.back();
for(int i = num.size() - 2, x; i >= 0; i --)
putchar(num[i] + '0');
}
}K, bigK, bigKK;
vec GetPow(vec aa, big k, int n){
int km = 0, ke = 0;
{/*ling : realK 看做真正的k， 前档0 的个数就是realK * num_zore
这里的big 是高精数，若非则不需要
km 就是 realK % mod，即的多项式的指数
ke 就是 realK % (eular(mod) == mod - 1) ,即AA[0]的指数
*/
if(bigK.num.size() >= 8) realK = 1e8;
else
for(int i = (int)bigK.num.size() - 1; i >= 0; i --)
realK = realK * 10 + bigK.num[i];// ling ： 这里的realK也非真k，它唯一的作用就是与n / num_zore 比较大小
{
for(int i = bigK.num.size() - 1; i >= 1; i --)
bigK.num[i - 1] = (bigK.num[i - 1] + (long long)bigK.num[i] % mod * 10) % mod;
bigK.num.pop_back();
km = bigK.num[0] % mod; // ling : km is 真k % mod;
bigK = bigKK;
for(int i = bigK.num.size() - 1; i >= 1; i --)
bigK.num[i - 1] = (bigK.num[i - 1] + (long long)bigK.num[i] % (mod - 1) * 10) % (mod - 1);
bigK.num.pop_back();
ke = bigK.num[0] % (mod - 1); // ling : ke is 真k % Eular(mod)
}
}
aa.resize(n);
aa = opR(aa);
while(aa.size() && aa[aa.size() - 1] == 0) aa.pop_back();
int num_zore = n - aa.size();
aa = opR(aa);
if(num_zore * realK >= n)
return vec(n);
int first = aa[0], invfirst = quickpow(aa[0], mod - 2, mod);
int first_ke = quickpow(first, ke, mod);
for(int i = 0; i <= aa.size() - 1; i ++)
aa[i] = 1ll * aa[i] * invfirst % mod;
aa = GetLin(n, aa);
for(int i = 0; i <= aa.size() - 1; i ++)
aa[i] = 1ll * aa[i] * km % mod;
aa = GetExp(n, aa);
for(int i = 0; i <= aa.size() - 1; i ++)
aa[i] = 1ll * aa[i] * first_ke % mod;
num_zore = min(1ll * num_zore * realK, (long long)n); // ling：计算答案的前档0的个数 应用真k，但是如果真k太大就在前边return了
aa = opR(aa);
for(int i = 1; i <= num_zore; i ++)
aa.push_back(0);
aa = opR(aa);
aa.resize(n);
return aa;
}
vec MulT(vec a, vec b){// ling : 玄学的转置乘法
int n = a.size(), m = b.size();
reverse(b.begin(), b.end()),
b = Multiply(a, b);
for(int i = 0; i <= n - 1; i ++)
a[i] = b[i + m - 1];
return a;
}
void Init(vec &A, int k, int l, int r){
if(l == r){
Q[k].resize(2);
Q[k][0] = 1, Q[k][1] = (mod - A[l]);
return ;
}
int mid = (l + r) >> 1;
Init(A, k << 1, l, mid), Init(A, k << 1 | 1, mid + 1, r);
Q[k] = Multiply(Q[k << 1], Q[k << 1 | 1]);
return ;
}
void Multipoints(int k, int l, int r, vec F, vec &g){
F.resize(r - l + 1);
if(l == r) return void(g[l] = F[0]);
int mid = (l + r) >> 1;
Multipoints(k << 1, l, mid, MulT(F, Q[k << 1 | 1]), g);
Multipoints(k << 1 | 1, mid + 1, r, MulT(F, Q[k << 1]), g);
}
void Multipoint(vec f, vec a, vec &v, int n){
f.resize(n + 1), a.resize(n);
// ling ；调用此函数进行多项式多点快速求值
// f : 系数数组， a ： 所求自变量数组 ， v 与自变量对应的最终答案数组
// n = max（多项式项数，所求点数量）
for(int i = 0; i <= a.size() - 1; i ++)
a[i] %= mod;
Init(a, 1, 0, n - 1), v.resize(n);
Multipoints(1, 0, n - 1, MulT(f, GetInv(n + 1, Q[1])), v);
return ;
}
void GetP(int p, int l, int r, const vec &xi){ // P(x) = PAI (0 <= i <= n - 1) (x - xi)
if(l == r){
P[p].resize(2);
P[p][0] = -xi[l], P[p][1] = 1; // x - xi[i]
return ;
}
int mid = l + r >> 1;
GetP(p << 1, l, mid, xi);
GetP(p << 1 | 1, mid + 1, r, xi);
P[p] = Multiply(P[p << 1], P[p << 1 | 1]);
}
vec InterPart(int p, int l, int r, const vec &xi, const vec &yi, const vec &InterPxi){
if(l == r){
vec ret;
ret.push_back(1ll * yi[l] * InterPxi[l] % mod);
return ret;
}
int mid = l + r >> 1;
return PolyAdd(Multiply(InterPart(p << 1, l, mid, xi, yi, InterPxi), P[p << 1 | 1]), Multiply(InterPart(p << 1 | 1, mid + 1, r, xi, yi, InterPxi), P[p << 1]));
}
void Interpolation(const vec &xi, const vec &yi, vec &f){
// ling 调用此函数进行多项式快速插值
// x, y 已知横纵坐标
// f 返回的系数数组
GetP(1, 0, xi.size() - 1, xi);
InterP = differ(P[1]);
InterP.pop_back(); // ling ：求导之后最高位为0 ,pop()一下
Multipoint(InterP, xi, InterPxi, xi.size());
InterPxi = GetInvSum(InterPxi); // InterPxi -> InvInterPxi
f = InterPart(1, 0, xi.size() - 1, xi, yi, InterPxi);
return ;
}
vec ChaMultiply(int n, vec A, vec B, int mod){
//多项式差卷积：
//Ret[k] = Sum(x = k~n) A[x - k] * B[x];
A.resize(n + 1);
B.resize(n + 1);
vec Ans = MTT::mtt(A, opR(B), mod);
Ans.resize(n + 1);
return opR(Ans);
}
}
多项式模板（jiangly版）
#define add(x, y) ((x) + (y) >= mod ? (x) + (y) - mod : (x) + (y))
#define mus(x, y) ((x) < (y) ? (x) - (y) + mod : (x) - (y))
#define div2(x) (((x) & 1) ? (((x) + mod) >> 1) : ((x) >> 1))
#define vec vector<int>
constexpr int mod = 998244353;
void out(vec x){
for(int i = 0; i <= x.size() - 1; i ++)
cout << x[i] << " ";
cout << endl;
}
namespace poly{
const int maxn = 1e6+10;
int rev[maxn], revv = 0;
int inv[maxn], invv = 1;
vector<int> none;
vector<int> roots{0, 1};
inline int quickpow(int x, int b, int mod){
int j = 1;
for(; b; b >>= 1, x = 1ll * x * x % mod)
if(b & 1)
j = 1ll * j * x % mod;
return j % mod;
}
void getrev(int n){ //计算蝴蝶变换的序列
if(revv == n) return ;
revv = n;
int k = __builtin_ctz(n) - 1;
for(int i = 0; i < n; i ++)
rev[i] = rev[i >> 1] >> 1 | (i & 1) << k;
}
void getinv(int n){
if(invv >= n) return ;
inv[1] = 1;
for(int i = invv + 1; i <= n; i ++)
inv[i] = 1ll * (mod - mod / i) * inv[mod % i] % mod;
invv = n;
}
void dft(vec &a) {
int n = a.size();
getrev(n);
for (int i = 0; i < n; i++)
if (rev[i] < i) swap(a[i], a[rev[i]]);
if (int(roots.size()) < n) {
int k = __builtin_ctz(roots.size());
roots.resize(n);
while ((1 << k) < n) {
int e = quickpow(3, (mod - 1) >> (k + 1), mod); // 这里的 3 是 mod 的原根
for (int i = 1 << (k - 1); i < (1 << k); i ++) {
roots[i << 1] = roots[i];
roots[i << 1 | 1] = 1ll * roots[i] * e % mod;
}
k ++;
}
}
for (int k = 1; k < n; k <<= 1) {
for (int i = 0; i < n; i += (k << 1)) {
for (int j = 0; j < k; j ++) {
int u = a[i + j];
int v = 1ll * a[i + j + k] * roots[k + j] % mod;
a[i + j] = add(u, v);
a[i + j + k] = mus(u, v);
}
}
}
}
void idft(vec &a) {
int n = a.size();
reverse(a.begin() + 1, a.end());
dft(a);
int inv = quickpow(n, mod - 2, mod);
for (int i = 0; i < n; i ++) {
a[i] = 1ll * a[i] * inv % mod;
}
}
vec operator * (vec a, vec b) {
if (a.size() == 0 || b.size() == 0) {
return none;
}
int sz = 1, tot = a.size() + b.size() - 1;
while (sz < tot) sz *= 2;
a.resize(sz);
b.resize(sz);
dft(a);
dft(b);
for (int i = 0; i < sz; i ++) {
a[i] = 1ll * a[i] * b[i] % mod;
}
idft(a);
a.resize(tot);
for(int i = 0; i < a.size(); i ++)
if(a[i] < 0) a[i] += mod;
return a;
}
vec Multiply(vec a, vec b, int len = -1){
if(len == -1) return a * b;
vec ret = a * b;
ret.resize(len);
return ret;
}
vec ChaMultiply(int n, vec A, vec B){
//多项式差卷积：
//Ret[k] = Sum(x = k~n) A[x - k] * B[x];
A.resize(n + 1);
B.resize(n + 1);
reverse(B.begin(), B.end());
vec Ans = Multiply(A, B);
Ans.resize(n + 1);
reverse(Ans.begin(), Ans.end());
return Ans;
}
vec operator - (vec a, const vec &b){
a.resize(max(a.size(), b.size()));
for(int i = 0; i < b.size(); i ++)
a[i] = mus(a[i], b[i]);
return a;
}
vec operator + (vec a, const vec & b){
a.resize(max(a.size(), b.size()));
for(int i = 0; i < b.size(); i ++)
a[i] = add(a[i], b[i]);
return a;
}
vec Deriv(vec a){//求导
if(a.empty()) return none;
for(int i = 1; i < a.size(); i ++)
a[i - 1] = 1ll * a[i] * i % mod;
a.pop_back();
return a;
}
vec Integ(vec a){//积分
if(a.empty()) return none;
a.push_back(0);
getinv(a.size() - 1);
for(int i = a.size() - 1; i > 0; i --)
a[i] = 1ll * a[i - 1] * inv[i] % mod;
a[0] = 0;
return a;
}
vec cut(vec x, int m){
if(x.size() > m)
x.resize(m);
return x;
}
vec GetInv(vec a, int len = -1){
if(len == -1) len = a.size();
vec ret(1, quickpow(a[0], mod - 2, mod));
for(int l = 4; (l >> 2) < len; l <<= 1 ){
getrev(l);
vec tmp(a);
tmp.resize(l >> 1);
tmp.resize(l); dft(tmp);
ret.resize(l); dft(ret);
for(int i = 0; i < l; i ++)
ret[i] = (2 + mod - 1ll * tmp[i] * ret[i] % mod) * ret[i] % mod;
idft(ret);
ret.resize(l >> 1);
}
ret.resize(len);
return ret;
}
vec GetLn(vec a, int len = -1){ //对数，需满足f(0) = 1
if(len == -1) len = a.size();
a = Integ(Deriv(a) * GetInv(a, len));
a.resize(len);
return a;
}
vec GetExp(vec a,int len = -1){ //指数，需满足f(0) = 0
if(len == -1) len = a.size();
vec b(1, 1), c;
int n = a.size();
for(int l = 2; (l >> 1) < len; l <<= 1){
c = GetLn(b, l);
for(int i = 0; i < l; i ++)
c[i] = (i < n) ? mus(a[i], c[i]) : mus(0, c[i]);
c[0] = add(c[0], 1);
b = b * c;
b.resize(l);
}
b.resize(len);
return b;
}
vec GetSqrt(vec a, int len = -1){ //多项式开根，这里只实现了f(0) = 1的情况，若不等于1需要先算二次剩余
if(len == -1) len = a.size();
vec ret(1, 1);
for(int l = 4; (l >> 2) < len; l <<= 1){
getrev(l);
vec tmp(a);
tmp.resize(l >> 1);
tmp.resize(l); dft(tmp);
vec tt(GetInv(ret, l >> 1));
tt.resize(l); dft(tt);
for(int i = 0; i < l; i ++)
tmp[i] = 1ll * tmp[i] * tt[i] % mod;
idft(tmp);
ret.resize(l >> 1);
for(int i = 0; i < (l >> 1); i ++){
ret[i] = add(ret[i], tmp[i]);
ret[i] = div2(ret[i]);
}
}
ret.resize(len);
return ret;
}
vec GetInvSum(vec a){
int n = a.size() - 1;
vec s(n + 1), sinv(n + 1);
s[0] = a[0];
for(int i = 1; i <= n; i ++)
s[i] = 1ll * s[i - 1] * a[i] % mod;
sinv[n] = quickpow(s[n], mod - 2, mod);
for(int i = n - 1; i >= 0; i --)
sinv[i] = 1ll * sinv[i + 1] * a[i + 1] % mod;
for(int i = 1; i <= n; i ++)
sinv[i] = 1ll * sinv[i] * s[i - 1] % mod;
return sinv;
}
vec InterP, InterPxi, Q[maxn], P[maxn];// ling : 多项式快速多点求值和快速插值
vec MulT(vec a, vec b){// ling : 玄学的转置乘法
int n = a.size(), m = b.size();
reverse(b.begin(), b.end());
b = a * b;
for(int i = 0; i <= n - 1; i ++)
a[i] = b[i + m - 1];
return a;
}
void Init(vec &A, int k, int l, int r){
if(l == r){
Q[k].resize(2);
Q[k][0] = 1, Q[k][1] = (mod - A[l]);
return ;
}
int mid = (l + r) >> 1;
Init(A, k << 1, l, mid), Init(A, k << 1 | 1, mid + 1, r);
Q[k] = Q[k << 1] * Q[k << 1 | 1];
return ;
}
void Multipoints(int k, int l, int r, vec F, vec &g){
F.resize(r - l + 1);
if(l == r) return void(g[l] = F[0]);
int mid = (l + r) >> 1;
Multipoints(k << 1, l, mid, MulT(F, Q[k << 1 | 1]), g);
Multipoints(k << 1 | 1, mid + 1, r, MulT(F, Q[k << 1]), g);
}
void Multipoint(vec f, vec a, vec &v, int n){
f.resize(n + 1), a.resize(n);
// ling ；调用此函数进行多项式多点快速求值
// f : 系数数组， a ： 所求自变量数组 ， v 与自变量对应的最终答案数组
// n = max（多项式项数，所求点数量）
// ↑项数 = 模指 + 1
for(int i = 0; i < a.size(); i ++)
a[i] %= mod;
Init(a, 1, 0, n - 1), v.resize(n);
Multipoints(1, 0, n - 1, MulT(f, GetInv(Q[1], n + 1)), v);
return ;
}
void GetP(int p, int l, int r, const vec &xi){ // P(x) = PAI (0 <= i <= n - 1) (x - xi)
if(l == r){
P[p].resize(2);
P[p][0] = -xi[l], P[p][1] = 1; // x - xi[i]
return ;
}
int mid = l + r >> 1;
GetP(p << 1, l, mid, xi);
GetP(p << 1 | 1, mid + 1, r, xi);
P[p] = P[p << 1] * P[p << 1 | 1];
}
vec InterPart(int p, int l, int r, const vec &xi, const vec &yi, const vec &InterPxi){
if(l == r){
vec ret;
ret.push_back(1ll * yi[l] * InterPxi[l] % mod);
return ret;
}
int mid = l + r >> 1;
return InterPart(p << 1, l, mid, xi, yi, InterPxi) * P[p << 1 | 1]
+ InterPart(p << 1 | 1, mid + 1, r, xi, yi, InterPxi) * P[p << 1];
}
void Interpolation(const vec &xi, const vec &yi, vec &f){
// ling 调用此函数进行多项式快速插值
// x, y 已知横纵坐标
// f 返回的系数数组
GetP(1, 0, xi.size() - 1, xi);
InterP = Deriv(P[1]);
Multipoint(InterP, xi, InterPxi, xi.size());
InterPxi = GetInvSum(InterPxi); // InterPxi -> InvInterPxi
f = InterPart(1, 0, xi.size() - 1, xi, yi, InterPxi);
return ;
}
void ContinInter(const vec &yi, vec &ans, int m){
// ling : n 次多项式 给定f(0 ~ n) and m, ask f(m + 0 ~ m + n)
int n = yi.size() - 1;
vec invfac(n + 1); invfac[0] = 1;
for(int i = 1; i <= n; i ++)
invfac[i] = 1ll * invfac[i - 1] * i % mod;
invfac = GetInvSum(invfac);
vec A(n + 1), B(n << 1 | 1);
for(int i = 0; i <= n; i ++)
A[i] = 1ll * yi[i]
* ((n - i & 1) ? (mod - 1) : 1) % mod
* invfac[i] % mod
* invfac[n - i] % mod;
for(int i = 0; i <= (n << 1); i ++)
B[i] = add(i + m - n, mod);
B = GetInvSum(B);
ans = A * B;
ans.resize(n << 1 | 1);
reverse(ans.begin(), ans.end());
for(int i = 0; i <= n - 1; i ++)
ans.pop_back();
reverse(ans.begin(), ans.end());
vec t(n + 1); t[0] = 1;
vec it(n + 1); it[0] = 1;
for(int i = 1; i <= n; i ++)
it[i] = m + i - n - 1;
it = GetInvSum(it);
for(int i = m - n; i <= m; i ++)
t[0] = 1ll * t[0] * i % mod;
for(int i = 1; i <= n; i ++)
t[i] = 1ll * t[i - 1]
* (m + i) % mod
* it[i] % mod;
for(int i = 0; i <= n; i ++)
ans[i] = 1ll * ans[i] * t[i] % mod;
}
}
快速阶乘计算
#include<bits/stdc++.h>
namespace FACT{
#define IMP(lim,act) for(int qwq=(lim),i=0;i^qwq;++i)act
typedef double db;
const int M=3e5+5;
const db Pi=acos(-1);
int n,P,F[M],ifac[M];
struct Barrett{
typedef unsigned long long ull;
typedef __uint128_t LL;
ull m,B;
Barrett(const ull&m=2):m(m),B((LL(1)<<64)/m){}
friend inline ull operator%(const ull&a,const Barrett&mod){
ull r=a-mod.m*(LL(mod.B)*a>>64);return r>=mod.m?r-mod.m:r;
}
}mod;
struct complex{
db x,y;
complex(const db&x=0,const db&y=0):x(x),y(y){}
inline complex operator+(const complex&it)const{
return complex(x+it.x,y+it.y);
}
inline complex operator-(const complex&it)const{
return complex(x-it.x,y-it.y);
}
inline complex operator*(const complex&it)const{
return complex(x*it.x-y*it.y,x*it.y+y*it.x);
}
}buf[M],*w[20];
inline int Getlen(const int&n){
int len(0);while((1<<len)<n)++len;return len;
}
inline void swap(complex&a,complex&b){
complex c=a;a=b;b=c;
}
inline int Get(const db&x){
return((long long)(x+.5))%mod;
}
inline int qpow(int a,int b=P-2){
int ans(1);for(;b;b>>=1,a=1ll*a*a%mod)if(b&1)ans=1ll*ans*a%mod;return ans;
}
inline void NTT_init(const int&n){
const int&m=Getlen(n)-1;complex*now=buf;w[m]=now;now+=1<<m;
IMP(1<<m,w[m][i]=complex(std::cos(i*Pi/(1<<m)),std::sin(i*Pi/(1<<m))));
for(int k=m-1;k>=0&&(w[k]=now,now+=1<<k);--k)IMP(1<<k,w[k][i]=w[k+1][i<<1]);
}
inline void DFT(complex*f,const int&M){
const int&n=1<<M;
for(int len=n>>1,d=M-1;d>=0;--d,len>>=1)for(int k=0;k^n;k+=len<<1){
complex*W=w[d],*L=f+(k),*R=f+(k|len),x,y;IMP(len,(x=*L,y=*R)),*L++=(x+y),*R++=*W++*(x-y);
}
}
inline void IDFT(complex*f,const int&M){
const int&n=1<<M;
for(int len=1,d=0;d^M;++d,len<<=1)for(int k=0;k^n;k+=len<<1){
complex*W=w[d],*L=f+(k),*R=f+(k|len),x,y;IMP(len,(x=*L,y=*W++**R)),*L++=(x+y),*R++=(x-y);
}
IMP(n,(f[i].x/=n,f[i].y/=n));for(int i=1;(i<<1)<n;++i)swap(f[i],f[n-i]);
}
inline void MTT(int*f,int*g,int*h,const int&n,const int&m,const int&t,const int&Len=-1){
static complex Q[M],P[M],T[M];const int&len=Getlen(!~Len?n+m-1:Len);
IMP(n,(Q[i].x=f[i]&32767,P[i].x=f[i]>>15));IMP(m,T[i]=complex(g[i]&32767,g[i]>>15));
DFT(Q,len);DFT(P,len);DFT(T,len);IMP(1<<len,(Q[i]=Q[i]*T[i],P[i]=P[i]*T[i]));IDFT(Q,len);IDFT(P,len);
IMP(t,h[i]=(Get(Q[i].x)+(1ll*Get(Q[i].y+P[i].x)<<15)+(1ll*Get(P[i].y)<<30))%mod);IMP(1<<len,Q[i]=P[i]=T[i]=0);
}
inline void Getinv(int*f,const int&n){
static int pre[M];pre[0]=f[0];for(int i=1;i<n;++i)pre[i]=1ll*pre[i-1]*f[i]%mod;
int t,c=qpow(pre[n-1]);for(int i=n-1;i>=1;--i)t=f[i],f[i]=1ll*pre[i-1]*c%mod,c=1ll*c*t%mod;f[0]=c;
}
inline void PT(int*f,int*g,const int&n,const int&m){
static int F[M],G[M],H[M];H[0]=1;IMP(n,H[0]=1ll*H[0]*(m-i)%mod);IMP(n+n,G[i]=m+i-n);G[0]=1;Getinv(G,n+n);
IMP(n,F[i]=1ll*ifac[i]*(n-i-1&1?P-ifac[n-i-1]:ifac[n-i-1])%mod*f[i]%mod);
const int&len=Getlen(n+n);for(int i=1;i<n;++i)H[i]=1ll*(m+i)*G[i]%mod*H[i-1]%mod;MTT(F,G,F,n,n+n,n+n,n+n);
IMP(n,g[i]=1ll*F[n+i]*H[i]%mod);IMP(1<<len,F[i]=G[i]=H[i]=0);
}
inline int Getfac(const int&n,const int&p){
static int F[M],G[M];const int&B=sqrt(n),m=Getlen(B)-1;
mod=Barrett(P=p);ifac[0]=ifac[1]=1;for(int i=2;i<=B;++i)ifac[i]=1ll*(P-P/i)*ifac[P%i]%mod;
for(int i=1;i<=B;++i)ifac[i]=1ll*ifac[i-1]*ifac[i]%mod;F[0]=1;F[1]=B+1;
for(int i=m;i>=0;--i){
const int&q=B>>i+1,p=B>>i;if(!q)continue;
PT(F,G,q+1,q+1);IMP(q,F[i+q+1]=G[i]),G[i]=0;G[q]=0;
PT(F,G,q*2+1,1ll*q*qpow(B)%mod);IMP(q*2+1,F[i]=1ll*F[i]*G[i]%mod),G[i]=0;
if(q*2+1==p){
IMP(q*2+1,F[i]=1ll*F[i]*(i*B+p)%mod);F[p]=1;IMP(p,F[p]=1ll*F[p]*(p*B+i+1)%mod);
}
}
int ans(1);IMP(B,ans=1ll*ans*F[i]%mod),F[i]=0;F[B]=0;for(int i=B*B+1;i<=n;++i)ans=1ll*ans*i%mod;return ans;
}
int Get(int N, int P){
NTT_init(1<<17); // 可以预处理
return Getfac(N,P);
}
}
int main(){
}
快速求\((1 + x + x ^ 2 + ... x ^ {L - 1}) ^ m\)的方法
\[f = (1 + x + x^2 + ... x^{L - 1})
\\
f = (1 - x^L)(1 + x + x ^ 2 + x ^ 3...)
\\
f^m =(1 - x^L) ^ m (1 + x + x^2 + x^3...) ^ m
\\ = \big ( \sum_{k = 0}^\infty C_m^k (-1) ^ kx^{Lk}\big)\big (\sum_{k = 0}^m C_{k + m - 1} ^ k\big)
\\
= \big ( \sum_{k = 0}^\infty [L | k] C_m^{k / L} (-1) ^ {k / L}x^{k}\big)\big (\sum_{k = 0}^m C_{k + m - 1} ^ k\big)
\\
\therefore \ \ A[i] =[L | i] C_m^{i / L} (-1) ^ {i / L}
\ \ \ \ \ \ \ \
B[i] = C_{i + m - 1} ^ i
f = A * B(直接NTT 在O (nlogn)卷过去)
\]
*2022/10/10 新增快速求\((1 + x + x ^ 2 + ... x ^ {L - 1}) ^ m\)的方法
*2022/10/12 新增了jiangly的多项式模板
*2022/10/12 新完善了jiangly的多项式模板，现支持多项式差卷积，多点求值，快速插值
【推荐】国内首个AI IDE，深度理解中文开发场景，立即下载体验Trae
【推荐】编程新体验，更懂你的AI，立即体验豆包MarsCode编程助手
【推荐】抖音旗下AI助手豆包，你的智能百科全书，全免费不限次数
【推荐】轻量又高性能的 SSH 工具 IShell：AI 加持，快人一步
· 被坑几百块钱后，我竟然真的恢复了删除的微信聊天记录！
· 没有Manus邀请码？试试免邀请码的MGX或者开源的OpenManus吧
· 【自荐】一款简洁、开源的在线白板工具 Drawnix
· 园子的第一款AI主题卫衣上架——"HELLO! HOW CAN I ASSIST YOU TODAY
· 无需6万激活码！GitHub神秘组织3小时极速复刻Manus，手把手教你使用OpenManus搭建本