CF542D Superhero's Job Solution

更好的阅读体验戳此进入

题面

存在 $J(x)=\sum _{k\mid x,gcd(x,\frac{x}{k})=1}k$$J(x) = \sum_{k \mid x, \gcd(x, \frac{x}{k}) = 1}k$，给定 $A$$A$ 求满足 $J(x)=A$$J(x) = A$ 的正整数解的个数。

Solution

首先 $J(x)$$J(x)$ 是个不完全积性函数，这里尝试（可能不严谨）详细地说一下：

首先考虑若我们将 $x$$x$ 质因数分解为 $x=\prod p_i^{k_i}$$x = \prod p_i^{k_i}$ 的形式，那么枚举其因数就相当于是枚举其所有质因子所有不同幂次能够组成的所有数 $k$$k$，而对应的 ${\displaystyle \frac{x}{k}}$$\dfrac{x}{k}$ 也就是剩下的质因数合起来，于是问题转化为将所有质因数划分为两个可空的集合，而考虑 $k\perp {\displaystyle \frac{x}{k}}$$k \perp \dfrac{x}{k}$ 的条件，其可以转化为两数质因数分解后无相同质因子，于是我们不难发现对于 $x$$x$ 分解出的 $p_i$$p_i$，每一个集合中要么有全部的 $p_i$$p_i$ 要么没有 $p_i$$p_i$，否则两个集合将会有公因子 $p_i$$p_i$。于是最终问题转化为将质因子种类数进行可空地划分。这时我们想到，如果 $a\perp b$$a \perp b$，那么 $a,b$$a, b$ 质因数分解出来的质因子一定完全不同，那么假设 $ a $ 有 $x$$x$ 种，$b$$b$ 有 $y$$y$ 种，$a\times b$$a \times b$ 就一定有 $x+y$$x + y$ 种，所以从 $x$$x$ 种里选择 $[0,x]$$[0, x]$ 的所有方案乘上 $y$$y$ 种里选择 $[0,y]$$[0, y]$ 的所有方案一定等价于从 $x+y$$x + y$ 种里选择 $[0,x+y]$$[0, x + y]$ 的方案，这个可以从类似组合意义的思想上去考虑，所以一定有 $J(a\times b)=J(a)\times J(b)(a\perp b)$$J(a \times b) = J(a) \times J(b) (a \perp b)$。

于是不难发现对于质数 $p$$p$，有 $J(p)=p+1$$J(p) = p + 1$，更具体地说应为 $J(p^k)=p^k+1$$J(p^k) = p^k + 1$，所以当我们对 $x$$x$ 分解为 $x=\prod p_i^{k_i}$$x = \prod p_i^{k_i}$ 时 $J(x)$$J(x)$ 就变为了 $J(x)=\prod (p_i^{k_i}+1)$$J(x) = \prod (p_i^{k_i} + 1)$。

至此不难想到，我们现在已知 $J(x)=A$$J(x) = A$，那么只需要找到所有 $A=\prod (p_i^{k_i}+1)$$A = \prod (p_i^{k_i} + 1)$ 的分解方案，就会一一对应着 $x$$x$ 的分解方案，也就是说分解方案的个数就是最终的答案。

考虑如何分解，显然对于 ${10}^{12}$$10^{12}$ 的范围根号分治是一个显而易见的想法，考虑枚举 ${10}^6$$10^6$ 以内的所有质数，及其范围内的若干次方丢到 unordered_map 里，同时需要记录基于的质数，而若其存在 ${10}^6$$10^6$ 以外的因子，那么一定最多为一次方，反之 $p^2+1$$p^2 + 1$ 就会超过 ${10}^{12}$$10^{12}$ 从而超过范围。

于是可以考虑根号枚举 $A$$A$ 的所有因数，${10}^6$$10^6$ 内的直接查一下之前预处理的里面有没有，反之判断其减 $1$$1$ 后是否为质数即可。

此时可以通过背包处理，即 $dp(i,j)$$dp(i, j)$ 表示考虑前 $i$$i$ 个质数（及其若干次方）能够凑成 $j$$j$ 的方案数，转移显然，同时注意整个 $dp$$dp$ 数组是需要离散下来的，也就是可以通过 unordered_map 维护，而记录基于某个质数便是为了在转移的时候不出现同时使用多个质数的情况，所以对每个质数开个 basic_string 记录其所有合法的 $p^k$$p^k$，具体来说这里也需要离散化，可以套个 unordered_map。

尝试分析复杂度，考虑 DP 时，显然枚举素数时最多有 $d(n)$$d(n)$，我们在 DP 过程中第二维枚举直接用前一维的状态和新的质数合并，这样合并之后的数一定亦为 $n$$n$ 的因子，故这一维也是 $d(n)$$d(n)$，则 DP 的复杂度为 $d^2(n)$$d^2(n)$，前面不论是线性筛亦或是枚举因数均为 $\sqrt{n}$$\sqrt{n}$ 的，对于 $x$$x$ 的 Miller-Rabin 朴素地进行 $k$$k$ 轮测试的复杂度为 $k\log x$$k\log x$，我们会对约 ${\displaystyle \frac{d(n)}{2}}$$\dfrac{d(n)}{2}$ 个数进行测试。则最终复杂度 $O(\sqrt{n}+d(n)(d(n)+k\log n))$$O(\sqrt{n} + d(n)(d(n) + k\log n))$。

