专栏首页申龙斌的程序人生通过欧拉计划学Rust编程(第650题)

通过欧拉计划学Rust编程(第650题)

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

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

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

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

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

这次解答的是第650题:

https://projecteuler.net/problem=650

题目描述:

这道题挺折磨人,花了我好几天时间才搞定。

如果你会其它编程语言,也可以一试。

解题过程:

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

第一步: 手工试算几个,找找规律

可以看到B(n)的求解与排列组合数的因子有关系。

B(2) = 2^1
B(3) = 3^2
B(4) = (2^5) * (3^1)
B(5) = (2^2) * (5^4)

如果一个数能够分解成几个素数乘积,利用等比数列的求和公式,可以推导出所有因子之和的公式:

公式中只有2个不同的质因子,可以推广到更多质因子的情况。

现在可以暴力计算了:

extern crate num_bigint;
use num_bigint::BigUint;

fn main() {
    // 先把一些阶乘计算好,保存起来
    let mut fact = vec![BigUint::from(1 as u64); 101];
    let mut a = BigUint::from(1 as u64);
    for n in 2..=100 {
        a *= BigUint::from(n as u64);
        fact[n] = a.clone();
    }
    //println!("{:?}", fact);

    let mut s = 0;
    for n in 1..=10 {
        let mut prod = BigUint::from(1_u64);
        for r in 1..=n {
            let comb = &fact[n] / &fact[r] / &fact[n - r];
            prod *= &comb;
            //println!("{} {} {}", n, r, comb.to_string());
        }
        let b = prod.to_string().parse::<u64>().unwrap();
        println!("B({}) = {}", n, b);

        let f_all = primes::factors(b);
        let f_uniq = primes::factors_uniq(b);
        //println!("{:?}", f_all);
        //println!("{:?}", f_uniq);

        let mut d = 1;
        for f in f_uniq {
            let c = f_all.iter().filter(|&n| *n == f).count();
            d *= (f.pow(c as u32 + 1) - 1) / (f-1);
        }
        //println!("D({}) = {}", n, d);

        s += d;
        println!("S({}) = {}", n, s);
    }
}

可以正确地计算出S(10),但在计算S(11)时就会出现溢出错误。

第二步:

溢出发生在pow()函数的计算上,求排列组合数和求阶乘的运算量太大,没必要把乘积计算出来,可以将因子保存在一个向量中,不断添加和删除相应的元素即可。

extern crate num_bigint;
use num_bigint::BigUint;

fn main() {
    let mut s = 0;
    for n in 1..=100 {
        let mut factors = vec![];
        for i in 1..=n {
            let mut f = comb_factors(n, i);
            factors.append(&mut f);
        }
        factors.sort();

        let d = factors_sum(&factors);
        println!("D({}) = {:?}", n, d);

        s = (s + d) % 1_000_000_007_u64;
        println!("S({}) = {}", n, s);
    }
}

fn factors_sum(v: &Vec<u64>) -> u64 {
    let mut uniq = v.clone();
    uniq.dedup();

    let mut prod = BigUint::from(1_u64);
    for p in uniq {
        let c = v.iter().filter(|&x| *x == p).count() as u64;
        let t = (big_pow(p, c + 1) - BigUint::from(1_u64)) / (BigUint::from(p - 1));
        //println!("{} {} {}", p, c, t);
        prod = prod * t % 1_000_000_007_u64;
    }
    let prod = prod % 1_000_000_007_u64;
    prod.to_string().parse::<u64>().unwrap()
}

fn big_pow(a: u64, b: u64) -> BigUint {
    let mut prod = BigUint::from(1 as u64);
    for _i in 0..b {
        prod *= BigUint::from(a as u64);
    }
    prod
}

// va中元素已经排好序
fn vec_remove(va: &mut Vec<u64>, vb: &Vec<u64>) {
    for item in vb {
        let index = va.binary_search(&item).unwrap();
        //println!("{:?} {:?} {}", va, vb, index);
        va.remove(index);
    }
}

