考虑设
f
[
x
]
[
y
]
f[x][y]
f[x][y] 为从
(
x
,
y
)
(x,y)
(x,y) 出发走到
(
0
,
0
)
(0,0)
(0,0) 的期望步数。有
f
[
x
]
[
y
]
=
p
∗
f
[
x
+
1
]
[
y
]
+
q
∗
f
[
x
]
[
y
+
1
]
+
1
f[x][y]=p*f[x+1][y]+q*f[x][y+1]+1
f[x][y]=p∗f[x+1][y]+q∗f[x][y+1]+1
推出
f
[
x
+
1
]
[
y
]
=
(
f
[
x
]
[
y
]
−
q
∗
f
[
x
]
[
y
+
1
]
−
1
)
/
p
f[x+1][y]=(f[x][y]-q*f[x][y+1]-1)/p
f[x+1][y]=(f[x][y]−q∗f[x][y+1]−1)/p 于是设第一行为主元,可以递推出第0行用这些主元表示的系数
g
(
x
,
y
)
g(x,y)
g(x,y) 。设
G
x
(
t
)
=
∑
i
=
0
m
−
1
g
(
x
,
i
)
t
i
G_x(t) = \sum_{i=0}^{m-1}g(x,i)t^i
Gx(t)=∑i=0m−1g(x,i)ti 。在循环卷积意义下有
G
x
(
t
)
=
G
x
−
1
(
t
)
1
−
q
t
p
G_x(t)=G_{x-1}(t) \frac{1-qt}{p}
Gx(t)=Gx−1(t)p1−qt ,于是通过循环卷积可以快速求出第0行的系数,再利用
f
[
x
]
[
y
]
=
p
∗
f
[
x
+
1
]
[
y
]
+
q
∗
f
[
x
]
[
y
+
1
]
+
1
f[x][y]=p*f[x+1][y]+q*f[x][y+1]+1
f[x][y]=p∗f[x+1][y]+q∗f[x][y+1]+1 和
f
[
0
]
[
0
]
=
0
f[0][0]=0
f[0][0]=0 列出方程组,高斯消元即可算出主元,接着对于每个询问,同样做循环卷积算出系数,再带入算出的主元即可
复杂度
O
(
m
3
+
t
m
log
m
log
n
)
O(m^3+tm \log m \log n)
O(m3+tmlogmlogn)
#include <bits/stdc++.h>
#define int long long
using std::array;
using std::cerr;
using std::cin;
using std::copy;
using std::cout;
using std::function;
using std::get;
using std::make_pair;
using std::make_tuple;
using std::map;
using std::mt19937;
using std::pair;
using std::reverse;
using std::set;
using std::sort;
using std::swap;
using std::tuple;
using std::uniform_int_distribution;
using std::unique;
using std::vector;
using ll = long long;
namespace qwq {
mt19937 eng;
void init(int Seed) { return eng.seed(Seed); }
int rnd(int l = 1, int r = 1000000000) { return uniform_int_distribution<int>(l, r)(eng); }
} // namespace qwq
template <typename T>
inline void chkmin(T &x, T y) {
if (x > y)
x = y;
}
template <typename T>
inline void chkmax(T &x, T y) {
if (x < y)
x = y;
}
template <typename T>
inline T min(const T &x, const T &y) {
return x < y ? x : y;
}
template <typename T>
inline T max(const T &x, const T &y) {
return x > y ? x : y;
}
char buf[100000], *bufs, *buft;
#define gc() \
((bufs == buft && (buft = (bufs = buf) + fread(buf, 1, 100000, stdin))), bufs == buft ? -1 : *bufs++)
#define O(x) cerr << #x << " : " << x << '\n'
const double Pi = acos(-1);
const int MAXN = 262144, MOD = 998244353, inv2 = (MOD + 1) / 2, I32_INF = 0x3f3f3f3f;
const long long I64_INF = 0x3f3f3f3f3f3f3f3f;
auto Ksm = [](int x, int y) -> int {
if (y < 0) {
y %= MOD - 1;
y += MOD - 1;
}
int ret = 1;
for (; y; y /= 2, x = (long long)x * x % MOD)
if (y & 1)
ret = (long long)ret * x % MOD;
return ret;
};
auto Mod = [](int x) -> int {
if (x >= MOD)
return x - MOD;
else if (x < 0)
return x + MOD;
else
return x;
};
template <const int N_num, const int M_num>
struct Graph {
int H[N_num];
struct Edge {
int to, lac;
} e[M_num];
inline void add_edge(int x, int y) {
e[*H] = { y, H[x] };
H[x] = (*H)++;
}
inline void init() {
memset(H, -1, sizeof H);
*H = 0;
}
};
#define go(x, y) for (int i = x.H[y], v; (v = x.e[i].to) && ~i; i = x.e[i].lac)
inline int ls(int k) { return k << 1; }
inline int rs(int k) { return k << 1 | 1; }
using ull = unsigned long long;
int fac[MAXN], ifac[MAXN];
namespace POLY {
int SZ, R[MAXN], W[MAXN], INV[MAXN + 1];
void POLYINIT() {
INV[1] = 1;
for (int i = 2; i <= MAXN; ++i) INV[i] = (long long)(MOD - MOD / i) * INV[MOD % i] % MOD;
}
void INIT(int len) {
if (SZ == len)
return;
SZ = len;
for (int i = 1; i < len; ++i) R[i] = (R[i >> 1] >> 1) | (i & 1 ? (len >> 1) : 0);
int wn = Ksm(3, (MOD - 1) / len);
W[len >> 1] = 1;
for (int i = (len >> 1) + 1; i < len; ++i) W[i] = (long long)W[i - 1] * wn % MOD;
for (int i = (len >> 1) - 1; i > 0; --i) W[i] = W[i << 1];
}
unsigned long long c[MAXN];
void Ntt(vector<int> &F, int limit, int type) {
copy(F.begin(), F.begin() + limit, c);
for (int i = 1; i < limit; ++i)
if (i < R[i])
swap(c[i], c[R[i]]);
for (int o = 2, j = 1; o <= limit; o <<= 1, j <<= 1) {
for (int i = 0; i < limit; i += o) {
for (int k = 0; k < j; ++k) {
unsigned long long OI = c[i + j + k] * W[k + j] % MOD;
c[i + j + k] = c[i + k] + MOD - OI;
c[i + k] += OI;
}
}
}
if (type == -1) {
reverse(c + 1, c + limit);
int inv = INV[limit];
for (int i = 0; i < limit; ++i) c[i] = c[i] % MOD * inv % MOD;
}
for (int i = 0; i < limit; ++i) F[i] = c[i] % MOD;
}
int w;
typedef std::pair<int, int> complex;
complex operator+(complex &x, complex &y) {
return std::make_pair((x.first + y.first) % MOD, (x.second + y.second) % MOD);
}
complex operator-(complex &x, complex &y) {
return std::make_pair((x.first - y.first + MOD) % MOD, (x.second - y.second + MOD) % MOD);
}
complex operator*(complex &x, complex &y) {
return std::make_pair(((long long)x.first * y.first + (long long)w * x.second % MOD * y.second) % MOD,
((long long)x.first * y.second + (long long)y.first * x.second) % MOD);
}
complex ksm(complex x, int y) {
complex res(1, 0);
for (; y; y /= 2, x = x * x)
if (y & 1)
res = res * x;
return res;
}
bool check(int x) { return Ksm(x, (MOD - 1) / 2) == 1; }
int Sqrt(int x) {
if (!x)
return x;
else {
int a = qwq::rnd() % MOD;
while (check(((long long)a * a + MOD - x) % MOD)) a = qwq::rnd() % MOD;
w = ((long long)a * a + MOD - x) % MOD;
complex b(a, 1);
int ans1(ksm(b, (MOD + 1) / 2).first);
return min(ans1, MOD - ans1);
}
}
int M, N;
struct Poly {
vector<int> v;
int &operator[](const int &pos) { return v[pos]; }
int len() { return v.size(); }
void set(int l) { return v.resize(l); }
void adjust() {
while (v.size() > 1 && !v.back()) v.pop_back();
}
void rev() { reverse(v.begin(), v.end()); }
void Ntt(int L, int type) {
int limit = 1 << L;
INIT(limit);
set(limit);
POLY::Ntt(v, limit, type);
}
void Squ() {
int L = ceil(log2(len())) + 1, limit = 1 << L;
Ntt(L, 1);
for (int i = 0; i < limit; ++i) v[i] = (long long)v[i] * v[i] % MOD;
Ntt(L, -1);
adjust();
}
void operator+=(Poly &x) {
if (len() < x.len())
set(x.len());
for (int i = 0; i < x.len(); ++i) v[i] = Mod(v[i] + x[i]);
adjust();
}
void operator-=(Poly &x) {
if (len() < x.len())
set(x.len());
for (int i = 0; i < x.len(); ++i) v[i] = Mod(v[i] - x[i]);
adjust();
}
Poly operator*(Poly &x) {
Poly ret, tmp0 = *this, tmp1 = x;
int L = ceil(log2(tmp0.len() + tmp1.len() - 1)), n = 1 << L;
tmp0.Ntt(L, 1);
tmp1.Ntt(L, 1);
ret.set(n);
for (int i = 0; i < n; ++i) ret[i] = (long long)tmp0[i] * tmp1[i] % MOD;
ret.Ntt(L, -1);
for (int i = 0; i < M; ++i) ret[i] = (ret[i] + ret[i + M]) % MOD;
ret.set(M);
return ret;
}
Poly operator-(Poly &x) {
Poly ret;
ret.set(max(len(), x.len()));
for (int i = 0; i < len(); ++i) ret[i] = v[i];
for (int i = 0; i < x.len(); ++i) ret[i] = Mod(ret[i] - x[i]);
return ret;
}
Poly operator+(Poly &x) {
Poly ret;
ret.set(max(len(), x.len()));
for (int i = 0; i < len(); ++i) ret[i] = v[i];
for (int i = 0; i < x.len(); ++i) ret[i] = Mod(ret[i] + x[i]);
return ret;
}
void operator*=(Poly &x) {
Poly tmp = x;
int L = ceil(log2(len() + x.len() - 1)), n = 1 << L;
Ntt(L, 1);
tmp.Ntt(L, 1);
for (int i = 0; i < n; ++i) v[i] = (long long)v[i] * tmp[i] % MOD;
Ntt(L, -1);
for (int i = 0; i < M; ++i) v[i] = (v[i] + v[i + M]) % MOD;
set(M);
}
Poly GetInv(int deg = -1) {
if (deg == 1)
return { { Ksm(v[0], MOD - 2) } };
Poly ret = GetInv((deg + 1) / 2), tmp;
int L = ceil(log2(deg)) + 1, n = 1 << L, mx = min(len(), deg);
tmp.set(deg);
for (int i = 0; i < mx; ++i) tmp[i] = v[i];
tmp.Ntt(L, 1);
ret.Ntt(L, 1);
for (int i = 0; i < n; ++i) ret[i] = (2 - (long long)tmp[i] * ret[i] % MOD + MOD) * ret[i] % MOD;
ret.Ntt(L, -1);
ret.set(deg);
return ret;
}
pair<Poly, Poly> operator%(Poly &x) {
if (x.len() > len())
return { { { 0 } }, *this };
Poly tmp0 = *this, tmp1 = x;
tmp0.rev();
tmp1.rev();
tmp1 = tmp1.GetInv(len() - x.len() + 1);
tmp0 *= tmp1;
tmp0.set(len() - x.len() + 1);
tmp0.rev();
tmp1 = tmp0 * x;
Poly ret = *this - tmp1;
ret.set(x.len() - 1);
return { tmp0, ret };
}
Poly Dif(int deg = -1) {
Poly tmp;
tmp.set(max(len() - 1, 1ll));
for (int i = 0; i < len() - 1; ++i) tmp[i] = v[i + 1] * (i + 1LL) % MOD;
if (~deg)
tmp.set(deg);
return tmp;
}
Poly GetSqrt(int deg = -1) {
if (deg == 1)
return { { POLY::Sqrt(v[0]) } };
Poly ret = GetSqrt((deg + 1) / 2), tmp0 = ret.GetInv(deg), tmp1;
int L = ceil(log2(deg)) + 1, mx = min(len(), deg);
tmp1.set(deg);
for (int i = 0; i < mx; ++i) tmp1[i] = v[i];
tmp0 *= tmp1;
tmp0.set(deg);
ret += tmp0;
for (auto &i : ret.v) i = (long long)i * inv2 % MOD;
return ret;
}
Poly Int(int deg = -1) {
Poly tmp;
tmp.set(len() + 1);
for (int i = 1; i < tmp.len(); ++i) tmp[i] = (long long)v[i - 1] * INV[i] % MOD;
if (~deg)
tmp.set(deg);
return tmp;
}
Poly GetLn(int deg = -1) {
Poly tmp0 = Dif(deg), tmp1 = GetInv(deg);
tmp0 *= tmp1;
tmp0.set(deg);
return tmp0.Int(deg);
}
Poly GetExp(int deg = -1) {
if (deg == 1)
return { { 1 } };
Poly tmp0 = GetExp((deg + 1) / 2), tmp1 = tmp0.GetLn(deg);
for (int i = 0; i < deg; ++i) tmp1[i] = Mod(v[i] - tmp1[i]);
tmp1[0] = Mod(tmp1[0] + 1);
tmp0 *= tmp1;
tmp0.set(deg);
return tmp0;
}
Poly GetKsm(int deg, int K) {
Poly tmp0 = GetLn(deg);
for (auto &i : tmp0.v) i = (long long)i * K % MOD;
return tmp0.GetExp(deg);
}
};
vector<int> getmulpointvalue(Poly &a, vector<int> &x) {
static Poly tmp[MAXN * 4];
function<void(int, int, int)> get = [&](int u, int l, int r) -> void {
if (l == r) {
tmp[u] = { { Mod(-x[l]), 1 } };
return;
}
int mid = (l + r) / 2;
get(ls(u), l, mid);
get(rs(u), mid + 1, r);
tmp[u] = tmp[ls(u)] * tmp[rs(u)];
};
get(1, 0, x.size() - 1);
vector<int> ret(x.size());
function<void(int, int, Poly, int)> solve = [&](int l, int r, Poly f, int u) -> void {
if (l == r) {
ret[l] = f[0];
return;
}
int mid = (l + r) / 2;
solve(l, mid, (f % tmp[ls(u)]).second, ls(u));
solve(mid + 1, r, (f % tmp[rs(u)]).second, rs(u));
};
solve(0, x.size() - 1, (a % tmp[1]).second, 1);
return ret;
}
vector<int> Interpolation(vector<pair<int, int>> &x) {
static Poly tmp[MAXN * 4];
vector<int> tmpx;
for (auto &i : x) tmpx.push_back(i.first);
function<void(int, int, int)> get = [&](int u, int l, int r) -> void {
if (l == r) {
tmp[u] = { { Mod(-tmpx[l]), 1 } };
return;
}
int mid = (l + r) / 2;
get(ls(u), l, mid);
get(rs(u), mid + 1, r);
tmp[u] = tmp[ls(u)] * tmp[rs(u)];
};
get(1, 0, tmpx.size() - 1);
Poly ret = tmp[1];
ret = ret.Dif();
vector<int> tmpy = getmulpointvalue(ret, tmpx);
for (int i = 0; i < x.size(); ++i) tmpy[i] = (long long)Ksm(tmpy[i], MOD - 2) * x[i].second % MOD;
function<Poly(int, int, int)> solve = [&](int u, int l, int r) -> Poly {
if (l == r)
return { { tmpy[l] } };
int mid = (l + r) / 2;
Poly tmp0 = solve(ls(u), l, mid) * tmp[rs(u)];
Poly tmp1 = solve(rs(u), mid + 1, r) * tmp[ls(u)];
return tmp0 + tmp1;
};
return solve(1, 0, x.size() - 1).v;
}
// Poly t[MAXN * 4];
} // namespace POLY
using namespace POLY;
Poly f, g;
Poly ksm(Poly x, int y) {
Poly res;
res.set(M);
res[0] = 1;
while (y) {
// cout<<res[0]<<' '<<x[0]<<'\n';
if (y & 1)
res *= x;
// cout<<res[0]<<'\n';
x *= x; //
y >>= 1;
}
return res;
}
ll val[MAXN], mat[2100][2100];
void solve() {
for (int i = 1; i <= M; i++) {
for (int j = i; j <= M; j++)
if (mat[j][i]) {
if (j != i)
swap(mat[i], mat[j]);
break;
}
ll inv = Ksm(mat[i][i], MOD - 2);
for (int j = i + 1; j <= M; j++) {
ll mul = inv * (MOD - mat[j][i]) % MOD;
(mat[j][0] += mat[i][0] * mul) %= MOD;
for (int k = i; k <= M; k++) (mat[j][k] += mat[i][k] * mul) %= MOD;
}
}
for (int i = 1; i <= M; i++) mat[i][0] = MOD - mat[i][0];
for (int i = M; i; i--) {
ll v = mat[i][0] * Ksm(mat[i][i], MOD - 2) % MOD;
val[i - 1] = v;
for (int j = i - 1; j > 0; j--) (mat[j][0] += MOD - v * mat[j][i] % MOD) %= MOD;
}
}
inline int read() {
int x = 0, f = 1;
char ch = getchar();
while (ch < '0' || ch > '9') {
if (ch == '-')
f = -1;
ch = getchar();
}
while (ch >= '0' && ch <= '9') {
x = x * 10 + ch - '0';
ch = getchar();
}
return x * f;
}
signed main() {
freopen("huawei.in", "r", stdin);
freopen("huawei.out", "w", stdout);
POLYINIT();
N = read(), M = read();
int p = read(), q = read();
int s = Ksm(p + q, MOD - 2);
(p *= s) %= MOD, (q *= s) %= MOD;
int ip = Ksm(p, MOD - 2), iq = Ksm(q, MOD - 2);
f.set(M);
f[0] = ip, f[1] = (MOD - q) * ip % MOD;
f = ksm(f, N - 1);
// cout<<f[0]<<" "<<f[1]<<'\n';
for (int i = 1; i <= M; i++) {
mat[i][0] = (MOD - ip) * (N - 1) % MOD;
for (int j = 1; j <= M; j++) mat[i][j] = f[(i - j + M) % M];
if (i > 1) {
mat[i][0] += MOD - 1 - (MOD - ip) * (N - 1) % MOD * q % MOD;
mat[i][0] %= MOD;
for (int j = 1; j <= M; j++) (mat[i][j] += MOD - f[(i - j + M - 1) % M] * q % MOD) %= MOD;
(mat[i][i] += MOD - p) %= MOD;
}
}
/* for(int i=1;i<=M;i++){
for(int j=0;j<=M;j++)
printf("%lld ",mat[i][j]);
puts("");
}*/
solve();
int t = read();
while (t--) {
int x = read(), y = read();
ll ans = (MOD - ip) * (N - 1 - x) % MOD;
f.set(0);
f.set(M);
f[0] = ip, f[1] = (MOD - q) * ip % MOD;
f = ksm(f, N - 1 - x);
for (int i = 0; i < M; i++) (ans += f[(y - i + M) % M] * val[i]) %= MOD;
printf("%lld\n", (ans % MOD + MOD) % MOD);
}
return 0;
}