Project Euler 362

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

Project Euler 362

题目

Squarefree factors

Consider the number \(54\).

\(54\) can be factored in \(7\) distinct ways into one or more factors larger than \(1\):

\(54, 2\times27, 3\times18, 6\times9, 3\times3\times6, 2\times3\times9\) and \(2\times3\times3\times3\).

If we require that the factors are all squarefree only two ways remain: \(3\times3\times6\) and \(2\times3\times3\times3\).

Let’s call \(\text{Fsf}(n)\) the number of ways \(n\) can be factored into one or more squarefree factors larger than \(1\), so \(\text{Fsf}(54)=2\).

Let \(S(n)\) be \(\sum \text{Fsf}(k)\) for \(k=2\) to \(n\).

\(S(100)=193\).

Find \(S(10 000 000 000)\).

解决方案

令\(N=10^{10}\)。注意：每一种“平方自由因子分解”\(k=a_1a_2\cdots a_t, 2\le a_1\le a_2\le\cdots\le a_t, a_i \text{ squarefree}\)都会给 \(\mathrm{Fsf}(k)\) 贡献 1；因此它也会给 \(S(N)\) 贡献 1 当且仅当 \(k\le N\)。

所以：\(S(N)\) 等于所有满足\(2\le a_1\le a_2\le \cdots\le a_t, a_i \text{ squarefree}, \prod a_i\le N\)的序列（长度 \(t\ge 1\)）的总数。

这就把“对 \(k\) 求和”变成了“数所有合法序列”。

定义 \(G(x,m)\)：所有满足

  • 因子都是 squarefree 且 \(\ge m\)，
  • 因子按非递减排列，
  • 乘积 \(\le x\)，

的序列数量，允许空序列（空序列代表“不再选因子”）。

那么我们要的答案是\(S(N)=G(N,2)-1\)（减去空序列那 1 个）。

如果选择第一个因子为 \(f\)（要求 \(f\ge m\)、squarefree、且 \(f\le x\)），后面还要选一个非递减序列，乘积上限变成 \(\left\lfloor x/f\right\rfloor\)，且最小因子变为 \(f\)。

可以写出关于\(G\)的递推式：

\[ G(x,m)=1+\sum_{\substack{f \text{ squarefree}\\m\le f\le x}} G\left(\left\lfloor\frac{x}{f}\right\rfloor,f\right). \]

其中前面的 “1” 就是“空序列”。

但是这个式子不能够直接使用，因为要遍历到 \(x\)（这里 \(x\) 可达 \(N\)）。接下来开始进行优化。

若 \(f>\sqrt{x}\)，则在子问题 \(G(\lfloor x/f\rfloor,f)\) 里，任何再选的因子都 \(\ge f>\sqrt{x}\)，两次相乘必定 \(>x\)，所以子问题只能是“空序列”——也就是说每个这样的 \(f\) 都只贡献 1。

因此当 \(m\le \lfloor\sqrt{x}\rfloor\) 时可改写：

