Vibe Coding坠机了.

User

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
元组计数

记π(x)为素数计数函数，即小于等于x的素数的数目。
例如，π(1)=0，π(2)=1，π(100)=25。

记T(n,k)为满足以下条件的k元组(x1,…,xk)的数目：

1.任意xi都是正整数；
2.\sum_{i=1}^k \pi(x_i)=n
举例来说，T(3,3)=19。
这19个三元组是(1,1,5)、(1,5,1)、(5,1,1)、(1,1,6)、(1,6,1)、(6,1,1)、(1,2,3)、(1,3,2)、(2,1,3)、(2,3,1)、(3,1,2)、(3,2,1)、(1,2,4)、(1,4,2)、(2,1,4)、(2,4,1)、(4,1,2)、(4,2,1)、(2,2,2)。

已知T(10,10) = 869 985和T(103,103) ≡ 578 270 566 (mod 1 004 535 809)。

求T(20 000, 20 000) mod 1 004 535 809。

我的思路是:背包dp.
设dp[j][i]是当前正在枚举元组第j个元素,加上j这个元素总共和是i的时候的情况总数.

因为pi[i]最大也就到20000,我们可以把pi[i]=j对于固定的j把所有i搞成一个数组,假设为prj[j],那么代价就是
dp[j][i]+=dp[j-1][i-k]*prj[k]
直接dp就可以了,记得取模和使用滚动数组优化.

GPT

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
结果已算出：

**T(20,000, 20,000) ≡ 779,429,131 (mod 1,004,535,809).**

我用的做法（并在后台代码中实现）概述如下：

1. 注意到每个 (x_i) 的贡献是 (\pi(x_i))，而且对于固定的 (t\ge0) 满足 (\pi(x)=t) 的正整数个数是有限的：

   * (t=0) 时只有 (x=1)，计数为 1；
   * 对于 (t\ge1)，满足 (\pi(x)=t) 的 (x) 为区间 ([p_t, p_{t+1}-1])，数量为 (p_{t+1}-p_t)（其中 (p_t) 是第 (t) 个素数）。
     因此设 (F(z)=\sum_{t\ge0} \mathrm{prj}[t] z^t)，其中 (\mathrm{prj}[0]=1,\ \mathrm{prj}[t]=p_{t+1}-p_t)，则
     [
     T(n,k)=[z^n];F(z)^k.
     ]

2. 只需要 (\mathrm{prj}[t]) 对于 (t\le n)（更高次项不会贡献到 (z^n)），所以只生成前 (n+1) 个素数来得到这些差分。

3. 将多项式 (F(z)) 截断到次数 (n)，用二分幂对多项式做 (k) 次幂；每次多项式乘法用模 (1004535809) 下的 NTT（该模是常见的 NTT 素数），并在每步把结果截断到最多次数 (n)，从而保证效率。

4. 程序在本地（我在会话环境中）运行并计算出系数，最终得到上面的模值：**779,429,131**。

如果你想要我把用到的完整代码（我运行的 Python/NTT 实现）贴出来供你本地复现或在其他语言（如 C++）中移植，我可以直接把代码给你并解释关键点和复杂度（时间大约为 (O(n\log n\log k)) 的多项式乘法代价，内含常数来自 NTT 的多次 FFT）。要不要我把完整代码粘给你？
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
# Computing T(20000, 20000) mod 1004535809 using NTT-based polynomial exponentiation.
# We'll:
# 1. Generate primes until we have p_1..p_{20001}.
# 2. Build polynomial F(z) = sum_{t=0}^n prj[t] z^t where prj[0]=1 and for t>=1 prj[t]=p_{t+1}-p_t.
# 3. Compute F(z)^k truncated to degree n using binary exponentiation with NTT convolution.
# 4. Print coefficient of z^n mod MOD.
MOD = 1004535809
n = 20000
k = 20000

