$$ a_n=\sum_{1\le i\le k}{f_ia_{n-i}} \\ \left\{ a_0,a_1….a_{k-1} \right\} \text{已知} \\ \\ \text{设}F\left( x \right) \text{表示从第}k\text{项开始的该数列的生成函数} \\ \sum_{i\ge k}{a_ix^i}=\sum_{i\ge k}{\sum_{1\le j\le k}{f_ja_{i-j}x^i}} \\ =\sum_{1\le j\le k}{f_j\sum_{i\ge k}{a_{i-j}x^i}} \\ =\sum_{1\le j\le k}{f_j\sum_{i\ge k-j}{a_ix^{i+j}}} \\ =\sum_{1\le j\le k}{f_jx^j\sum_{i\ge k-j}{a_ix^i}} \\ =\sum_{1\le j\le k}{f_jx^j\left( F\left( x \right) +\sum_{k-j\le i\le k-1}{a_ix^i} \right)} \\ \\ F\left( x \right) =F\left( x \right) \sum_{1\le j\le k}{f_jx^j}+\sum_{1\le j\le k}{f_jx^j\sum_{k-j\le i\le k-1}{a_ix^i}} \\ F\left( x \right) =\frac{\sum_{1\le j\le k}{f_jx^j\sum_{k-j\le i\le k-1}{a_ix^i}}}{1-\sum_{1\le j\le k}{f_jx^j}} \\ =\frac{\sum_{1\le j\le k}{f_jx^j\sum_{k-j\le i-j\le k-1}{a_{i-j}x^{i-j}}}}{1-\sum_{1\le j\le k}{f_jx^j}} \\ =\frac{\sum_{1\le j\le k}{f_j\sum_{k\le i\le k-1+j}{a_{i-j}x^i}}}{1-\sum_{1\le j\le k}{f_jx^j}} \\ =\frac{\sum_{k\le i\le 2k-1}{\begin{array}{c} x^i\sum_{i-k+1\le j\le k}{a_{i-j}f_j}\\ \end{array}}}{1-\sum_{1\le j\le k}{f_jx^j}} \\ F\left( x \right) +\sum_{0\le i\le k-1}{a_ix^i}=\frac{\left( 1-\sum_{1\le j\le k}{f_jx^j} \right) \left( \sum_{0\le i\le k-1}{a_ix^i} \right) +\sum_{1\le j\le k}{f_j\sum_{k\le i\le k-1+j}{a_{i-j}x^i}}}{1-\sum_{1\le j\le k}{f_jx^j}} \\ \frac{\sum_{0\le i\le k-1}{a_ix^i}-\sum_{1\le j\le k}{f_j\sum_{0\le i\le k-1}{a_ix^{i+j}}}+\sum_{1\le j\le k}{f_j\sum_{k\le i\le k-1+j}{a_{i-j}x^i}}}{1-\sum_{1\le j\le k}{f_jx^j}} \\ =\frac{\sum_{0\le i\le k-1}{a_ix^i}-\sum_{1\le j\le k}{f_j\sum_{j\le i\le j+k-1}{a_{i-j}x^i}}+\sum_{1\le j\le k}{f_j\sum_{k\le i\le k-1+j}{a_{i-j}x^i}}}{1-\sum_{1\le j\le k}{f_jx^j}} \\ \frac{\sum_{0\le i\le k-1}{a_ix^i}+\sum_{1\le j\le k}{f_j\left( \sum_{k\le i\le k-1+j}{a_{i-j}x^i}-\sum_{j\le i\le j+k-1}{a_{i-j}x^i} \right)}}{1-\sum_{1\le j\le k}{f_jx^j}} \\ \frac{\sum_{0\le i\le k-1}{a_ix^i}+\sum_{1\le j\le k}{f_j\left( \sum_{k\le i\le k-1+j}{a_{i-j}x^i}-\sum_{k\le i\le j+k-1}{a_{i-j}x^i}-\sum_{j\le i\le k-1}{a_{i-j}x^i} \right)}}{1-\sum_{1\le j\le k}{f_jx^j}} \\ \frac{\sum_{0\le i\le k-1}{a_ix^i}+\sum_{1\le j\le k}{f_j\left( -\sum_{j\le i\le k-1}{a_{i-j}x^i} \right)}}{1-\sum_{1\le j\le k}{f_jx^j}} \\ =\frac{\sum_{0\le i\le k-1}{a_ix^i}-\sum_{1\le j\le k}{f_j\left( \sum_{j\le i\le k-1}{a_{i-j}x^i} \right)}}{1-\sum_{1\le j\le k}{f_jx^j}} \\ =\frac{\sum_{0\le i\le k-1}{a_ix^i}-\sum_{1\le j\le k-1}{f_j\left( \sum_{j\le i\le k-1}{a_{i-j}x^i} \right)}}{1-\sum_{1\le j\le k}{f_jx^j}} \\ =\frac{\sum_{0\le i\le k-1}{a_ix^i}-\sum_{1\le i\le k-1}{x^i\sum_{1\le j\le i}{a_{i-j}f_j}}}{1-\sum_{1\le j\le k}{f_jx^j}} \\ =\frac{a_0+\sum_{1\le i\le k-1}{x^ia_i}-\sum_{1\le i\le k-1}{x^i\sum_{1\le j\le i}{a_{i-j}f_j}}}{1-\sum_{1\le j\le k}{f_jx^j}} \\ =\frac{a_0+\sum_{1\le i\le k-1}{x^i\left( a_i-\sum_{1\le j\le i}{a_{i-j}f_j} \right)}}{1-\sum_{1\le j\le k}{f_jx^j}} \\ \text{我们得到了原数列}\left\{ a_i \right\} \text{的生成函数}G\left( x \right) =\frac{P\left( x \right)}{Q\left( x \right)} \\ \text{考虑计算这个有理函数的}k\text{次项} \\ \text{令}G_k\left( x \right) =\frac{P\left( x \right)}{Q\left( x \right)}=\frac{P\left( x \right) P\left( -x \right)}{Q\left( x \right) Q\left( -x \right)}=\frac{xA\left( x^2 \right) +B\left( x^2 \right)}{C\left( x^2 \right)} \\ \text{即分母只有偶数次方项},\text{分子的奇数次方项和偶数次方项拆开} \\ \text{如果}k\text{为奇数},\text{那么}x^k\text{之可能存在于}x\frac{A\left( x^2 \right)}{C\left( x^2 \right)}\text{中}\left( \text{只有这边才有奇数次方} \right) \\ \text{同理},\text{如果}k\text{为偶数，那么}x^k\text{只能出现在}\frac{B\left( x^2 \right)}{C\left( x^2 \right)}\text{中，} \\ \text{然后根据}k\text{的奇偶性，继续递归} \\ \text{如果}k\text{是奇数},\text{那么计算}G_{\frac{k-1}{2}}\left( x \right) =\frac{A\left( x \right)}{C\left( x \right)}\left( \text{项的次数除以}2 \right) ,\text{答案为}\left[ x^{\frac{k-1}{2}} \right] G_{\frac{k-1}{2}}\left( x \right) \\ \text{如果}k\text{是偶数},\text{那么计算}G_{\frac{k}{2}}\left( x \right) =\frac{B\left( x \right)}{C\left( x \right)}\text{，答案为}\left[ x^k \right] G_{\frac{k}{2}}\left( x \right) \\ \text{然后有两种选择}:\text{递归到}k=0\text{时，计算}\frac{P\left( 0 \right)}{G\left( 0 \right)}\mathrm{mod}p\text{，此时不需要多项式求逆运算} \\ \text{但常数较大}\left( \text{递归次数多} \right) \\ \text{递归到}k<n\text{时，直接求逆计算。} \\ \\ \\ \\ \\ $$
代码1: 递归到$k=0$时执行整数运算，较慢
#include <cstring>
#include <iostream>
using int_t = long long int;
using std::cin;
using std::cout;
using std::endl;
const int_t LARGE = 2e6;
const int_t mod = 998244353;
const int_t g = 3;
int_t power(int_t b, int_t i) {
int_t r = 1;
if (i < 0)
i = (i % (mod - 1) + mod - 1) % (mod - 1);
while (i) {
if (i & 1)
r = r * b % mod;
b = b * b % mod;
i >>= 1;
}
return r;
}
void makeflip(int_t* arr, int_t size2) {
int_t len = (1 << size2);
arr[0] = 0;
for (int_t i = 1; i < len; i++) {
arr[i] = (arr[i >> 1] >> 1) | ((i & 1) << (size2 - 1));
}
}
int_t upper2n(int_t x) {
int_t r = 0;
while ((1 << r) < x)
r++;
return r;
}
template <int_t arg = 1>
void transform(int_t* A, int_t size2, int_t* flip) {
int_t len = (1 << size2);
for (int_t i = 0; i < len; i++) {
int_t r = flip[i];
if (r > i)
std::swap(A[i], A[r]);
}
for (int_t i = 2; i <= len; i *= 2) {
int_t mr = power(g, arg * (mod - 1) / i);
for (int_t j = 0; j < len; j += i) {
int_t curr = 1;
for (int_t k = 0; k < i / 2; k++) {
int_t u = A[j + k], v = curr * A[j + k + i / 2] % mod;
A[j + k] = (u + v) % mod;
A[j + k + i / 2] = (u - v + mod) % mod;
curr = curr * mr % mod;
}
}
}
}
void poly_inv(const int_t* A, int_t n, int_t* result) {
/*
计算B(x)*A(x) = 1 mod x^n, 其中A(x)已知
假设已知A(x)*B(x) = 1 mod x^{ceil(n/2)}
假设C(x)*A(x) = 1 mod x^n
(A(x)B(x)-1)^2 = A^2(x)B^2(x)-2A(x)B(x)+1= 0
A(x)B^2(x)-2B(x)+C(x) = 0
C(x) = 2B(x)-A(x)B^2(x)
*/
if (n == 1) {
result[0] = power(A[0], -1);
return;
}
int_t next = n / 2 + n % 2;
poly_inv(A, next, result);
//次数不要选错了，应该用n次的A和B去卷
int_t size2 = upper2n(n + 2 * next + 1);
static int_t X[LARGE];
static int_t Y[LARGE];
int_t len = (1 << size2);
memset(X + next, 0, sizeof(X[0]) * (len - n));
memset(Y + next, 0, sizeof(Y[0]) * (len - next));
memcpy(X, A, sizeof(A[0]) * n);
memcpy(Y, result, sizeof(result[0]) * next);
static int_t fliparr[LARGE];
makeflip(fliparr, size2);
transform<1>(X, size2, fliparr);
transform<1>(Y, size2, fliparr);
for (int_t i = 0; i < len; i++) {
X[i] = (2 * Y[i] - X[i] * Y[i] % mod * Y[i] % mod + mod) % mod;
}
transform<-1>(X, size2, fliparr);
const int_t inv = power(len, -1);
for (int_t i = 0; i < n; i++)
result[i] = X[i] * inv % mod;
}
int_t poly_divat(const int_t* F, const int_t* G, int_t n, int_t k) {
/*
n次多项式F和G
计算F(x)/G(x)的k次项前系数
考虑F(x)*G(-x)/G(x)*G(-x)，分母只有偶数次项，写为C(x^2)；分子写成xA(x^2)+B(x^2)，如果k是奇数，那么递归(A,C,n,k/2)，如果k是偶数，那么递归(B,C,n,k/2)
到k<=n时直接计算
*/
int_t size2 = upper2n(2 * n + 1);
int_t len = 1 << size2;
static int_t fliparr[LARGE];
makeflip(fliparr, size2);
static int_t T1[LARGE], T2[LARGE], T3[LARGE];
for (int_t i = 0; i < len; i++) {
if (i <= n)
T1[i] = F[i], T2[i] = G[i];
else
T1[i] = T2[i] = T3[i] = 0;
}
const int_t inv = power(len, -1);
while (k != 0) {
for (int_t i = 0; i < len; i++) {
if (i <= n) {
T3[i] = T2[i] * (i % 2 ? (mod - 1) : 1);
} else
T3[i] = 0;
}
transform(T1, size2, fliparr);
transform(T2, size2, fliparr);
transform(T3, size2, fliparr);
for (int_t i = 0; i < len; i++) {
T1[i] = T1[i] * T3[i] % mod;
T2[i] = T2[i] * T3[i] % mod;
}
transform<-1>(T1, size2, fliparr);
transform<-1>(T2, size2, fliparr);
for (int_t i = 0; i < len; i++) {
if (i * 2 < len) {
T2[i] = T2[i * 2] * inv % mod;
} else
T2[i] = 0;
}
int_t b = k % 2;
for (int_t i = 0; i < len; i++) {
if (i % 2 == b) {
T1[i / 2] = T1[i] * inv % mod;
}
if (i > 0) //防止把T1[0]改为0
T1[i] = 0;
}
k >>= 1;
}
return T1[0] * power(T2[0], -1) % mod;
}
void poly_mul(const int_t* A, int_t n, const int_t* B, int_t m, int_t* C) {
int_t size2 = upper2n(n + m + 1);
int_t len = 1 << size2;
static int_t T1[LARGE], T2[LARGE];
for (int_t i = 0; i < len; i++) {
if (i <= n)
T1[i] = A[i];
else
T1[i] = 0;
if (i <= m)
T2[i] = B[i];
else
T2[i] = 0;
}
static int_t fliparr[LARGE];
makeflip(fliparr, size2);
transform(T1, size2, fliparr);
transform(T2, size2, fliparr);
for (int_t i = 0; i < len; i++)
T1[i] = T1[i] * T2[i] % mod;
transform<-1>(T1, size2, fliparr);
int_t inv = power(len, -1);
for (int_t i = 0; i <= n + m; i++)
C[i] = T1[i] * inv % mod;
}
int main() {
std::ios::sync_with_stdio(false);
static int_t A[LARGE], F[LARGE];
static int_t T1[LARGE], T2[LARGE];
int_t n, k;
cin >> n >> k;
for (int_t i = 1; i <= k; i++) {
cin >> F[i];
F[i] = (F[i] % mod + mod) % mod;
}
F[0] = 0;
for (int_t i = 0; i < k; i++) {
cin >> A[i];
A[i] = (A[i] % mod + mod) % mod;
}
A[k] = 0;
poly_mul(A, k, F, k, T1);
T1[0] = A[0];
for (int_t i = 1; i <= k - 1; i++) {
T1[i] = (A[i] - T1[i] + mod) % mod;
}
for (int_t i = k; i <= 2 * k; i++)
T1[i] = 0;
T2[0] = 1;
for (int_t i = 1; i <= k; i++)
T2[i] = (mod - F[i]) % mod;
int_t r = poly_divat(T1, T2, k, n);
cout << r << endl;
return 0;
}
/*
(1+2*x)/(1+x+x^2)
2 5
1 2 0
1 1 1
ans = -2 = 998244351
*/
代码2: 递归到$k<n$时直接求逆运算，较快，大约是代码1的一半时间
#include <cstring>
#include <iostream>
using int_t = long long int;
using std::cin;
using std::cout;
using std::endl;
const int_t LARGE = 2e6;
const int_t mod = 998244353;
const int_t g = 3;
int_t power(int_t b, int_t i) {
int_t r = 1;
if (i < 0)
i = (i % (mod - 1) + mod - 1) % (mod - 1);
while (i) {
if (i & 1)
r = r * b % mod;
b = b * b % mod;
i >>= 1;
}
return r;
}
void makeflip(int_t* arr, int_t size2) {
int_t len = (1 << size2);
arr[0] = 0;
for (int_t i = 1; i < len; i++) {
arr[i] = (arr[i >> 1] >> 1) | ((i & 1) << (size2 - 1));
}
}
int_t upper2n(int_t x) {
int_t r = 0;
while ((1 << r) < x)
r++;
return r;
}
template <int_t arg = 1>
void transform(int_t* A, int_t size2, int_t* flip) {
int_t len = (1 << size2);
for (int_t i = 0; i < len; i++) {
int_t r = flip[i];
if (r > i)
std::swap(A[i], A[r]);
}
for (int_t i = 2; i <= len; i *= 2) {
int_t mr = power(g, arg * (mod - 1) / i);
for (int_t j = 0; j < len; j += i) {
int_t curr = 1;
for (int_t k = 0; k < i / 2; k++) {
int_t u = A[j + k], v = curr * A[j + k + i / 2] % mod;
A[j + k] = (u + v) % mod;
A[j + k + i / 2] = (u - v + mod) % mod;
curr = curr * mr % mod;
}
}
}
}
void poly_inv(const int_t* A, int_t n, int_t* result) {
/*
计算B(x)*A(x) = 1 mod x^n, 其中A(x)已知
假设已知A(x)*B(x) = 1 mod x^{ceil(n/2)}
假设C(x)*A(x) = 1 mod x^n
(A(x)B(x)-1)^2 = A^2(x)B^2(x)-2A(x)B(x)+1= 0
A(x)B^2(x)-2B(x)+C(x) = 0
C(x) = 2B(x)-A(x)B^2(x)
*/
if (n == 1) {
result[0] = power(A[0], -1);
return;
}
int_t next = n / 2 + n % 2;
poly_inv(A, next, result);
//次数不要选错了，应该用n次的A和B去卷
int_t size2 = upper2n(n + 2 * next + 1);
static int_t X[LARGE];
static int_t Y[LARGE];
int_t len = (1 << size2);
memset(X + next, 0, sizeof(X[0]) * (len - n));
memset(Y + next, 0, sizeof(Y[0]) * (len - next));
memcpy(X, A, sizeof(A[0]) * n);
memcpy(Y, result, sizeof(result[0]) * next);
static int_t fliparr[LARGE];
makeflip(fliparr, size2);
transform<1>(X, size2, fliparr);
transform<1>(Y, size2, fliparr);
for (int_t i = 0; i < len; i++) {
X[i] = (2 * Y[i] - X[i] * Y[i] % mod * Y[i] % mod + mod) % mod;
}
transform<-1>(X, size2, fliparr);
const int_t inv = power(len, -1);
for (int_t i = 0; i < n; i++)
result[i] = X[i] * inv % mod;
}
int_t poly_divat(const int_t* F, const int_t* G, int_t n, int_t k) {
/*
n次多项式F和G
计算F(x)/G(x)的k次项前系数
考虑F(x)*G(-x)/G(x)*G(-x)，分母只有偶数次项，写为C(x^2)；分子写成xA(x^2)+B(x^2)，如果k是奇数，那么递归(A,C,n,k/2)，如果k是偶数，那么递归(B,C,n,k/2)
到k<=n时直接计算
*/
int_t size2 = upper2n(2 * n + 1);
int_t len = 1 << size2;
static int_t fliparr[LARGE];
makeflip(fliparr, size2);
static int_t T1[LARGE], T2[LARGE], T3[LARGE];
for (int_t i = 0; i < len; i++) {
if (i <= n)
T1[i] = F[i], T2[i] = G[i];
else
T1[i] = T2[i] = T3[i] = 0;
}
const int_t inv = power(len, -1);
while (k >= n) {
for (int_t i = 0; i < len; i++) {
if (i <= n) {
T3[i] = T2[i] * (i % 2 ? (mod - 1) : 1);
} else
T3[i] = 0;
}
transform(T1, size2, fliparr);
transform(T2, size2, fliparr);
transform(T3, size2, fliparr);
for (int_t i = 0; i < len; i++) {
T1[i] = T1[i] * T3[i] % mod;
T2[i] = T2[i] * T3[i] % mod;
}
transform<-1>(T1, size2, fliparr);
transform<-1>(T2, size2, fliparr);
for (int_t i = 0; i < len; i++) {
if (i * 2 < len) {
T2[i] = T2[i * 2] * inv % mod;
} else
T2[i] = 0;
}
int_t b = k % 2;
for (int_t i = 0; i < len; i++) {
if (i % 2 == b) {
T1[i / 2] = T1[i] * inv % mod;
}
if (i > 0) //防止把T1[0]改为0
T1[i] = 0;
}
k >>= 1;
}
poly_inv(T2, k + 1, T3);
int_t result = 0;
//计算结果的k次项
for (int_t i = 0; i <= k; i++) {
result = (result + T1[i] * T3[k - i] % mod) % mod;
}
return result;
}
void poly_mul(const int_t* A, int_t n, const int_t* B, int_t m, int_t* C) {
int_t size2 = upper2n(n + m + 1);
int_t len = 1 << size2;
static int_t T1[LARGE], T2[LARGE];
for (int_t i = 0; i < len; i++) {
if (i <= n)
T1[i] = A[i];
else
T1[i] = 0;
if (i <= m)
T2[i] = B[i];
else
T2[i] = 0;
}
static int_t fliparr[LARGE];
makeflip(fliparr, size2);
transform(T1, size2, fliparr);
transform(T2, size2, fliparr);
for (int_t i = 0; i < len; i++)
T1[i] = T1[i] * T2[i] % mod;
transform<-1>(T1, size2, fliparr);
int_t inv = power(len, -1);
for (int_t i = 0; i <= n + m; i++)
C[i] = T1[i] * inv % mod;
}
int main() {
std::ios::sync_with_stdio(false);
static int_t A[LARGE], F[LARGE];
static int_t T1[LARGE], T2[LARGE];
int_t n, k;
cin >> n >> k;
for (int_t i = 1; i <= k; i++) {
cin >> F[i];
F[i] = (F[i] % mod + mod) % mod;
}
F[0] = 0;
for (int_t i = 0; i < k; i++) {
cin >> A[i];
A[i] = (A[i] % mod + mod) % mod;
}
A[k] = 0;
poly_mul(A, k, F, k, T1);
T1[0] = A[0];
for (int_t i = 1; i <= k - 1; i++) {
T1[i] = (A[i] - T1[i] + mod) % mod;
}
for (int_t i = k; i <= 2 * k; i++)
T1[i] = 0;
T2[0] = 1;
for (int_t i = 1; i <= k; i++)
T2[i] = (mod - F[i]) % mod;
int_t r = poly_divat(T1, T2, k, n);
cout << r << endl;
return 0;
}
/*
(1+2*x)/(1+x+x^2)
2 5
1 2 0
1 1 1
ans = -2 = 998244351
*/
发表回复