题目要求

$\frac{1}{\pi} \int_{0}^{\infty} \frac{1}{ \prod_{i = 1}^{n} (a_i^2 + x^2)} {\rm d}x$

这里我们需要进行列项，列项的方法是部分分式展开法。

部分分式展开法

这部分内容可以参考《信号与系统》这门课当中的复频域分析部分。（这题据说也可参考《复变函数》当中的留数定理）

假设分式 $F(s) = \frac{N(s)}{D(s)}$ ， $N(s)$ 和 $D(s)$ 分别是关于 $s$ 的多项式。

分母 $D(s)$ 可以进行因式分解：

$D(s) = c(s - s_1)(s - s2) \cdots (s - s_n)$

这里 $s_i$ 是 $D(s)$ 的了零点，也是 $F(s)$ 的极点，当它们互不相等的时候， $F(s)$ 可以表示为

$\begin{aligned} \frac{N(s)}{D(s)} = & \frac{N(s)}{c(s - s_1)(s - s2) \cdots (s - s_n)} \\ = & \frac{k_1}{s - s_1} + \frac{k_2}{s - s_2} + \cdots + \frac{k_n}{s - s_n} \end{aligned}$

式子里， $k_i$ 为待定系数，在等式两边同时乘以因子 $(s-s_i)$ ，再令 $s=s_i$ ，于是右边只留下了 $k_i$ 项：

$k_i = \left. \left[ (s-s_i) \frac{N(s)}{D(s)} \right] \right|_{s=s_i} = \left. \left[ \frac{N(s)}{\prod_{j \neq i} (s - s_j) } \right] \right|_{s = s_i}$

本题题解

观察我们要列项的式子，可以认为 $s = x^2$ ， $s_i = -a_i^2$ ，接着利用部分分式展开得到

$k_i = \frac{1}{\prod_{j \neq i} [(-a_i^2) - (-a_j^2)] } =\frac{1}{\prod_{j \neq i} (a_j^2-a_i^2) }$

接下来对于展开的每一项分别积分后求和，每一项分别是

$\int_{0}^{\infty} \frac{k_i}{a_i^2 +x^2} {\rm d}x = \frac{k_i}{a_i} \int_{0}^{\infty} \frac{1}{1 + \left( \frac{x}{a} \right) ^ 2 } {\rm d}{\frac{x}{a_i}} =\frac{k_i \pi}{2a_i}$

注意对 $k_i$ 求逆元的时候先把分母全部相乘再求逆元，这样只需要求一次。

整个算法的渐进时间复杂度为 $O(n^2 + n\log P)$ 。

#include <bits/stdc++.h>

using namespace std;

typedef long long LL;

const int MAX_N = 1000 + 5;
const LL P = LL(1E9) + 7;

int n;
LL a[MAX_N];

inline LL getPow(LL x, LL y, LL mod) {
    LL ret = 1;
    while (y) {
        if (y & 1) ret = ret * x % mod;
        x = x * x % mod;
        y >>= 1;
    }
    return ret;
}

inline LL sub(LL u, LL v) {
    return ((u - v) % P + P) % P;
}

inline LL getC(int x) {
    LL ret = 1;
    for (int i = 1; i <= n; ++i) {
        if (i == x) continue;
        ret = ret * sub(a[i] * a[i], a[x] * a[x]) % P;
    }
    ret = getPow(ret, P - 2, P);
    return ret;
}

int main() {
    while (scanf("%d", &n) != EOF) {
        for (int i = 1; i <= n; ++i) scanf("%lld", a + i);
        LL sum = 0;
        for (int i = 1; i <= n; ++i) {
            LL ans = 1;
            ans = ans * getC(i) % P;
            ans = ans * getPow(a[i] << 1, P - 2, P) % P;
            sum = (sum + ans) % P;
        }
        printf("%lld\n", sum);
    }
    return 0;
}

我的博客

www.cfzhao.com
(欢迎dalao们py友链接啊！)

所以对于每个 i 要找到满足 p(j) <= sum - p(i) 的 j 有几个
于是可以把所有的 p(i) 和 sum - p(i) 离散化，树状数组，

2019牛客多校第一场B - Integration（题解）

部分分式展开法

本题题解

我的博客