Project Euler 861

发表于 分类于 阅读次数： 本文字数： 2.7k 阅读时长 ≈ 2 分钟

Project Euler 861

题目

Products of Bi-Unitary Divisors

A unitary divisor of a positive integer \(n\) is a divisor \(d\) of \(n\) such that \(\gcd\left(d,\dfrac{n}{d}\right)=1\).

A bi-unitary divisor of \(n\) is a divisor \(d\) for which 1 is the only unitary divisor of \(d\) that is also a unitary divisor of \(\dfrac{n}{d}\).

For example, \(2\) is a bi-unitary divisor of \(8\), because the unitary divisors of \(2\) are \(\{1,2\}\), and the unitary divisors of \(8/2\) are \(\{1,4\}\), with \(1\) being the only unitary divisor in common.

The bi-unitary divisors of \(240\) are \(\{1,2,3,5,6,8,10,15,16,24,30,40,48,80,120,240\}\).

Let \(P(n)\) be the product of all bi-unitary divisors of \(n\). Define \(Q_k(N)\) as the number of positive integers \(1 \lt n \leq N\) such that \(P(n)=n^k\). For example, \(Q_2\left(10^2\right)=51\) and \(Q_6\left(10^6\right)=6189\).

Find \(\sum_{k=2}^{10}Q_k\left(10^{12}\right)\).

解决方案

设\(n\)的质因子分解为\(\displaystyle n=\prod_{i=1}^r p_i^{e_i},d=\prod_{i=1}^r p_i^{a_i}\mid n.\)

先看单个素幂 \(p^e\)。对任意 \(0\le a\le e\)：

  • \(p^a\) 的 unitary divisors 仅有 \({1,p^a}\)；
  • \(p^{e-a}\) 的 unitary divisors 仅有 \({1,p^{e-a}}\)。

因此，\(p^a\) 与 \(p^{e-a}\) 的公共 unitary divisor 除了 \(1\) 以外，还会出现额外公共项，当且仅当\(p^a=p^{e-a}\iff 2a=e.\)这就得到：

  • 若 \(e\) 为奇数，\(a=e/2\) 不存在，故全部 \(e+1\) 个约数都 bi-unitary；
  • 若 \(e\) 为偶数，恰好剔除 \(a=e/2\)，故有 \(e\) 个 bi-unitary 约数。

记 \(\tau^{**}(n)\) 为 bi-unitary 约数个数，则\(\tau^{**}(p^e)=e+(e\bmod 2).\) 又因各素因子独立，\(\tau^{**}\) 乘法性成立，于是 \[ \tau^{**}(n)=\prod_{i=1}^r\bigl(e_i+(e_i\bmod 2)\bigr). \]

题设 \(P(n)\) 是所有 bi-unitary 约数之积。若 \(d\) 是 bi-unitary 约数，则 \(n/d\) 也是，且二者配对乘积恒为 \(n\)。该配对没有不动点：\(d=n/d\) 要求每个指数取一半，而这正对应偶指数时被禁止的中间指数。

