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

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

https://projecteuler.net/problem=650

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

```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);
}
}
```

```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()
}

```

```{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 => {}
}
}
}

```

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

```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
}
```

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

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

```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_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 + c ≡ b + d (mod m)
a − c ≡ b − d (mod m)
a × c ≡ b × d (mod m)
```

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

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

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

```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_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 ---

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应用程序存档。

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

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

• TW洞见 | 郑晔：测试时间

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

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

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

• Java管理扩展特殊MBean之MXBean学习

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