通过欧拉计划学Rust编程（第684题）

申龙斌

发布于 2020-01-15 16:10:18

54000

代码可运行

文章被收录于专栏：申龙斌的程序人生申龙斌的程序人生

运行总次数：0

代码可运行

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

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

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

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

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

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

题目描述：

定义s(n)是数字和为n的最小整数，例如s(10)=19。

记S(k) = s(1) + s(2) + ... + s(k)，可以知道 S(20)=1074。

设斐波那契数列f(n)按如下方式定义：

求：

请

先

不

要

直

接

看

答

案

，

最

好

自

己

先

尝

试

一

下

。

解题过程：

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

第一步：

题目中知道S(20)=1074，那就先求S(20)。

手算前20个，发现一个规律。每9个为一组，后面的数字是9，前面的数字从0到8。

可以推导出一个公式：s(n) = (a+1) * 10 ^ b - 1

S(20)很容易求出来：

let mut ss = 0;
for n in 1..=20 {
    let a = n % 9;
    let b = n / 9;
    let s = (a + 1) * 10_u32.pow(b) - 1;
    println!("{}", s);
    ss += s;
}
println!("S(20): {}", ss);

但求S(200)时就会溢出，因为 10 ^ b是一个超出u64范围的大整数。

第二步

把90个斐波那契数求出来，以后要用。

let mut fib = [0_u64; 91];
fib[1] = 1;
for i in 2..=90 {
    fib[i] = fib[i - 1] + fib[i - 2];
    println!("fib({}): {}", i, fib[i]);
}

运行一下，可以发现第90个数非常非常大。

fib(88): 1100087778366101931
fib(89): 1779979416004714189
fib(90): 2880067194370816120

第三步

首先能够想到的是用大整数库num-bigint优化算法。

fn fs(n: u64) -> u64 {
    let a = n % 9;
    let b = n / 9;
    let mut s = BigUint::from(a + 1);
    for i in 0..b {
        s = s * BigUint::from(10_u64);
    }
    s = s - BigUint::from(1_u64);
    let result = s % BigUint::from(1_000_000_007_u64);
    result.to_string().parse::<u64>().unwrap()
}

现在可以比较快地计算出s(20000)和S(20000)，但离目标2880067194370816120还有相当大的距离，bigint这条路不通，还得改进算法。

第四步

看看S(n)的计算是否还有其它规律。每9个为一组，计算9个数之和，可以找到规律。

可以发现，当n = 9 * m - 1时，有公式：

S(n) = 5 * 10 ^ m - 5 - 9 * m

有了这个公式，在计算S(20)时，可以先快速计算出S(17)，再加上s(18)+s(19)+s(20)就可以得到最终结果，算法复杂度相当于计算四次s(n)。

现在仍有一个关键问题没有解决，10 ^ m 是一个非常大的数，必须找到快速计算10 ^ m mod 1_000_000_007_u64 的办法。

这里要利用费马小定理，如果p是一个质数，而整数a不是p的倍数，则有a^(p-1) mod p = 1。

看一个特殊的例子，a=10, p=7时，有助于大致理解费马小定理的含义。

有个这个定理，10 ^ m mod 1_000_000_007_u64 = 10 ^ (m mod 1_000_000_006_u64) mod 1_000_000_007_u64，可以大大加速计算过程，25秒计算出结果。

最后的源代码：

use std::time::SystemTime;

const PRIME: u64 = 1_000_000_007_u64;

#[macro_use]
extern crate lazy_static;
lazy_static! {
    static ref ARRAY: Vec<u64> = {
        println!("initializing ARRAY ...");
        let mut arr = vec![1];
        let mut x = 1;
        for _i in 1..PRIME - 1 {
            x = x * 10 % PRIME;
            arr.push(x as u64);
        }
        arr
    };
}

fn main() {
    let start = SystemTime::now();
    let mut result = 0;
    let mut fib = vec![0_u64, 1];
    for i in 2..=90 {
        let n = fib[i - 1] + fib[i - 2];
        fib.push(n);
        let ss = fss(n);
        result = (result + ss) % PRIME;
        println!("n:{} S:{} result: {}", n, ss, result);
    }
    println!("{:?}", start.elapsed());
}

fn ten_power_mod(n: u64) -> u64 {
    let m = n % (PRIME - 1);
    ARRAY[m as usize]
}

fn fs(n: u64) -> u64 {
    let a = n % 9;
    let b = n / 9;
    let s = (a + 1) * ten_power_mod(b) - 1;
    s % PRIME
}

fn sum_group(m: u64) -> u64 {
    let temp = (9 * m) % PRIME;
    let s = 5 * ten_power_mod(m) + PRIME - temp - 5;
    s % PRIME
}

fn fss(n: u64) -> u64 {
    let m = n / 9;
    let mut s = sum_group(m);
    for i in 9 * m..=n {
        s += fs(i);
    }
    s % PRIME
}

你可能会错过：