Tips：注意计算过程中包括但不限于 Miller-Rabin 有多处会超过 long long 范围，需要开 __int128_t 或用 long double 快速乘，可以通过 Sanitizer 检测 signed-integer-overflow 并输入最大的数据进行测试寻找报错。

Upd：注意在 Miller-Rabin 中若我们使用固定的基底，那么验证时需要对于基底进行 $a\leftarrow amodn$$a \leftarrow a \bmod{n}$，并且若 $a\equiv 0(\mathrm{mod}n)$$a \equiv 0 \pmod{n}$，那么应该直接认为其通过本轮测试，否则在进行检测固定基底的时候会概率出错。

Code

x
1
#define _USE_MATH_DEFINES
2
#include <bits/stdc++.h>
3

4
#define PI M_PI
5
#define E M_E
6
#define npt nullptr
7
#define SON i->to
8
#define OPNEW void* operator new(size_t)
9
#define ROPNEW void* Edge::operator new(size_t){static Edge* P = ed; return P++;}
10
#define ROPNEW_NODE void* Node::operator new(size_t){static Node* P = nd; return P++;}
11

12
using namespace std;
13

14
mt19937 rnd(random_device{}());
15
int rndd(int l, int r){return rnd() % (r - l + 1) + l;}
16
bool rnddd(int x){return rndd(1, 100) <= x;}
17

18
typedef unsigned int uint;
19
typedef unsigned long long unll;
20
typedef long long ll;
21
typedef long double ld;
22

23
#define LIM (1100000)
24

25
template < typename T = int >
26
inline T read(void);
27

28
ll N;
29
unordered_map < ll, ll > dp[2];
30
basic_string < ll > Prime;
31
unordered_map < ll, ll > pows;
32
bitset < LIM + 100 > notPrime;
33
unordered_map < ll, basic_string < ll > > tpow;
34
basic_string < ll > tot;
35

36
namespace MillerRabin{
37
    int radix[10] = {0,  2, 3, 5, 7, 11, 13, 19, 17};
38
    ll qpow(ll a, ll b, ll MOD){
39
        ll ret(1), mul(a);
40
        while(b){
41
            if(b & 1)ret = (__int128_t)ret * mul % MOD;
42
            b >>= 1;
43
            mul = (__int128_t)mul * mul % MOD;
44
        }return ret;
45
    }
46
    bool Check(ll N){
47
        if(N <= 2 || !(N & 1))return N == 2;
48
        ll base(N - 1); int cpow(0);
49
        while(base % 2 == 0)base >>= 1, ++cpow;
50
        for(int t = 1; t <= 8; ++t){
51
            if(radix[t] % N == 0)continue;
52
            ll cur = qpow(radix[t] % N, base, N);
53
            if(cur == 1)continue;
54
            for(int i = 1; i <= cpow; ++i){
55
                if(cur == N - 1)break;
56
                cur = (__int128_t)cur * cur % N;
57
                if(i == cpow)return false;
58
            }
59
        }return true;
60
    }
61
};
62

63
void Sieve(void){
64
    for(ll i = 2; i <= LIM; ++i){
65
        if(!notPrime[i])Prime += i;
66
        for(auto p : Prime){
67
            if(p * i > LIM)break;
68
            notPrime[p * i] = true;
69
            if(i % p == 0)break;
70
        }
71
    }
72
}
73
void Init(void){
74
    Sieve();
75
    for(auto p : Prime){
76
        ll base(p);
77
        while(base <= N)pows[base] = p, base *= p;
78
    }
79
    for(ll i = 1; i * i <= N; ++i){
80
        if(N % i)continue;
81
        if(pows.find(i - 1) != pows.end())
82
            tpow[pows[i - 1]] += i, tot += pows[i - 1];
83
        if(i * i != N){
84
            if(pows.find(N / i - 1) != pows.end())
85
                tpow[pows[N / i - 1]] += N / i, tot += pows[(N / i) - 1];
86
            else if(MillerRabin::Check(N / i - 1))
87
                tpow[N / i - 1] += N / i, tot += N / i - 1;
88
        }
89
    }
90
    sort(tot.begin(), tot.end());
91
    tot.erase(unique(tot.begin(), tot.end()), tot.end());
92
}
93

94
int main(){
95
    N = read < ll >();
96
    Init();
97
    dp[0][1] = 1;
98
    int base(0);
99
    for(auto p : tot){
100
        base ^= 1, dp[base] = dp[base ^ 1];
101
        for(auto lst : dp[base ^ 1])
102
            for(auto pri : tpow[p])
103
                if((__int128_t)pri * lst.first <= N && N % (pri * lst.first) == 0)
104
                    dp[base][pri * lst.first] += lst.second;
105
    }printf("%lld\n", dp[base][N]);
106
    fprintf(stderr, "Time: %.6lf\n", (double)clock() / CLOCKS_PER_SEC);
107
    return 0;
108
}
109

110
template < typename T >
111
inline T read(void){
112
    T ret(0);
113
    int flag(1);
114
    char c = getchar();
115
    while(c != '-' && !isdigit(c))c = getchar();
116
    if(c == '-')flag = -1, c = getchar();
117
    while(isdigit(c)){
118
        ret *= 10;
119
        ret += int(c - '0');
120
        c = getchar();
121
    }
122
    ret *= flag;
123
    return ret;
124
}

UPD

update-2023_02_10 初稿

update-2023_02_10 修正一些不符合规范的错误

update-2023_02_10 修正 Miller-Rabin 中的错误