【题解】51nod-1965-奇怪的式子

Posted on 2018-10-09 Edited on 2018-10-10 In ACM Disqus:

题目

http://www.51nod.com/onlineJudge/problemSolution.html#!problemId=1965

题意

$\prod_{i=1}^n \sigma_0(i)^{\mu(i)+i}$

题解

设 $S_1=\prod_{i=1}^n \sigma_0(i)^i,S_2=\prod_{i=1}^n \sigma_0(i)^{\mu(i)}$。

若 x = ∏p_i^α_i，那么有 σ₀(x) = ∏(α_i + 1)。

枚举素数 p 的指数 k 的贡献，有 $$ \begin{eqnarray} S_1 &=& \prod_p \prod_k (k+1)^{\sum_{i=1}^n i[p^k|i][p^{k+1}\nmid i]} \\ &=& \prod_p \prod_k (k+1)^{\sum_{i=1}^n i[p^k|i] -\sum_{i=1}^ni[p^k|i]} \\ &=& \prod_p \prod_k (k+1)^{p^ks(\lfloor \frac n {p^k}\rfloor)-p^{k+1}s(\lfloor \frac n {p^{k+1}}\rfloor)} ,s(x)=\frac{x(x+1)}2 \\ S_2 &=& \prod_p 2^{\sum_{i=1}^n \mu(i)[p|i][p^{2}\nmid i]} \\ f(n) &=& \sum_{i=1}^n \mu(i)[p|i][p^{2}\nmid i] \\ &=& \sum_{i=1}^{\lfloor \frac np \rfloor} \mu(ip) [i \perp p] \\ &=& -\sum_{i=1}^{\lfloor \frac np \rfloor} \mu(i) [i \perp p] \\ &=& -\sum_{k\ge 1}\sum_{i=1}^{\lfloor \frac n{p^k} \rfloor} \mu(i) \end{eqnarray} $$ 对于 S₁ 和 S₂，$p \le \sqrt n$ 的部分直接暴力计算，设剩下的部分是 S_1^′ 和 S_2^′。 $$ \begin{eqnarray} S_1'&=& 2^{\sum_{p>\sqrt n} p\cdot s(\lfloor \frac n p \rfloor)} \\ S_2'&=& 2^{-\sum_{p>\sqrt n}\sum_{i=1}^{\lfloor \frac np \rfloor}\mu(i)} \end{eqnarray} $$ 还需要完成的部分有：

杜教筛预处理 $\sum_{i=1}^{\lfloor \frac n x \rfloor} \mu(i)$。
Min_25 筛 / 洲阁筛中前半部分预处理 $\sum_{p \le \lfloor \frac nx \rfloor} p$ 以及 $\sum_{p \le \lfloor \frac n x \rfloor}$

代码（在 51nod 上会 CE，但不给 CE 信息，没办法。保证不 TLE，测过大数据）

#include <bits/stdc++.h>
using namespace std;
//using LL = long long;
#define FOR(i, x, y) for (decay<decltype(y)>::type i = (x), _##i = (y); i < _##i; ++i)
#define FORD(i, x, y) for (decay<decltype(x)>::type i = (x), _##i = (y); i > _##i; --i)
#ifdef zerol
#define dbg(args...) do { cout << "\033[32;1m" << #args << " -> "; err(args); } while (0)
#else
#define dbg(...)
#endif
void err() { cout << "\033[39;0m" << endl; }
template<template<typename...> class T, typename t, typename... Args>
void err(T<t> a, Args... args) { for (auto x: a) cout << x << ' '; err(args...); }
template<typename T, typename... Args>
void err(T a, Args... args) { cout << a << ' '; err(args...); }
// -----------------------------------------------------------------------------
using LL = __int128;
using ll = long long;
const LL MOD = 1000000000039, MOD1 = MOD - 1, INV2 = (MOD + 1) / 2;
namespace dujiao {
    const int M = 2E7;
    LL f[M] = {0, 1};
    ll pr[2000000], p_sz, d;
    void init() {
        static bool vis[M];
        FOR (i, 2, M) {
            if (!vis[i]) { pr[p_sz++] = i; f[i] = -1; }
            FOR (j, 0, p_sz) {
                if ((d = pr[j] * i) >= M) break;
                vis[d] = 1;
                if (i % pr[j] == 0) {
                    f[d] = 0;
                    break;
                } else f[d] = -f[i];
            }
        }
        FOR (i, 2, M) f[i] += f[i - 1];
    }
    inline LL s_fg(LL n) { return 1; }
    inline LL s_g(LL n) { return n; }

