BZOJ 4652 - [Noi2016]循环之美

#数学 #莫比乌斯反演 #数论 #NOI #杜教筛

Published on 2016-08-17

描述

分析

记 $x\perp y$ 为 $x$ 与 $y$ 互质，如果 $\frac x y, x\perp y$ 是纯循环小数，那么必定存在一个正整数 $l$ 使得

\begin{aligned} x\cdot k^l &\equiv x \pmod y\\ k^l &\equiv 1 \pmod y \end{aligned}

根据欧拉定理，如果 $k\perp y$ ，那么有 $k^{\varphi(y)}\equiv 1 \pmod y$ ，由于逆元是唯一的，所以存在 $l$ 的充分必要条件是 $k\perp y$ 。

则问题转化为求

\sum_{i=1}^n\sum_{j = 1}^m[j\perp k][i\perp j]

其中 $n, m \le {10}^9, k\le 2000$ 。
我们利用莫比乌斯函数化简一下：

\begin{aligned} &\sum_{i=1}^n\sum_{j = 1}^m[j\perp k][i\perp j]\\ =&\sum_{i=1}^n\sum_{j = 1, j\perp k}^m[i\perp j]\\ =&\sum_{i=1}^n\sum_{j = 1, j\perp k}^m\sum_{x\mid\gcd(i, j)}\mu(x)\\ =&\sum_{i=1}^n\sum_{j = 1, j\perp k}^m\sum_{x\mid i, x\mid j}\mu(x)\\ =&\sum_{x = 1, x\perp k}\mu(x)\sum_{u=1}^{\left\lfloor\frac n x\right\rfloor}\sum_{v = 1}^{\left\lfloor\frac m x\right\rfloor}[v\perp k]\\ =&\sum_{x = 1, x\perp k}\mu(x)\left\lfloor\frac n x\right\rfloor\sum_{v = 1}^{\left\lfloor\frac m x\right\rfloor}[v\perp k] \end{aligned}

这是一个分块的形式，我们需要快速求两个东西的前缀和

\begin{aligned} f(n) &= \sum_{i = 1}^{n}[i\perp k]\\ g(n) &= \sum_{i = 1}^n\mu(i)[i\perp k] \end{aligned}

对于 $f(n)$ ，因为有 $\gcd(a, b) = \gcd(b, a\bmod b)$ ，所以 $f(n) = \left\lfloor\frac n k\right\rfloor f(k) + f(n\bmod k)$ ，只需要预处理 $\le k$ 的 $f(n)$ 便可以 $O(1)$ 计算。

重点是求 $g(n)$ ，而且 $n$ 高达 ${10}^9$ ，怎么办？

一种方法

这个方法是我一年前学了解到的，当时还不会证复杂度，学了洲阁筛以后顿时觉得这个方法十分 naive。

设 $s(n, k) = \sum_{i = 1, i\perp k}^n \mu(i)$ ，则有：

\begin{aligned} s(n, k)&= \sum_{i = 1, i\perp k}^n \mu(i)\\ &= \sum_{i = 1}^n \mu(i)\sum_{x\mid k, x\mid i}\mu(x)\\&= \sum_{x \mid k}\mu(x)\sum_{u = 1}^{\left\lfloor\frac n x\right\rfloor}{\mu(u\cdot x)}\\ &= \sum_{x \mid k}\mu^2(x)\sum_{u = 1, u\perp x}^{\left\lfloor\frac n x\right\rfloor}{\mu(u)}\\ &= \sum_{x \mid k}\mu^2(x)s\left(\left\lfloor\frac n x\right\rfloor, x\right) \end{aligned}

其中倒数第二个等号是非常妙的，根据 $\mu(n)$ 的定义，如果 $u\perp x$ ，那么 $\mu(u\cdot x) =\mu(u)\mu(x)$ ，否则就是 $0$ 。这样我们就可以递归下去了，边界为 $k = 1$ ，此时需要求 $\mu$ 的前缀和，直接使用杜教筛即可。

对于复杂度分析，显然所有 $s(n, k)$ 中用到的 $n$ 只有 $O(\sqrt n)$ 个，当 $k = 1$ 时杜教筛复杂度为 $O(n^{2/3})$ 。其余的 $k$ 至多只有 $\sigma_0(k)$ 种，转移也只有 $\sigma_0(k)$ 个，复杂度可以粗略估计为 $O(n^{2/3} + \sigma_0^2(k)\sqrt n)$ ，事实上由于我们每次递归到 $s(\left\lfloor\frac n x\right\rfloor, x)$ ，两维总有一个会比较小，所以这个上界是远远达不到的。

杜教筛学习：浅谈一类积性函数的前缀和、BZOJ 3944、BZOJ 4176。

另一种方法

观察要求的式子，由于我们要求 $i\perp k$ ，这等价于 $i$ 中不能包含任何 $k$ 中拥有的质因子，我们可以使用类似洲阁筛的思想，依次「筛去」 $k$ 中包含的质因子。