fn comb_factors(m: u64, n: u64) -> Vec<u64> {
    let mut factors = vec![];
    let mut x = m;
    for i in 0..n {
        let mut f = primes::factors(x);
        factors.append(&mut f);
        x -= 1;
    }
    factors.sort();
    //println!("{:?}", factors);
    for i in 2..=n {
        let f = primes::factors(i);
        //println!("{} {:?}", n, f);
        vec_remove(&mut factors, &f);
    }
    factors.to_vec()
}

这次可以很容易地计算出S(100),但要计算S(1000)则要花相当长的时间,主要是因为数组越来越大,数组排序太花时间。

还得优化。

第三步 用HashMap

数组中元素的排序太慢,尝试换成字典来存储因子,在Rust中用HashMap来实现。比如,B(10)中含有12个2,10个3,8个5,3个7。

而B(100)中含有如下因子,字典中key是质因子,value是个数。

{97: 93, 19: 65, 83: 65, 41: 44, 59: 17, 79: 57, 61: 21, 67: 33, 2: 335, 3: 192, 37: 20, 23: 56, 31: 69, 47: 80, 13: 21, 71: 41, 7: 148, 73: 45, 11: 81,
43: 56, 17: 5, 53: 5, 5: 176, 89: 77, 29: 45}

修改后的代码:

extern crate num_bigint;
use num_bigint::BigUint;

use std::collections::HashMap;
fn main() {
    let mut s = 0;
    for n in 1..=10000 {
        let map = comb_factors_hash_map(n);
        //println!("{:?}", map);

        let temp = factors_hash_map_sum(&map);
        //println!("D({}) = {:?}", n, temp);

        s = (s + temp) % 1_000_000_007_u64;
        println!("S({}) = {}", n, s);
    }
}

fn big_pow(a: u64, b: u64) -> BigUint {
    let mut prod = BigUint::from(1 as u64);
    for _i in 0..b {
        prod *= BigUint::from(a as u64);
    }
    prod
}

fn factors_hash_map_sum(map: &HashMap<u64, u64>) -> u64 {
    let mut prod = BigUint::from(1_u64);
    for (&f, count) in map {
        let t = (big_pow(f, count + 1) - BigUint::from(1_u64)) / (BigUint::from(f - 1));
        prod = prod * t % 1_000_000_007_u64;
    }
    let prod = prod % 1_000_000_007_u64;
    prod.to_string().parse::<u64>().unwrap()
}

fn comb_factors_hash_map(x: u64) -> HashMap<u64, u64> {
    let mut map = HashMap::new();

    let mut count = x as i64 - 1;
    for n in (2..=x).rev() {
        let f = primes::factors(n);
        let a = factors_to_hash_map(&f);
        if count >= 0 {
            hash_map_add_count(&mut map, &a, count as u64);
        } else {
            hash_map_substract_count(&mut map, &a, count.abs() as u64);
        }
        count -= 2;
    }
    map
}

fn factors_to_hash_map(factors: &Vec<u64>) -> HashMap<u64, u64> {
    let mut map = HashMap::new();
    for f in factors {
        let v = map.get(f).cloned(); // 如果不写cloned(),有警告,不理解原因
        match v {
            Some(x) => {
                map.insert(*f, x + 1);
            }
            None => {
                map.insert(*f, 1);
            }
        }
    }
    map
}

fn hash_map_add_count(map: &mut HashMap<u64, u64>, a: &HashMap<u64, u64>, times: u64) {
    for (f, count) in a {
        let v = map.get(f).cloned();
        match v {
            Some(x) => {
                map.insert(*f, x + count * times);
            }
            None => {
                map.insert(*f, *count * times);
            }
        }
    }
}

fn hash_map_substract_count(map: &mut HashMap<u64, u64>, a: &HashMap<u64, u64>, times: u64) {
    for (f, count) in a {
        let v = map.get(f).cloned();
        match v {
            Some(x) => {
                map.insert(*f, x - count * times);
            }
            None => {}
        }
    }
}

这次可以轻松地计算到S(1000),但后面越来越慢。还得继续优化。

第四步 加速pow()运算