    LL N, rd[M];
    bool vis[M];
    LL go(LL n) {
        if (n < M) return f[n];
        LL id = N / n;
        if (vis[id]) return rd[id];
        vis[id] = true;
        LL& ret = rd[id] = s_fg(n);
        for (LL l = 2, v, r; l <= n; l = r + 1) {
            v = n / l; r = n / v;
            ret -= (s_g(r) - s_g(l - 1)) * go(v);
        }
        return ret;
    }
    LL solve(LL n) {
        N = n;
        memset(vis, 0, sizeof vis);
        return go(n);
    }
}

namespace min25 {
    const int M = 1E6 + 100;
    LL B, N;

    // g(x)
    inline LL pg(LL x) { return 1; }
    inline LL ph(LL x) { return x % MOD1; }
    // Sum[g(i),{x,2,x}]
    inline LL psg(LL x) { return x % MOD1 - 1; }
    inline LL psh(LL x) {
        return x * (x + 1) / 2 % MOD1 - 1;
    }

    ll pr[M], pc;
    LL sg[M], sh[M];
    void get_prime(ll n) {
        static bool vis[M]; pc = 0;
        FOR (i, 2, n + 1) {
            if (!vis[i]) {
                pr[pc++] = i;
                sg[pc] = (sg[pc - 1] + pg(i)) % MOD1;
                sh[pc] = (sh[pc - 1] + ph(i)) % MOD1;
            }
            FOR (j, 0, pc) {
                if (pr[j] * i > n) break;
                vis[pr[j] * i] = 1;
                if (i % pr[j] == 0) break;
            }
        }
    }

    LL w[M];
    LL id1[M], id2[M], h[M], g[M];
    inline LL id(LL x) { return x <= B ? id1[x] : id2[N / x]; }

    void init() { get_prime(M); }
    void solve(LL _N) {
        N = _N;
        B = sqrt(N + 0.5);
        int sz = 0;
        for (LL l = 1, v, r; l <= N; l = r + 1) {
            v = N / l; r = N / v;
            w[sz] = v; g[sz] = psg(v); h[sz] = psh(v);
            if (v <= B) id1[v] = sz; else id2[r] = sz;
            sz++;
        }
        FOR (k, 0, pc) {
            LL p = pr[k];
            FOR (i, 0, sz) {
                LL v = w[i]; if (p * p > v) break;
                LL t = id(v / p);
                g[i] = (g[i] - (g[t] - sg[k]) * pg(p)) % MOD1;
                h[i] = (h[i] - (h[t] - sh[k]) * ph(p)) % MOD1;
            }
        }
    }
}

inline LL s(LL x) {
    return x * (x + 1) / 2 % MOD1;
}
LL bin(LL x, LL n) {
    LL ret = 1;
    for (n = (n % MOD1 + MOD1) % MOD1; n; n >>= 1, x = x * x % MOD)
        if (n & 1) ret = ret * x % MOD;
    return ret;
}

int main() {
#ifdef zerol
    freopen("in", "r", stdin);
#endif
    int T; cin >> T;
    dujiao::init(); min25::init();
    while (T--) {
        long long _n; cin >> _n;
        LL n = _n;
        dujiao::solve(n);
        min25::solve(n);

        LL b = sqrt(n), bb = (n / (n / b) + 1);

        LL ans = 1;
        FOR (i, 0, dujiao::p_sz) {
            LL p = dujiao::pr[i];
            if (p >= bb) break;
            for (LL e = 1, pp = p; pp <= n; ++e, pp *= p) {
                ans = ans *
                      bin(e + 1, pp * s(n / pp) - pp * p * s(n / pp / p))
                      % MOD;
                ans = ans * bin(INV2, dujiao::go(n / pp)) % MOD;
            }
        }
        for (LL l = bb, v, r; l <= n; l = r + 1) {
            using namespace min25;
            v = n / l; r = n / v;
            ans = ans * bin(2, s(v) * (h[id(r)] - h[id(l - 1)])) % MOD;
            ans = ans * bin(2, -dujiao::go(v) * (g[id(r)] - g[id(l - 1)])) % MOD;
        }

        long long Ans = (ans % MOD + MOD) % MOD;
        cout << Ans << endl;
    }
}