Project Euler 432

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

<– more –>

Project Euler 432

题目

Totient sum

Let \(S(n,m) = \sum \varphi(n \times i)\) for \(1 \le i \le m\). (\(\varphi\) is Euler’s totient function)

You are given that \(S(510510,10^6 )= 45480596821125120\).

Find \(S(510510,10^{11})\).

Give the last \(9\) digits of your answer.

解决方案

首先介绍两个基本恒等式：

  1. 根据欧拉函数的定义，有\(\displaystyle{\sum_{d\mid x}\varphi(d)=x}\)
  2. \(\varphi(ab)\)的按 \(\gcd\)分解形式为\(\displaystyle{\varphi(ab)=\varphi(b)\sum_{d\mid\gcd(a,b)}\varphi\left(\frac{a}{d}\right)}\)。

接下来证明第二条恒等式，只需在素数幂上验证，再利用积性推广即可。令 \(a=p^\alpha,b=p^\beta\)，那么 \(\gcd(a,b)=p^{\min(\alpha,\beta)}\)。

第二条恒等式的右侧为\(\displaystyle{\varphi(p^\beta)\sum_{t=0}^{\min(\alpha,\beta)}\varphi(p^{\alpha-t})}\)。注意 \(\varphi(p^k)=p^{k-1}(p-1)(k\ge1)\)，并且 \(\varphi(1)=1\)。于是\(\displaystyle{\sum_{t=0}^{\min(\alpha,\beta)}\varphi(p^{\alpha-t})= \varphi(p^\alpha)+\varphi(p^{\alpha-1})+\cdots+\varphi(p^{\alpha-\min(\alpha,\beta)})}.\)

而这串求和恰好等于 \(p^\alpha\)（因为 \(\sum_{d\mid p^\alpha}\varphi(d)=p^\alpha\)，取其中连续一段即可与 \(\gcd\) 对齐），最终化简得到：\(\varphi(p^\beta)\cdot p^{\min(\alpha,\beta)}=\varphi(p^{\alpha+\beta})\)与左侧一致。由于两边都对 \(a,b\) 具有积性，因此结论对一般整数成立。

从恒等式出发：\(\displaystyle{\varphi(ni)=\varphi(i)\sum_{d\mid(n,i)}\varphi\left(\frac{n}{d}\right)}\)，于是\(\displaystyle{S(n,m)=\sum_{i=1}^{m}\varphi(i)\sum_{d\mid(n,i)}\varphi\left(\frac{n}{d}\right)}\)。

交换求和顺序，对每个 \(d\mid n\)，条件 \(d\mid(n,i)\) 等价于 \(d\mid i\)，得到：\(\displaystyle{S(n,m)=\sum_{d\mid n}\varphi\left(\frac{n}{d}\right)\sum_{\substack{1\le i\le m\\d\mid i}}\varphi(i)}\)，令 \(i=d\cdot t\)，内层和变为：\(\displaystyle{\sum_{t=1}^{\lfloor m/d\rfloor}\varphi(dt)=S(d,\lfloor m/d\rfloor)}\)。因此得到一个非常干净的递推式：

\[ S(n,m)=\sum_{d\mid n}\varphi\left(\frac{n}{d}\right)S\left(d,\left\lfloor\frac{m}{d}\right\rfloor\right) \]

本题\(N=510510=2\cdot 3\cdot 5\cdot 7\cdot 11\cdot 13\cdot 17\)是无平方因子（每个素因子指数都是 1）。对无平方因子更适合使用下面的等价递归（更容易展开成枚举）。

