Project Euler 468

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

Project Euler 468

题目

Smooth divisors of binomial coefficients

An integer is called B-smooth if none of its prime factors is greater than \(B\).

Let \(S_B(n)\) be the largest \(B\)-smooth divisor of \(n\).

Examples:

\(\begin{aligned} &S_1(10)=1\\ &S_4(2100) = 12\\ &S_{17}(2496144) = 5712 \end{aligned}\)

Define \(\displaystyle F(n)=\sum_{B=1}^n \sum_{r=0}^n S_B\left(\binom n r\right)\). Here, \(\displaystyle \binom n r\) denotes the binomial coefficient.

Examples:

\(\begin{aligned} &F(11) = 3132\\ &F(1111) \bmod 1\,000\,000\,993 = 706036312\\ &F(111\,111) \bmod 1\,000\,000\,993 = 22156169 \end{aligned}\)

Find \(F(11\,111\,111) \bmod 1\,000\,000\,993\).

解决方案

令\(\displaystyle{T_n(k)=\sum_{B=1}^n S_B(k),}\)则交换求和次序立即得到\(\displaystyle{F(n)=\sum_{r=0}^n T_n\left(\binom nr\right).}\)因此核心变成：对每个二项式系数 \(k=\dbinom nr\)，高效计算 \(T_n(k)\)。

对任意整数 \(k\)，有标准分解\(\displaystyle{k=\prod_{p}p^{v_p(k)}.}\)注意到\(\displaystyle{S_B(k)=\prod_{p\le B}p^{v_p(k)}.}\)因此当 \(B\) 在相邻质数之间变化时，\(S_B(k)\) 不会变。

令 \(n\) 内的质数为\(2=p_1<p_2<\cdots<p_t\le n,\)并补上\(p_0=1, p_{t+1}=n+1.\)定义前缀乘积\(\displaystyle{P_i(k)=\prod_{j=1}^{i}p_j^{v_{p_j}(k)}, P_0(k)=1.}\)

那么当 \(B\in[p_i,p_{i+1}-1]\) 时，\(S_B(k)=P_i(k).\)于是\(\displaystyle{T_n(k)=\sum_{B=1}^n S_B(k)=\sum_{i=0}^{t}(p_{i+1}-p_i)P_i(k).}\)设权重\(w_i=p_{i+1}-p_i,\)则上式是一个“权重 \(\times\) 前缀积”的形式：\(\displaystyle{T_n(k)=\sum_{i=0}^{t}w_iP_i(k)}\)

对 \(k=1\)，显然 \(P_i(1)=1\)，故\(\displaystyle{T_n(1)=\sum_{i=0}^t w_i=n.}\)。

利用经典递推：\(\dbinom n{r+1}=\dbinom nr\cdot \dfrac{n-r}{r+1}.\)从 \(r=0\) 的 \(k=\dbinom n0=1\) 出发，每一步用分子 \(n-r\) 乘上、分母 \(r+1\) 除掉。

若当前 \(k\) 被乘上某个质数幂 \(p^e\)，那么对所有 \(B\ge p\) 的 smooth 因子都会多出 \(p^e\)，即对所有 \(i\) 满足 \(p_i\ge p\) 的前缀积都要乘上 \(p^e\)，等价于进行更行：\(P_i(k)\leftarrow P_i(k)\cdot p^e, \forall i\ge I(p).\)其中\(I(p)\)表示质数的序号（\(I(2)=1,I(3)=2\)）。

因此对表达式\(\displaystyle{T_n(k)=\sum_{i=0}^{t}w_i,P_i(k)}\) 来说，乘上 \(p^e\) 就是对所有 \(i\ge I(p)\) 的项整体乘 \(p^e\)，即一次 后缀乘法更新。同理，除以 \(p^e\) 就是乘上 \((p^e)^{-1}\bmod M\)。

直接维护 \(P_i(k)\) 并做后缀乘法会很麻烦。关键是引入“乘法差分”：定义一列因子 \(g_i\) 满足\(\displaystyle{P_i=\prod_{j=0}^{i}g_j,g_0=1.}\)

当我们要让所有 \(P_i\)（\(i\ge s\)）统一乘上一个因子 \(a\) 时，只需令\(g_s\leftarrow g_s\cdot a,\)因为对 \(i\ge s\) 的前缀积都会包含 \(g_s\)，而对 \(i<s\) 的前缀积不包含它。

因此可以我们可以考虑用线段树一些信息，单点更新时沿路径回溯更新即可。

于是\(\displaystyle{T_n(k)=\sum_{i=0}^{t}w_i\prod_{j=0}^{i}g_j.}\)现在每次乘/除一个质数幂，只对应 对某个位置 \(g_s\) 的单点乘法更新。

我们需要支持两种操作：

  1. 单点乘法：\(g_s\leftarrow g_s\cdot a\)
  2. 查询：\(\displaystyle{\sum_{i=0}^{t}w_i\prod_{j=0}^{i}g_j.}\)