计算一个数的n次方,有一个常用的算法可以加速,称为Exponentiation by squaring算法,相关网站:

https://en.wikipedia.org/wiki/Exponentiation_by_squaring

修改后程序运行速度提高,多运行一段时间,可以计算出S(5000)。

fn big_pow(base: u64, exp: u64) -> BigUint {
    let mut result = BigUint::from(1 as u64);
    let mut b = BigUint::from(base);
    let mut e = exp;
    while e > 0 {
        if e % 2 == 1 {
            result = result * &b;
        }
        e = e >> 1;
        b = &b * &b;
    }
    result
}

第五步 递推公式

程序运行几个小时,应该能够计算出S(20000),在它计算的过程中,我还要尝试进一步的优化。但在这之后,我开始走弯路了,尝试缓存一些中间的计算结果来进行加速,效果都不理想。

google "euler project problem 650",发现一篇文章。

http://goatleaps.xyz/euler/maths/Project-Euler-650.html

里面提到一个B(n-1)递推出B(n)的公式,可以进一步优化。

我也用笨办法推导出了这条公式:

use std::collections::HashMap;

const MODULUS: u64 = 1_000_000_007_u64;
extern crate num_bigint;
use num_bigint::BigUint;

fn main() {
    let mut factorial_factors = HashMap::new();
    let mut map = HashMap::new();
    let mut s = 1;
    for n in 2..=20000 {
        let factor_n = factors_map(n);
        map_add_count(&mut map, &factor_n, n);

        // 缓存 n! 阶乘的因子
        map_add(&mut factorial_factors, &factor_n);

        map_substract(&mut map, &factorial_factors);

        //println!("{:?}", &map);
        let d = map_sum(&map);
        
        s = (s + d) % MODULUS;
        if n == 10 {
            assert_eq!(s, 141740594713218418 % MODULUS);
        }
        if n == 100 {
            assert_eq!(s, 332792866_u64);
        }
        if n % 100 == 0 {
            println!("D({}) = {:?} \t S({}) = {}", n, d, n, s);
        }
    }
    println!("{}", s);

}

fn map_sum(map: &HashMap<u64, u64>) -> u64 {
    let mut prod = BigUint::from(1_u64);
    for (&f, count) in map {
        let t = (big_pow(f, count + 1) - BigUint::from(1_u64)) / (BigUint::from(f - 1));
        prod = prod * t % 1_000_000_007_u64;
    }
    let prod = prod % 1_000_000_007_u64;
    prod.to_string().parse::<u64>().unwrap()
}

fn big_pow(base: u64, exp: u64) -> BigUint {
    let mut result = BigUint::from(1 as u64);
    let mut b = BigUint::from(base);
    let mut e = exp;
    while e > 0 {
        if e % 2 == 1 {
            result = result * &b;
        }
        e = e >> 1;
        b = &b * &b;
    }
    result
}

fn factors_map(n:u64) -> HashMap<u64, u64> {
    let mut map = HashMap::new();
    let all_factors = primes::factors(n);
    for f in &all_factors {
        let v = map.get(f).cloned(); // 如果不写cloned(),有警告,不理解原因
        match v {
            Some(x) => {
                map.insert(*f, x + 1);
            }
            None => {
                map.insert(*f, 1);
            }
        }
    }
    map
}

fn map_add(map: &mut HashMap<u64, u64>, a: &HashMap<u64, u64>) {
    for (f, count) in a {
        let v = map.get(f).cloned();
        match v {
            Some(x) => {
                map.insert(*f, x + count);
            }
            None => {
                map.insert(*f, *count);
            }
        }
    }
}

fn map_add_count(map: &mut HashMap<u64, u64>, a: &HashMap<u64, u64>, times: u64) {
    for (f, count) in a {
        let v = map.get(f).cloned();
        match v {
            Some(x) => {
                map.insert(*f, x + count * times);
            }
            None => {
                map.insert(*f, *count * times);
            }
        }
    }
}

fn map_substract(map: &mut HashMap<u64, u64>, a: &HashMap<u64, u64>) {
    for (f, count) in a {
        let v = map.get(f).cloned();
        match v {
            Some(x) => {
                map.insert(*f, x - count);
            }
            None => {}
        }
    }
}

