Project Euler 447

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

Project Euler 447

题目

Retractions C

For every integer \(n>1\), the family of functions \(f_{n,a,b}\) is defined by

\(f_{n,a,b}(x)\equiv a x + b \bmod n\) for \(a,b,x\) integer and \(0< a <n, 0 \le b < n,0 \le x < n\).

We will call \(f_{n,a,b}\) a retraction if \(f_{n,a,b}(f_{n,a,b}(x)) \equiv f_{n,a,b}(x) \bmod n\) for every \(0 \le x < n\).

Let \(R(n)\) be the number of retractions for \(n\).

\(\displaystyle F(N)=\sum_{n=2}^N R(n)\).

\(F(10^7)\equiv 638042271 \bmod 1\,000\,000\,007\).

Find \(F(10^{14})\).

Give your answer modulo \(1\,000\,000\,007\).

解决方案

第445题直接给出了如下关于\(R(n)\)的式子，将\(n\)写成质因数分解的形式\(n=\prod_{i=1}^k p_i^{e_i}\)。那么有：

\[R(n)=\prod_{i=1}^k(p_i^{e_i}+1)-n\]

定义\(\displaystyle{U(n)=\prod_{i=1}^k(1+p_i^{e_i}).}\)则\(R(n)=U(n)-n.\)因此

\[ F(N)=\sum_{n=2}^N (U(n)-n)=\left(\sum_{n=1}^N U(n)\right)-\frac{N(N+1)}{2}, \]

因为 \(U(1)=1\) 且 \(\displaystyle{\sum_{n=2}^N n=\dfrac{N(N+1)}{2}-1}\)，两边的 \(-1\) 抵消。于是核心变成计算\(U(n)\)。

注意 \(U(n)\) 也是单位因子的和：\(\displaystyle{U(n)=\sum_{\substack{d\mid n\\\gcd(d,n/d)=1}} d.}\)

用指示函数写为\(\displaystyle{U(n)=\sum_{d\mid n} d\cdot [\gcd(d,n/d)=1].}\)利用经典恒等式 \(\displaystyle{\mathbb{1}\{\gcd(x,y)=1\}=\sum_{k\mid \gcd(x,y)} \mu(k),}\) 其中 \(\mu\) 为 Möbius 函数。令 \(x=d,y=n/d\)，构造Möbius反演，得到\(\displaystyle{U(n)=\sum_{d\mid n} d \sum_{k\mid \gcd(d,n/d)} \mu(k).}\)

把满足 \(k\mid d\) 且 \(k\mid n/d\) 的 \(k\) 提取出来：设 \(d=k t\)，则 \(n/d = n/(kt)\) 也需被 \(k\) 整除，等价于 \(k^2\mid n\)。写 \(n=k^2 m\)，则 \(t\mid m\)。代入上式：

\[ U(n)=\sum_{k^2\mid n}\mu(k)\sum_{t\mid n/k^2} (kt) =\sum_{k^2\mid n}\mu(k)k\sum_{t\mid n/k^2} t =\sum_{k^2\mid n}\mu(k)k\sigma\left(\frac{n}{k^2}\right), \]

其中 \(\displaystyle{\sigma(m)=\sum_{d\mid m} d}\) 是约数和函数。

令 \[ S(X)=\sum_{m=1}^X\sigma(m). \] 则 \[ \begin{aligned} U(N)&=\sum_{n\le N}\sum_{k^2\mid n}\mu(k)k\sigma(n/k^2) \\ &=\sum_{k\le \sqrt N}\mu(k)k\sum_{m\le N/k^2}\sigma(m) \\ &=\sum_{k\le \sqrt N}\mu(k)kS\left(\left\lfloor\frac{N}{k^2}\right\rfloor\right). \end{aligned} \]

因此，\(S(X)\)可以通过数论分块求解。由\(\displaystyle{\sigma(m)=\sum_{d\mid m} d,}\)可交换求和得到\(\displaystyle{S(X)=\sum_{m\le X}\sum_{d\mid m} d=\sum_{d\le X} d\left\lfloor\frac{X}{d}\right\rfloor.}\)

用整除分块：当 \(v=\left\lfloor X/i\right\rfloor\) 固定时，\(i\) 落在区间\(i\in\left[L,R\right],L=i,R=\left\lfloor\dfrac{X}{v}\right\rfloor.\)