\[ G(x,m)=1+\sum_{\substack{m\le f\le \lfloor\sqrt{x}\rfloor \\f \text{ squarefree}}} G\left(\left\lfloor\frac{x}{f}\right\rfloor,f\right) +\#\{f: \lfloor\sqrt{x}\rfloor<f\le x, f \text{ squarefree}\}. \]

上式最后一项只需要“数区间内 squarefree 的个数”。

考虑子问题 \(G(x',f)\) 其中 \(x'=\lfloor x/f\rfloor\)。若 \(f>\sqrt[3]{x}\)，则 \(f^2>x/f\ge x'\)（忽略整除误差不影响结论），于是 \(f>\sqrt{x'}\)。这意味着在子问题中最多只能再选一个因子（因为两个 \(\ge f\) 的因子相乘就超过上限了）。

因此：

  • 空序列：1 种；
  • 选 1 个因子 \(g\)（squarefree 且 \(f\le g\le x'\)）：每个 \(g\) 贡献 1。

所以\(G(x',f)=1+\#\{g: f\le g\le x', g \text{ squarefree}\}.\)这让大量递归调用变成了一个区间计数。

设 \(Q(n)\) 表示 \(\le n\) 的 squarefree 数量（包含 1）。那么有\(\displaystyle{Q(n)=\sum_{k=1}^{\lfloor\sqrt{n}\rfloor}\mu(k)\left\lfloor\frac{n}{k^2}\right\rfloor,}\)其中 \(\mu\) 是 Möbius 函数。

对本题 \(n\le N\)，\(\sqrt{n}\le \lfloor\sqrt{N}\rfloor\)，我们可以先线性筛出 \(\mu(1\dots\lfloor\sqrt{N}\rfloor)\)，并做前缀和\(\displaystyle{M(t)=\sum_{k\le t}\mu(k).}\)

然后用分块：当 \(v=\left\lfloor n/k^2\right\rfloor\) 固定时，\(k\) 落在一个连续区间\(k\in\left[L,R\right], R=\left\lfloor\sqrt{\dfrac{n}{v}}\right\rfloor.\)于是\(\displaystyle{\sum_{k=L}^{R}\mu(k)\left\lfloor\frac{n}{k^2}\right\rfloor= v\cdot\left(M(R)-M(L-1)\right)}.\)

并且这种不同的 \(v\) 的个数大约是 \(O(N^{1/3})\)，所以一次 \(Q(n)\) 很快。因此，区间 \([A,B]\) 内 squarefree 个数就是 \(Q(B)-Q(A-1)\)。

最终，令 \(\ell=\lfloor\sqrt{x}\rfloor\)，\(c=\lfloor\sqrt[3]{x}\rfloor\)。

  • 若 \(m>\ell\)：\(G(x,m)=1+Q(x)-Q(m-1).\)
  • 否则：\(\displaystyle{G(x,m)=1+\bigl(Q(x)-Q(\ell)\bigr)+\sum_{\substack{m\le f\le \min(\ell,c)\\ f \text{ squarefree}}}G(\lfloor x/f\rfloor,f)+\sum_{\substack{\max(m,c+1)\le f\le \ell\\ f \text{ squarefree}}}\left(1+Q(\lfloor x/f\rfloor)-Q(f-1)\right).}\)

最后一大段用按 \(\lfloor x/f\rfloor\) 分块把 \(\sum Q(\lfloor x/f\rfloor)\) 聚合，避免逐个 \(f\) 递归计算。

代码

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
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
#include <bits/stdc++.h>
#include "tools.h"
using namespace std;

typedef long long ll;
ll N = 10000000000LL;

ll icbrtll(ll n)
{
    ll x = (ll)llround(powl((long double)n, 1.0L / 3.0L));
    while ((x + 1) * (x + 1) * (x + 1) <= n)
        ++x;
    while (x * x * x > n)
        --x;
    return x;
}

// pack key for memo: (x << 17) | m, since m <= 1e5 < 2^17 and x <= 1e10 < 2^34
ll pack_key(ll x, int m)
{
    return (ll(x) << 17) | (ll)m;
}

// ---------- globals ----------
ll GN;   // N in problem
int LIM; // floor(sqrt(N))

vector<int> mu;                // mu[1..LIM]
vector<int> mu_pref;           // prefix sum of mu
vector<char> is_sqfree;        // squarefree up to LIM
vector<int> sf_pref;           // Q_small prefix counts including 1
vector<int> sq_list;           // squarefree numbers >=2 up to LIM
vector<ll> sq_Qm1_pref;        // prefix sum over squarefree f: Q(f-1) for f in sq_list (small only)
unordered_map<ll, ll> Q_cache; // Q(x) for x > LIM
unordered_map<ll, ll> G_cache; // G(x,m)

// ---------- build mobius up to LIM ----------
void build_mobius(int n)
{
    mu.assign(n + 1, 0);
    vector<int> primes;
    vector<int> spf(n + 1, 0);
    mu[1] = 1;
    for (int i = 2; i <= n; ++i)
    {
        if (spf[i] == 0)
        {
            spf[i] = i;
            primes.push_back(i);
            mu[i] = -1;
        }
        for (int p : primes)
        {
            ll v = 1LL * i * p;
            if (v > n)
                break;
            spf[(int)v] = p;
            if (i % p == 0)
            {
                mu[(int)v] = 0;
                break;
            }
            mu[(int)v] = -mu[i];
        }
    }
    mu_pref.assign(n + 1, 0);
    int s = 0;
    for (int i = 1; i <= n; ++i)
    {
        s += mu[i];
        mu_pref[i] = s;
    }
}

// ---------- build squarefree prefix up to LIM ----------
void build_squarefree_prefix(int n)
{
    is_sqfree.assign(n + 1, true);
    is_sqfree[0] = false;
    for (int p = 2; 1LL * p * p <= n; ++p)
    {
        int pp = p * p;
        for (int k = pp; k <= n; k += pp)
            is_sqfree[k] = false;
    }
    sf_pref.assign(n + 1, 0);
    int s = 0;
    for (int i = 0; i <= n; ++i)
    {
        if (is_sqfree[i])
            ++s;
        sf_pref[i] = s;
    }
    sq_list.clear();
    for (int i = 2; i <= n; ++i)
        if (is_sqfree[i])
            sq_list.push_back(i);

    sq_Qm1_pref.assign(sq_list.size() + 1, 0);
    for (size_t i = 0; i < sq_list.size(); ++i)
    {
        int f = sq_list[i];
        // for f <= LIM, Q(f-1) = sf_pref[f-1]
        sq_Qm1_pref[i + 1] = sq_Qm1_pref[i] + sf_pref[f - 1];
    }
}

// ---------- Q(x): count squarefree <= x INCLUDING 1 ----------
ll Q(ll x)
{
    if (x <= LIM)
        return sf_pref[(int)x];
    auto it = Q_cache.find(x);
    if (it != Q_cache.end())
        return it->second;

    ll res = 0;
    ll k = 1;
    while (k * k <= x)
    {
        ll v = x / (k * k);
        ll kmax = int_square(x / v);
        res += v * ((ll)mu_pref[(int)kmax] - (ll)mu_pref[(int)k - 1]);
        k = kmax + 1;
    }
    Q_cache.emplace(x, res);
    return res;
}

int count_sqfree_interval_small(int a, int b)
{
    if (a > b)
        return 0;
    return sf_pref[b] - sf_pref[a - 1];
}

ll sum_Qm1_interval_small(int a, int b)
{
    if (a > b)
        return 0;
    auto l = lower_bound(sq_list.begin(), sq_list.end(), a) - sq_list.begin();
    auto r = upper_bound(sq_list.begin(), sq_list.end(), b) - sq_list.begin();
    return sq_Qm1_pref[r] - sq_Qm1_pref[l];
}

// ---------- G(x,m): # nondecreasing sequences of squarefree factors >=m, product<=x, allowing empty ----------
ll G(ll x, int m)
{
    ll key = pack_key(x, m);
    auto it = G_cache.find(key);
    if (it != G_cache.end())
        return it->second;

    ll lim = int_square(x);
    ll res;

    if (m > lim)
    {
        res = 1 + (Q(x) - Q((ll)m - 1));
        G_cache.emplace(key, res);
        return res;
    }

    // empty + choose a single factor > sqrt(x)
    res = 1 + (Q(x) - Q(lim));

    ll cbr = icbrtll(x);
    int upper = (int)min(lim, cbr);

    // recursive part: m <= f <= upper (squarefree)
    if (m <= upper)
    {
        auto li = lower_bound(sq_list.begin(), sq_list.end(), m) - sq_list.begin();
        auto ri = upper_bound(sq_list.begin(), sq_list.end(), upper) - sq_list.begin();
        for (size_t idx = li; idx < ri; ++idx)
        {
            int f = sq_list[idx];
            res += G(x / f, f);
        }
    }

    // compressed part: upper < f <= lim  (here lim <= LIM always)
    int A = max(m, upper + 1);
    int B = (int)lim;
    if (A <= B)
    {
        int cnt = count_sqfree_interval_small(A, B);
        ll sum_qm1 = sum_Qm1_interval_small(A, B);

        // sum_{squarefree f in [A,B]} Q(x/f), grouped by q = floor(x/f)
        ll sum_qq = 0;
        int t = A;
        while (t <= B)
        {
            ll q = x / t;
            int tmax = (int)min<ll>(B, x / q);
            int c = count_sqfree_interval_small(t, tmax);
            if (c)
                sum_qq += Q(q) * (ll)c;
            t = tmax + 1;
        }

        // add sum_f (1 + Q(x/f) - Q(f-1))
        res += (ll)cnt + sum_qq - sum_qm1;
    }

    G_cache.emplace(key, res);
    return res;
}

// S(N) = sum_{k=2..N} Fsf(k) = G(N,2) - 1
ll S(ll N)
{
    GN = N;
    LIM = (int)int_square(GN);
    build_mobius(LIM);
    build_squarefree_prefix(LIM);

    return G(GN, 2) - 1;
}

int main()
{
    cout << S(N) << "\n";
    return 0;
}