程序并没有想像中提高运行速度,还得想办法。

第六步 倒数的模运算

问题仍出在计算因子和的公式上,计算乘方时用大整数库,当数字越来越大时,速度越来越慢。

这时得利用模运算的数学知识了,模运算有一些基本性质:如果 a ≡ b (mod m) 且有 c ≡ d (mod m),那么下面的模运算律成立:

a + c ≡ b + d (mod m)
a − c ≡ b − d (mod m)
a × c ≡ b × d (mod m)

这里可以看到加、减、乘都满足这种同余定理,但没有除法!

问题没完,在数学运算中:a / b == a * (1 / b),其中1 / b叫b的倒数。那么在模运算中,有没有这种类似于倒数的存在呢?答案是有!它就是乘法逆元!

相关文章链接: http://conw.net/archives/6/

有人已经写好了Rust的代码,直接可以拿来用。

http://rosettacode.org/wiki/Modular_inverse#Rust

pow()乘方函数也有带模运算的Rust代码,直接拿来用。

https://rob.co.bb/posts/2019-02-10-modular-exponentiation-in-rust/

这次优化彻底去掉了大整数运算,速度提升了不知多少倍。

在我的笔记本电脑上,8秒钟计算出S(20000)!

最终的源代码:

use std::collections::HashMap;

const MODULUS: u64 = 1_000_000_007_u64;

fn main() {
    let mut factorial_factors = HashMap::new();
    let mut map = HashMap::new();
    let mut s = 1;
    for n in 2..=20000 {
        let factor_n = factors_map(n);
        map_add_count(&mut map, &factor_n, n);

        // 缓存 n! 阶乘的因子
        map_add(&mut factorial_factors, &factor_n);

        map_substract(&mut map, &factorial_factors);

        //println!("{:?}", &map);
        let d = map_sum(&map);
        
        s = (s + d) % MODULUS;
        if n == 10 {
            assert_eq!(s, 141740594713218418 % MODULUS);
        }
        if n == 100 {
            assert_eq!(s, 332792866_u64);
        }
        if n % 100 == 0 {
            println!("D({}) = {:?} \t S({}) = {}", n, d, n, s);
        }
    }
    println!("{}", s);

}

// http://rosettacode.org/wiki/Modular_inverse#Rust
// 求乘法逆元
fn mod_inv(a: isize, module: isize) -> isize {
    let mut mn = (module, a);
    let mut xy = (0, 1);

    while mn.1 != 0 {
        xy = (xy.1, xy.0 - (mn.0 / mn.1) * xy.1);
        mn = (mn.1, mn.0 % mn.1);
    }

    while xy.0 < 0 {
        xy.0 += module;
    }
    xy.0
}

// https://rob.co.bb/posts/2019-02-10-modular-exponentiation-in-rust/
fn mod_pow(mut base: u64, mut exp: u64, modulus: u64) -> u64 {
    let mut result = 1;
    base = base % modulus;
    while exp > 0 {
        if exp % 2 == 1 {
            result = result * base % modulus;
        }
        exp = exp >> 1;
        base = base * base % modulus
    }
    result
}


fn map_sum(map: &HashMap<u64, u64>) -> u64 {
    let mut prod = 1;
    for (&f, count) in map {
        if *count > 0 {
            // 计算 f^(count+1) / (f-1)
            let temp = mod_pow(f, count + 1, MODULUS) - 1;
            // 计算(f-1)的乘法逆元
            let inv = mod_inv(f as isize - 1, MODULUS as isize) as u64;
            //println!("inv:{}", inv);
            let temp = temp * inv % MODULUS;
            prod = prod * temp % MODULUS;
            //println!("f:{} count+1:{} prod: {}", f, count+1, prod);
        }
    }
    prod
}

