题解 accoders::NOI 5515【Precise】

题解 accoders::NOI 5515【Precise】

这个题的名字就叫 Precise /prɪˈsaɪs/，我说的。

我的天，这是 EI 题。

problem

定义 \(a_k=(\sum_{i = 0}^{m - 1}k ^ i)\bmod M\)。对于所有 \(1\leq k\leq n\) 求出 \(a_k\)，输出异或和。（异或和没有性质，所以还是需要求全部答案）

然而数据范围是：\(n\leq 10^7, m\leq 10^{10^6}, M\leq 10^9\)，不保证 \(M\) 是质数。

solution 0

如果 \(M\) 是个质数，我们直接求 \(\dfrac{k^m-1}{k-1}\) 即可。

如果允许 \(O(n\log m)\)，我们可以对原式分治，分成两半，并发现两半的答案实际上只用算一边，做到 \(O(\log m)\) 求一个值。

若 \(m\) 是高精度数字，根据扩展欧拉定理：

\[a^b\equiv \begin{cases} a^b, &b<\varphi(p),\\ a^{b\bmod \varphi(p)+\varphi(p)}, &b \geq \varphi(p). \end{cases}\pmod p \]

对于质数，我们可以将 \(k^m\) 降幂；对于 \(O(\log m)\)，我们还是可以扩展欧拉定理降幂，但是需要计算跳过的环的贡献（实际上我们的扩展欧拉定理是跳环，然而环的答案在这里并不是 \(0\)）。

获得 76 分的高分。

solution 1

考虑到如果将 \(M\) 分解为 \(\prod p_i^{\alpha_i}\)，且对于一个 \(k\) 我们能对所有 \(p^\alpha\) 求出答案，则我们可以 CRT 合并。若 \(k-1\) 在 \(\pmod {p^\alpha}\) 意义下有逆元（即 \(\gcd(k-1, p)=1\)，也就是 \(p\nmid (k-1)\)，那么直接跑路；求逆元可以线性筛或者用什么 前后缀积 + 全部积的逆元 做到线性，平凡。<- 笑死，一点也不平凡，最好使用 exgcd 求逆元

否则设 \(k=rp+1\)，得到

\[\dfrac{k^m-1}{k-1}=\dfrac{(rp+1)^m-1}{rp}=\dfrac{\sum_{i=0}^m\binom m i (rp)^i-1}{rp}=\sum_{i=0}^{\min(m, \alpha-1)}\binom m {i+1} (rp)^i \]

先不管组合数并记 \(f_i=\binom m {i+1}\)，一共有 \(n/p\) 个 \(k\) 要计算，每次都是 \(O(\alpha)\)，时间复杂度 \(O(\sum_pn\alpha/p)\)，给个 \(2^{30}\) 就原地爆炸了。

考虑 BSGS，记 \(rp=u+vp^\beta\)。

\[f(rp)=\sum_{i=0}^{\alpha-1}f_i(rp)^i=\sum_{i=0}^{\alpha-1}f_i(u+vp^\beta)^i \]

对后面做二项式展开后交换求和顺序，原来的 \(i\) 变为 \(i+j\)。

\[f(rp)=\sum_{i+j<\alpha}f_{i+j}u^iv^jp^{\beta j}\binom {i+j}{i} \]

此时可以使 \(\beta j<\alpha\)。若预处理 \(h(u,j)=\sum_{i+j<\alpha}f_{i+j}u^i\binom {i+j}{i}\)（对一个 \(p\) 的预处理复杂度为 \(O(p^{\beta-1}\alpha^2/\beta)\)），则 \(f(rp)=\sum_{j<\alpha/\beta}h(u,j)v^jp^{\beta j}\)。

取 \(\beta=\left\lfloor\alpha/2\right\rfloor\) 得到复杂度为预处理 \(O(p^{\alpha/2})\)，询问 \(O(1)\)。胜利。

solution 2

还没有结束，\(f_i=\dbinom m {i+1}\) 仍是未解之谜。我们需要对每个质因子 \(p\) 算 \(\alpha\) 个。模 \(p^\alpha\)。

我们的做法是，观察到 \(m\) 与 \(m\bmod \varphi(M)M\) 的答案是相同的，因为这样是一次跳过 \(M\) 个长为 \(\varphi(M)\) 的环，一次跳过的总和是 \(0\)，可以忽略。

然后我们有 \((a\bmod bc)\bmod c=a\bmod c\)。

所以可以算了。

复杂度 \(O(n\omega(M))\)。

code

