diff --git a/solutions/0053.Combinatoric/README.md b/solutions/0053.Combinatoric/README.md new file mode 100644 index 0000000..383db19 --- /dev/null +++ b/solutions/0053.Combinatoric/README.md @@ -0,0 +1,115 @@ +# 组合数的思考 + +我最惊讶的地方在于,直接计算比重新整理计算公式更快。这是不是因为搜索空间不大,所以才如此? + +那么另一个角度来说,是不是基于数学原理来直接确认数量是更快的方式?显然,是这样的。 + +简单来说,组合数 $\displaystyle\binom{n}{r}$ 相对于 $r$ 具有单峰对称性, +所以只有找到最小的 $r^*$ 就可以快速得到具体数量了,同时根据单峰对称性,计算量也能有大幅减少, +至少不用计算大数的连乘。 + +## 核心数学观察 + +对于固定的 $n$,组合数 $\displaystyle\binom{n}{r}$ 关于 $r$ 具有**单峰对称性**: + +$$\binom{n}{0} < \binom{n}{1} < \cdots < \binom{n}{\lfloor n/2\rfloor} = \binom{n}{\lceil n/2\rceil} > \cdots > \binom{n}{n}$$ + +因此,只要找到**最小的** $r^*$ 使得 $\displaystyle\binom{n}{r^*} > 10^6$,则该 $n$ 下所有满足条件的 $r$ 构成连续区间: + +$$r^* \le r \le n - r^*$$ + +满足条件的个数为: + +$$\boxed{n - 2r^* + 1}$$ + +这已经把问题从"枚举所有 $r$"转化为"对每个 $n$ 找临界点 $r^*$"。 + +--- + +## 方法一:递推截断法(最实用) + +利用递推关系,从 $r=0$ 开始逐项计算,一旦超过 $10^6$ 立即停止: + +$$\binom{n}{r} = \binom{n}{r-1} \cdot \frac{n-r+1}{r}$$ + +**优点**: +- 不需要计算阶乘 +- 不需要处理大整数(因为中途就停了) +- 不需要计算对称右侧的所有值 + +例如 $n=100$ 时,只需算到 $r=4$ 就发现 $\binom{100}{4}=3{,}921{,}225 > 10^6$,立刻得知该 $n$ 下有 $100-2\times4+1=93$ 个满足条件的 $r$。 + +--- + +## 方法二:固定 $r$,反向求阈值 $N_r$ + +这是一个更"数学化"的思路。对于每个固定的 $r$,由于 $\displaystyle\binom{n}{r}$ 关于 $n$ 严格递增,存在唯一的最小整数 $N_r$ 使得: + +$$\binom{N_r}{r} > 10^6 \quad\text{且}\quad \binom{N_r-1}{r} \le 10^6$$ + +计算得到(只需对 $r=1$ 到 $50$ 各跑一个递增序列): + +| $r$ | $N_r$(最小 $n$) | +|:---:|:---:| +| 4 | 72 | +| 5 | 44 | +| 6 | 33 | +| 7 | 28 | +| 8 | 25 | +| 9 | 24 | +| 10 | 23 | +| 11 | 23 | +| 12 | 23 | +| 13 | 23 | +| 14 | 24 | +| ... | ... | + +**关键推论**:对于给定的 $n$,临界点恰好是 + +$$r^*(n) = \min\{r \mid N_r \le n\}$$ + +即:**找到第一个阈值不超过 $n$ 的 $r$**。 + +例如: +- $n=25$:$N_8=25 \le 25$,而 $N_7=28 > 25$,故 $r^*=8$,个数为 $25-16+1=10$ +- $n=72$:$N_4=72 \le 72$,而 $N_3=183 > 72$,故 $r^*=4$,个数为 $72-8+1=65$ + +然后对 $n=23$ 到 $100$ 求和 $\sum(n-2r^*+1)$ 即得 **4075**。 + +这种方法把"对每个 $(n,r)$ 算组合数"变成了"对每个 $r$ 找一次阈值",计算量从 $O(n^2)$ 降到约 $O(n \cdot r_{\max})$,且每个序列找到阈值即停。 + +--- + +## 方法三:对数法(避免所有大整数) + +预先计算 $\ln(k)$ 的前缀和,则: + +$$\ln\binom{n}{r} = \sum_{i=1}^{n}\ln i - \sum_{i=1}^{r}\ln i - \sum_{i=1}^{n-r}\ln i$$ + +只需判断该值是否大于 $\ln(10^6) \approx 13.8155$。全程只涉及 double 类型的加减,**完全不出现大整数**。 + +--- + +## 方法四:多项式求根(针对小 $r$) + +对于固定的 $r \le 4$,$\displaystyle\binom{n}{r}$ 是关于 $n$ 的低次多项式,可以直接解方程: + +- $r=1$: $n = 10^6$(超出范围) +- $r=2$: $n(n-1) = 2\times 10^6$,解得 $n \approx 1414$(超出范围) +- $r=3$: $n(n-1)(n-2) = 6\times 10^6$,解得 $n \approx 182$(超出范围) +- $r=4$: $n(n-1)(n-2)(n-3) = 24\times 10^6$,正实根约为 $71.5$,故 $N_4=72$ + +对于 $r \ge 5$,虽然是高次方程无解析解,但可以用牛顿迭代等数值方法快速求根,从而确定 $N_r$。 + +--- + +## 总结 + +| 方法 | 是否避免大整数 | 是否减少计算次数 | 核心数学工具 | +|:---:|:---:|:---:|:---:| +| 递推截断 | ✅ | ✅ | 单调性 + 递推公式 | +| 阈值法 ($N_r$) | ✅ | ✅ | 单调性 + 固定 $r$ 求根 | +| 对数法 | ✅ | ⚠️(仍需逐对判断) | $\ln$ 前缀和 | +| 多项式求根 | ✅ | ✅ | 代数方程 / 数值分析 | + +**不存在一个封闭公式能直接写出答案 4075 而完全不做任何数值运算**——因为问题的本质是离散的不等式判断。但利用组合数的单调性和对称性,可以把计算量压缩到极小,且全程避开巨大整数的运算。 diff --git a/solutions/0053.Combinatoric/euler_53_math.py b/solutions/0053.Combinatoric/euler_53_math.py new file mode 100644 index 0000000..bb8b7fc --- /dev/null +++ b/solutions/0053.Combinatoric/euler_53_math.py @@ -0,0 +1,43 @@ +import time + + +def timer(func): + def wrapper(*args, **kwargs): + start_time = time.time() + result = func(*args, **kwargs) + end_time = time.time() + elapsed_time = end_time - start_time + print(f"{func.__name__} time: {elapsed_time:.6f} seconds") + return result + + return wrapper + + +@timer +def count_combinations_above_million(limit_n=100, threshold=1_000_000): + """ + 对于 1 ≤ n ≤ limit_n,统计组合数 C(n, r) 大于 threshold 的个数。 + 利用递推 C(n, r) = C(n, r-1) * (n - r + 1) / r 以及对称性 C(n, r) = C(n, n-r)。 + """ + count = 0 + + for n in range(1, limit_n + 1): + comb = 1 # C(n, 0) + # 只需检查 r 从 1 到 n//2 + for r in range(1, n // 2 + 1): + # 递推更新组合数 + comb = comb * (n - r + 1) // r + + if comb > threshold: + # 一旦找到最小的 r 使得 C(n, r) > threshold, + # 则 r 到 n-r 之间的所有值都满足条件。 + # 符合条件的个数 = (n - r) - r + 1 = n - 2*r + 1 + count += n - 2 * r + 1 + break # 该 n 的剩余情况不必再检查 + + return count + + +if __name__ == "__main__": + result = count_combinations_above_million() + print(result)