[2022 Guangdong saim] Lagrange interpolation (multivariate function extreme value divide and conquer NTT)

The question

Seeking satisfaction $\sum _{i=1}^k\frac{x_i^2}{a_i^2}=1$ Under the condition of , From the length to $m$ Array of $b$ Middle selection $k$ It's made up of numbers $a_1,a_2,\cdots ,a_k$,$\prod _{i=1}^k{x}_i$ The expectation of the maximum value of ,$k$ For the even .

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

analysis ：

First, the Lagrange multiplier method of conditional extremum of multivariate functions in advanced mathematics is used to solve the maximum value , set up
$$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)$$
Find the partial derivative of each variable , Let the partial derivative be $0$ have to
$$\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}$$
Then simplify it a little , about $1\le i\le k$ There are
$$\prod _{i=1}^k{x}_i=\frac{-2\lambda x_i^2}{a_i^2}$$
Through any two formulas $1\le i,j\le k$ Simultaneous elimination $\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}$$
Simplify to
$$\frac{x_i}{a_i}=\frac{x_j}{a_j}$$
So if and only if $\frac{x_1}{a_1}=\frac{x_2}{a_2}=\cdots =\frac{x_k}{a_k}$ Maximum value when , And $\sum _{i=1}^k\frac{x_i^2}{a_i^2}=1$, So for any $1\le i\le k$ There are $\frac{x_i}{a_i}=\pm \sqrt{\frac{1}{k}}$, that $\prod _{i=1}^k{x}_i=k^{-\frac{k}{2}}\prod _{i=1}^k{a}_i$, because $k$ For the even , So it must be positive , And $\frac{k}{2}$ It must be an integer .

Seek from $b$ Select... From the array $k$ Sum of all products of numbers , Consider constructing a generating function
$$F(x)=\prod _{i=1}^k(1+b_{i}x)$$
that $[x^k]F(x)$ Is to choose $k$ Sum of all products of numbers , All in all $\left(\genfrac{}{}{0ex}{}{m}{k}\right)$ Seed selection , So the expectation is
$$k^{-\frac{k}{2}}\times \frac{[x^k]F(x)}{\left(\genfrac{}{}{0ex}{}{m}{k}\right)}$$
$F(x)$ Divide and conquer available $\text{NTT}$ Calculation , Total time complexity $O(n{\log }^{2}n)$

Code ：

#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;
}