🌑

小羊儿的心情天空

2019CCPC湘潭邀请赛 F. Neko and function

Jun 1, 2019 由 小羊

题目链接:HDU-6537

题目有 f(n,k)f(n,k) 的定义,但是显然不够好直接计算。不妨先将限定条件改为 ai1a_i \ge 1,考虑 g(n,k)g(n,k)

此时就可以发现,每个素数和其幂次对 g(n,k)g(n,k) 的贡献是独立的,也就是说,g(n,k)g(n,k)是个积性函数。所以,我们来考虑 g(pe,k)g(p^e, k)。这个过程中,相当于将 ee 个小球放进 kk 个盒子,并且盒子有序可空,球无序。

g(n,k)=penCe+k1k1g(n, k) = \prod_{p^e || n} C_{e+k-1}^{k-1}

用容斥原理去除重复项。考虑长度为 ss{an}\lbrace a_n \rbrace 存在 11 的方案对 g(n,x)g(n,x) 的贡献。那么,在长度为 kk 的数组里选择 xx 个放 11,剩下的将之前方案排列进去,也就是 CkxC_k^x

f(n,k)=x=0k(1)kxg(n,x)Ckxf(n, k) = \sum_{x=0}^k (-1)^{k-x} g(n, x) C_k^x

所以

i=1nf(n,k)=x=0k(1)kxCkxi=1ng(n,x)\sum_{i=1}^n f(n, k) = \sum_{x=0}^k (-1)^{k-x} C_k^x \sum_{i=1}^n g(n, x)

那么用 min_25 筛就可以解决 g(n,k)g(n,k) 的前缀和啦,再特判一些可能浪费时间的样例。

吐槽一下,这个时限卡的太狠了,用long long过不了。本身就是个人类智慧题,竟然还要卡常。

#include <bits/stdc++.h>
using namespace std;
const int MAXN = 5e6+5;
const int MOD = 1e9+7;
bool isnp[MAXN];
int pri[MAXN], Pcn, id1[MAXN], id2[MAXN], m;
int w[MAXN], n, Sqrt, h[MAXN]; // h:sum 1
int inv[MAXN], invs[MAXN], fac[MAXN], pcn, K, kk;
inline int mul(int a, int b) { return int(1ll * a * b % MOD); }
inline int mul(int a, int b, int c) { return int(1ll * a * b % MOD * c % MOD); }
inline void modu(int &a) { while (a >= MOD) a -= MOD; while (a < 0) a += MOD; }
inline int add(int a, int b) { modu(a); modu(b); modu(a += b); return a; }
inline int add(int a, int b, int c) { return add(a, add(b, c)); }
inline int C(int N, int M) { return N < 0 || M > N || M < 0 ? 0 : mul(fac[N], invs[M], invs[N-M]); }

void sieve()
{
    inv[1] = invs[0] = fac[0] = invs[1] = fac[1] = 1;

    for (int i = 2; i < MAXN; i++)
    {
        if (!isnp[i]) pri[++Pcn] = i;
        inv[i] = mul(MOD-MOD/i, inv[MOD%i]);
        fac[i] = mul(fac[i-1], i);
        invs[i] = mul(invs[i-1], inv[i]);
        for (int j = 1; j <= Pcn && i * pri[j] < MAXN; j++)
        {
            isnp[i * pri[j]] = true;
            if (i % pri[j] == 0) break;
        }
    }
}

int S(int x, int y)
{
    if (x <= 1 || pri[y] > x) return 0;
    int k = (x <= Sqrt) ? id1[x] : id2[n/x];
    int ans = mul(add(h[k], MOD - y + 1), K);

    for (int i = y; i <= pcn && pri[i] <= x / pri[i]; i++)
    {
        int p1 = pri[i], p2 = p1 * p1;
        for (int e = 1; p1 <= x / pri[i]; ++e, p1 = p2, p2 *= pri[i])
            ans = add(ans, mul(S(x/p1, i+1), C(e+K-1, K-1)), C(e+K, K-1));
    }

    return ans;
}

inline int g(int x)
{
    if (K == 0) return 1;
    else if (K == 1) return x % MOD;
    else return S(n,1) + 1;
}

inline int calc()
{
    if (kk == 1) return (n-1) % MOD;
    else if (n < (1<<kk)) return 0;
    else if (n < (1<<kk)+(1<<(kk-1))) return 1;

    Sqrt = sqrt(n);
    for (pcn = 1, m = 0; pri[pcn] <= Sqrt; pcn++);

    for (int l = 1, r; l <= n; l = r+1)
    {
        r = n / (n / l); w[++m] = n / l;
        if (w[m] <= Sqrt) id1[w[m]] = m; else id2[n/w[m]] = m;
        h[m] = add(w[m], MOD - 1);
    }

    for (int j = 1; j <= pcn; j++)
    {
        for (int i = 1; i <= m && pri[j] <= w[i] / pri[j]; i++)
        {
            int k = (w[i]/pri[j]<=Sqrt)?id1[w[i]/pri[j]]:id2[n/(w[i]/pri[j])];
            h[i] = add(h[i], MOD - h[k], j - 1);
        }
    }

    int ans = 0;
    for (K = 0; K <= kk; K++)
        ans = add(ans, mul(g(n), C(kk, K), ((kk-K)&1) ? MOD-1 : 1));
    return (ans % MOD + MOD) % MOD;
}

int main()
{
    sieve();
    while (scanf("%d %d", &n, &kk) == 2)
        printf("%d\n", calc());
    return 0;
}

vjudge 的 Leaderboard 上还有个神仙做法,看懂后更新。