题解 P9091 【「SvR-2」Let's Meet at a Higher Place】

Leasier · 2023-02-18 20:40:07 · 题解

令人遗憾的是本题唯一赛时通过者 @Silver187 写的是经过大力卡常后的算法三。

算法一

容易发现前缀 \gcd 的限制相当于 \forall 2 \leq i \leq m, b_i \mid b_{i - 1}。

枚举 b_1 后记忆化搜索即可。时间复杂度为 O(nm^2 \log n)，可以通过 \text{Subtask } 1。

算法二

考虑 f(n, m, k) 中相邻项相等出现的次数的限制实际上相当于钦定某些项，使之与后一项相等，于是我们设 g(n, m) 表示以 n 结尾的合法的 m 项的方案数，则 f(n, m, 0) = \displaystyle\sum_{i = 1}^n g(i, m)。

接下来考虑 g(n, m)，显然可以容斥。我们希望求的是恰好有 0 个 1 \leq i < m, b_i = b_{i + 1} 的方案数，通过容斥我们可以将其转化为至少有 k 个 1 \leq i < m, b_i = b_{i + 1} 的方案数。

考虑钦定其中 k 项与后一项相同，方案数是 C_{m - 1}^k（因为不能钦定最后一项）；对于余下的项，其可以选择与后一项相同，也可以不同，且每一项一定是后一项的因数，方案数是 I_{m - k}(b_m)。于是我们可以得到 g(n, m) = \displaystyle\sum_{i = 0}^{m - 1} (-1)^i C_{m - 1}^i I_{m - i}(n)。

注意到 f(n, m, k) = \displaystyle\sum_{i = 0}^k C_{m - 1}^i f(n, m - i, 0)，令 S_f(n) = \displaystyle\sum_{i = 1}^n f(i)，则 f(n, m, k) = \displaystyle\sum_{i = 0}^k C_{m - 1}^i \sum_{j = 1}^n \sum_{x = 0}^{m - i - 1} (-1)^x C_{m - i - 1}^x I_{m - i - x}(j) = \sum_{i = 0}^k C_{m - 1}^i \sum_{x = 0}^{m - i - 1} (-1)^x C_{m - i - 1}^x S_{I_{m - i - x}}(n)。

原式 = \displaystyle\sum_{i = 1}^n \sum_{j = 1}^m \sum_{k = 0}^{j - 1} \sum_{l = 0}^k C_{j - 1}^l \sum_{x = 0}^{j - l - 1} (-1)^x C_{j - l - 1}^x S_{I_{j - l - x}}(\lfloor \frac{n}{i} \rfloor)

注意到我们可以直接枚举 T = l + x，有：

原式 = \displaystyle\sum_{i = 1}^n \sum_{j = 1}^m \sum_{T = 0}^{j - 1} S_{I_{j - T}}(\lfloor \frac{n}{i} \rfloor) \sum_{l = 0}^T (-1)^{T - l} C_{j - 1}^l C_{j - l - 1}^{T - l} \sum_{k = l}^{j - 1} 1

注意到最后一个求和符号的结果只与 j, T 有关且 m 很小，于是令 g_{n, m} = \displaystyle\sum_{i = 0}^m (-1)^{m - i} (n - i) C_{n - 1}^i C_{n - i - 1}^{m - i}，则可以在 O(m^3) 的时间复杂度内预处理出。

原式 = \displaystyle\sum_{i = 1}^n \sum_{j = 1}^m \sum_{T = 0}^{j - 1} g_{j, T} S_{I_{j - T}}(\lfloor \frac{n}{i} \rfloor)

注意到我们可以直接枚举 U = j - T，有：

注意到最后一个求和符号的值只与 T 有关，于是令 h_n = \displaystyle\sum_{i = 0}^{m - n} g_{n + i, i}，则可以在 O(m^2) 的时间复杂度内预处理出。

原式 = \displaystyle\sum_{i = 1}^n \sum_{U = 1}^m h_U S_{I_U}(\lfloor \frac{n}{i} \rfloor)

