大致题意:在二维坐标系中,每个整点位置都有一垂直于y=−xy=-xy=−x的玻璃,光线可折射和反射。折射率为aaa,反射率为bbb。现从(−0.5,0)(-0.5,0)(−0.5,0)朝xxx轴正方向射出光强为111的光线,问(n+0.5,m)(n+0.5,m)(n+0.5,m)处的光强。
解:
令x=n+1,y=mx=n+1,y=mx=n+1,y=m,期间共经过了x+yx+yx+y次玻璃。假设与玻璃接触后射向xxx轴正方向记为000,射向yyy轴正方向记为111,那么就可以得到一个以000结尾,xxx个000和yyy个111组成的010101字符串。其中,相邻两个字符相同会贡献×a\times a×a,否则贡献×b\times b×b。
现转化为:所有可能的字符串权值之和。
假设有k1k_1k1个111段,k0k_0k0个000段。
当起始位置为111的时候,形如1…01\dots 01…0,此时k0=k1k_0=k_1k0=k1。另,xxx个000分成k0k_0k0段,会出现x−k0x-k_0x−k0个000000。所以此时权值为:ax−k1+y−k1×b2k1−1a^{x-k_1+y-k_1}\times b^{2k_1-1}ax−k1+y−k1×b2k1−1,除此,由于是朝xxx轴正方向发射,故还会多出一个010101,即:ax+y−2k1×b2k1a^{x+y-2k_1}\times b^{2k_1}ax+y−2k1×b2k1。
当起始位置为000的时候,形如0…00\dots 00…0,此时k0=k1+1k_0=k_1+1k0=k1+1。算上初始发射时候的000000,权值为同样为ax+y−2k1×b2k1a^{x+y-2k_1}\times b^{2k_1}ax+y−2k1×b2k1。
从这里,我们可以开始考虑枚举k1k_1k1。
如果起始位置为111,总的答案为:Cx−1k1−1×Cy−1k1−1×ax+y−2k1×b2k1C_{x-1}^{k_1-1}\times C_{y-1}^{k_1-1}\times a^{x+y-2k_1}\times b^{2k_1}Cx−1k1−1×Cy−1k1−1×ax+y−2k1×b2k1。因为xxx个000有x−1x-1x−1个空位,从这些空位选取k1−1k_1-1k1−1个截断会得到k1k_1k1个段。
如果起始位置为000,那么答案就是Cx−1k1×Cy−1k1−1×ax+y−2k1×b2k1C_{x-1}^{k_1}\times C_{y-1}^{k_1-1}\times a^{x+y-2k_1}\times b^{2k_1}Cx−1k1×Cy−1k1−1×ax+y−2k1×b2k1。
二者加起来得到最终式子:Cxk1×Cy−1k1−1×ax+y−2k1×b2k1C_{x}^{k_1}\times C_{y-1}^{k_1-1}\times a^{x+y-2k_1}\times b^{2k_1}Cxk1×Cy−1k1−1×ax+y−2k1×b2k1。接下来就枚举k∈[1,min(x,y)]k\in[1,min(x,y)]k∈[1,min(x,y)]。
此时有一个需要注意的,如果每次都logloglog计算会tletletle。需要预处理出facfacfac和invfacinvfacinvfac及invinvinv。同时:记fi=Cxi×Cy−1i−1×ax+y−2i×b2if_i=C_{x}^{i}\times C_{y-1}^{i-1}\times a^{x+y-2i}\times b^{2i}fi=Cxi×Cy−1i−1×ax+y−2i×b2i,那么fi+1=fi×(y−i)×(x−i)×inv(i)×inv(i+1)f_{i+1}=f_i\times (y-i)\times (x-i)\times inv(i)\times inv(i+1)fi+1=fi×(y−i)×(x−i)×inv(i)×inv(i+1),用此式子递推可以避免每次都暴力用快速幂计算ax+y−2ia^{x+y-2i}ax+y−2i和b2ib^{2i}b2i。
#include<bits/stdc++.h>
#define endl '\n'
#define pii pair<int,int>
#define int long long
#define uint unsigned long long
using namespace std;
const double eps = 1e-10;
const int mod = 1e9 + 7;
const int inf = 1e18;
mt19937_64 rng(chrono::steady_clock::now().time_since_epoch().count());
inline int sub(int a, int b) {
a -= b;
a %= mod;
if (a < 0) return a + mod;
else return a;
}
inline int add(int a, int b) {
return (a + b) % mod;
}
inline int mul(int a, int b) {
return a * b % mod;
}
inline int lowbit(int x) {
return x & (-x);
}
inline int get(int x, int pos) {
return (x >> pos) & 1;
}
inline int qp(int a, int b, int p) {
if (b == 0) return 1;
if (a == 0) return 0;
int res = 1;
while (b) {
if (b & 1) {
res = res % p * (a % p) % p;
}
b >>= 1;
a = a % p * (a % p) % p;
}
return res % p;
}
inline uint rdm() {
return rng();
}
const int N = 63;
int dx[] = {0, 1, -1, 0};
int dy[] = {1, 0, 0, -1};
const int NN = 1000000;
vector<int> fac(NN + 4);
vector<int> invfac(NN + 4);
vector<int> inv(NN + 4);
void init() {
fac[0] = 1;
for (int i = 1; i <= NN; i++) {
fac[i] = mul(fac[i - 1], i);
inv[i] = qp(i, mod - 2, mod);
}
invfac[NN] = qp(fac[NN], mod - 2, mod);
for (int i = NN - 1; i >= 0; i--) {
invfac[i] = mul(invfac[i + 1], i + 1);
}
}
int cnm(int n, int m) {
if (m < 0 || m > n) return 0;
return mul(fac[n], mul(invfac[m], invfac[n - m]));
}
void solve() {
int n, m, c, d;
cin >> n >> m >> c >> d;
int x = n + 1, y = m;
int ans = 0;
int t = qp(c + d, mod - 2, mod);
int a = mul(c, t), b = mul(d, t);
int inva = qp(a, mod - 2, mod);
if (y == 0) {
cout << qp(a, x, mod) << endl;
return;
}
if (b == 0) {
cout << 0 << endl;
return;
}
if (a == 0) {
if (x != y) cout << 0 << endl;
else cout << 1 << endl;
return;
}
int ab = mul(mul(b, b), mul(inva, inva));
int l = mul(x, mul(qp(a, x + y - 2, mod), qp(b, 2, mod)));
for (int i = 1; i <= min(x, y); i++) {
ans = add(ans, l);
l = mul(l, sub(x, i));
l = mul(l, sub(y, i));
l = mul(l, mul(inv[i], inv[i + 1]));
l = mul(l, ab);
}
cout << ans << endl;
}
signed main() {
ios::sync_with_stdio(false);
cin.tie(nullptr);
cout.tie(nullptr);
int T = 1;
cin >> T;
init();
while (T--) {
solve();
}
}