取 \(p\mid n\)，写 \(n=p\cdot n'\)。将求和按 \(i\) 是否被 \(p\) 整除分两类。

欧拉函数在乘一个素数时满足：\(\varphi(p x)=\begin{cases}(p-1)\varphi(x), & p\nmid x\\p\varphi(x), & p\mid x\end{cases}\)。于是：

  • 若 \(p\nmid i\)，则 \(\varphi(ni)=\varphi(p n' i)=(p-1)\varphi(n'i)\)
  • 若 \(p\mid i\)，写 \(i=p t\)，则\(\varphi(ni)=\varphi(p n' \cdot p t)=p\varphi(n'\cdot p t)=\varphi(n\cdot t)\)

因此\(\displaystyle{S(n,m)=\sum_{p\nmid i}(p-1)\varphi(n'i)+\sum_{t\le \lfloor m/p\rfloor}\varphi(nt)}\)，即

\[ S(n,m)=(p-1)S(n',m)+S\left(n,\left\lfloor\frac{m}{p}\right\rfloor\right) \]

这条式子非常关键：右边第二项会把 \(m\) 除以 \(p\)，不断递归就会把参数迅速降下去。

对 \(n=\prod_{j=1}^{7}p_j\)（\(p_j\)互不相同），反复应用上面的递归，每次遇到“第二项”就相当于把当前的 \(m\) 多除一个 \(p_j\)，最终会产生形如\(q=\prod_{j=1}^{7}p_j^{e_j}\le m\)的所有数（它们的素因子只来自 \({2,3,5,7,11,13,17}\)）。递归完全展开后，会把问题化成只需要 \(S(1,x)\)：\(\displaystyle{S(1,x)=\sum_{k=1}^{x}\varphi(k)}\)

记\(\displaystyle{F(x)=\sum_{k=1}^{x}\varphi(k)}\)并注意每个素数 \(p_j\) 在展开时都会贡献固定因子 \((p_j-1)\)，因此总系数是：\(\varphi(n)\)。

最终得到：

\[ S(n,m)=\varphi(n)\sum_{q} F\left(\left\lfloor\frac{m}{q}\right\rfloor\right) \]

其中\(n\)是无平方因子数，\(q\) 枚举所有满足 \(q\le m\) 且素因子只来自 \({2,3,5,7,11,13,17}\) 的数。

接下来计算\(F(n)\)，由题目415的结论，可知

\[ F(n)=\frac{n(n+1)}{2}-\sum_{t=2}^{n}F\left(\left\lfloor\frac{n}{t}\right\rfloor\right) \]

这就是可递归计算 \(F(n)\) 的公式，由此我们可以计算数论分块完成。

代码

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

typedef long long ll;
const int MOD = 1000000000;
int K = 510510;
const ll M = 100000000000LL;
const int PR[7] = {2, 3, 5, 7, 11, 13, 17};
const int LIMIT = 10000000;

vector<int> pref;
unordered_map<ll, int> memo;

void build_prefix_phi()
{
    vector<int> phi(LIMIT + 1);
    vector<int> primes;
    vector<int> is_comp(LIMIT + 1, 0);
    phi[1] = 1;
    for (int i = 2; i <= LIMIT; i++)
    {
        if (!is_comp[i])
        {
            primes.push_back(i);
            phi[i] = i - 1;
        }
        for (int p : primes)
        {
            ll v = 1LL * i * p;
            if (v > LIMIT)
                break;
            is_comp[v] = 1;
            if (i % p == 0)
            {
                phi[v] = phi[i] * p;
                break;
            }
            else
            {
                phi[v] = phi[i] * (p - 1);
            }
        }
    }
    pref.assign(LIMIT + 1, 0);
    for (int i = 1; i <= LIMIT; i++)
    {
        pref[i] = (pref[i - 1] + phi[i]) % MOD;
    }
}

int phi_sum(ll n)
{
    if (n <= LIMIT)
        return pref[n];
    auto it = memo.find(n);
    if (it != memo.end())
        return it->second;

    ll a = n % MOD;
    ll b = (n + 1) % MOD;
    if ((n & 1) == 0)
        a = (n / 2) % MOD;
    else
        b = ((n + 1) / 2) % MOD;
    ll res = (a * b) % MOD;

    ll l = 2;
    while (l <= n)
    {
        ll v = n / l;
        ll r = n / v;
        ll cnt = (r - l + 1) % MOD;
        res -= cnt * phi_sum(v) % MOD;
        res %= MOD;
        l = r + 1;
    }
    int ans = (res + MOD) % MOD;
    memo[n] = ans;
    return ans;
}

ll acc = 0;

void dfs(int i, ll cur)
{

    if (i == 7)
    {
        acc += phi_sum(M / cur);
        acc %= MOD;
        return;
    }
    int p = PR[i];
    ll x = cur;
    while (x <= M)
    {
        dfs(i + 1, x);
        if (x > M / p)
            break;
        x *= p;
    }
}

int main()
{
    build_prefix_phi();
    dfs(0, 1);
    ll ans = acc * phi(K) % MOD;
    cout << setw(9) << setfill('0') << ans << "\n";
    return 0;
}