[2022 广东省赛M] 拉格朗日插值 (多元函数极值 分治NTT)

题意

求在满足 $\sum _{i=1}^k\frac{x_i^2}{a_i^2}=1$ 的条件下，从长度为 $m$ 的数组 $b$ 中选 $k$ 个数组成 $a_1,a_2,\cdots ,a_k$，$\prod _{i=1}^k{x}_i$ 的最大值的期望，$k$ 为偶数。

$(1\le k\le m\le 10^5,0<b_i<10^9)$

分析：

首先求解最大值需要用到高等数学中多元函数条件极值的拉格朗日乘数法，设
$$L(x_1,x_2,\cdots ,x_k,\lambda )=\prod _{i=1}^k{x}_i+\lambda (\sum _{i=1}^k\frac{x_i^2}{a_i^2}-1)$$
对每个变量求偏导数，令偏导数为 $0$ 得
$$\begin{gathered} \frac{\partial L}{\partial x_1}=\frac{\prod _{i=1}^k{x}_i}{x_1}+\frac{2\lambda x_1}{a_1^2}=0 \\ \frac{\partial L}{\partial x_2}=\frac{\prod _{i=1}^k{x}_i}{x_2}+\frac{2\lambda x_2}{a_2^2}=0 \\ \cdots \\ \frac{\partial L}{\partial x_k}=\frac{\prod _{i=1}^k{x}_i}{x_k}+\frac{2\lambda x_k}{a_k^2}=0 \\ \frac{\partial L}{\partial \lambda }=\sum _{i=1}^k\frac{x_i^2}{a_i^2}-1=0 \end{gathered}$$
那么稍微化简一下，对于 $1\le i\le k$ 都有
$$\prod _{i=1}^k{x}_i=\frac{-2\lambda x_i^2}{a_i^2}$$
通过任意两式 $1\le i,j\le k$ 联立消掉 $\lambda$
$$\frac{a_i^2\prod _{i=1}^k{x}_i}{-2x_i^2}=\frac{a_j^2\prod _{i=1}^k{x}_i}{-2x_j^2}$$
化简得
$$\frac{x_i}{a_i}=\frac{x_j}{a_j}$$
所以当且仅当 $\frac{x_1}{a_1}=\frac{x_2}{a_2}=\cdots =\frac{x_k}{a_k}$ 时取得最大值，且 $\sum _{i=1}^k\frac{x_i^2}{a_i^2}=1$，所以对任意 $1\le i\le k$ 都有 $\frac{x_i}{a_i}=\pm \sqrt{\frac{1}{k}}$，那么 $\prod _{i=1}^k{x}_i=k^{-\frac{k}{2}}\prod _{i=1}^k{a}_i$，因为 $k$ 为偶数，所以一定为正，且 $\frac{k}{2}$ 一定是整数。

求从 $b$ 数组中选出 $k$ 个数的所有乘积之和，考虑构造生成函数
$$F(x)=\prod _{i=1}^k(1+b_{i}x)$$
那么 $[x^k]F(x)$ 就是选出 $k$ 个数的所有乘积之和，总共有 $\left(\genfrac{}{}{0ex}{}{m}{k}\right)$ 种选法，所以期望就为
$$k^{-\frac{k}{2}}\times \frac{[x^k]F(x)}{\left(\genfrac{}{}{0ex}{}{m}{k}\right)}$$
$F(x)$ 可用分治 $\text{NTT}$ 计算，总时间复杂度 $O(n{\log }^{2}n)$

代码：

#include <bits/stdc++.h>
#define int long long
#define poly vector<int>
#define len(x) ((int)x.size())
using namespace std;
const int N = 3e5 + 5, G = 3, Ginv = 332748118, mod = 998244353;
int rev[N], lim;
int qmi(int a, int b) {
    
    int res = 1;
    while (b) {
    
        if (b & 1) res = res * a % mod;
        a = a * a % mod;
        b >>= 1;
    }
    return res;
}
void polyinit(int n) {
    
    for (lim = 1; lim < n; lim <<= 1);
    for (int i = 0; i < lim; i ++) rev[i] = (rev[i >> 1] >> 1) | (i & 1 ? lim >> 1 : 0);
}
void NTT(poly &f, int op) {
    
    for (int i = 0; i < lim; i ++) {
    
        if (i < rev[i]) swap(f[i], f[rev[i]]);
    }
    for (int mid = 1; mid < lim; mid <<= 1) {
    
        int Gn = qmi(op == 1 ? G : Ginv, (mod - 1) / (mid << 1));
        for (int i = 0; i < lim; i += mid * 2) {
    
            for (int j = 0, G0 = 1; j < mid; j ++, G0 = G0 * Gn % mod) {
    
                int x = f[i + j], y = G0 * f[i + j + mid] % mod;
                f[i + j] = (x + y) % mod, f[i + j + mid] = (x - y + mod) % mod;
            }
        }
    }
    if (op == -1) {
    
        int inv = qmi(lim, mod - 2);
        for (int i = 0; i < lim; i ++) f[i] = f[i] * inv % mod;
    }
}
poly operator * (poly f, poly g) {
    
    int n = len(f) + len(g) - 1;
    polyinit(n), f.resize(lim), g.resize(lim);
    NTT(f, 1), NTT(g, 1);
    for (int i = 0; i < lim; i ++) f[i] = f[i] * g[i] % mod;
    NTT(f, -1), f.resize(n);
    return f;
}
vector<int> fact, infact;
void init(int n) {
    
    fact.resize(n + 1), infact.resize(n + 1);
    fact[0] = infact[0] = 1;
    for (int i = 1; i <= n; i ++) {
    
        fact[i] = fact[i - 1] * i % mod;
    }
    infact[n] = qmi(fact[n], mod - 2);
    for (int i = n; i; i --) {
    
        infact[i - 1] = infact[i] * i % mod;
    }
}
int C(int n, int m) {
    
    if (n < 0 || m < 0 || n < m) return 0ll;
    return fact[n] * infact[n - m] % mod * infact[m] % mod;
}
signed main() {
    
    cin.tie(0) -> sync_with_stdio(0);
    init(1e5);
    int n, k;
    cin >> n >> k;
    vector<int> b(n + 1);
    vector<poly> f(n + 1, poly(2));
    for (int i = 1; i <= n; i ++) {
    
        cin >> b[i];
        f[i][0] = 1, f[i][1] = b[i];
    }
    function<poly(int, int)> dc = [&](int l, int r) {
    
        if (l == r) return f[l];
        int mid = l + r >> 1;
        return dc(l, mid) * dc(mid + 1, r);
    };
    poly ans = dc(1, n);
    int res = 1;
    cout << qmi(qmi(k, k / 2), mod - 2) * ans[k] % mod * qmi(C(n, k), mod - 2) % mod << endl;
}