区间内 \(\left\lfloor X/t\right\rfloor=v\) 恒定，因此\(\displaystyle{\sum_{t=L}^{R} t\left\lfloor\frac{X}{t}\right\rfloor= v\sum_{t=L}^{R} t= v\cdot \frac{(L+R)(R-L+1)}{2}.}\)

为了避免计算\(\mu\)，我们构建一个递推式。设\(q=\left\lfloor\sqrt N\right\rfloor,A(d)=S\left(\left\lfloor\dfrac{N}{d^2}\right\rfloor\right).\)

定义新函数

\[g_q(d)=\sum_{k=1}^{\lfloor q/d\rfloor} \mu(k)kA(kd).\]

显然目标就是\(U(N)=g_q(1).\)接下来证明一个关键恒等式：

\[ \begin{aligned} \sum_{k=1}^{\lfloor q/d\rfloor} kg_q(kd) &=\sum_{k} k \sum_{t} \mu(t)tA(tkd) \\ &=\sum_{m} A(md)\sum_{t\mid m} \mu(t)t\cdot \frac{m}{t} \\ &=\sum_{m} A(md)m \sum_{t\mid m}\mu(t). \end{aligned} \]

而可见，

\[ \sum_{t\mid m}\mu(t)= \begin{cases} 1,& m=1,\\ 0,& m>1, \end{cases} \]

所以只剩 \(m=1\) 项：\(\displaystyle{\sum_{k} kg_q(kd)=A(d).}\)因此得到了反演形式的递推：

\[ g_q(d)= \left \{\begin{aligned} &A(d) & & \text{if}\quad 2d>q\\ &A(d)-\sum_{k=2}^{\left\lfloor q/d\right\rfloor} k\cdot g_q(kd) & & \text{if}\quad 2d\le q \end{aligned}\right. \]

因此，按 \(d=q,q-1,\dots,1\) 递减计算。因为右侧用到的都是 \(kd>d\)，已在更大的 \(d\) 处算出。

最后计算 \[ F(N) = g(1)-\frac{N(N+1)}{2}. \]

代码

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

typedef long long ll;
const ll N = 100000000000000LL; // 1e14
const ll MOD = 1000000007LL;
const ll INV2 = (MOD + 1) / 2;

ll tri_sum(ll L, ll R)
{
    // sum_{t=L}^R t mod MOD
    // = (L+R)*(R-L+1)/2
    ll a = (L + R) % MOD;
    ll b = (R - L + 1) % MOD;
    ll prod = a * b % MOD * INV2 % MOD;
    return prod;
}

ll sum_sigma_prefix(ll X)
{
    // S_sigma(X) = sum_{d<=X} d * floor(X/d)
    ll res = 0;
    for (ll i = 1; i <= X;)
    {
        ll v = X / i;
        ll j = X / v; // max index with same quotient
        ll sum_i = tri_sum(i, j);
        ll add = (ll)(v % MOD) * sum_i % MOD;
        res += (ll)add;
        if (res >= MOD)
            res -= MOD;
        i = j + 1;
    }
    return res;
}

unordered_map<ll, int> cache;

ll getA(int d)
{
    ll dd = (ll)d * d;
    ll K = (ll)(N / (ll)dd);
    auto it = cache.find(K);
    if (it != cache.end())
        return it->second;
    ll val = sum_sigma_prefix(K);
    cache.emplace(K, (int)val);
    return val;
};

int main()
{
    ios::sync_with_stdio(false);
    cin.tie(nullptr);
    const int q = (int)llround(floor(sqrt((long double)N)));
    vector<int> g(q + 1, 0);

    // Cache A(d) = S_sigma(N / d^2) by the key K = floor(N / d^2).

    // DP / inversion recursion: g(d) = A(d) - sum_{k>=2} k*g(kd)
    for (int d = q; d >= 1; --d)
    {
        ll cur = getA(d);
        // iterate multiples m = k*d
        for (int m = d + d, k = 2; m <= q; m += d, ++k)
        {
            cur -= 1ll * k * g[m] % MOD;
            if (cur < 0)
                cur += MOD;
        }
        g[d] = cur;
    }
    ll triN = N % MOD * ((N + 1) % MOD) % MOD * INV2 % MOD;
    ll ans = 1ll * g[1] - triN;
    ans = (ans % MOD + MOD) % MOD;
    cout << ans << "\n";
    return 0;
}