通过欧拉计划学Rust编程（第69题）
由于研究Libra等数字货币编程技术的需要，学习了一段时间的Rust编程，一不小心刷题上瘾。
“欧拉计划”的网址： https://projecteuler.net
英文如果不过关，可以到中文翻译的网站： http://pe-cn.github.io/
这个网站提供了几百道由易到难的数学问题，你可以用任何办法去解决它，当然主要还得靠编程，编程语言不限，论坛里已经有Java、C#、Python、Lisp、Haskell等各种解法，当然如果你直接用google搜索答案就没任何乐趣了。
这次解答的是第69题：
https://projecteuler.net/problem=69
题目描述：
欧拉总计函数与最大值
在小于n的数中，与n互质的数的数目记为欧拉总计函数φ(n)（有时也称为φ函数）。例如，因为1、2、4、5、7和8均小于9且与9互质，故φ(9)=6。
n | 互质的数 | φ(n) | n/φ(n) |
|---|---|---|---|
2 | 1 | 1 | 2 |
3 | 1, 2 | 2 | 1.5 |
4 | 1, 3 | 2 | 2 |
5 | 1, 2, 3, 4 | 4 | 1.25 |
6 | 1, 5 | 2 | 3 |
7 | 1, 2, 3, 4, 5, 6 | 6 | 1.1666… |
8 | 1, 3, 5, 7 | 4 | 2 |
9 | 1, 2, 4, 5, 7, 8 | 6 | 1.5 |
10 | 1, 3, 7, 9 | 4 | 2.5 |
可以看出，对于n ≤ 10，当n=6时n/φ(n)取得最大值。
当n ≤ 1,000,000时，求使得n/φ(n)取得最大值的n。
请
先
不
要
直
接
看
答
案
，
最
好
自
己
先
尝
试
一
下
解题过程：
遇到一个复杂的问题，可以先从简单的情况入手，寻找一些规律，再慢慢逼近最终的问题。
第一步： 暴力求解
先根据题意暴力求解。
fn main() {
let mut max_ratio = 0_f64;
for n in 2..=1_000_000 {
// 最大公约数为1，表示互质
let phi = (1..n).filter(|&x| gcd(x, n) == 1).count();
let ratio = (n as f64) / (phi as f64);
if ratio > max_ratio {
println!("n= {:6} phi={:6} n/phi= {:.4}", n, phi, ratio);
max_ratio = ratio;
}
}
}
// 最大公约数
fn gcd(a: u64, b: u64) -> u64 {
if b == 0 {
a
} else {
gcd(b, a % b)
}
}
可以得到部分结果，计算到2310非常快，计算到30030则要1分钟，计算到1百万估计得几天。
n= 2 phi= 1 n/phi= 2.0000
n= 6 phi= 2 n/phi= 3.0000
n= 30 phi= 8 n/phi= 3.7500
n= 210 phi= 48 n/phi= 4.3750
n= 2310 phi= 480 n/phi= 4.8125
n= 30030 phi= 5760 n/phi= 5.2135
第二步： 找规律
通过前面几个结果，可以发现是连续几个素数的乘积，所以问题可以变为求哪些连续素数的乘积小于1百万。
// 2 * 3 * 5 * 7 * 11 * ... 找到最多的素数，乘积在1百万以内
let mut pset = PrimeSet::new();
let mut prod = 1;
for p in pset.iter() {
if prod * p > 1_000_000 {
break;
}
prod *= p;
println!("prime={} {}", p, prod);
}
println!("{}", prod);
这里可以练习一下函数式编程，好像可读性很差。
let mut some_primes = (2..).filter(|&x| primes::is_prime(x));
let mut prod = 1;
let some_products = std::iter::repeat_with(|| {
let tmp = prod;
prod *= some_primes.next().unwrap();
tmp}).take_while(|&x| x <= 1_000_000)
.collect::<Vec<_>>();
println!("{:?}", some_products);
println!("{:?}", some_products.last().unwrap());
用itertools写一下，这样好一些。
let mut some_primes = (2..).filter(|&x| primes::is_prime(x));
let result = itertools::iterate(1, |&prod| prod * some_primes.next().unwrap())
.take_while(|&x| x <= 1_000_000)
.last()
.unwrap();
println!("{:?}", result);
第三步： 补一下数学知识
这个欧拉总计函数，可以推导出数学公式，网上很容易搜到一大堆相关内容。
let mut max_ratio = 0_f64;
for n in 2..=1_000_000 {
let all_factors = primes::factors(n);
let uniq_factors = primes::factors_uniq(n);
let mut phi = 1;
for p in uniq_factors {
let k = all_factors.iter().filter(|&x| *x == p).count();
phi *= p.pow(k as u32 - 1) * (p-1);
}
let ratio = (n as f64) / (phi as f64);
if ratio > max_ratio {
println!("n= {:6} phi={:6} n/phi= {:.4}", n, phi, ratio);
max_ratio = ratio;
}
}
看看后一种写法的例子：
用这个公式，程序写起来更简单。
let mut max_ratio = 0_f64;
for n in 2..=1_000_000 {
let uniq_factors = primes::factors_uniq(n);
let mut phi = n;
for p in uniq_factors {
phi = phi * (p-1) / p;
}
let ratio = (n as f64) / (phi as f64);
if ratio > max_ratio {
println!("n= {:6} phi={:6} n/phi= {:.4}", n, phi, ratio);
max_ratio = ratio;
}
}
--- END ---
我把解决这些问题的过程记录了下来，写成了一本《用欧拉计划学 Rust 编程》PDF电子书，请随意下载。
链接：https://pan.baidu.com/s/1NRfTwAcUFH-QS8jMwo6pqw
提取码：qfha