fn factors_map(n:u64) -> HashMap<u64, u64> {
    let mut map = HashMap::new();
    let all_factors = primes::factors(n);
    for f in &all_factors {
        let v = map.get(f).cloned(); // 如果不写cloned(),有警告,不理解原因
        match v {
            Some(x) => {
                map.insert(*f, x + 1);
            }
            None => {
                map.insert(*f, 1);
            }
        }
    }
    map
}

fn map_add(map: &mut HashMap<u64, u64>, a: &HashMap<u64, u64>) {
    for (f, count) in a {
        let v = map.get(f).cloned();
        match v {
            Some(x) => {
                map.insert(*f, x + count);
            }
            None => {
                map.insert(*f, *count);
            }
        }
    }
}

fn map_add_count(map: &mut HashMap<u64, u64>, a: &HashMap<u64, u64>, times: u64) {
    for (f, count) in a {
        let v = map.get(f).cloned();
        match v {
            Some(x) => {
                map.insert(*f, x + count * times);
            }
            None => {
                map.insert(*f, *count * times);
            }
        }
    }
}

fn map_substract(map: &mut HashMap<u64, u64>, a: &HashMap<u64, u64>) {
    for (f, count) in a {
        let v = map.get(f).cloned();
        match v {
            Some(x) => {
                map.insert(*f, x - count);
            }
            None => {}
        }
    }
}

优化要点:

1)不要轻易用大整数运算库

2)质因子分解

3)同余定理

4)Exponentiation by squaring乘方运算加速

5)找到递推公式

6)乘法逆元

卡了我好几天的优化算法是因为以前不知道“乘法逆元”的存在。

--- END ---

本文分享自微信公众号 - 申龙斌的程序人生(slbGTD)

原文出处及转载信息见文内详细说明,如有侵权,请联系 yunjia_community@tencent.com 删除。

原始发表时间:2020-01-26

本文参与腾讯云自媒体分享计划,欢迎正在阅读的你也加入,一起分享。

我来说两句

0 条评论
登录 后参与评论

相关文章

  • 用欧拉计划学Rust编程(第92题)

    将一个数的所有数字的平方相加得到一个新的数,不断重复直到新的数已经出现过为止,这构成了一条数字链。

    申龙斌
  • 云币网及KYC【区块链生存训练】

    李笑来在7月5日发布了Press.One即将 ICO 的消息,大批小白开始在云币网注册开户,我的“区块链生存训练”饭团也在一夜之间加入了40多人。有位新人在饭团...

    申龙斌
  • 通过欧拉计划学Rust编程(第684题)

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

    申龙斌
  • python的map()函数

    原文链接:https://www.runoob.com/python/python-func-map.html

    于小勇
  • 高可用数据库主从复制延时的解决

    如果主从复制之间出现延时,就会影响主从数据的一致性。 此时发生容灾切换,且在新的主库写入了数据,那么从业务角度上,会产生意想不到的严重后果。 复制延时问题,...

    山山仙人
  • Apache Maven WAR Plugin

    Apache Maven WAR插件负责收集所有工件的依赖性,web应用程序的类和资源,包装成一个web应用程序存档。

    用户3004328
  • macos系统占用闪存过多的解决方案

    随着长期使用(通常也没有关机习惯),mac的缓存垃圾越堆越多,最终系统占到了80g以上,严重挤占了本应留给其他文件的闪存空间,这里谈一下“系统”中可能存在的垃圾...

    gojam
  • TW洞见 | 郑晔:测试时间

    在单元测试中,与时间相关的测试总是让人很头疼。举个例子,我们希望做一个定期过期缓存,比如30分钟过期,这该怎么测试呢?等30分钟?那要是过期时间是3天,你打算把...

    ThoughtWorks
  • php7.0添加curl,mbstring,pdo,openssl扩展

    进入php安装源码,若安装源码在/var/local/lnmp1/php-7.0,安装的途径在/var/local/lnmp/php-7.0

    93年的老男孩
  • Java管理扩展特殊MBean之MXBean学习

    MXBean是一种引用预定义数据类型的MBean。通过这种方式,您可以确保任何客户机(包括远程客户机)都可以使用您的MBean,而不需要客户机访问代表MBean...

    用户1289394

扫码关注云+社区

领取腾讯云代金券