敬启，致那时的我

题面

题目描述

实乃理给你两个整数 \(S, k\)，你需要帮她求出以下式子的值对 \(1,000,000,007\) 取模的结果：

\[\sum_{X = 0}^S [\mathrm{popc}(X) = k]F(X) \]

其中 \(F\) 为斐波那契数列，即 \(F(0) = F(1) = 1, F(n) = F(n - 1) + F(n - 2)\)。\(\mathrm{popc}(X)\) 表示 \(X\) 的二进制表示中 \(1\) 的个数。

输入格式

第一行两个整数 \(n, k\)。

第二行一个长度为 \(n\) 的 \(01\) 串 \(S_{2}\)，代表 \(S\) 的二进制表示。保证 \(S_2\) 没有前导零。

输出格式

输出一行一个整数表示答案对 \(1,000,000,007\) 取模的结果。

输入输出样例 #1

输入 #1

3 0
110

输出 #1

1

输入输出样例 #2

输入 #2

3 1
110

输出 #2

8

输入输出样例 #3

输入 #3

3 2
110

输出 #3

24

输入输出样例 #4

输入 #4

17 8
11011111101010010

输出 #4

894224396

说明/提示

样例解释

对于样例一，答案为 \(F(0) = 1\)；

对于样例二，答案为 \(F(1) + F(2) + F(4) = 1 + 2 + 5 = 8\)；

对于样例三，答案为 \(F(3) + F(5) + F(6) = 3 + 8 + 13 = 24\)。

数据范围

对于 \(100\%\) 的数据，\(0 \le k \le n \le 2 \times 10^3\)，\(n, S \ge 1\)。

子任务编号	特殊性质	分数
\(1\)	A	\(20\)
\(2\)	B	\(30\)
\(3\)	C	\(30\)
\(4\)	无	\(20\)

特殊性质 A：\(S \le 10^6\)。

特殊性质 B：\(k \le 2\)。

特殊性质 C：\(S\) 可以表示为 \(2^y - 1\) 的形式，其中 \(y\) 是一个正整数。

题解

这道题发现是FIB序列，那这个FIB序列没啥好的性质，只能用矩阵乘法快速求出。

这道题开始通过特殊性质B，发现如果枚举两个一在哪里然后再去算矩阵快速幂这样是 \(O(n^3)\) 会超时，那我们发现快速幂这个东西本质不就是在每个二进制为一的时候去乘一下当前的阶乘吗？那我们就可以把快速幂数组预处理出来。那这样我们就不用每次去重新做一遍快速幂了。

那这个部分分启示我们要预处理快速幂数组，那我们考虑预处理完这个题目变成了什么别的形式，那现在这个题目变成了在这n个预处理出的快速幂元素中，选择k个相乘的和，那这个是好做的：设 \(f_{i,j}\) 表示前i个元素选了j个的乘积和。

那本题还要考虑是否超过了上界，那这个就写一个类似数位dp的东西，记录每一位是否已经达到上界即可。

最终状态设计：\(f_{i,j,0/1}\) 表示前i个元素选了j个是否到达上界的乘积和。

代码

#include <bits/stdc++.h>
using namespace std;
using ll = long long;
const ll mod = 1e9 + 7;
struct MAT
{
    ll a[2][2];
    MAT operator*(const MAT &x) const
    {
        MAT c;
        c.a[0][0] = c.a[0][1] = c.a[1][0] = c.a[1][1] = 0;
        for (int k = 0; k < 2; k++)
        {
            for (int i = 0; i < 2; i++)
            {
                for (int j = 0; j < 2; j++)
                {
                    c.a[i][j] += a[i][k] * x.a[k][j] % mod;
                    c.a[i][j] %= mod;
                }
            }
        }
        return c;
    }
    MAT operator+(const MAT &x) const
    {
        MAT c;
        c.a[0][0] = c.a[0][1] = c.a[1][0] = c.a[1][1] = 0;
        for (int i = 0; i < 2; i++)
        {
            for (int j = 0; j < 2; j++)
            {
                c.a[i][j] = (a[i][j] + x.a[i][j]) % mod;
            }
        }
        return c;
    }
    void print()
    {
        for (int i = 0; i < 2; i++)
        {
            for (int j = 0; j < 2; j++)
            {
                cout << a[i][j] << " \n"[j == 1];
            }
        }
    }
} base, bse, ksm[2005], f[2][2005][2], zro, fib;
int main()
{
    ios::sync_with_stdio(0), cin.tie(0);
    int n, K;
    cin >> n >> K;
    string s;
    cin >> s;
    // reverse(s.begin(), s.end());
    fib.a[0][0] = fib.a[0][1] = 1;
    fib.a[1][0] = fib.a[1][1] = 0;
    bse.a[0][0] = bse.a[1][1] = 1;
    bse.a[0][1] = bse.a[1][0] = 0;
    base.a[0][0] = base.a[0][1] = base.a[1][0] = 1;
    base.a[1][1] = 0;
    ksm[n - 1] = base;
    for (int i = n - 2; i >= 0; i--)
    {
        base = base * base;
        ksm[i] = base;
    }
    int cur = 0, nxt = 1;
    f[cur][0][1] = bse;
    for (int i = 0; i < n; i++)
    {
        for (int j = 0; j <= K; j++)
        {
            for (int tt = 0; tt <= 1; tt++)
            {
                f[nxt][j][tt] = zro;
            }
        }
        int bit = s[i] - '0';
        MAT MT = ksm[i];
        for (int j = 0; j <= K; j++)
        {
            for (int tt = 0; tt <= 1; tt++)
            {
                if (f[cur][j][tt].a[0][0] == 0 && f[cur][j][tt].a[0][1] == 0 && f[cur][j][tt].a[1][0] == 0 && f[cur][j][tt].a[1][1] == 0)
                {
                    continue;
                }
                int zd = 1;
                if (tt)
                {
                    zd = bit;
                }
                for (int tb = 0; tb <= zd; tb++)
                {
                    int ntt = (tt & (tb == zd));
                    int nj = j + tb;
                    if (nj > K)
                    {
                        continue;
                    }
                    MAT nw = f[cur][j][tt];
                    MAT nMT = nw * (tb ? MT : bse);
                    f[nxt][nj][ntt] = f[nxt][nj][ntt] + nMT;
                }
            }
        }
        swap(cur, nxt);
    }
    cout << ((f[cur][K][0]).a[0][0] + (f[cur][K][1]).a[0][0]) % mod;
    return 0;
}