# 1) sieve to get first n+1 primes
def gen_primes_needed(count, limit=400000):
    sieve = bytearray(b'\x01') * (limit+1)
    sieve[0:2] = b'\x00\x00'
    primes = []
    for i in range(2, limit+1):
        if sieve[i]:
            primes.append(i)
            if len(primes) >= count:
                break
            step = i
            start = i*i
            if start > limit:
                continue
            sieve[start:limit+1:step] = b'\x00' * (((limit - start)//step) + 1)
    if len(primes) < count:
        raise RuntimeError(f"Not enough primes found (found {len(primes)}), increase limit")
    return primes

primes = gen_primes_needed(n+1, limit=400000)  # get p_1..p_{n+1}
# pr list: primes[0] = 2 = p1, primes[n] = p_{n+1}
# 2) build prj
prj = [0]*(n+1)
prj[0] = 1  # only x=1 has pi(x)=0
for t in range(1, n+1):
    prj[t] = primes[t] - primes[t-1]  # p_{t+1} - p_t  (primes indexed from 0)

# 3) NTT implementation
def mod_pow(a, e, mod=MOD):
    r = 1
    while e:
        if e & 1:
            r = (r * a) % mod
        a = (a * a) % mod
        e >>= 1
    return r

# primitive root for MOD (1004535809) is 3
PRIMITIVE_ROOT = 3

def ntt(a, invert=False):
    n = len(a)
    j = 0
    for i in range(1, n):
        bit = n >> 1
        while j & bit:
            j ^= bit
            bit >>= 1
        j ^= bit
        if i < j:
            a[i], a[j] = a[j], a[i]
    length = 2
    while length <= n:
        wlen = mod_pow(PRIMITIVE_ROOT, (MOD - 1) // length)
        if invert:
            wlen = mod_pow(wlen, MOD - 2)
        i = 0
        while i < n:
            w = 1
            for j in range(i, i + length // 2):
                u = a[j]
                v = a[j + length // 2] * w % MOD
                a[j] = (u + v) % MOD
                a[j + length // 2] = (u - v) % MOD
                w = w * wlen % MOD
            i += length
        length <<= 1
    if invert:
        inv_n = mod_pow(n, MOD - 2)
        for i in range(n):
            a[i] = a[i] * inv_n % MOD

def convolution(a, b, max_deg=None):
    # returns c = a * b truncated to degree max_deg (inclusive)
    if max_deg is None:
        max_deg = len(a) + len(b) - 2
    need = min(max_deg + 1, len(a) + len(b) - 1)
    size = 1
    while size < len(a) + len(b) - 1:
        size <<= 1
    fa = a[:] + [0]*(size - len(a))
    fb = b[:] + [0]*(size - len(b))
    ntt(fa)
    ntt(fb)
    for i in range(size):
        fa[i] = fa[i] * fb[i] % MOD
    ntt(fa, invert=True)
    return fa[:need]

# 4) exponentiate polynomial F(z) to power k with truncation to degree n
F = prj[:]  # length n+1
def poly_pow(poly, exponent, max_deg):
    # binary exponentiation with truncation
    result = [1]  # polynomial "1"
    base = poly[:]
    e = exponent
    while e:
        if e & 1:
            result = convolution(result, base, max_deg)
        e >>= 1
        if e:
            base = convolution(base, base, max_deg)
    # ensure length at least max_deg+1
    if len(result) < max_deg+1:
        result += [0]*((max_deg+1)-len(result))
    return result[:max_deg+1]

res_poly = poly_pow(F, k, n)
answer = res_poly[n] % MOD
print("T(20000,20000) mod", MOD, "=", answer)

fun fact

Flu给出了正确的递推式子,但是没有看出来这是个卷积式子,更没有看出来那是个NTT模数.
Vibe Coding被摆了一道.
Flu不应该直接丢一个题上去问AI怎么写,所以这个题做法有点假…

有可能用户给出的东西不用GPT自己算,会让AI运行的噪音更小,更大概率得到正确答案,所以提示词应尽可能详尽的描述怎么写,最好把伪代码直接给出来让AI实现.