// accoders
#include <cstdio>
#include <vector>
#include <cstring>
#include <cassert>
#include <utility>
#include <tuple>
#include <algorithm>
#include <functional>
using namespace std;
#ifdef LOCAL
#define debug(...) fprintf(stderr, ##__VA_ARGS__)
#else
#define debug(...) void(0)
#endif
typedef long long LL;
template <int id>
struct dynamic_modint {
    static unsigned P;
    unsigned v;
    dynamic_modint() : v(0) {}
    template <class T>
    dynamic_modint(T x) {
        x %= (int)P, v = x < 0 ? x + P : x;
    }
    dynamic_modint operator+() const { return *this; }
    dynamic_modint operator-() const { return dynamic_modint(0) - *this; }
    dynamic_modint inv() const { return assert(v), qpow(*this, P - 2); }
    friend int raw(const dynamic_modint &self) { return self.v; }
    template <class T>
    friend dynamic_modint qpow(dynamic_modint a, T b) {
        dynamic_modint r = 1;
        for (; b; b >>= 1, a *= a)
            if (b & 1)
                r *= a;
        return r;
    }
    dynamic_modint &operator+=(const dynamic_modint &rhs) {
        if (v += rhs.v, v >= P)
            v -= P;
        return *this;
    }
    dynamic_modint &operator-=(const dynamic_modint &rhs) {
        if (v -= rhs.v, v >= P)
            v += P;
        return *this;
    }
    dynamic_modint &operator*=(const dynamic_modint &rhs) {
        v = 1ull * v * rhs.v % P;
        return *this;
    }
    dynamic_modint &operator/=(const dynamic_modint &rhs) { return *this *= rhs.inv(); }
    friend dynamic_modint operator+(dynamic_modint lhs, const dynamic_modint &rhs) { return lhs += rhs; }
    friend dynamic_modint operator-(dynamic_modint lhs, const dynamic_modint &rhs) { return lhs -= rhs; }
    friend dynamic_modint operator*(dynamic_modint lhs, const dynamic_modint &rhs) { return lhs *= rhs; }
    friend dynamic_modint operator/(dynamic_modint lhs, const dynamic_modint &rhs) { return lhs /= rhs; }
    friend bool operator==(const dynamic_modint &lhs, const dynamic_modint &rhs) { return lhs.v == rhs.v; }
    friend bool operator!=(const dynamic_modint &lhs, const dynamic_modint &rhs) { return lhs.v != rhs.v; }
};
template <int id>
unsigned dynamic_modint<id>::P;
typedef dynamic_modint<-1> mint;
template <int N>
struct siever {
    int pri[N + 10], cnt, div[N + 10];
    siever() : cnt(0) {
        memset(pri, 1, sizeof pri);
        pri[0] = pri[1] = 0;
        div[1] = 1;
        for (int i = 1; i <= N; i++) {
            if (pri[i])
                pri[++cnt] = i, div[i] = i;
            for (int j = 1; j <= cnt && i * pri[j] <= N; j++) {
                pri[i * pri[j]] = 0;
                div[i * pri[j]] = pri[j];
                if (i % pri[j] == 0)
                    break;
            }
        }
    }
};
template <class T>
tuple<T, T, T> exgcd(T a, T b, T c) {
    if (!b)
        return { c / a, 0, a };
    auto [x, y, d] = exgcd(b, a % b, c);
    return { y, x - a / b * y, d };
}
template <class T>
T gcd(T x, T y) {
    return !y ? x : gcd(y, x % y);
}
template <class T>
T solveEquation(T a, T b, T c) {  // get the vaild solution x of ax + by = c
    auto [x, y, d] = exgcd(a, b, c);
    auto mod = [](T x, T y) { return (x % y + y) % y; };
    return c % d == 0 ? mod(x, b / d) : -1;
}
template <class T>
T getInv(T a, T p) {
    return solveEquation(a, p, 1);
}
LL modolus(char buf[], LL mod) {
    int len = strlen(buf);
    unsigned long long x = 0;
    bool flag = 0;
    for (int i = 0; i < len; i++) {
        x = x * 10 + buf[i] - '0';
        if (x >= mod)
            flag = 1, x %= mod;
    }
    if (flag)
        x += mod;
    return x;
}
pair<int, vector<pair<int, int>>> divide(int n) {
    int ans = n;
    vector<pair<int, int>> res = {};
    for (int i = 2; i * i <= n; i++) {
        if (n % i == 0) {
            int c = 0;
            while (n % i == 0) ++c, n /= i;
            res.emplace_back(i, c);
            ans = ans / i * (i - 1);
        }
    }
    if (n > 1) {
        res.emplace_back(n, 1);
        ans = ans / n * (n - 1);
    }
    return { ans, res };
}
mint binom(LL n, int m) {
    mint ans = 1;
    LL s1 = 1, s2 = 1;
    for (int i = 1; i <= min(20, m); i++) s1 *= i;
    for (int i = 21; i <= m; i++) s2 *= i;
    for (int i = 0; i < m; i++) {
        LL now = n - i, d;
        d = gcd(now, s1);
        now /= d, s1 /= d;
        d = gcd(now, s2);
        now /= d, s2 /= d;
        ans *= mint(now);
    }
    assert(s1 == 1 && s2 == 1);
    return ans;
}
template <int N>
struct solverInv {
    typedef dynamic_modint<0> modi;
    modi pre[N + 10], ipre[N + 10], a[N + 10];
    int cnt;
    void set_mod(int p) { cnt = 0, modi::P = p, pre[0] = ipre[0] = 1; }
    void add(modi x) { a[++cnt] = x, pre[cnt] = pre[cnt - 1] * x; }
    void build() {
        ipre[cnt] = getInv(raw(pre[cnt]), int(modi::P));
        for (int i = cnt; i >= 1; i--) ipre[i - 1] = ipre[i] * a[i];
    }
    mint query(int i) { return raw(ipre[i] * pre[i - 1]); }
};
int n;
LL m;
mint answer[10000010], tmp[10000010], F[50], C[50][50], qpn[10000010];
siever<10000010> S;
mint inv[10000010];
void part1(int p, int alpha, function<void(int, mint)> insert) {
    int pa = 1;
    for (int i = 0; i < alpha; i++) pa *= p;
    int now = 0;
    for (int i = 1; i <= pa && i <= n; i++) {
        if (S.div[i] == i)
            inv[i] = getInv(i, pa);
        else
            inv[i] = inv[i / S.div[i]] * inv[S.div[i]];
        if (i % p == 1 && alpha == 1)
            insert(i, tmp[i] = m);
        else if ((i - 1) % p != 0)
            insert(i, tmp[i] = (qpn[i] - 1) * inv[i - 1]);
    }
}
void part2(int p, int alpha, function<void(int, mint)> insert) {
    int beta = (alpha + 1) / 2, lim = 1;
    for (int i = 0; i < beta; i++) lim *= p;
    int pa = 1;
    for (int i = 0; i < alpha; i++) pa *= p;
    static mint H[1 << 15][2];
    for (int u = 0; u < lim; u++) {
        for (int j = 0; j < 2; j++) {
            mint &res = H[u][j] = 0, now = 1;
            for (int i = 0; i + j < alpha; i++, now *= u) res += F[i + j] * now * C[i + j][i];
        }
    }
    auto query = [&](int k) { return H[k % lim][0] + H[k % lim][1] * (k / lim) * lim; };
    for (int i = 1; i <= pa && i <= n; i += p) insert(i, tmp[i] = query(i - 1));
}
void solve(int p, int alpha) {
    debug("solve(%d, %d)\n", p, alpha);
    int _coe = 1;
    for (int i = 0; i < alpha; i++) _coe *= p;
    mint coe = mint(mint::P / _coe) * getInv(int(mint::P / _coe), _coe);
    debug("coe = %d\n", raw(coe));
    auto insert = [&](int i, mint ans) { answer[i] += mint(raw(ans) % _coe) * coe; };
    part1(p, alpha, insert);
    if (alpha > 1)
        part2(p, alpha, insert);
    int pa = 1;
    for (int i = 0; i < alpha; i++) pa *= p;
    if (pa <= n)
        tmp[0] = tmp[pa];
    for (int i = pa + 1; i <= n; i++) insert(i, tmp[i] = tmp[i - pa]);//注意这里有巨大优化
}
char buf[1 << 20];
int main() {
#ifdef ONLINE_JUDGE
    freopen("precise.in", "r", stdin);
    freopen("precise.out", "w", stdout);
#endif
    scanf("%d%s%u", &n, buf, &mint::P);
    auto [phiP, ps] = divide(mint::P);
    m = modolus(buf, 1ll * mint::P * phiP);
    debug("m = %d\n", m);
    for (int i = 0; i <= 30 && i < m; i++) F[i] = binom(m, i + 1);
    for (int i = 0; i <= 30; i++) {
        C[i][0] = 1;
        for (int j = 1; j <= i; j++) C[i][j] = C[i - 1][j] + C[i - 1][j - 1];
    }
    int rm = m < phiP ? m : m % phiP + phiP;
    for (int i = 1; i <= n; i++) {
        if (S.div[i] == i)
            qpn[i] = qpow(mint(i), rm);
        else
            qpn[i] = qpn[i / S.div[i]] * qpn[S.div[i]];
    }
    fprintf(stderr, "hjhakccpc\n");
    for (auto &&[p, alpha] : ps) solve(p, alpha);
    int output = 0;
    for (int i = 1; i <= n; i++) output ^= raw(answer[i]);
    printf("%d\n", output);
    return 0;
}