<think>题目要求计算n秒后，Cirno被击中的次数（记为C）小于Aya被击中的次数（记为A）的概率，即P(C < A)。注意，每一秒是独立的，且每一秒：
- Aya被击中的概率为p = xp / yp
- Cirno被击中的概率为q = xq / yq
同时，每一秒有四种可能的情况：
1. 两者都不中：概率为 (1-p)*(1-q)
2. 仅Aya中：概率为 p*(1-q)
3. 仅Cirno中：概率为 (1-p)*q
4. 两者都中：概率为 p*q
我们关心的是在n秒后，满足C < A的概率。
由于事件独立，我们可以考虑每一秒对A和C的贡献。设第i秒：
- 定义随机变量X_i: 如果Aya被击中而Cirno未被击中，则X_i=1（对A有贡献）；否则为0。
- 定义随机变量Y_i: 如果Cirno被击中而Aya未被击中，则Y_i=1（对C有贡献）；否则为0。
- 注意：当两者都中或都不中时，对A和C的贡献均为0（因为都中时，两者都加1，抵消了比较？不对，实际上我们比较的是总次数，所以两者都中时，A和C都加1，这样A和C的差值不变。而都不中时，两者都不变）。
但是，我们也可以从另一个角度考虑：每一秒可以认为对A和C的差值有影响。然而，题目要求的是C < A，即A - C >= 1。
实际上，我们可以将问题转化为：在n秒后，A - C >= 1的概率？不对，注意是C < A，即A - C > 0，也就是A - C >= 1。但是，由于每一秒的贡献可能为负（当Cirno中而Aya不中时，A-C会减少1；当Aya中而Cirno不中时，A-C增加1；两者都中或都不中时，A-C不变）。因此，我们可以定义每一秒对差值D（初始为0）的贡献：
- 仅Aya中：D增加1（概率为p*(1-q)）
- 仅Cirno中：D减少1（概率为(1-p)*q）
- 其他两种情况：D不变（概率为p*q + (1-p)*(1-q)）
那么，n秒后，要求D>0（即A-C>0）的概率。
因此，问题转化为：一个随机游走，每一步以概率a增加1，以概率b减少1，以概率c不变（a+b+c=1），求n步后位置大于0的概率。
但是，注意：我们要求的是严格大于0，即D>0。
然而，由于n可能很大（最大10^7），我们需要一个高效的方法。
设a = p*(1-q), b = (1-p)*q, c = 1 - a - b = p*q + (1-p)*(1-q)
我们要求的是n步后位置D>0的概率。
注意：D的初始值为0，第i步的变动为δ_i，δ_i的取值为+1（概率a）、-1（概率b）、0（概率c）。则D = δ_1+δ_2+...+δ_n。
我们要求P(D>0)。
由于每一步独立，我们可以考虑用生成函数或者动态规划，但是n最大10^7，动态规划O(n^2)不可行。
另一种思路：枚举非零步数。设总共有k步是非零的（即发生了增加1或减少1），那么这k步中，有i步是+1，k-i步是-1，则最后位置为i - (k-i) = 2i-k。要求2i-k>0，即i>k/2。
同时，非零步数k的概率为C(n,k) * (a+b)^k * c^(n-k)。而对于固定的k，非零步中+1的步数i服从二项分布，概率为C(k,i) * [a/(a+b)]^i * [b/(a+b)]^{k-i}。
因此，总概率为：
P = Σ_{k=0}^{n} [ C(n,k) * (a+b)^k * c^(n-k) * Σ_{i=floor(k/2)+1}^{k} C(k,i) * (a/(a+b))^i * (b/(a+b))^{k-i} ]
注意：当k=0时，位置为0，不满足>0，所以k>=1。
但是，这个式子需要双重循环，k从0到n，i从floor(k/2)+1到k，而n最大10^7，所以直接计算不可行。
我们需要优化这个表达式。
注意到内层求和是关于i的：Σ_{i=floor(k/2)+1}^{k} C(k,i) * r^i * (1-r)^{k-i}，其中r = a/(a+b)。这是二项分布中成功概率为r时，成功次数大于k/2的概率。这个部分没有闭式解，而且k很大时计算困难。
另一种思路：考虑用生成函数。设F(x) = (c + a*x + b*x^{-1})^n，则我们要求的是所有指数>0的项的系数之和。但是，由于有负指数，我们可以乘以x^n，转化为求多项式(x^0 * c + a*x^1 + b*x^{-1}*x)^n = (c*x + a*x^2 + b)^n 中指数>n的项的系数？不对，重新调整：
原生成函数：G(x) = (b * x^{-1} + c + a*x)^n
要求的是：G(x)中所有指数大于0的项的系数之和。
那么，总概率为：
P = Σ_{d=1}^{n} [x^d] G(x)
但是，由于d的范围很大，而且我们无法展开一个n次多项式（n=10^7），所以需要其他方法。
注意到：整个概率的生成函数为G(x)= (b*x^{-1} + c + a*x)^n，那么我们可以用卷积或者FFT？但是n=10^7，FFT也不可行（O(n log n)在10^7下也很慢，且常数大，且模数不是NTT模数？模数为10^9+7，不是质数？不对，10^9+7是质数，但是NTT需要模数有原根，10^9+7有原根，但是n=10^7，做NTT需要O(n log n)，在10^7下可能超时，而且题目有多个subtask，n最大10^7，T最多多少？题目没有说T，但输入格式中n是直接给的，然后输入概率，所以T=1？因为输入格式第一行一个正整数n，然后两行概率。而题目说“第一行一个正整数n”，所以每组数据只有一组？但是题目没有说T，看样例输入输出：input1只有一组，input2也是一组。所以T=1？再看数据范围：n最大10^7，所以只有一组数据。
但是，题目中并没有说测试组数，而输入格式第一行是n，所以应该只有一组数据。
因此，我们只需要处理一组数据，但是n最大10^7。
现在，我们要求的是多项式G(x)= (a*x + c + b/x)^n 中所有正次幂的系数和。注意，这里x是形式变量。
我们可以将G(x)写成：
G(x) = x^{-n} * (a*x^2 + c*x + b)^n
那么，原来的多项式G(x)中x^k的系数等于多项式P(x)= (a*x^2 + c*x + b)^n 中x^{k+n}的系数。
因此，G(x)中指数大于0的项，即k>0，那么
对应P(x)中指数大于n的项（因为k>0 => k+n>n）。
所以，我们要求的是多项式P(x)= (a*x^2 + c*x + b)^n 中，次数大于n的项的系数之和。
设P(x) = Σ_{i=0}^{2n} coeff[i] * x^i，则我们要求的是：Σ_{i=n+1}^{2n} coeff[i]
而整个多项式P(x)的系数和等于P(1) = (a+b+c)^n = (a+b+c)^n，但注意a+b+c=1，所以P(1)=1。
那么，我们只需要知道Σ_{i=0}^{n} coeff[i] 和 Σ_{i=n+1}^{2n} coeff[i] 的和为1，但是我们需要的是后半部分。
然而，我们还可以考虑用另一种对称性：注意P(x)是一个多项式，而我们可以用多项式取模技术？或者利用生成函数的性质。
注意：我们要求的概率为P(D>0)，由对称性，我们还可以考虑P(D<0)和P(D=0)。但是这里a和b不一定相等，所以不对称。
另一种思路：考虑用动态规划结合生成函数，但n太大。
我们回到之前的双重求和，并尝试用组合恒等式化简。
设r = a/(a+b)，则内层求和为：Σ_{i=floor(k/2)+1}^{k} C(k,i) * r^i * (1-r)^{k-i} = 1 - F_{r,k}(floor(k/2))，其中F_{r,k}是二项分布累积分布函数。这个值没有简单的闭式，但是我们可以利用递推？或者用中心极限定理近似？但题目要求精确值取模。
因此，我们需要一个O(n)或O(n log n)的方法。
我们考虑直接递推：设dp[i][j]表示前i秒，差值为j的概率。但是j的范围是[-i,i]，空间和时间都是O(n^2)，n=10^7显然不行。
我们考虑生成函数的系数提取。我们要求的是：
S = Σ_{d=1}^{n} [x^d] (b*x^{-1} + c + a*x)^n
= Σ_{d=1}^{n} [x^{n+d}] (b + c*x + a*x^2)^n （因为乘以了x^n）
所以，S = Σ_{d=n+1}^{2n} [x^d] (b + c*x + a*x^2)^n
而多项式Q(x)=b + c*x + a*x^2是一个二次多项式，它的n次幂我们可以用生成函数的方法，但直接展开系数有公式吗？
我们可以用二项式定理展开，但是这里是三项，所以用多项式定理：
Q(x)^n = (b + c*x + a*x^2)^n = Σ_{i+j+k=n} [n!/(i!j!k!)] * b^i * c^j * a^k * x^{j+2k}
其中，i,j,k>=0，且i+j+k=n。
那么，x^m的系数为：Σ_{j+2k=m, i+j+k=n} [n!/(i!j!k!)] * b^i * c^j * a^k
这里，i = n - j - k，所以系数为：
coeff[m] = Σ_{j,k: j+2k=m, j>=0, k>=0, j+k<=n} [n!/( (n-j-k)! j! k! )] * b^{n-j-k} * c^j * a^k
那么，我们要求的是m从n+1到2n的系数和。
但是，m=j+2k，而m的范围是n+1到2n，即j+2k>=n+1。
同时，j,k满足：j>=0, k>=0, j+k<=n。
所以，S = Σ_{m=n+1}^{2n} coeff[m] = Σ_{j,k: j+2k>=n+1, j>=0,k>=0,j+k<=n} [n!/( (n-j-k)! j! k! )] * b^{n-j-k} * c^j * a^k
这个式子仍然需要枚举j,k，而j,k的范围都是0~n，所以O(n^2)不可行。
我们需要优化这个求和。注意，这个求和可以写成：
S = Σ_{j=0}^{n} Σ_{k=0}^{n} [条件：j+2k>=n+1, j+k<=n] [n!/( (n-j-k)! j! k! )] * b^{n-j-k} * c^j * a^k
改变枚举顺序：先枚举k，再枚举j，但仍然是O(n^2)。
另一种思路：考虑用容斥，用总概率1减去D<=0的部分，即P(D>0)=1-P(D<=0)=1-P(D<0)-P(D=0)。
P(D=0)是比较好求的：就是上述三重求和（多项式系数）中m=0的部分？不对，注意我们调整后的多项式是P(x)= (a*x^2+c*x+b)^n，而D>0
对应的是x^{n+d}（d>=1），而D=0
对应的是x^n（因为调整后指数为n+d，d=0时指数为n）。所以P(D=0) = [x^n] (a*x^2+c*x+b)^n。
而P(D<0)
对应的是d<0，即调整后指数小于n，即m<n（在P(x)中指数为m，m<n）。
所以，P(D>0) = 1 - P(D<0) - P(D=0)
那么，问题转化为求P(D<=0) = P(D<0)+P(D=0) = Σ_{m=0}^{n} [x^m] (b + c*x + a*x^2)^n
即求多项式Q(x)=b+c*x+a*x^2的n次幂的前n+1项系数和。
现在，我们要求一个多项式的前n+1项系数和。注意，这个多项式是二次的，我们可以用递推的方法来求系数。
设F(n) = Q(x)^n = (b + c*x + a*x^2)^n，我们想求F(n)的前n+1项系数和。
考虑生成函数的递推关系。但是，我们也可以利用线性递推：因为Q(x)是二次的，所以F(n)的系数满足一个线性递推关系。
具体地，设F(n) = Σ_{m=0}^{2n} f(m) * x^m。
对F(n)求导：F(n)' = n * (c + 2a*x) * (b + c*x + a*x^2)^{n-1} = n * (c + 2a*x) * F(n-1)
所以，F(n)' = n * c * F(n-1) + 2n a x F(n-1)
即：Σ_{m=1}^{2n} m * f(m) * x^{m-1} = n*c * Σ_{m=0}^{2n-2} f_{n-1}(m) x^m + 2n a * Σ_{m=0}^{2n-2} f_{n-1}(m) x^{m+1}
调整指数：
左边 = Σ_{m=0}^{2n-1} (m+1) * f_{n}(m+1) * x^m
右边 = n*c * Σ_{m=0}^{2n-2} f_{n-1}(m) x^m + 2n a * Σ_{m=1}^{2n-1} f_{n-1}(m-1) x^m
比较x^m的系数，对于m>=0：
(m+1) * f_n(m+1) = n*c * f_{n-1}(m) + 2n a * f_{n-1}(m-1) [当m>=1时]
当m=0时：
1 * f_n(1) = n*c * f_{n-1}(0)
但是，这个递推需要f_{n-1}，而且n很大，我们无法递推所有系数。
我们要求的是F(n)的前n+1项系数和，即S(n) = Σ_{m=0}^{n} f(m)。
考虑S(n) = Σ_{m=0}^{n} [x^m] F(n)(x)
有没有生成函数的方法？比如用积分或赋值。
我们赋值：如果我们能求出一个封闭形式就好了。
注意：我们要求的是多项式在x=1处的前缀和？即Σ_{m=0}^{n} f(m) = ? 但这不是生成函数在x=1的值，而是生成函数在x=1时，用某种算子作用的结果。
一种方法是利用生成函数的积分：∫(F(n)(x)) from 0 to 1? 但这和系数和无关。
另一种思路：构造一个辅助生成函数：G(x) = F(n)(x) * (1+x+x^2+...+x^n) = F(n)(x) * (1-x^{n+1})/(1-x)
则S(n) = G(1)？不对，当x=1时，G(1) = F(n)(1) * (n+1) = (a+b+c)^n * (n+1) = 1 * (n+1) = n+1，但这显然不是我们想要的，因为F(n)(1)=1，但S(n)是前n+1项系数和，而F(n)(x)的系数和是1，所以S(n)<=1，而n+1>1，矛盾。
我们要求的是系数和，不是值的和。
所以，我们需要用其他方法。
我们考虑用离散积分：设H(x) = Σ_{m=0}^{n} (Σ_{k=0}^{m} f(k)) x^m，那么我们要求H(1) = S(n)？不对，H(1) = Σ_{m=0}^{n} (Σ_{k=0}^{m} f(k)) = Σ_{m=0}^{n} (n-m+1) f(m) ? 不是S(n)。
S(n) = Σ_{m=0}^{n} f(m) = [x^n] (F(n)(x) * (1+x+...+x^n)) = [x^n] (F(n)(x) * (1-x^{n+1})/(1-x))
因此，S(n) = [x^n] [F(n)(x) / (1-x)] - [x^n] [F(n)(x) * x^{n+1}/(1-x)] = [x^n] [F(n)(x) / (1-x)] 因为第二项x^{n+1}/(1-x)展开后x^n的系数为0。
所以，S(n) = [x^n] [F(n)(x) / (1-x)]
现在，[x^n] [F(n)(x) / (1-x)] = Σ_{k=0}^{n} f(k) [ 1/(1-x) = Σ_{i=0}^{infty} x^i, 所以F(n)(x)/(1-x)的x^n系数为Σ_{k=0}^{n} f(k) ]
因此，S(n) = [x^n] [ (b+c*x+a*x^2)^n / (1-x) ]
现在，我们要求一个生成函数的x^n系数，这个生成函数是两个函数的乘积：U(x)= (b+c*x+a*x^2)^n, V(x)=1/(1-x)。
由卷积，[x^n] (U(x)*V(x)) = Σ_{k=0}^{n} [x^k] U(x) * [x^{n-k}] V(x) = Σ_{k=0}^{n} f(k) * 1 = S(n)， tautology.
这并没有给出新的表达式。
我们考虑用生成函数的封闭形式。 (b+c*x+a*x^2)^n / (1-x) = (b+c*x+a*x^2)^n * Σ_{i=0}^{infty} x^i
所以，[x^n] = Σ_{i=0}^{n} [x^{n-i}] (b+c*x+a*x^2)^n = Σ_{i=0}^{n} f(n-i) [f(m) is [x^m] U(x)]
也就是S(n) = Σ_{i=0}^{n} f(i) = S(n) again.
所以，我们回到原点。
我们可能需要一个线性的递推
for S(n) as a whole.
notice that the polynomial U(x) = (b+c*x+a*x^2) is of degree 2, and we are convolving it with 1/(1-x) and taking the n-th coefficient. This is the same as the sum of the first n+1 coefficients of U(x), which is what we want.
then, we can use the following: consider the generating function
for the partial sums.
in fact, we have the matrix
for the
linear recurrence of the coefficients of U(x) is known, but the summation also has a recurrence.
we can use the following identity: if we let S(n,k) = sum_{m=0}^{k} [x^m] (b+c*x+a*x^2)^n
for fixed n, then we are to compute S(n,n).
and
for the generating function, we have a recurrence on n.
but n is up to 10^7, and k is up to n, so we cannot do a 2D state.
we need to find a recurrence
for the sum.
after some research, we might use the
method of
for the sum of coefficients of a polynomial.
in fact, we have:
(b+c*x+a*x^2)^n = F(n)(x)
then, the sum of the first n+1 coefficients can be written as:
S(n) = (1/(2πi)) ∫_C F(n)(z) * (1-z^{n+1}) / (1-z) dz / z^{0} ? not helpful.
we can try to use generating function in another way.
we can use the representation by the roots. let the roots of the quadratic equation a*x^2+c*x+b=0 be r1, r2.
then, F(n)(x) = (a*x^2+c*x+b)^n = a^n (x-r1)^n (x-r2)^n.
then, the sum of the coefficients from 0 to n might be expressed by the roots, but it is not trivial.
but we can use the following: the sum of the first n+1 coefficients of a polynomial P(x) is P(1) if we set the higher coefficients to 0, which is not the case.
a better idea: we can use the fact that the polynomial (a*x^2+c*x+b)^n is a
linear combination of (x-r1)^{n} and (x-r2)^{n} if we use partial fractions, but it is a product.
we write:
(a*x^2+c*x+b)^n = a^n (x-r1)^n (x-r2)^n
then, we want to compute the sum_{k=0}^{n} [x^k] (a*x^2+c*x+b)^n = a^n * [x^n] ( (x-r1)^n (x-r2)^n * (1+x+...+x^n) ) ? not exactly.
actually, we want the sum of the coefficients of x^0 to x^n, which is the same as the constant term of (a*x^2+c*x+b)^n * (1+x+...+x^n) * x^{-n} (if we use the residue and.in.
a more direct
method: use that the summation of the first n+1 coefficients can be done by
(1/2) [F(n)(1) + something] but not.
we can use the following integral from -1 to 0 of F(n)(x) dx or something, but it's not discrete.
after thinking, a known approach
for this kind of problem is to use the and then use the prefix sum by a
linear recurrence.
in fact, the coefficients of (a*x^2+c*x+b)^n satisfy a
linear recurrence of length 2.
specifically,
for a polynomial of degree 2, the coefficients of its n-th power satisfy a recurrence of order 2.
indeed,
for a polynomial with two roots, we have a recurrence.
assume the roots are r1, r2, then we have the generating function as a^n (x-r1)^n (x-r2)^n.
then, the coefficient of x^m in this is a^n times the coefficient of x^m in (x^2 - (r1+r2)x + r1r2)^n.
and
for a polynomial of the
form (x^2 + p x + q)^n, the coefficients satisfy a recurrence.
in fact,
for the polynomial (x^2 + p x + q)^n, we have:
(m+1)u_{m+1} = n p u_m + 2n q u_{m-1} - (n+1) (p u_{m+1} + 2q u_m) ? not exactly.
we can use:
(x^2 + p x + q) * U'(x) = n * (2x+p) * U(x)
then, as be
fore, we have a recurrence
for u_m.
from the earlier differential equation:
U'(x) = n (2x+p) U(x) / (x^2+p x+q)
then, (x^2+p x+q) U'(x) = n (2x+p) U(x)
then, in coefficients:
(x^2 U') = n (2x+p) U(x) - p x U'(x) - q U'(x)
and then take the coefficient of x^{m} on both left and right.
left: (x^2 U') = Σ_{m} m u_m x^{m+1} = Σ_{m} (m-1) u_{m-1} x^m
right: n (2x+p) U(x) = n (2 Σ_{m} u_m x^{m+1} + p Σ_{m} u_m x^m) = n (2 Σ_{m} u_{m-1} x^m + p Σ_{m} u_m x^m)
and -p x U'(x) = -p x * (Σ_{m} m u_m x^{m-1}) = -p Σ_{m} m u_m x^m
and -q U'(x) = -q Σ_{m} m u_m x^{m-1} = -q Σ_{m} (m+1) u_{m+1} x^m
so, in total,
for x^m (m>=0) on the right:
n (2 u_{m-1} + p u_m) - p m u_m - q (m+1) u_{m+1}
and on the left: (m-1) u_{m-1} (
for m>=1)
so,
for m>=1:
(m-1) u_{m-1} = n (2 u_{m-1} + p u_m) - p m u_m - q (m+1) u_{m+1}
then, we can solve
for u_{m+1}:
q (m+1) u_{m+1} = n (2 u_{m-1} + p u_m) - p m u_m - (m-1) u_{m-1}
= (2n - (m-1)) u_{m-1} + (n p - p m) u_m
so, u_{m+1} = [ (2n - m+1) u_{m-1} + p (n - m) u_m ] / [ q (m+1) ]
and this recurrence allows us to compute the coefficients in O(n) time, provided we have initial conditions.
initial conditions
for U(x) = (x^2+p x+q)^n:
u_0 = [x^0] = (q)^n ( when n>=0)
u_1 = [x^1] = n * p * q^{n-1} ( because the only way to get x is to take the p*x term in one factor and q in the others)
and
for m>2n, u_m=0.
so, we can compute the coefficients from m=0 to m=2n in O(n) time.
then, the sum S = sum_{m=0}^{n} u_m can be computed by a simple loop from 0 to n.
the whole algorithm
for the one data case would be:
read n
read xp, yp -> p = xp * inv(yp) mod mod
read xq, yq -> q = xq * inv(yq) mod mod
then compute:
a = p * (1-q) % mod [ probability only Aya ]
b = (1-p) * q % mod [ probability only Cirno ]
c = 1 - a - b [ probability either both or none, but note: a+b+c=1 ]
but wait, we defined earlier:
a = p*(1-q) [
for only Aya hit ]
b = (1-p)*q [
for only Cirno hit ]
c = 1 - a - b [
for both or none ]
and then we built the polynomial: U(x) = a*x^2 + c*x + in the probability
for the difference?
however, in the polynomial
for the difference, we had the generating function per second: b*x^{-1} + c + a*x, and then we adjusted to (b + c*x + a*x^2) * x^{-1}?
but in our later adjustment, we defined the polynomial
for the adjusted generating function as: b + c*x + a*x^2.
so, the polynomial is: P(x) = a*x^2 + c*x + b.
so, in the recurrence
for the coefficients of (a*x^2+c*x+b)^n, we use the recurrence with:
p = c / a ? -> wait, in the recurrence above, we wrote the polynomial as (x^2 + p x + q) which is not the case.
actually, we should scale the polynomial to have the leading coefficient 1.
in the recurrence, we assumed the polynomial is x^2 + p x + q, but here it is a*x^2 + c*x + b.
so, we can write: a*x^2 + c*x + b = a (x^2 + (c/a) x + (b/a)) = a (x^2 + p x + q) where p = c/a, q = b/a.
then, (a*x^2+c*x+b)^n = a^n (x^2 + p x + q)^n.
then, the recurrence is
for the coefficients of (x^2 + p x +)^n, say v_m = [x^m] (x^2+p x+q)^n, then the coefficients of the original are u_m = a^n * v_m.
and the recurrence
for v_m is:
v_{m} = [ (2n - (m-2) +1) v_{m-2} + p (n - (m-1)) v_{m-1} ] / [ q * (m) ] ? from the recurrence we had
for u_{m+1} in terms of u_{m-1} and u_{m}:
in the recurrence
for the coefficients of (x^2+p x+q)^n, we had:
q * (m+1) * v_{m+1} = (2n - m+1) * v_{m-1} + p * (n - m) * v_{m}
and this is
for m>=1.
initial conditions
for (x^2+p x+q)^n:
v_0 = q^n
v_1 = n * p * q^{n-1}
and
for m>2n, v_m=0.
so, we can compute the array v
for m from 0 to 2n in O(n) time.
then, the coefficients of the original polynomial are u_m = a^n * v_m.
then, the sum of the first n+1 coefficients is:
S = sum_{m=0}^{n} u_m = a^n * sum_{m=0}^{n} v_m
then, the probability we want is P(D>0) = 1 - (S + something) wait, recall:
we wanted P(D<=0) = the sum of the coefficients of x^m
for m in [0, n] in the polynomial (a*x^2+c*x+b)^n.
and then the answer
for the probability that C<A is P(D>0) = 1 - P(D<=0) = 1 - S.
but wait, is that correct?
remember: D = A - C, and we want P(D>0) = 1 - P(D<=0).
so, if we let S = sum_{m=0}^{n} u_m = sum_{m=0}^{n} [x^m] (a*x^2+c*x+b)^n = P(D<=0)
then, answer = 1 - S.
however, note that the polynomial (a*x^2+c*x+b)^n has degrees from 0 to 2n, and the sum of all coefficients is (a+b+c)^n = 1^n = 1.
and the sum
for m>n is the rest, so P(D>0) = 1 - S.
but is that right?
P(D<=0) = sum_{m=0}^{n} [x^m] (a*x^2+c*x+b)^n = S.
then P(D>0) = 1 - S.
however, note that in our earlier adjustment, the original generating function
for the difference was
for the exponent of the difference in the range [-n, n], and we have a polynomial of degree 2n, and the coefficients sum to 1.
so, indeed, the probability of every possible difference is accounted
for.
there
fore, the answer is (1 - S) mod (10^9+7).
steps
for the algorithm:
mod = 10^9+7
step 1: read n.
step 2: read xp, yp; then p = xp * inv(yp) mod mod.
read xq, yq; then q = xq * inv(yq) on mod.
step 3:
a = p * (1 - q) % mod; // probability only Aya
b = (1 - p) * q % mod; // probability only Cirno
// note: both or neither: c = 1 - a - b, but we also need to use it in the polynomial: the middle coefficient is c.
// and the polynomial is: a*x^2 + c*x + b, where c = 1 - a - b.
c = (1 - a - b) % mod;
if (c<0) c += mod;
step 4: if a==0, then the polynomial is c*x+b, then (c*x+b)^n is easy: but we can still use the recurrence or direct. to be safe, we use the recurrence
for the coefficients of the polynomial (x^2 + p x + q) scaled by a.
let A = a, B = c, C = b. // the polynomial is A*x^2+B*x+C.
if A==0:
// then the polynomial is B*x+C, then the n-th power is: (B*x+C)^n = Σ_{k=0}^{n} C(n,k) * (B)^k * (C)^{n-k} * x^k.
// then, the sum of coefficients from x^0 to x^n is the entire polynomial since degree=n, so S = (B+C)^n = (c+b)^n = (1-a)^n? but wait, a might be not zero even if A=0?
// actually, if a=0 then A=0, and then the polynomial is B*x+C = c*x+b.
// and the sum of the coefficients from x^0 to x^n is the entire sum: (c+b)^n, because there are no terms beyond x^n.
// then, answer = 1 - (c+b)^n.
// but note: a=0 means p*(1-q)=0, then either p=0 or q=1.
// then, the entire probability difference: only Aya hit never happens.
// then, the only way to have a positive difference is only Cirno hit, but wait, if a=0, then only Aya hit never happens, so the difference can only be non-positive.
// in fact, if a=0, then the only way to increase the difference is gone, so D can only be<=0, so P(D>0)=0.
// so answer=0.
// but also, (c+b) = (1-a) = 1, so (c+b)^n=1, then 1-1=0. so it holds.
similarly, if A==0 and B==0, then polynomial is C, then (C)^n, then S=C^n, answer=1-C^n.
but we can handle in recurrence.
generally, we let:
if A==0:
then the polynomial is of degree 1, and we can avoid the recurrence
for per
formance? but n can be 10^7, and the recurrence is O(n), so we can do the recurrence even
for a
linear polynomial? in the recurrence
for a quadratic, if A=0, then we have a
linear polynomial, and the recurrence might not hold because we assumed A!=0.
so we need to handle the case A=0 separately.
if A==0:
// then (B*x+C)^n, the sum_{m=0}^{n} [x^m] = (B+C)^n, because it is the evaluation at x=1 of the polynomial, and there are only terms from 0 to n.
S = pow(B+C, n, mod)
answer = (1 - S) % mod
else:
// let in the scaled polynomial: we factor out A: = A^n * (x^2 + (B/A) x + (C/A))^n.
// then, the coefficients of the scaled polynomial (x^2 + p x + q)^n, where p = B/A, q = C/A, in the ring mod.
p_val = B * inv(A, mod) % mod
q_val = C * inv(A, mod) % mod
// then, we want to compute the coefficients v[0..2n]
for (x^2+p_val*x+q_val)^n.
// we only need from 0 to n, and then sum_{m=0}^{n} v_m.
// recurrence
for v_m:
// v_0 = q_val^n
// v_1 = n * p_val * q_val^{n-}
//
for m from 2 to 2n:
// v_m = [ (2n - m+1) * v_{m-2} + p_val*(n - m+1) * v_{m-1} ] * inv(q_val * m, mod) % mod
// but wait, our recurrence was:
// q_val * (m) * v_{m} = (2n - (m-2)+1) * v_{m-2} + p_val * (n - (m-1)) * v_{m-1}
// =>
for m>=2:
// v_m = [ (2n - m+3) * v_{m-2} + p_val * (n - m+1) * v_{m-1} ] * inv(q_val * m, mod) % mod
// let me re-express the recurrence:
// from:
// q_val * (m) * v_{m} = (2n - (m-2) +1) * v_{m-2} + p_val * (n - (m-1)) * v_{m-1}
// = (2n - m+3) * v_{m-2} + p_val * (n - m+1) * v_{m-1}
// so,
for m>=2:
// v_m = [ (2n-m+3) * v_{m-2} + p_val*(n-m+1)*v_{m-1} ] * inv(q_val * m, mod) % mod
// however, note that the recurrence we derived was
for v_{m+1} in terms of v_{m} and v_{m-1}:
// q_val * (m+1) * v_{m+} = (2n - m+1) * v_{m-1} + p_val * (n-m) * v_{m}
// so, to compute v_{m}
for m from 0 to 2n, we can do:
// v[0] = pow(q_val, n, mod)
// if n==0: then v[1]=0, else v[1] = n * p_val * pow(q_val, n-1, mod) % mod
//
for m from 2 to 2n:
// v[m] = ( (2n - (m-2)+1) * v[m-2] + p_val * (n - (m-1)) * v[m-1] ) * inv(q_val * m, mod) % mod
// = ( (2n - m+3) * v[m-2] + p_val * (n-m+1) * v[m-1] ) * inv(q_val * m, mod) % mod
// but note: when m>2n, v_m=0, and we only need up to m=n.
// then, the sum of v[0..n] is let's call it T.
// then, the sum
for the original polynomial is: S = pow(A, n, mod) * T % mod.
// then, answer = (1 - S) % mod.
step 5: output answer.
but note: the recurrence
for m from 0 to n requires O(n) time, and n=10^7, which might be borderline in C++ but in Python might be slow, but the problem says data range up to 10^7, and we are in C++? the user asked
for C++ with very short variable names? but we are in Python.
however, the problem does not specify language, so we try to write in Python and hope that it passes in PyPy or C++ would be required
for 10^7.
but the user said "用c++变量名长度不超过2", so we are to write in C++.
but the user's first sentence is "请使用中文回答问题", and then the problem is in Chinese, and then the last constraint is "用c++变量名长度不超过2", so it requires C++ code with very short variable names (<=2 characters).
so we will write in C++.
however, the recurrence is O(n) and n=10^7, which is acceptable in C++.
but we must use modular arithmetic and precompute the modular inverses
for the
numbers in the recurrence? the recurrence requires division by q_val * m, and m is from 1 to 2n, so we need the modular inverse
for each number in the range [1, 2n] (which is 2e7) and also
for q_val.
we can precompute an array of inverses
for numbers from 1 to 2n, but 2n=20e6, which is acceptable in C++.
steps in C++:
const int mod = 1000000007;
int n;
cin >> n;
long long xp, yp, xq, yq;
cin >> xp >> yp >> xq >> yq;
// compute p = xp/yp mod mod, q=xq/yq mod mod.
auto inv = [mod](long long a) {
return pow(a, mod-2, mod); // but we need to do modular exponentiation
};
long long p = xp * inv(yp) % mod;
long long q = xq * inv(yq) % mod;
long long a = p * (1 - q + mod) % mod; // probability only Aya
long long b = (1 - p + mod) * q % mod; // probability only Cirno
long long c = (1 - a - b) % mod;
if (c<0) c += mod;
// the polynomial is: a*x^2 + c*x + b.
if (a==0) {
// then the polynomial is c*x + b.
long long S = pow(b+c, n, mod); // because the entire polynomial is (b+c*x), then the sum of coefficients from x^0 to x^n is (b+c)^n.
long long ans = (1 - S + mod) % mod;
cout << ans << endl;
} else {
// compute p_val = c / a, q_val = b / a, in mod.
long long ia = inv(a);
long long p_val = c * ia % mod;
long long q_val = b * ia % mod;
// we will compute the coefficients v[0] to v[2n]? but we only need up to n, and then sum_{m=0}^{n} v_m.
// but the recurrence
for v_{m} requires v_{m-1} and v_{m-2}
for m>=2, and we only need m from 0 to n, so we only need to compute v[0] to v[n].
// Precompute an array of inverses
for integers from 1 to n (since in the recurrence we have m from 1 to n, and we need inv(m)
for each m, and also inv(q_val) is known.
// But note: the recurrence has division by (q_val * m), so we need inv(q_val * m) = inv(q_val) * inv(m) mod mod.
// Precompute an array inv_m
for m in 1..n (n<=10^7) -> we can do a
linear sieve for inverses, or use a loop with 1..n and use the
formula: inv[i] = (mod - mod/i) * inv[mod%i] % mod.
// However, we only need up to n, not 2n.
// Steps
for recurrence:
vector<long long> v(n+1, 0); // v[0] to v[n]
// v[0] = q_val^n
v[0] = pow(q_val, n, mod);
if (n>=1) {
// v[1] = n * p_val * q_val^(n-1) % mod
v[1] = n % mod * p_val % mod * pow(q_val, n-1, mod) % mod;
}
// Precompute the array of inverses
for integers from 1 to n:
vector<long long> inv_i(n+1);
inv_i[1] = 1;
for (int i=2; i<=n; i++) {
inv_i[i] = (mod - mod/i) * inv_i[mod % i] % mod;
}
long long inv_q = inv(q_val);
// Now,
for m from 2 to n:
for (int m=2; m<=n; m++) {
// recurrence:
// v[m] = [ (2n - m+3) * v[m-2] + p_val * (n-m+1) * v[m-1} ] * inv(q_val * m) % mod
// = [ (2n - m+3) * v[m-2] + p_val * (n-m+1) * v[m-1} ] * inv_q * inv_i[m] % mod
long long term1 = (2*n - m + 3) % mod * v[m-2] % mod;
long long term2 = p_val * ((n - m + 1) % mod) % mod * v[m-1] % mod;
v[m] = (term1 + term2) % mod;
v[m] = v[m] * inv_q % mod;
v[m] = v[m] * inv_i[m] % mod;
}
// Now, T = sum_{m=0}^{n} v[m] mod mod.
long long T = 0;
for (int i=0; i<=n; i++) {
T = (T + v[i]) % mod;
}
// Then, the sum
for the original polynomial is:
long long A_n = pow(a, n, mod);
long long S = A_n * T % mod;
long long ans = (1 - S + mod) % mod;
cout << ans << endl;
}
however, note: the recurrence may be numerically unstable
for large n and mod, but since we are in modular arithmetic with a
prime mod, and the recurrence is
linear, it should be fine.
but we must be cautious about negative
numbers in modular arithmetic.
also, we have to take care of the range: n=10^7, and we are creating a vector of size n+1 (10^7+1) which is about 80 MB in C++, which is acceptable.
but the problem says n up to 10^7, and we are in the worst case, so we must be efficient.
we can avoid storing the whole array if we only need the sum T. we only need v[m-1] and v[m-2] to compute v[m], so we can use two variables to store the last two and accumulate T.
so, we can do:
v0 = q_val^n
T = v0
if n>=1:
v1 = n * p_val * q_val^{n-1} % mod
T = (T + v1) % mod
then
for m from 2 to n:
v2 = [ (2n-m+3) * v0 + p_val*(n-m+1)*v1 ] * inv(q_val * m) % mod
T = (T + v2) % mod
v0 = v1
v1 = v2
then we don't need the array.
so, we can change to:
long long v0 = pow(q_val, n, mod);
long long T = v0;
long long v1 = 0;
if (n>=1) {
v1 = n % mod * p_val % mod;
v1 = v1 * pow(q_val, n-1, mod) % mod;
T = (T + v1) % mod;
}
// Precompute the array of inverses
for integers from 1 to n? but we only need
for m from 2 to n, so we can compute the inverses on the fly using the
formula:
for a given m, inv_m = pow(m, mod-2, mod) is too slow.
// But we need to do it
for each m, and there are 10^7 ms, and pow(m, mod-2, mod) is O(log mod) which is about 30 iterations, so total 10^7 * 30 = 300e6, which might be borderline in C++.
// Alternatively, we can precompute the inverses
for 1..n in a separate array, but that is 80 MB, and then we can use a pointer to it? or we can do a
sieve for inverses in a loop without storing the whole array? but we need to store them
for the recurrence? no, we can compute on the fly with the recursive
formula:
// inv[i] = (mod - mod/i) * inv[mod%i] % mod, and we only need one at a time.
// We can do:
// long long inv_m = 1;
//
for (int m=2; m<=n; m++) {
// inv_m = (mod - mod/m) * inv_m % mod; // no, the
formula is
for the inverse of m, but it requires the inverse of mod%i, which we haven't computed.
// Instead, we can precompute an array of size n+1
for the inverses, but we want to avoid the memory? 80 MB is acceptable in C++.
// But the problem says n up to 10^7, and in C++ we can allocate 80 MB.
// However, we can avoid precomputation by using the
formula with a temporary variable
for the current inverse, but the recurrence
for the inverses
for integers 1..n can be computed in a loop without storing all if we only need the current one. But the
formula:
// inv_i = (mod - mod / i) * inv_i[mod % i] % mod;
// requires the inverse of mod%i, which is less than i, and we have computed it? but we are not storing them.
// So we must either store an array of size n+1
for the inverses, or compute each with Fermat's little theorem (which is O(log mod)).
// Since mod is fixed (10^9+7), and we need n=10^7, we can do:
//
for m from 2 to n:
// long long inv_m = pow(m, mod-2, mod); // this is O(log mod) = about 30 iterations per call, so total 10^7 * 30 = 300e6 iterations, which is about 3 seconds in C++ if optimized.
// or we can precompute an array of inverses in a single pass in O(n) time and O(n) space.
// Given the constraints, we choose to precompute the array of inverses
for 1..n.
// But then we also need the inverse of q_val once.
// So, we do:
if (n>=2) {
// Precompute inv_arr
for 1 to n:
vector<long long> inv_arr(n+1);
inv_arr[1] = 1;
for (int i=2; i<=n; i++) {
inv_arr[i] = (mod - mod/i) * inv_arr[mod % i] % mod;
if (inv_arr[i] < 0) inv_arr[i] += mod;
}
long long inv_q_val = inv(q_val); // only once
long long v2 = 0;
for (int m=2; m<=n; m++) {
// v2 = [ (2n - m + 3) * v0 + p_val * (n - m + 1) * v1 ] * inv_q_val * inv_arr[m] % mod;
long long term1 = (2*n - m + 3) % mod;
term1 = term1 * v0 % mod;
long long term2 = p_val * ( (n - m + 1) % mod ) % mod;
term2 = term2 * v1 % mod;
v2 = (term1 + term2) % mod;
v2 = v2 * inv_q_val % mod;
v2 = v2 * inv_arr[m] % mod;
T = (T + v2) % mod;
// update v0, v1
for next iteration
v0 = v1;
v1 = v2;
}
}
// then continue as be
fore.
but note: the terms (2n-m+3) and (n-m+1) might be negative? we do % mod, but we should take modulo and then adjust negative.
and also, the recurrence might not be 100% stable
for large n, but mathematically it is.
finally, we compute S = A^n * T, and then answer = 1 - S.
let's test with the sample: n=2, xp=1, yp=2, xq=1, yq=2.
then p = 1/2, q=1/2.
a = p*(1-q) = (1/2)*(1/2)=1/4.
b = (1-p)*q = (1/2)*(1/2)=1/4.
c = 1 - a - b = 1/2.
then the polynomial = (1/4) x^2 + (1/2) x + 1/4.
then, (1/4 x^2 + 1/2 x + 1/4)^2 = (1/16) (x^2 + 2x+1)^2 = (1/16)(x+1)^4.
= (1/16)(x^4 + 4x^3+6x^2+4x+1)
then, the sum of coefficients
for x^0 to x^2 (m=0,1,2) is: 1/16 * (1 + 4 + 6) = 11/16.
then, P(D>0) = 1 - 11/16 = 5/16.
but the sample output is 812500006.
and 5/16 = 0.3125, and 0.3125 * 1000000007 mod is not 812500006.
wait, 5/16 mod 10^9+7 = 5 * inv(16) mod (10^9+7) = 5 * 562500004 % mod = 2812500020 % mod = 2812500020 % 1000000007 = 2812500020 - 2*1000000007 = 2812500020 - 2000000014 = 812500006.
so it is 5/16.
but in our calculation, we have S = 11/16, then 1-S=5/16.
so the sample is correct.
now, in the code
for the sample:
a = 1/4 mod mod = (1 * inv(4)) % mod = 250000002
similarly, b = 250000002, c=500000004.
then in the else branch:
p_val = c / a = 500000004 * inv(250000002) % mod = 2.
q_val = b / a = 250000002 * inv(250000002) = 1.
then v0 = q_val^n = 1^2 = 1.
v1 = 2 * 2 * 1^(2-1) = 4. [n=2, p_val=2, so v1= n * p_val * q_val^{n-1} = 2*2*1=4]
then
for m=2:
v2 = [ (2*2 - 2 + 3) * v0 + p_val*(2-2+1)*v1 ] * inv(q_val * 2)
= [ (4-2+3)=5 * 1 + 2*(1)*4 ] * inv(2) [because q_val=1, so 1*2=2]
= (5+8) * inv(2) = 13 * inv(2) = 13 * 500000004 = 6500000052 mod mod = 6500000052 % 1000000007 = 6500000052 - 6*1000000007 = 6500000052 - 6000000042 = 500000010, then mod = 500000010.
but the expected v2 should be 6 (from the expansion: 1+4+6), but wait, in the scaled polynomial we have (x^2+2x+1)^2 = x^4+4x^3+6x^2+4x+1, so v0=1, v1=4, v2=6.
but we have v2=500000010, which is not 6.
what is inv(2) in mod? 500000004.
then 13 * 500000004 = 6500000052, then 6500000052 % 1000000007 = 6500000052 - 6*1000000007 = 6500000052-6000000042=500000010, which is not 6.
why? because the recurrence is
for the scaled polynomial (x^2+2x+1)^2, and we are not scaling by a^n at the end.
in the recurrence, we are computing the coefficients of (x^2+2x+1)^2, which is 1,4,6,4,1.
then we sum v0, v1, v2: 1+4+6=11.
then S = a^n * T = (1/4)^2 * 11 = (1/16)*11 = 11/16.
then answer = 1-11/16=5/16.
so we don't need the coefficients to be integers in the ordinary sense, because we are in modular arithmetic.
and 11/16 mod mod = 11 * inv(16) = 11 * 562500004 = 6187500044 % mod = 687500033?
actually, 16 * 562500004 = 9000000064 % mod = 9000000064 % 1000000007 = 9000000064 - 9*1000000007 = 9000000064-9000000063=1, so inv(16)=562500004.
then 11 * 562500004 = 6187500044 % 1000000007 = 6187500044 % 1000000007 = 6187500044 - 6*1000000007 = 6187500044-6000000042=187500002.
then 1 - 187500002 = 812500005? but expected 812500006.
wait, 1 - 187500002 = 1000000007 - 187500002 = 812500005, but the sample answer is 812500006.
let me compute 5/16: 5 * 562500004 = 2812500020 % 1000000007:
1000000007 * 2 = 2000000014, 2812500020 - 2000000014 = 812500006.
so why is our S = 11/16 = 11 * 562500004 = 6187500044.
then 6187500044 / 1000000007 = quotient 6, remainder 6187500044 - 6*1000000007 = 6187500044-6000000042=187500002.
then 1 - 187500002 = 812500005, but expected 812500006.
what is 1 in mod? it is 1, so 1 - 187500002 = 1 - 187500002 = negative, so we add mod: 1 - 187500002 + 1000000007 = 1000000008 - 187500002 = 812500006.
so we must do: (1 - S + mod) % mod.
in the code, we do: ans = (1 - S + mod) % mod.
then with S=187500002, we have 1 - 187500002 = -187500001, then + mod: 1000000007 - 187500001 = 812500006.
so it is correct.
but why is v2=500000010? that is not 6.
in the recurrence
for the scaled polynomial (x^2+2x+1)^2, the coefficient
for x^2 is 6, but we computed v2=500000010.
but 500000010 is 6 in the ordinary sense? no, 500000010 is not 6.
what is the recurrence
for m=2:
v2 = [ (2*2 - 2 + 3) * v0 + 2 * (2-2+1) * v1 ] / (1 * 2)
= [ (4-2+3)=5 * 1 + 2 * 1 * 4 ] / 2 = (5+8)/2 = 13/2.
but in the polynomial (x^2+2x+1)^2, the coefficient
for x^2 is 6, not 13/2.
so there is a mistake in the recurrence.
let me re-derive the recurrence
for the scaled polynomial (x^2 + p x + q)^n.
we had:
q * (m+1) * v_{m+1} = (2n - m+1) * v_{m-1} + p * (n - m) * v_{m}
for m=1:
left: q * 2 * v_2
right: (2n-1+1) * v_{0} + p * (n-1) * v_{1} = (2n) * v0 + p*(n-1)*v1
so, v_2 = [ 2n * v0 + p*(n-1)*v1 ] / (2*q)
for our example: n=2, p=2, q=1, v0=1, v1=4.
v2 = [ 4*1 + 2*(1)*4 ] / (2*1) = (4+8)/2 = 12/2=6.
so it should be 6.
in the recurrence
for m=2 in the code, we used:
v2 = [ (2n - m+3) * v0 + p_val*(n-m+1)*v1 ] * inv(q_val * m)
= [ (4-2+3)*1 + 2*(2-2+1)*4 ] * inv(1*2)
= [5 + 2*1*4] * inv(2) = (5+8)*inv(2)=13/2.
so the recurrence
formula in the code is not the same as the derived recurrence.
why in the code we used (2n-m+3) and (n-m+1)
for the indices?
in the derived recurrence,
for v_{m} (which is the coefficient of x^m), the recurrence is:
for the current v_{m} (we are about to compute), we have an equation in terms of v_{m-2} and v_{m-1}?
but the recurrence we derived is
for v_{m+1} in terms of v_{m} and v_{m-1} (
for a given m).
specifically, to compute v_{m}
for a given m, we use the recurrence with m replaced by m-1:
to compute v_{m} (which is at index m), set in the equation:
let k = m-1, then the equation is:
q * (k+1) * v_{k+1} = (2n - k+1) * v_{k-1} + p * (n - k) * v_{k}
which is:
q * m * v_{m} = (2n - (m-1)+1) * v_{m-2} + p * (n - (m-1)) * v_{m-1}
= (2n - m+2) * v_{m-2} + p * (n - m+1) * v_{m-1}
so the recurrence
for v_{m} is:
v_{m} = [ (2n - m+2) * v_{m-2} + p * (n - m+1) * v_{m-1} ] * inv(q * m)
in our example
for m=2:
v2 = [ (4-2+2)*v0 + 2*(2-2+1)*v1 ] * inv(1*2) = [4*1 + 2*1*4] * inv(2) = (4+8)/2 = 12/2=6.
so the code should be:
term1 = (2n - m + 2) * v0
term2 = p_val * (n - m + 1) * v1
v2 = (term1 + term2) * inv(q_val * m) % mod
but in the code we had (2n - m+3) and (n-m+1) [which is the same as above] but the first term is (2n-m+3) instead of (2n-m+2).
so the mistake is in the first term: we used (2n-m+3) but should be (2n-m+2).
there
fore, we should change the recurrence to:
v2 = [ (2n - m + 2) * v0 + p_val * (n - m + 1) * v1 ] * inv(q_val * m) % mod
then
for the sample:
v2 = [ (4-2+2)*1 + 2*(2-2+1)*4 ] * inv(2) = (4+8)/2=6.
then T = v0+v1+v2 = 1+4+6=11.
then S = a^n * T = (1/4)^2 * 11 = 11/16.
then answer = 1 - 11/16 = 5/16 = 812500006.
so we must change the recurrence to (2n-m+2).
and the derivation is: from the equation
for v_{m} (which is the coefficient of x^m) we have:
q_val * m * v_{m} = (2n - (m-1) + 1) * v_{m-2} + p_val * (n - (m-1)) * v_{m-1}
= (2n - m+2) * v_{m-2} + p_val * (n - m+1) * v_{m-1}
so the code is fixed.
summary
for the else branch:
p_val = (c * inv(a)) % mod;
q_val = (b * inv(a)) % mod;
v0 = pow(q_val, n, mod); // coefficient
for x^0
T = v0;
if (n>=1) {
v1 = n * p_val % mod * pow(q_val, n-1, mod) % mod;
T = (T + v1) % mod;
}
if (n>=2) {
// Precompute the array of inverses
for 1..n (
for m in [2, n])
vector<long long> inv_arr(n+1);
inv_arr[1] = 1;
for (int i=2; i<=n; i++) {
inv_arr[i] = (mod - mod / i) * inv_arr[mod % i] % mod;
if (inv_arr[i] < 0) inv_arr[i] += mod;
}
long long inv_q_val = inv(q_val, mod); // modular inverse of q_val
long long v2;
for (int m=2; m<=n; m++) {
long long term1 = (2*n - m + 2) % mod;
term1 = term1 * v0 % mod;
long long term2 = p_val * ((n - m + 1) % mod) % mod;
term2 = term2 * v1 % mod;
v2 = (term1 + term2) % mod;
v2 = v2 * inv_q_val % mod;
v2 = v2 * inv_arr[m] % mod;
T = (T + v2) % mod;
// update
for next iteration
v0 = v1;
v1 = v2;
}
}
long long A_n = pow(a, n, mod);
long long S_val = A_n * T % mod;
long long ans = (1 - S_val + mod) % mod;
cout << ans << endl;
we must be cautious about the range of integers: n up to 10^7, and mod=10^9+7, so we use long long.
and also, the terms (2n-m+2) and (n-m+1) might be negative? we do % mod, but it is better to keep them non-negative.
we can do: (2n-m+2) % mod, but if 2n-m+2 is negative, we add mod. but 2n-m+2 is at least 2n-n+2 = n+2>=0.
similarly, (n-m+1) is non-negative
for m<=n.
so we are safe.
finally, we must consider the case n=0, but the data range says n>=1.
and the case n=1.
for n=1, the polynomial is a*x^2 + c*x + b.
then the sum of coefficients
for x^0 and x^1:
v0 = q_val^1 = b/a?
then T = v0 + v1 = b/a + 1 * p_val * (b/a)^0 = b/a + p_val.
then S_val = a^1 * T = a * (b/a + p_val) = b + a * p_val = b + a*(c/a) = b+c.
then answer = 1 - (b+c) = 1 - (b+c) = a?
but in fact,
for n=1, the difference can be:
only Aya: D=1 -> not in [0,1] (which is the sum
for m=0,1) -> we include it?
we want P(D<=0) = P(D=0 or negative) = 1 - P(D=1) = 1 - a.
and in our result: 1 - (b+c) = a, which is not 1-a.
what is the probability of D<=0
for n=1?
D = A-C, where A and C are either 0 or 1.
the possibilities:
both not hit: D=0 -> included.
only Aya: D=1 -> not included.
only Cirno: D=-1 -> included.
both hit: D=0 -> included.
so P(D<=0) = 1 - a.
but our result is a.
so there is a mistake.
in the polynomial
for the adjusted generating function, we have the polynomial (a*x^2 + c*x + b) which is
for the adjusted generating function of the difference per second?
recall: the original generating function per second is:
b * x^{-1} + c + a * x.
then
for n seconds, the generating function is (b*x^{-1}+c+a*x)^n.
then we adjusted by multiplying by x^n: (b + c*x + a*x^2)^n.
then, the coefficient of x^{n+d} in the adjusted polynomial corresponds to difference d.
then,
for n=1, the adjusted polynomial is (b + c*x + a*x^2).
then, the coefficient
for x^0: b (d = 0 - n = -1) -> difference -1.
the coefficient
for x^1: c (d = 1 - n = 0).
the coefficient
for x^2: a (d = 2 - n = 1).
then, the difference d<=0 corresponds to the coefficients
for x^{n+d} with d<=0, i.e.,
for the adjusted polynomial: the coefficients
for x^0 and x^1.
so the sum
for m=0 and m=1 is b+c.
then P(D<=0) = b+c.
and then the answer = 1 - (b+c) = a.
but in reality, the probability
for D<=0 is 1 - a = b+c.
and we have 1 - (b+c) = a, which is not the same as the desired answer.
we want to output P(D>0) = a.
and indeed, in n=1, the only way to have D>0 is only Aya hit, which is a.
so it is correct.
and the probability
for D<=0 is b+c, which is 1-a.
then our answer = 1 - (b+c) = a.
so it is correct.
there
fore, the code is fixed.
let's run the sample n=2, and we get 812500006.
we will write the code in C++ with very short variable names (<=2 characters).
but the user said "变量名长度不超过2", so we use at most 2 characters
for variable names.
we will use:
T:
for the number of test cases? but the input only has one test case, so we don't need an outer loop.
n: the first integer.
xp, yp, xq, yq:
for the input.
p = xp * inv(yp) % mod, q = xq * inv(yq) % mod.
a = p*(1-q) % mod, b = (1-p)*q % mod, c = (1-a-b) % mod.
if a==0:
// then S = (b+c)^n, ans = 1 - S.
else:
// p_val = c * inv(a) % mod, q_val = b * inv(a) % mod.
// then compute the sum T
for the scaled polynomial.
we must use long long.
and we will use:
v0, v1, v2:
for the coefficients.
T:
for the sum.
inv_arr:
for the array of inverses, but we can name it as iv or something.
but the variable names must be at most 2 characters.
we can use:
int n;
long long xp, yp, xq, yq;
long long a, b, c, p, q;
long long v0, v1, v2, T;
vector<long long> iv; //
for inv_arr
for 1..n
and
for the function inv, we can write a function.
let's code accordingly.
note: the exponentiation must be done with modular exponentiation.
we will write a fast exponentiation function.
but the standard pow in c++ is not
for modular, so we write one.
however, we can use recursion or iterative.
also, the mod is 1000000007.
we will write:
long long pw(long long a, long long b) {
long long r = 1;
while (b) {
if (b&1) r = r*a % mod;
a = a*a % mod;
b >>=1;
}
return r;
}
then, the modular inverse of a is pw(a, mod-2).
but the recurrence
for the inverses
for integers from 1 to n we do with the
linear recurrence.
code structure:
const int mod = 1000000007;
int n;
cin >> n;
long long xp, yp, xq, yq;
cin >> xp >> yp >> xq >> yq;
auto inv = [mod](long long a) { return pw(a, mod-2); };
long long p = xp * inv(yp) % mod;
long long q = xq * inv(yq) % mod;
long long a = p * (1 - q + mod) % mod;
long long b = (1 - p + mod) * q % mod;
long long c = (1 - a - b) % mod;
if (c<0) c += mod;
if (a==0) {
long long S = pw((b+c)%mod, n);
long long ans = (1 - S + mod) % mod;
cout << ans << endl;
} else {
long long ia = inv(a);
long long p_val = c * ia % mod;
long long q_val = b * ia % mod;
long long v0 = pw(q_val, n);
long long T = v0;
long long v1 = 0;
if (n>=1) {
v1 = n % mod * p_val % mod;
v1 = v1 * pw(q_val, n-1) % mod;
T = (T + v1) % mod;
}
if (n>=2) {
// precompute iv
for 1..n
vector<long long> iv(n+1);
iv[1] = 1;
for (int i=2; i<=n; i++) {
iv[i] = (mod - mod/i) * iv[mod%i] % mod;
if (iv[i]<0) iv[i] += mod;
}
long long iq = inv(q_val);
for (int m=2; m<=n; m++) {
long long t1 = (2*n - m + 2) % mod;
t1 = t1 * v0 % mod;
long long t2 = p_val * ((n - m + 1) % mod) % mod;
t2 = t2 * v1 % mod;
v2 = (t1 + t2) % mod;
v2 = v2 * iq % mod;
v2 = v2 * iv[m] % mod;
T = (T + v2) % mod;
v0 = v1;
v1 = v2;
}
}
long long a_n = pw(a, n);
long long S_val = a_n * T % mod;
long long ans = (1 - S_val + mod) % mod;
cout << ans << endl;
}
we must be cautious about the memory: the vector iv has size n+1, about 80 MB
for n=10^7.
and the time is O(n), which is acceptable.
let me test with n=1.
input: 1
1 2 1 2
then p=1/2, q=1/2, a=1/4, b=1/4, c=1/2.
then in the else branch:
p_val = (1/2) * inv(1/4) = (1/2)*4 = 2.
q_val = (1/4) * inv(1/4) = 1.
v0 = 1^1 = 1.
T=1.
then n>=1: v1 = 1 * 2 * 1^0 = 2.
T=1+2=3.
then n<2, so skip the loop.
a_n = (1/4)^1 = 250000002.
S_val = 250000002 * 3 = 750000006 % mod.
ans = 1 - 750000006 = 1000000007-750000006=250000001, which is 1/4.
and the only way to have D>0 is only Aya hit, which is a=1/4, so it is correct.
now, the sample input2, input3, input4 are large, so we trust the recurrence.
but note: the recurrence
for the coefficients is
for the scaled polynomial (x^2 + p_val x + q_val)^n, and we are summing the coefficients
for x^0 to x^n.
we output the code in C++ with short variable names.
note: the variable names must be at most 2 characters.
we have used: n, xp, yp, xq, yq, a, b, c, p, q, v0, v1, v2, T, iv, ia, iq, p_val, q_val, a_n, S_val, ans, t1, t2, m.
all are at most 2 characters.
we will write accordingly.
we must include headers.
the entire code:
note: the exponentiation function: we write a function.
and the mod is 1000000007.
also, use vector.
and note:
for n=0, the problem says n>=1.
code:
#include <iostream>
#include <vector>
#include <cmath>
using namespace std;
const int M = 1000000007;
long long pw(long long a, long long b) {
a %= M;
long long r = 1;
while (b) {
if (b&1) r = r*a % M;
a = a*a % M;
b /= 2;
}
return r;
}
int main() {
int n;
cin >> n;
long long xp, yp, xq, yq;
cin >> xp >> yp >> xq >> yq;
// compute p = xp/yp, q=xq/yq in mod M
long long p = xp * pw(yp, M-2) % M;
long long q = xq * pw(yq, M-2) % M;
// a = p*(1-q), b = (1-p)*q, c=1-a-b
long long a = p * (1 - q + M) % M;
long long b = (1 - p + M) * q % M;
long long c = (1 - a - b) % M;
if (c<0) c += M;
if (a == 0) {
long long S = pw((b+c) % M, n);
long long ans = (1 - S + M) % M;
cout << ans << endl;
} else {
long long ia = pw(a, M-2); // inv(a)
long long p_val = c * ia % M; // c/a
long long q_val = b * ia % M; // b/a
// v0 = (q_val)^n
long long v0 = pw(q_val, n);
long long T = v0;
long long v1 = 0;
if (n >= 1) {
v1 = (n % M) * p_val % M;
v1 = v1 * pw(q_val, n-1) % M;
T = (T + v1) % M;
}
if (n >= 2) {
// Precompute inverses
for 1..n
vector<long long> iv(n+1);
iv[1] = 1;
for (int i=2; i<=n; i++) {
iv[i] = (M - M/i) * iv[M%i] % M;
if (iv[i] < 0) iv[i] += M;
}
long long iq = pw(q_val, M-2); // inv(q_val)
long long v2;
for (int m=2; m<=n; m++) {
// term1 = (2n - m + 2) * v0
long long t1 = (2*(long long)n - m + 2) % M;
t1 = t1 * v0 % M;
// term2 = p_val * (n - m + 1) * v1
long long t2 = p_val * ((n - m + 1) % M) % M;
t2 = t2 * v1 % M;
v2 = (t1 + t2) % M;
v2 = v2 * iq % M;
v2 = v2 * iv[m] % M; // division by m
T = (T + v2) % M;
// shift
v0 = v1;
v1 = v2;
}
}
long long a_n = pw(a, n); // a^n
long long S_val = a_n * T % M; // the sum
for the original polynomial
long long ans = (1 - S_val + M) % M;
cout << ans << endl;
}
return 0;
}
let me test with n=2, and the sample input.
input: 2
1 2
1 2
then p = 1/2, q=1/2, a=1/4, b=1/4, c=1/2.
then in the else branch:
p_val = (1/2) * 4 = 2, q_val= (1/4)*4=1.
v0 = 1^2 = 1.
T=1.
v1 = 2 * 2 * 1^(1) = 4.
T=5.
then
for m=2:
t1 = (2*2-2+2)=4 * v0=4*1=4.
t2 = 2*(2-2+1)=2*1=2 * v1=2*4=8.
v2 = (4+8)=12 * inv(1) (which is 1) * inv(2) (which is 500000004) = 12 * 500000004 = 6000000048 % M = 6000000048 % 1000000007 = 6000000048 - 6*1000000007 = 6000000048 - 6000000042 = 6.
T = 5+6 = 11.
a_n = (1/4)^2 = 250000002^2 = 250000002*250000002 % M = ... we know 250000002^2 = (1/4)^2 = 1/16, and 1/16 mod M is 562500004.
but wait, we have a=1/4, so a_n = (1/4)^2 = 1/16.
then S_val = (1/16)*11 = 11/16, which is 11 * 562500004 % M = 6187500044 % M = 187500002.
then ans = 1 - 187500002 = 812500005? -> but then we add M: 1000000007+1-187500002 = 812500006.
so it is correct.
we output the code.
note: the data range says n up to 10^7, and we are using a vector of size n+1 (iv) and a loop of n, so it is O(n) in time and space.
and 10^7 integers is about 80 MB, which is acceptable.
we hope it passes.</think>### 问题解答
Cirno 和 Aya 的弹幕游戏中，需要计算 Cirno 被击中的次数（记为 $C$）小于 Aya 被击中的次数（记为 $A$）的概率，即 $P(C < A)$。游戏持续 $n$ 秒，每秒事件独立：
- Aya 被击中的概率 $p = \frac{x_p}{y_p}$
- Cirno 被击中的概率 $q = \frac{x_q}{y_q}$
定义每秒钟的差值变化：
- $A - C$ 增加 1 的概率 $a = p(1-q)$
- $A - C$ 减少 1 的概率 $b = (1-p)q$
- $A - C$ 不变的概率 $c = 1 - a - b$
最终需求 $P(A - C > 0)$，等价于 $1 - P(A - C \leq 0)$。通过生成函数和线性递推计算多项式 $(ax^2 + cx + b)^n$ 的前 $n+1$ 项系数和 $S$，则 $P(A - C > 0) = 1 - S$。
#### C++ 实现（变量名长度 ≤ 2）
```cpp
#include <iostream>
#include <vector>
using namespace std;
const int M = 1000000007;
long long pw(long long a, long long b) {
long long r = 1;
for (a %= M; b; b /= 2, a = a * a % M)
if (b & 1) r = r * a % M;
return r;
}
int main() {
ios::sync_with_stdio(0);
cin.tie(0);
int n;
cin >> n;
long long x, y, u, v;
cin >> x >> y >> u >> v;
long long p = x * pw(y, M-2) % M;
long long q = u * pw(v, M-2) % M;
long long a = p * (M + 1 - q) % M;
long long b = (M + 1 - p) * q % M;
long long c = (M*2 + 1 - a - b) % M;
if (a == 0) {
long long S = pw((b + c) % M, n);
cout << (M + 1 - S) % M << endl;
} else {
long long ia = pw(a, M-2);
long long pv = c * ia % M;
long long qv = b * ia % M;
long long v0 = pw(qv, n), T = v0;
long long v1 = 0;
if (n >= 1) {
v1 = (n % M) * pv % M * pw(qv, n-1) % M;
T = (T + v1) % M;
}
if (n >= 2) {
vector<long long> iv(n + 1);
iv[1] = 1;
for (int i = 2; i <= n; ++i)
iv[i] = (M - M / i) * iv[M % i] % M;
long long iq = pw(qv, M - 2);
for (int m = 2; m <= n; ++m) {
long long t1 = (2LL * n - m + 2) % M * v0 % M;
long long t2 = pv * ((n - m + 1) % M) % M * v1 % M;
long long v2 = (t1 + t2) * iq % M * iv[m] % M;
T = (T + v2) % M;
v0 = v1;
v1 = v2;
}
}
long long an = pw(a, n);
long long S = an * T % M;
cout << (M + 1 - S) % M << endl;
}
return 0;
}
```
### 算法分析
1. **概率计算**：
- 计算 $p = x_p/y_p \mod 10^9+7$ 和 $q = x_q/y_q \mod 10^9+7$
- 计算 $a = p(1-q)$, $b = (1-p)q$, $c = 1 - a - b$（模意义下）
2. **特殊情况处理**：
- 若 $a = 0$，则 $P(C < A) = 0$（无正差值可能）
3. **多项式系数和计算**：
- 构造多项式 $ax^2 + cx + b$
- 递推计算 $(ax^2 + cx + b)^n$ 的前 $n+1$ 项系数和 $S$
- 递推公式：$v_m = \frac{(2n - m + 2)v_{m-2} + p_v(n - m + 1)v_{m-1}}{q_v \cdot m}$（$v_m$ 为系数）
4. **最终概率**：
- $P(C < A) = 1 - S$
### 复杂度
- **时间复杂度**：$O(n)$（递推计算系数）
- **空间复杂度**：$O(n)$（存储递推数组）
本文介绍了一种优化的素数筛选算法，通过避免重复删除合数来提高效率。利用每个合数都能分解为两个素数乘积的特性，仅需删除特定形式的合数即可找出所有素数。
扫一扫
172
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
