通过欧拉计划学Rust编程：第100题

由于研究Libra等数字货币编程技术的需要，学习了一段时间的Rust编程，一不小心刷题上瘾。

• 刷完欧拉计划中的63道基础题，能学会Rust编程吗？

“欧拉计划”的网址：https://projecteuler.net

英文如果不过关，可以到中文翻译的网站：http://pe-cn.github.io/

这个网站提供了几百道由易到难的数学问题，你可以用任何办法去解决它，当然主要还得靠编程，编程语言不限，论坛里已经有Java、C#、Python、Lisp、Haskell等各种解法，当然如果你直接用google搜索答案就没任何乐趣了。

这次解答的是第100题：https://projecteuler.net/problem=100

题目描述：

在一个盒子中装有21个彩色碟子，其中15个是蓝的，6个是红的。如果随机地从盒子中取出两个碟子，取出两个蓝色碟子的概率是P(BB) = (15/21) × (14/20) = 1/2。

下一组使得取出两个蓝色盘子的概率恰好为50%的情况，有120个碟子，其中85个蓝色碟子和35个红色碟子。

当盒子中装有超过10^12 = 1,000,000,000,000个碟子时，找出第一组满足上述要求的安排，并求此时盒子中蓝色碟子的数量。

解题过程：

遇到一个复杂的问题，首先可以尝试先解决简单的情况，然后慢慢逼近最终的问题。

第一步：

假设碟子数量不超过1000个，计算满足条件的蓝色碟子的数量。

先运用一下数学知识，假设碟子总数为n，蓝碟数量为b，则：P(BB) = (b / n) * ((b-1) / (n-1)) = 0.5 即：n*(n-1) = 2 * b * (b-1)

代码的逻辑很简单：

for n in 2..1000 {
    for b in 2..n {
        if n * (n-1) == 2 * b * (b - 1) {
            println!("n:{} b:{}", n, b);
        }
    }
}

运行之后，可以找到几组答案：

n:4 b:3
n:21 b:15
n:120 b:85
n:697 b:493

第二步，尝试将n放大，寻找满足条件的解

将n的上限改到10万，注意用u64数据类型，即100_000_u64，程序仍可运行，快速找到几组解，但后面的求解速度会非常慢。

n:4 b:3
n:21 b:15
n:120 b:85
n:697 b:493
n:4060 b:2871
n:23661 b:16731

第三步，尝试优化算法性能，缩小b的寻找范围

这次利用不等式的性质，可以迅速缩小b的寻找范围，减小一层循环，只需浮点计算并取整即可。

let sq = 0.5_f64.sqrt();
for n in 2..100_000_000_u64 {
    let b = (sq * (n as f64) + (1.0 - sq)) as u64;
    if n * (n-1) == 2 * b * (b - 1) {
        println!("n:{} b:{}", n, b);
    }
}

现在可以轻松求出1亿之内的所有答案：

n:4 b:3
n:21 b:15
n:120 b:85
n:697 b:493
n:4060 b:2871
n:23661 b:16731
n:137904 b:97513
n:803761 b:568345
n:4684660 b:3312555
n:27304197 b:19306983

第四步，暴力求解

由于原题只需要求出1万亿以上的一组解，可以暴力求解。

let sq = 0.5_f64.sqrt();
for n in 1000_000_000_000_u64.. {
    let b = (sq * (n as f64) + (1.0 - sq)) as u64;
    if n * (n-1) == 2 * b * (b - 1) {
        println!("n:{} b:{}", n, b);
        break;
    }
}

程序运行了4分钟，最后求出了答案。

第五步，更深入地研究

在projecteuler网站上提交了正确答案之后，可以参与论坛的讨论，看看别人的解法和源代码。

如果你会使用Mathematica这款软件，可以用一行代码解决这道题：

Reduce[ 0 < b < 10^12 < n && 2 b(b-1) == n(n-1),{b,n},Integers ]

还有人发现下一组解n与上一组解(n-1)的比为5.8284的规律，可以大幅缩小n的搜索范围，30毫秒之内找到答案。

let sqrt_1_2 = 0.5_f64.sqrt();
let mut n = 21_u64;
loop {
    let b = (sqrt_1_2 * (n as f64) + (1.0 - sqrt_1_2)) as u64;
    if n * (n - 1) == 2 * b * (b - 1) {
        println!("n:{} b:{}", n, b);
        if n > 1000_000_000_000_u64 {
            break;
        }
        n = (n as f64 * 5.8284) as u64;
    }
    n += 1;
}

原来数学家已经给这个序列起了名字：Diophantine pairs，计算公式为：

a(0) = 1

a(1) = 3

a(n) = 6 * a(n-1) - a(n-2) - 2

感兴趣的朋友，可以在这个网站中找到有关这个序列的详细介绍：https://oeis.org/A011900

--- END ---