故所有 bi-unitary 约数可以两两配对，得到\(P(n)=n^{\tau^{**}(n)/2}.\)于是\(P(n)=n^k\iff \tau^{**}(n)=2k.\)所以\(Q_k(N)=\#\{1<n\le N:\tau^{**}(n)=2k\}.\)

问题即化为统计 \(\tau^{**}(n)\in\{4,6,\dots,20\}\) 的整数个数。把任意 \(n\) 唯一写成\(n=u\cdot s, \gcd(u,s)=1,\)其中\(u\) 的每个素因子指数都 \(\ge2\)（powerful 部分），\(s\) 是无平方因子数。

设 \(\omega(s)=m\)，则 \(s\) 的每个素因子指数为 \(1\)，对应局部因子都是 \(2\)，所以\(\tau^{**}(s)=2^m.\)因此\(\tau^{**}(n)=\tau^{**}(u)2^m.\)对固定 \(k\)，约束是\(\tau^{**}(u)2^m=2k.\)因为 \(2k\le20\)，所以 \(m\) 只能很小（至多 \(4\)），计数规模因此大幅收缩。

记\(c=\tau^{**}(u), L=\left\lfloor\dfrac N u\right\rfloor, U\)是\(u\)的质因子集合。则要求\(c\cdot2^m=2k.\)若 \(\dfrac{2k}{c}\) 不是 \(2\) 的幂，则该 \(u\) 对 \(Q_k\) 贡献为 \(0\)；若 \(\dfrac{2k}{c}=2^m\)：\(m=0\)：只可能 \(s=1\)，贡献 \(1\)；\(m\ge1\)：贡献为\(\displaystyle\#\left\{q_1<\cdots<q_m:\ q_i\notin U,\prod_{i=1}^m q_i\le L\right\}.\)本质上是在计数：从不整除 \(u\) 的素数中选出 \(m\) 个互异素数，其乘积不超过 \(L\)。

在上述计数中，递归会反复用到区间内素数个数\(\#\{p:p\in \Omega,a<p\le b\}=\pi(b)-\pi(a),\)其中\(\Omega\)是质数集合。

因此核心在于高效查询 \(\pi(x)\)。为此，先对所有可能出现的查询点做一次离线预处理：\(\{1,2,\dots,\lfloor\sqrt N\rfloor\}\cup\left\{\left\lfloor\dfrac Ni\right\rfloor:1\le i\le\lfloor\sqrt N\rfloor\right\}.\) 利用按 \(\lfloor N/i\rfloor\) 分块的素数计数递推（hyperbola 分块法），即可一次性得到该集合上的全部 \(\pi\) 值，从而将后续查询降为 \(O(1)\)。

又因为本题中 \(m\) 很小（至多为 \(4\)），组合计数只需浅层递归；其边界项可写为\(\pi(r)-\pi(l)-\#(U\cap(l,r]).\)

最终处理好后反向计算出对于每个\(k\)，\(\tau^{\ast\ast}(n)=2k\)的数量。最后把数量累计即可。

代码

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
#include <bits/stdc++.h>
#include "tools.h"

using namespace std;
typedef long long ll;

ll N = 1000000000000LL;
int K1, K2, R, LIMCB;
vector<int> primes;
vector<int> small_pi;
vector<ll> large_pi;
ll ans_sum;

void init(ll n)
{
    N = n;
    R = (int)floor(sqrt((long double)N));
    small_pi.assign(R + 1, 0);
    large_pi.assign(R + 1, 0);
    for (int i = 1; i <= R; i++)
    {
        small_pi[i] = i - 1;
        large_pi[i] = N / i - 1;
    }
    for (int p = 2; p <= R; p++)
    {
        if (small_pi[p] == small_pi[p - 1])
            continue;
        int sp = small_pi[p - 1];
        ll p2 = 1LL * p * p;
        int end = (int)min((ll)R, N / p2);
        for (int i = 1; i <= end; i++)
        {
            ll d = 1LL * i * p;
            ll v = (d <= R ? large_pi[(int)d] : small_pi[(int)(N / d)]);
            large_pi[i] -= (v - sp);
        }
        if (p2 <= R)
        {
            for (int i = R; i >= (int)p2; i--)
            {
                small_pi[i] -= (small_pi[i / p] - sp);
            }
        }
    }
}

inline ll pi_count(ll x)
{
    if (x <= 1)
        return 0;
    if (x >= N)
        return large_pi[1];
    if (x <= R)
        return small_pi[(int)x];
    return large_pi[(int)(N / x)];
}

inline bool used_has(const vector<int> &used, int p)
{
    for (int q : used)
        if (q == p)
            return true;
    return false;
}

ll count_sqfree(const vector<int> &used, int last, ll limit, int m)
{
    if (limit <= 1)
        return 0;
    if (m == 1)
    {
        ll res = pi_count(limit) - pi_count(last);
        if (res <= 0)
            return 0;
        for (int p : used)
        {
            if ((ll)p > last && (ll)p <= limit)
                res--;
        }
        return max(0ll, res);
    }
    ll res = 0;
    auto it = upper_bound(primes.begin(), primes.end(), last);
    for (; it != primes.end(); ++it)
    {
        int p = *it;
        if (1LL * p * p > limit)
            break;
        if (used_has(used, p))
            continue;
        res += count_sqfree(used, p, limit / p, m - 1);
    }
    return res;
}

inline int log2_if_pow2(ll x)
{
    if (x <= 0 || (x & (x - 1)) != 0)
        return -1;
    return __builtin_ctzll((unsigned long long)x);
}

void process_powerful(ll pw, ll cb, const vector<int> &used)
{
    if (cb <= 0 || cb > LIMCB)
        return;
    ll lim = N / pw;
    for (int k = K1; k <= K2; k++)
    {
        ll t = 2LL * k;
        if (t % cb)
            continue;
        ll x = t / cb;
        int m = log2_if_pow2(x);
        if (m < 0)
            continue;
        if (m == 0)
            ans_sum += 1;
        else
            ans_sum += count_sqfree(used, 1, lim, m);
    }
}

void dfs_powerful(int idx, ll nRem, ll curVal, ll curCb, vector<int> &used)
{
    if (curCb > LIMCB)
        return;
    process_powerful(curVal, curCb, used);
    for (int i = idx; i < (int)primes.size(); i++)
    {
        int p = primes[i];
        ll pp = 1LL * p * p;
        if (pp > nRem)
            break;
        ll nn = nRem / pp;
        ll cc = curVal * pp;
        int e = 2;
        ll cb = curCb;
        used.push_back(p);
        while (true)
        {
            int mult = e + (e & 1);
            if (cb > LIMCB / mult)
                break;
            cb *= mult;
            dfs_powerful(i + 1, nn, cc, cb, used);
            cb /= mult;
            if (nn < p)
                break;
            nn /= p;
            cc *= p;
            e++;
        }
        used.pop_back();
    }
}

ll solve_sum(ll n, int k1, int k2)
{
    K1 = k1;
    K2 = k2;
    LIMCB = 2 * K2;
    ans_sum = 0;
    init(n);
    Factorizer fac(R);
    primes = fac.primes;
    vector<int> used;
    dfs_powerful(0, n, 1, 1, used);
    return ans_sum;
}

int main()
{
    cout << solve_sum(N, 2, 10) << '\n';
    return 0;
}