对任意区间 \([L,R]\) 定义结点信息为一对：\((p,s)\)，其中\(\displaystyle{p=\prod_{i=L}^{R} g_i,s=\sum_{i=L}^{R} w_i\prod_{j=L}^{i} g_j.}\)若区间拆为左右子区间 \([L,M]\) 与 \([M+1,R]\)，则：

  • 乘积显然：\(p=p_L\cdot p_R.\)
  • 对加权和，右边的前缀积要先乘上左区间整体乘积：\(s=s_L+p_L\cdot s_R.\)

可见，根结点的 \(s\) 就是所求的 \(T_n(k)\)。

由于\(\dbinom nr=\dbinom n{n-r},\)从而\(T_n\left(\binom nr\right)=T_n\left(\binom n{n-r}\right).\)

令 \(h=\left\lfloor \dfrac n2\right\rfloor\)，则

  • 若 \(n\) 为奇数：\(\displaystyle{F(n)=2\sum_{r=0}^{h}T_n\left(\binom nr\right).}\)
  • 若 \(n\) 为偶数：\(\displaystyle{F(n)=2\sum_{r=0}^{h}T_n\left(\binom nr\right)-T_n\left(\binom nh\right)}.\)

所以只需遍历 \(r=0\sim h\)。

最终维护好\(T_n(\cdot)\)的值并相加即可。

代码

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

typedef long long ll;
const int N = 11111111;
const int MOD = 1000000993;

struct SegTree
{
    int n;
    vector<int> prod;
    vector<int> sum;

    SegTree(const vector<int> &w)
    {
        n = w.size();
        prod.assign(4 * n + 5, 1);
        sum.assign(4 * n + 5, 0);
        build(1, 0, n, w);
    }
    inline void pull(int u)
    {
        int lc = u << 1;
        int rc = u << 1 | 1;
        prod[u] = (ll)prod[lc] * prod[rc] % MOD;
        sum[u] = (sum[lc] + (ll)prod[lc] * sum[rc]) % MOD;
    }

    void build(int u, int l, int r, const vector<int> &w)
    {
        if (l + 1 == r)
        {
            prod[u] = 1;
            sum[u] = w[l] % MOD;
            return;
        }
        int m = (l + r) >> 1;
        build(u << 1, l, m, w);
        build(u << 1 | 1, m, r, w);
        pull(u);
    }

    void update(int u, int l, int r, int pos, int factor)
    {
        if (l + 1 == r)
        {
            prod[u] = (ll)prod[u] * factor % MOD;
            sum[u] = (ll)sum[u] * factor % MOD;
            return;
        }
        int m = (l + r) >> 1;
        if (pos < m)
            update(u << 1, l, m, pos, factor);
        else
            update(u << 1 | 1, m, r, pos, factor);
        pull(u);
    }

    int total() const
    {
        return sum[1];
    }
};

int solve(int n)
{
    Factorizer factorizer(n);
    vector<int> &spf = factorizer.spf, &primes = factorizer.primes;

    int t = (int)primes.size();
    int L = t + 1;

    vector<int> primeIndex(n + 1, 0);
    for (int i = 0; i < t; i++)
        primeIndex[primes[i]] = i + 1;

    vector<int> invIdx(L, 1);
    for (int i = 1; i <= t; i++)
        invIdx[i] = qpow(primes[i - 1], MOD - 2, MOD);

    vector<int> w(L, 0);
    int prev = 1;
    for (int i = 0; i < t; i++)
    {
        w[i] = primes[i] - prev;
        prev = primes[i];
    }
    w[t] = (n + 1) - prev;

    SegTree st(w);

    int half = n >> 1;
    ll total = st.total();

    for (int r = 0; r < half; r++)
    {
        int num = n - r;
        int den = r + 1;
        int g = __gcd(num, den);
        num /= g;
        den /= g;

        int x = num;
        while (x > 1)
        {
            int p = spf[x];
            int c = 0;
            while (x % p == 0)
            {
                x /= p;
                c++;
            }
            int mul = 1;
            for (int i = 0; i < c; i++)
                mul = (ll)mul * p % MOD;
            int idx = primeIndex[p];
            if (mul != 1)
                st.update(1, 0, L, idx, mul);
        }

        x = den;
        while (x > 1)
        {
            int p = spf[x];
            int c = 0;
            while (x % p == 0)
            {
                x /= p;
                c++;
            }
            int ip = invIdx[primeIndex[p]];
            int mul = 1;
            for (int i = 0; i < c; i++)
                mul = (ll)mul * ip % MOD;
            int idx = primeIndex[p];
            if (mul != 1)
                st.update(1, 0, L, idx, mul);
        }

        total += st.total();
        total %= MOD;
    }

    int mid = st.total();
    int ans = (2ll * total) % MOD;
    if (!(n & 1))
        ans = (ans + MOD - mid) % MOD;
    return ans;
}

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