通过欧拉计划学Rust编程：第100题
由于研究Libra等数字货币编程技术的需要，学习了一段时间的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 ---