将 $k$ 的质因子写为 $p_1, p_2, \ldots, p_m$ ，设 $s_i(n)$ 为 $\mu(j)$ 的和，其中 $j\in [1, n]$ 且 $j$ 中不包含前 $i$ 个质因子。 $s_0(n)$ 就是 $\mu$ 的前缀和，同样使用杜教筛求出。

转移就是一个不断筛数的过程

s_i(n) = s_{i - 1}(n) - \mu(p_i)s_{i - 1}\left(\left\lfloor\frac n {p_i}\right\rfloor\right)

所有需要用到的 $s_i(n)$ 中 $n$ 的值只会有 $O(\sqrt n)$ 个，所以复杂度为 $O(n^{2/3} + \omega(k)\sqrt n)$ ，其中 $\omega(k)$ 为 $k$ 的质因子个数。

由于 $\omega(k)$ 是非常小的，这种方法可以将 $k$ 出到 $10^9$ 甚至 $10^{18}$ ，不过此时 $f(n)$ 就不能预处理前 $k$ 项了，可以使用刚刚的方法做到一样的复杂度。

所以本题的最终复杂度为 $O(n^{2/3} + \omega(k)\sqrt n)$ 。

代码

//  Created by Sengxian on 8/16/16.
//  Copyright (c) 2016年 Sengxian. All rights reserved.
//  BZOJ 4652 莫比乌斯反演 + 杜教筛 NOI2016 D1T3
# pragma GCC optimize("O3")
#include <bits/stdc++.h>
using namespace std;

typedef long long ll;
const int maxn = 1000000 + 3, maxk = 2000 + 3;
int primes[maxn], mu[maxn], s[maxn], num;
bool isPrime[maxn];
int K;

void sieve() {
    fill(isPrime, isPrime + maxn, 1);
    mu[1] = 1, num = 0;
    for (int i = 2; i < maxn; ++i) {
        if (isPrime[i]) primes[num++] = i, mu[i] = -1;
        static int d;
        for (int j = 0; j < num && (d = i * primes[j]) < maxn; ++j) {
            isPrime[d] = false;
            if (i % primes[j] == 0) {
                mu[d] = 0;
                break;
            } else mu[d] = -mu[i];
        }
    }
    for (int i = 1; i < maxn; ++i) s[i] = s[i - 1] + mu[i];
}

int pre[maxk];
inline ll G(int n) {
    return (n / K) * pre[K - 1] + pre[n % K];
}

namespace hashTable {
    const int hashBase = 999997;
    typedef pair<int, int> state;
    typedef unsigned long long ull;
    struct node {
        node *next;
        ull val;
        int ans;
        node(node *next, ull val, ll ans): next(next), val(val), ans(ans) {}
        node() {}
    }pool[maxn], *first[hashBase];
    inline ull getHash(const state &val) {
        return (ull)val.first * maxk + val.second;
    }
    inline node* find(const ull &val) {
        int pos = val % hashBase;
        for (node *cur = first[pos]; cur; cur = cur->next)
            if (cur->val == val) return cur;
        return NULL;
    }
    inline void insert(const ull &val, const int &ans) {
        static node *pit;
        if (!pit) pit = pool;
        int pos = val % hashBase;
        first[pos] = new(pit++) node(first[pos], val, ans);
    }
}

using namespace hashTable;
inline int solve(int n, int k) {
    if (n == 0 || (k == 1 && n < maxn)) return s[n];
    ll val = getHash(state(n, k));
    node *cur = find(val);
    if (cur) return cur->ans;
    int ans = 0;
    if (k == 1) {
        ans = 1;
        for (int i = 2, j = 2; i <= n; i = j + 1) {
            j = n / (n / i);
            ans -= (j - i + 1) * solve(n / i, k);
        }
    } else {
        ans = 0;
        for (int i = 1; i * i <= k; ++i) if (k % i == 0) {
            if (mu[i]) ans += solve(n / i, i);
            if (i * i != k && mu[k / i])
                ans += solve(n / (k / i), k / i);
        }
    }
    insert(val, ans);
    return ans;
}

inline ll F(int n, int m) {
    ll ans = 0, last = 0, now = 0;
    for (int i = 1, j = 1; i <= min(n, m); i = j + 1) {
        j = min(n / (n / i), m / (m / i));
        now = solve(j, K);
        ans += (now - last) * (n / i) * G(m / i);
        last = now;
    }
    return ans;
}

int main() {
    int n, m;
    scanf("%d%d%d", &n, &m, &K);
    sieve();
    for (int i = 1; i <= K; ++i)
        pre[i] = __gcd(i, K) == 1, pre[i] += pre[i - 1];
    printf("%lld\n", F(n, m));
    return 0;
}

BZOJ 4174 - tty的求助 BZOJ 3930 - [CQOI2015]选数

Sengxian's Blog