线性筛所有 I_i 即可。时间复杂度为 O(nm + m^3)，可以通过 \text{Subtask } 1, 2。

算法三

事实上我们只需要对于 \forall 1 \leq i \leq m，求出 I_i 的块筛。

数论分块 + 根号分治或递推版 Min_25 筛求解即可。时间复杂度为 O(mn^{\frac{2}{3}} + m^3) 或 O(\frac{mn^{\frac{3}{4}}}{\ln n} + m^3)，可以通过 \text{Subtask } 1 \sim 3。实现常数足够优秀还可以通过 \text{Subtask } 4。

算法四

考虑 I_i 在质数处的表现，我们发现 \forall p \in \operatorname{Prime}，I_i(p) 的值与 p 具体是多少无关。

于是不难想到 Powerful Number 筛。构造 I_i'(p^k) = [k \leq 1] I_i(p^k)，现在问题转变为求出 I_i' 的块筛。

显然可以算贡献，可得：S_{I_i'}(n) = \displaystyle\sum_{i = 0} I_i'(p)^i \sum_{j = 1}^n [\omega(j) = i] \mu^2(j)，其中 p 为任意质数。后半部分显然可以用类似 Min25 筛 Part 1 的方法求出，因为 $\displaystyle\max{i = 1}^n w(i) = O(\frac{\ln n}{\ln \ln n})，则时间复杂度为 O(\frac{n^{\frac{3}{4}}}{\ln \ln n})，进而可以在 O(\frac{m \sqrt{n} \ln n}{\ln \ln n})$ 的时间复杂度内求出块筛。

但如果此时直接像杜教筛一样根号分治，块筛部分的时间复杂度为 O(\frac{n^{\frac{3}{4}}}{\ln \ln n} + mn^{\frac{2}{3}})，不能通过。

考虑 Powerful Number 的性质。对于 \lfloor \frac{n}{i} \rfloor，当 i \leq \sqrt[3]{n}，显然只有 O(\sqrt[3]{n}) 种值；当 i > \sqrt[3]{n}，\lfloor \frac{n}{i} \rfloor < n^{\frac{2}{3}}，所以只有 \sqrt{n^{\frac{2}{3}}} = O(\sqrt[3]{n}) 种值。据此我们可以通过根号分治将单次求块筛的时间复杂度优化至 O(n^{\frac{3}{5}})。

时间复杂度为 O(\frac{n^{\frac{3}{4}}}{\ln \ln n} + mn^{\frac{3}{5}} + m^3)，可以通过所有数据，若实现常数过大将无法通过 \text{Subtask } 5。

算法五

感谢 @渐变色 提供本部分做法。

注意到 h_i = [i = m] m（具体证明略去），则原式 = m \displaystyle\sum_{i = 1}^n S_{I_m}(\lfloor \frac{n}{i} \rfloor) = m S_{I_{m + 1}}(n)。

于是我们只需要求 S_{I_{m + 1}}(n)。可以直接采用上述方法，并还可以在 O(m \sqrt{n}) 的时间复杂度对于 1 \sim m 求解答案。

时间复杂度为 O(m \log n + \frac{n^{\frac{3}{4}}}{\ln \ln n})。

其实可以出加强版，但是懒得出了（

代码：

#include <stdio.h>
#include <math.h>

typedef unsigned int uint;
typedef long long ll;
typedef unsigned long long ull;
typedef __int128 lll;
typedef __uint128_t ulll;

typedef struct Division_tag {
    ulll a;
    Division_tag(){}
    Division_tag(ull mod){
        a = (((ulll)1 << 64) + mod - 1) / mod;
    }
} Division;

const int N = 1e5 + 7, M = 2e5 + 7, K = 67 + 7, P = 33 + 7, Q = 10 + 7;
int sqrt_n, prime_cnt = 0, id = 0;
ll cur_n;
int prime[N], max_prime_cnt[M];
uint c[K][P], g[M], h[Q][M], prime_power[P];
ll number[M];
bool p[N];
Division inv_prime[N];

ull operator /(const ull &a, const Division &b){
    return a * b.a >> 64;
}

inline int min(int a, int b){
    return a < b ? a : b;
}

inline int log2(ll n){
    int ans = log2((double)n);
    while ((1ll << ans) <= n) ans++;
    return ans - 1;
}

inline int get_max_prime_cnt(ll n, int cnt){
    int ans = 0;
    for (register ll i = 1; i <= n && ans <= cnt; i *= prime[ans]) ans++;
    return ans - 1;
}

inline int divide(ll x, ll y){
    return 1.0 * x / y;
}

inline int get_id(ll n){
    return n <= sqrt_n ? id - n + 1 : divide(cur_n, n);
}

inline void init(ll n, int m){
    int up1 = log2(n), up2 = up1 + m - 1;
    sqrt_n = sqrt(n);
    cur_n = n;
    c[0][0] = 1;
    for (register int i = 1; i <= up2; i++){
        int t = min(i, up1);
        c[i][0] = 1;
        for (register int j = 1; j <= t; j++){
            c[i][j] = c[i - 1][j] + c[i - 1][j - 1];
        }
    }
    p[0] = p[1] = true;
    for (register int i = 2; i <= sqrt_n; i++){
        if (!p[i]){
            prime_cnt++;
            prime[prime_cnt] = i;
            inv_prime[prime_cnt] = Division(i);
        }
        for (register int j = 1; j <= prime_cnt && i * prime[j] <= sqrt_n; j++){
            p[i * prime[j]] = true;
            if (i % prime[j] == 0) break;
        }
    }
    for (register ll i = 1, j; i <= n; i = j + 1){
        ll tn = n / i;
        j = n / tn;
        id++;
        number[id] = tn;
        max_prime_cnt[id] = get_max_prime_cnt(tn, prime_cnt);
        g[id] = tn - 1;
    }
    for (register int i = 1; i <= prime_cnt; i++){
        ll t = 1ll * prime[i] * prime[i];
        for (register int j = 1; j <= id && number[j] >= t; j++){
            g[j] -= g[get_id(number[j] / inv_prime[i])] - (i - 1);
        }
    }
    for (register int i = 1; i <= id; i++){
        h[0][i] = 1;
        h[1][i] = g[i];
    }
    for (register int i = prime_cnt; i >= 1; i--){
        ll t1 = 1ll * prime[i] * prime[i];
        for (register int j = 1; j <= id && number[j] >= t1; j++){
            int cur_id = get_id(number[j] / inv_prime[i]), t2 = max_prime_cnt[j];
            h[2][j] -= i;
            for (register int k = 2; k <= t2; k++){
                h[k][j] += h[k - 1][cur_id];
            }
        }
    }
}

inline uint get_f__sum(ll n, int m){
    int id = get_id(n);
    uint ans = 0;
    for (register int i = max_prime_cnt[id]; i >= 0; i--){
        ans = ans * m + h[i][id];
    }
    return ans;
}

uint get_f_sum(int cur, int m, ll val1, uint val2){
    uint ans = val2 * get_f__sum(cur_n / val1, m);
    for (register int i = cur + 1; i <= prime_cnt && (lll)val1 * prime[i] * prime[i] <= cur_n; i++){
        int j = 2;
        for (register ll x = val1 * prime[i] * prime[i]; x <= cur_n; j++, x *= prime[i]){
            ans += get_f_sum(i, m, x, val2 * prime_power[j]);
        }
    }
    return ans;
}

int main(){
    int m, mi;
    ll n;
    scanf("%lld %d", &n, &m);
    mi = m + 1;
    init(n, mi);
    prime_power[0] = 1;
    for (register int i = 1; (1ll << i) <= n; i++){
        prime_power[i] = c[i + m][i] - mi * prime_power[i - 1];
    }
    printf("%u", m * get_f_sum(0, mi, 1, 1));
    return 0;
}