user

1
2
3
4
5
6
7
8
9
10
11
巨大的二项式系数

二项式系数C(10^{18},10^9)是一个超过90亿（9×10^9）位的数。

记M(n,k,m)为二项式系数C(n,k)除以m的余数。

对于素数1000<p<q<r<5000，试计算∑M(10^{18},10^9,pqr)。

我的思路:首先筛出所有的p,然后用lucas定理求出这个组合数%p的值,然后用crt合并.

请给我c++代码,记得开long long.

gpt

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
下面我按**可实现、可通过时限**的思路给出**完整算法说明 + 可直接使用的 C++ 实现**。整体严格遵循你给出的技术路线（Lucas + CRT），并针对 (n=10^{18},k=10^9) 的极端规模做了必要的工程化处理。

---

## 一、问题结构拆解

要求计算

[
\sum_{\substack{1000 < p < q < r < 5000 \ p,q,r\ \text{prime}}}
M(10^{18},10^9,pqr)
]

其中

[
M(n,k,m) = \binom{n}{k} \bmod m
]

### 核心观察

* (p,q,r) **均为不同素数**
* (pqr) 互素 ⇒ 可用 **CRT**
* 只需计算
  [
  \binom{n}{k} \bmod p,\quad \binom{n}{k} \bmod q,\quad \binom{n}{k} \bmod r
  ]

---

## 二、Lucas 定理的可行性

由于：

* (k=10^9)
* 所有模数 (p < 5000)

可以用 **Lucas 定理**

[
\binom{n}{k} \bmod p
====================

\prod \binom{n_i}{k_i} \bmod p
]

其中 (n_i,k_i) 是 (p) 进制展开。

### 关键点

* (k < p^5) 对所有 (p>1000) 成立
* (n=10^{18}) 的 (p) 进制位数 (\le 6)
* 每次 Lucas 分解只需 **常数级小组合数**

---

## 三、CRT 合并（三模）

对于

[
x \equiv a \pmod p
]
[
x \equiv b \pmod q
]
[
x \equiv c \pmod r
]

使用标准 CRT（逐步合并）

---

## 四、整体算法流程

### Step 1：筛素数

* 筛出所有 (1000 < p < 5000) 的素数
* 数量约 **~500**

### Step 2：预处理组合数模 p

* 对每个素数 (p)

  * 预计算 `fac[0..p-1]` 与逆元
  * 以便 O(1) 计算 (\binom{a}{b} \bmod p)

### Step 3：Lucas(n,k,p)

* 将 (n,k) 转为 (p) 进制
* 按位计算

### Step 4：三重循环枚举 (p,q,r)

* 对每个三元组

  * 查表得到：
    [
    a_p = C(n,k) \bmod p
    ]
    [
    a_q = C(n,k) \bmod q
    ]
    [
    a_r = C(n,k) \bmod r
    ]
  * CRT 合并为 ( \bmod (pqr) )
  * 累加

---

## 五、完整 C++ 实现（`long long`）
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
#include <bits/stdc++.h>
using namespace std;
using ll = long long;

/* 快速幂 */
ll modpow(ll a, ll e, ll mod) {
    ll r = 1;
    while (e) {
        if (e & 1) r = r * a % mod;
        a = a * a % mod;
        e >>= 1;
    }
    return r;
}

/* 扩展欧几里得 */
ll exgcd(ll a, ll b, ll &x, ll &y) {
    if (!b) {
        x = 1; y = 0;
        return a;
    }
    ll x1, y1;
    ll g = exgcd(b, a % b, x1, y1);
    x = y1;
    y = x1 - a / b * y1;
    return g;
}

/* CRT: 合并 x ≡ a (mod m), x ≡ b (mod n) */
ll crt(ll a, ll m, ll b, ll n) {
    ll x, y;
    exgcd(m, n, x, y);
    ll mod = m * n;
    ll res = (a + (b - a + n) % n * x % n * m) % mod;
    if (res < 0) res += mod;
    return res;
}

/* Lucas 支持结构 */
struct Lucas {
    int p;
    vector<ll> fac, invfac;

    Lucas(int _p): p(_p), fac(p), invfac(p) {
        fac[0] = 1;
        for (int i = 1; i < p; i++) fac[i] = fac[i-1] * i % p;
        invfac[p-1] = modpow(fac[p-1], p-2, p);
        for (int i = p-2; i >= 0; i--)
            invfac[i] = invfac[i+1] * (i+1) % p;
    }

    ll Csmall(int n, int k) {
        if (k < 0 || k > n) return 0;
        return fac[n] * invfac[k] % p * invfac[n-k] % p;
    }

    ll C(ll n, ll k) {
        ll res = 1;
        while (n || k) {
            int ni = n % p;
            int ki = k % p;
            if (ki > ni) return 0;
            res = res * Csmall(ni, ki) % p;
            n /= p;
            k /= p;
        }
        return res;
    }
};

int main() {
    const ll N = 1000000000000000000LL;
    const ll K = 1000000000LL;

    /* 筛素数 */
    vector<int> primes;
    vector<bool> isprime(5000, true);
    isprime[0] = isprime[1] = false;
    for (int i = 2; i < 5000; i++) {
        if (isprime[i]) {
            if (i > 1000) primes.push_back(i);
            for (ll j = 1LL * i * i; j < 5000; j += i)
                isprime[j] = false;
        }
    }

    /* 预处理 Lucas */
    unordered_map<int, ll> Cmod;
    for (int p : primes) {
        Lucas luc(p);
        Cmod[p] = luc.C(N, K);
    }

    ll ans = 0;

    int sz = primes.size();
    for (int i = 0; i < sz; i++) {
        for (int j = i+1; j < sz; j++) {
            for (int k = j+1; k < sz; k++) {
                ll p = primes[i];
                ll q = primes[j];
                ll r = primes[k];

                ll a = Cmod[p];
                ll b = Cmod[q];
                ll c = Cmod[r];

                ll x = crt(a, p, b, q);
                x = crt(x, p*q, c, r);

                ans += x;
            }
        }
    }

    cout << ans << "\n";
    return 0;
}
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
---

## 六、复杂度评估

| 部分        | 复杂度                                    |
| --------- | -------------------------------------- |
| 筛素数       | (O(n\log\log n))                       |
| Lucas 预处理 | (O(\sum p))                            |
| 三重枚举      | (\binom{500}{3} \approx 2 \times 10^7) |
| CRT       | 常数                                     |

**整体在 C++ O2 下可跑**

---

## 七、可进一步优化（如需要）

* 使用 `int` 替代 `ll` 存模值
* 提前剪枝 `C(n,k) ≡ 0 mod p`
* OpenMP 并行三重循环

---

如果你希望：

* **进一步数学化简（减少三重循环）**
* **并行 / SIMD 优化**
* **改写为多模 Lucas 一次性处理**

可以继续深入。