Project Euler 484

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

Project Euler 484

题目

Arithmetic Derivative

The arithmetic derivative is defined by

  • \(p' = 1\) for any prime \(p\)

  • \((ab)' = a'b + ab'\) for all integers \(a, b\) (Leibniz rule)

For example, \(20' = 24\).

Find \(\sum \mathbf{gcd}(k,k')\) for \(1 < k \le 5\times 10^{15}\).

Note: \(\mathbf{gcd}(x,y)\) denotes the greatest common divisor of \(x\) and \(y\).

解决方案

先从素数幂开始。由 Leibniz 规则可归纳得到：\((p^e)' = ep^{e-1}.\)令\(\displaystyle{g(n)=\gcd(n,n'),S(n)=\sum_{i=2}^n g(n)}\)。

因此\(g(p^e)=\gcd(p^e,(p^e)')=\gcd(p^e,e,p^{e-1})= p^{e-1}\gcd(p,e).\)

而 \(\gcd(p,e)\) 只可能是 \(1\) 或 \(p\)，所以

\[ g(p^e)= \begin{cases} p^{e-1}, & p\nmid e,\\ p^{e}, & p\mid e. \end{cases} \]

现在推广到一般形式。令 \(n=p^e m\)，其中 \(p\nmid m\)。则\(n'=(p^e m)'=(p^e)'m+p^e m' = e p^{e-1}m+p^e m'= p^{e-1}(e m + p m').\)

观察括号内对 \(p\) 的可整除性：

  • 若 \(p\nmid e\)，则 \(em\not\equiv 0\pmod p\)（因为 \(p\nmid m\)），所以\(e m + p m' \not\equiv 0\pmod p,\)从而 \(v_p(n')=e-1\)，于是 \(\gcd(n,n')\) 在 \(p\) 上的指数为 \(e-1\)。

  • 若 \(p\mid e\)，则 \(em\) 本身就含一个 \(p\)，而 \(pm'\) 也含一个 \(p\)，所以括号整体至少含一个 \(p\)，从而 \(v_p(n')\ge e\)，于是 \(\gcd(n,n')\) 在 \(p\) 上的指数为 \(e\)（但不会超过 \(e\)）。

因此对于每个素因子 \(p^e|n\)，\(\gcd(n,n')\) 的 \(p\)-部分只取决于 \((p,e)\)，并与其它素因子相互独立；于是可以知道，\(g(n)\)是乘法函数，且\(g(p^e)=p^{e-1}\gcd(p,e).\)

既然 \(g\) 是乘法函数，一个常用技巧是对每个素数 \(p\) 定义素数幂上的差分：\(f(p^0)=1, f(p^e)=g(p^e)-g(p^{e-1}), (e\ge 1).\)

那么对素数幂 \(p^e\) 之和：\(\displaystyle{\sum_{j=0}^{e} f(p^j)=g(p^e).}\)

由于 \(f,g\) 都是乘法函数，推广到一般 \(n\) 就得到 Dirichlet 卷积形式：\(\displaystyle{g(n)=\sum_{d\mid n} f(d).}\)

于是目标和可以交换求和次序：\(\displaystyle{\sum_{k\le N} g(k)=\sum_{k\le N}\sum_{d\mid k} f(d)=\sum_{d\le N} f(d)\left\lfloor \frac{N}{d}\right\rfloor.}\)

最后题目排除 \(k=1\)，所以答案是

\[ S(N)=\left(\sum_{d\le N} f(d)\left\lfloor \frac{N}{d}\right\rfloor\right)-1. \]

关键在于，什么情况 \(d\) 会满足非零的 \(f(d)\)。注意到\(g(1)=\gcd(1,0)=1, g(p)=\gcd(p,1)=1,\)所以\(f(p)=g(p)-g(1)=1-1=0.\)

也就是说：只要 \(d\) 的分解里出现某个素数的指数恰好为 \(1\)，就会让 \(f(d)=0\)。因此对 \(f(d)\ne 0\) 的 \(d\)，每个素因子指数只能是 \(0\) 或 \(\ge 2\)——这类数正是经典的 幂数：\(d\) 是幂数当且仅当对任意素数 \(p\)，若 \(p\mid d\) 则 \(p^2\mid d\)。

所以我们只需枚举 \(d\le N\) 的 幂数，并累加\(f(d)\left\lfloor \dfrac{N}{d}\right\rfloor.\)

由于幂数的结构是\(\displaystyle{d=\prod p_i^{e_i}(e_i\ge 2),}\)并且 \(p_i^{e_i}\le N\)。由于至少平方，若 \(p>\sqrt N\) 则 \(p^2>N\)，不可能出现在 \(d\) 中；因此只需要筛出所有 \(p\le \sqrt N\) 的素数。

枚举时用 DFS 逐个加入素数（按升序保证不重不漏）。为了避免反复做大整数除法，采用一个等价参数：

  • 维护 \(r=\left\lfloor \dfrac{N}{d}\right\rfloor\)；
  • 若要乘上 \(p^e\)（\(e\ge 2\)），只需检查 \(p^e\le r\)，并递归到\(r'=\left\lfloor\dfrac{r}{p^e}\right\rfloor=\left\lfloor\dfrac{N}{d\cdot p^e}\right\rfloor.\)

同时维护系数 \(c=f(d)\)。当选择乘入 \(p^e\) 时，新系数是\(c' = c \cdot f(p^e).\)DFS 到某个节点（对应某个 \(d\)）时立即贡献一项：\(c\cdot r = f(d)\left\lfloor \dfrac{N}{d}\right\rfloor.\)并继续尝试加入更大的素数。

代码

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
#include <bits/stdc++.h>
using namespace std;

typedef long long ll;

// --------- global state ----------
const ll N = 5000000000000000LL; // 5e15
vector<int> primes_global;
ll ans_global = 0;

// odd-only sieve: primes up to limit
vector<int> primes_upto(int limit)
{
    vector<int> primes;
    if (limit < 2)
        return primes;
    primes.push_back(2);

    int m = (limit >> 1) + 1; // i -> odd = 2*i+1
    vector<uint8_t> comp(m, 0);
    int r = (int)floor(sqrt((long double)limit));

    for (int odd = 3; odd <= r; odd += 2)
    {
        if (comp[odd >> 1])
            continue;
        ll step = 2LL * odd;
        ll start = 1LL * odd * odd;
        for (ll x = start; x <= limit; x += step)
        {
            comp[(int)(x >> 1)] = 1;
        }
    }
    for (int odd = 3; odd <= limit; odd += 2)
    {
        if (!comp[odd >> 1])
            primes.push_back(odd);
    }
    return primes;
}

// dfs over 幂数 d, encoded via r = floor(N/d)
// c = f(d)
void dfs(int idx, ll r, ll c)
{
    // contribution of this d:
    ans_global += c * r;

    for (int i = idx; i < primes_global.size(); ++i)
    {
        ll p = primes_global[i];
        ll p2 = p * p;
        if (p2 > r)
            break; // need p^2 <= r

        ll pw = p2; // p^e, start e=2
        ll e = 2;

        // g(p^1)=1
        ll g_prev = 1;

        while (pw <= r)
        {
            // g(p^e) = p^{e-1}*gcd(p,e)
            // since gcd(p,e) is either 1 or p:
            ll g_e = pw / p; // p^{e-1}
            if (e % p == 0)
                g_e *= p; // => p^e

            ll delta = g_e - g_prev; // f(p^e)=g(p^e)-g(p^{e-1})
            if (delta != 0)
            {
                dfs(i + 1, r / pw, c * delta);
            }
            g_prev = g_e;

            if (pw > r / p)
                break;
            pw *= p;
            ++e;
        }
    }
}

ll solve(ll N)
{
    ans_global = 0;

    int lim = (int)floor(sqrt((long double)N));
    primes_global = primes_upto(lim);

    // start with d=1, f(1)=1, r=N
    dfs(0, N, 1);

    // subtract gcd(1,1')=1 because problem asks 1<k<=N
    return ans_global - 1;
}

int main()
{
    ios::sync_with_stdio(false);
    cin.tie(nullptr);

    ll res = solve(N);
    cout << res << "\n";
    return 0;
}