一、题意理解

我们有：

• 一个集合 $S$，包含 $n$ 个正整数 $a_1, a_2, \dots, a_n$
• 一个线性递推数列 $f$： $f(i) = 2 f(i-1) + 3 f(i-2) \quad (i > 1)$ 已知 $f(0)$ 和 $f(1)$
• 我们要计算： $ \sum_{\substack{S' \subseteq S \\ |S'|=k}} f\left( \sum_{x \in S'} x \right) $ 即所有大小为 $k$ 的子集的元素和作为 $f$ 的输入，然后求和。

---

二、样例分析

样例1：

n=4, k=2
a = [1,2,1,3]
f0=1, f1=1

递推： $f_0=1, f_1=1, f_2=5, f_3=13, f_4=41, f_5=121, \dots$

子集和：

• {1,2} → 和=3 → $f(3)=13$
• {1,1} → 和=2 → $f(2)=5$
• {1,3} → 和=4 → $f(4)=41$
• {2,1} → 和=3 → $f(3)=13$
• {2,3} → 和=5 → $f(5)=121$
• {1,3} → 和=4 → $f(4)=41$

总和 = $13+5+41+13+121+41 = 234$

---

三、问题难点

• $n$ 最多 $10^5$，不能枚举所有 $\binom{n}{k}$ 个子集
• $a[i]$ 最大 $10^9$，但 $f$ 只对 $99991$ 取模，且递推是线性的
• 需要高效计算所有 $k$-子集的 $f(\text{sum})$

---

四、关键思路：生成函数与线性递推

设 $A(x) = \prod_{i=1}^n (1 + x^{a_i})$，则 $[x^s]A(x)$ 表示和为 $s$ 的子集个数（不考虑大小限制）。

我们要求的是： $ \sum_{s \ge 0} (\text{大小为 }k\text{ 且和为 }s\text{ 的子集数}) \cdot f(s) $ 即： $ \sum_{s \ge 0} [x^s y^k] \prod_{i=1}^n (1 + y x^{a_i}) \cdot f(s) $

---

五、利用线性递推的性质

$f$ 是线性递推，可以用矩阵快速幂计算 $f(s)$，但 $s$ 可能很大（最大 $k \cdot \max a_i \approx 10^{14}$），不能直接枚举 $s$。

但注意模数 $M = 99991$ 很小，$f(s)$ 在模 $M$ 下是周期性的（因为递推是线性的，模 $M$ 后状态数有限）。

---

1. 矩阵形式

递推： $ \begin{bmatrix} f(i) \\ f(i-1) \end{bmatrix} = \begin{bmatrix} 2 & 3 \\ 1 & 0 \end{bmatrix} \cdot \begin{bmatrix} f(i-1) \\ f(i-2) \end{bmatrix} $ 所以： $ \begin{bmatrix} f(s) \\ f(s-1) \end{bmatrix} = \begin{bmatrix} 2 & 3 \\ 1 & 0 \end{bmatrix}^s \cdot \begin{bmatrix} f(1) \\ f(0) \end{bmatrix} $ 记 $M = \begin{bmatrix} 2 & 3 \\ 1 & 0 \end{bmatrix}$，在模 $99991$ 下，$M$ 的阶是有限的（至多 $99991^2$ 量级，实际更小）。

---

2. 循环节

由于 $M$ 在 $GL_2(\mathbb{F}_{99991})$ 中，其乘法阶 $T$ 整除 $|\mathrm{GL}_2(\mathbb{F}_{99991})| = (99991^2-1)(99991^2-99991)$，所以 $T$ 不会太大（相比 $10^{14}$ 小很多）。

我们可以找到 $T$ 使得 $M^T \equiv I \pmod{99991}$，那么 $f(s+T) \equiv f(s) \pmod{99991}$。

---

六、问题转化

我们要求： $ \sum_{\substack{S' \subseteq S \\ |S'|=k}} f\left( \sum_{x \in S'} x \right) \pmod{99991} $ 由于 $f$ 在模 $99991$ 下有周期 $T$，我们可以将 $\sum x$ 模 $T$ 考虑。

定义： $ cnt[r] = \text{大小为 }k\text{ 的子集且和模 }T\text{ 为 }r\text{ 的数量} $ 那么： $ \text{答案} = \sum_{r=0}^{T-1} cnt[r] \cdot f(r) \pmod{99991} $ 注意：这里 $f(r)$ 是 $f$ 在 $r$ 处的值（模 $99991$），不是 $f(r \bmod T)$，因为 $f$ 本身在模 $99991$ 下周期为 $T$，所以 $f(r) \equiv f(r \bmod T)$。

---

七、计算 $cnt[r]$

我们要计算： $ cnt[r] = [x^r y^k] \prod_{i=1}^n (1 + y x^{a_i \bmod T}) $ 因为 $a_i$ 模 $T$ 后范围在 $[0,T-1]$，且 $T$ 不太大（可接受范围内），我们可以用生成函数+快速幂的方法。

---

设多项式 $P(x) = \sum_{j=0}^{T-1} c_j x^j$，其中 $c_j$ 是 $a_i \bmod T = j$ 的个数。

那么： $ \prod_{i=1}^n (1 + y x^{a_i \bmod T}) = \prod_{j=0}^{T-1} (1 + y x^j)^{c_j} $ 我们要求 $[y^k]$ 的这个乘积，这可以用对数-指数方法或分治NTT计算（但 $T$ 可能很大？需要估计 $T$ 的大小）。

---

八、$T$ 的估计

$M = \begin{bmatrix} 2 & 3 \\ 1 & 0 \end{bmatrix}$ 在 $\mathbb{F}_{99991}$ 上的特征多项式为： $\lambda^2 - 2\lambda - 3 = 0$ 判别式 $\Delta = 4 + 12 = 16$，有整数根 $3$ 和 $-1$，所以 $M$ 可对角化，且 $M^T = I$ 当且仅当 $3^T \equiv 1$ 且 $(-1)^T \equiv 1$ 模 $99991$。

$(-1)^T \equiv 1$ 意味着 $T$ 是偶数。

$3^T \equiv 1 \pmod{99991}$，$99991$ 是质数吗？检查：$99991 = 99991$，是质数。那么 $3$ 在 $\mathbb{F}_{99991}^\times$ 中的阶整除 $99990$。

所以 $T$ 是 $99990$ 的某个偶数因子。$99990 = 2 \times 3^3 \times 5 \times 7 \times 23 \times 37$，所以 $T$ 不会太大，最大 $99990$。

---

九、算法步骤

1. 找到 $T$ 使得 $M^T \equiv I \pmod{99991}$（可以取 $T=99990$ 或更小的偶数因子测试）。
2. 预处理 $f(0..T-1)$ 模 $99991$。
3. 计算 $a_i \bmod T$，得到计数 $c_0..c_{T-1}$。
4. 计算生成函数 $F(y) = \prod_{j=0}^{T-1} (1 + y x^j)^{c_j}$ 的 $[y^k]$ 系数（这是一个长度为 $T$ 的多项式），即 $cnt[r]$。
5. 答案 $= \sum_{r=0}^{T-1} cnt[r] \cdot f(r) \pmod{99991}$。

---

十、生成函数计算优化

$F(y)$ 的 $[y^k]$ 系数计算： $[y^k] \prod_{j=0}^{T-1} (1 + y x^j)^{c_j}$ 这相当于：有 $T$ 种物品，第 $j$ 种有 $c_j$ 个，每个重量为 $j$，选 $k$ 个物品的所有方案的生成函数。

可以用动态规划： 设 $dp[t][i]$ 表示前 $t$ 种物品选了 $i$ 个的所有方案的生成函数（多项式模 $x^T$）。

转移： $dp[t] = dp[t-1] * (1 + x^j)^{c_j}$ 其中 $*$ 是多项式卷积（模 $x^T$），可以用快速幂计算 $(1+x^j)^{c_j}$。

复杂度 $O(T^2 \log n)$，$T \approx 10^5$ 时太大。

实际上可以用分治NTT，但这里 $T$ 最大 $99990$，$n=10^5$，可能勉强可过。

---

十一、代码框架（C++）

#include <bits/stdc++.h>
#include <ext/pb_ds/tree_policy.hpp>
#include <ext/pb_ds/trie_policy.hpp>
#include <ext/pb_ds/assoc_container.hpp>
#include <ext/pb_ds/priority_queue.hpp>

using namespace std;
using namespace __gnu_pbds;
using namespace __gnu_cxx;
typedef long long ll;
//tree<int, null_type, less<int>, rb_tree_tag, tree_order_statistics_node_update> tr;
//__gnu_pbds::priority_queue<int, greater<int>, pairing_heap_tag> qu;
//typedef trie<string, null_type, trie_string_access_traits<>, pat_trie_tag, trie_prefix_search_node_update> pref_trie;
const int N = 100005;
int n, m, k;
int a[N];
namespace FFTComplex {
const double PI = acos(-1.0);

struct C {
    double x, y;

    C(double x = 0, double y = 0) : x(x), y(y) {}

    C operator!() const {
        return C(x, -y);    //共轭复数
    }
};

inline C operator*(const C &a, const C &b) {
    return C(a.x * b.x - a.y * b.y, a.x * b.y + a.y * b.x);
}

inline C operator+(const C &a, const C &b) {
    return C(a.x + b.x, a.y + b.y);
}

inline C operator-(const C &a, const C &b) {
    return C(a.x - b.x, a.y - b.y);
}

constexpr int L(1 << 20);
C w[L];
int _ = [] {
    for (int i = 0; i < L / 2; i++)
        w[i + L / 2] = C(cos(2 * PI *i / L), sin(2 * PI *i / L));

    for (int i = L / 2 - 1; i; i--)
        w[i] = w[i << 1];

    return 0;
}();

void dft(C *a, int n) {
    for (int k = n >> 1; k; k >>= 1)
        for (int i = 0; i < n; i += k << 1)
            for (int j = 0; j < k; j++) {
                C x = a[i + j], y = a[i + j + k];
                a[i + j + k] = (x - y) * w[j + k];
                a[i + j] = x + y;
            }
}

void idft(C *a, int n) {
    for (int k = 1; k < n; k <<= 1)
        for (int i = 0; i < n; i += k << 1)
            for (int j = 0; j < k; j++) {
                C &x = a[i + j], y = a[i + j + k] * w[j + k];
                a[i + j + k] = x - y, x = x + y;
            }

    for (int i = 0; i < n; i++)
        a[i].x /= n, a[i].y /= n;

    reverse(a + 1, a + n);
}
}
namespace polybase {//需要保证输入项均为正!!!
constexpr ll mod(99991);

ll fpow(ll x, ll r) {
    ll result = 1;

    while (r) {
        if (r & 1)
            result = result * x % mod;

        r >>= 1;
        x = x * x % mod;
    }

    return result;
}

using namespace FFTComplex;

int norm(int n) {
    return 1 << __lg(n * 2 - 1);
}

struct poly : public vector<ll> {
    using vector<ll>::vector;
#define T (*this)

    poly modxk(int k) const {
        k = min(k, (int) size());
        return poly(begin(), begin() + k);
    }

    poly rev() const {
        return poly(rbegin(), rend());
    }

    friend poly operator*(const poly &x, const poly &y) {
        if (x.empty() || y.empty())
            return poly();

        poly a(x), b(y);
        int len = a.size() + b.size() - 1, n = max(2, norm(len));//n=1时该板子卷积部分会有问题
        vector<C> f(n), g(n), p(n), q(n);

        for (int i = 0; i < a.size(); i++)
            f[i] = C(a[i] >> 15, a[i] & 0x7fff);

        for (int i = 0; i < b.size(); i++)
            g[i] = C(b[i] >> 15, b[i] & 0x7fff);

        dft(f.data(), n), dft(g.data(), n);

        for (int k = 1, i = 0, j; k < n; k <<= 1)
            for (; i < k * 2; i++) {
                j = i ^ (k - 1);
                p[i] = C(0.5, 0) * g[i] * (f[i] + !f[j]);
                q[i] = C(0, -0.5) * g[i] * (f[i] - !f[j]);
            }

        idft(p.data(), n), idft(q.data(), n);
        a.resize(len);

        for (int i = 0; i < len; i++) {
            ll x, y, z, w;
            x = p[i].x + 0.5;
            y = p[i].y + 0.5;
            z = q[i].x + 0.5;
            w = q[i].y + 0.5;
            a[i] = ((x % mod << 30) + ((y + z) % mod << 15) + w) % mod;
        }

        return a;
    }

    poly operator+(const poly &b) {
        poly a(T);

        if (a.size() < b.size())
            a.resize(b.size());

        for (int i = 0; i < b.size(); i++) { //用b.size()防止越界
            a[i] += b[i];

            if (a[i] >= mod)
                a[i] -= mod;
        }

        return a;
    }

    poly operator-(const poly &b) {
        poly a(T);

        if (a.size() < b.size())
            a.resize(b.size());

        for (int i = 0; i < b.size(); i++) {
            a[i] -= b[i];

            if (a[i] < 0)
                a[i] += mod;
        }

        return a;
    }

    poly operator*(const ll p) {
        poly a(T);

        for (auto &x : a)
            x = x * p % mod;

        return a;
    }

    poly &operator<<=(int r) {
        return insert(begin(), r, 0), T;    //注意逗号,F(x)*(x^r)
    }

    poly operator<<(int r) const {
        return poly(T) <<= r;
    }

    poly operator>>(int r) const {
        return r >= size() ? poly() : poly(begin() + r, end());
    }

    poly &operator>>=(int r) {
        return T = T >> r;    //F[x]/(x^r)
    }

#undef T
};
}
using namespace polybase;

template<const int x>
poly solve(int l, int r) {
    if (l == r)
        return poly{1, fpow(x, a[l])};

    int m = (l + r) >> 1;

    return solve<x>(l, m) * solve<x>(m + 1, r);
}

int main() {
    int p, q, u, v, w, x, y, z, T;
    cin >> n >> k;

    for (int i = 1; i <= n; i++)
        scanf("%d", &a[i]);

    ll f0, f1;
    cin >> f0 >> f1;
    ll A = (f0 + f1) * fpow(4, mod - 2) % mod;
    ll B = (f0 - A + mod) % mod;
    poly F = solve<3>(1, n);
    poly G = solve < mod - 1 > (1, n);
    cout << (A * F[k] + B * G[k]) % mod;
    return 0;
}

---

十二、总结

本题的关键在于：

1. 发现 $f$ 是线性递推，在模 $99991$ 下有周期 $T$。
2. 将问题转化为计算模 $T$ 下和为 $r$ 的大小为 $k$ 的子集数。
3. 使用生成函数和多项式卷积计算 $cnt[r]$。
4. 利用矩阵快速幂或直接递推预处理 $f(r)$。