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

申龙斌

发布于 2020-02-25 15:32:56

4110

发布于 2020-02-25 15:32:56

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

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

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

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

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

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

这次解答的是第345题：

https://projecteuler.net/problem=345

题目描述：

矩阵和

在一个矩阵中选择部分元素，使得每个元素都是所在的行和列中唯一被选中的元素，这些元素的和的最大值定义为矩阵的矩阵和。例如，如下矩阵的矩阵和为3315 ( = 863 + 383 + 343 + 959 + 767)：

7 53 183 439 863

497 383 563 79 973

287 63 343 169 583

627 343 773 959 943

767 473 103 699 303

求如下矩阵的矩阵和：

  7  53 183 439 863 497 383 563  79 973 287  63 343 169 583
627 343 773 959 943 767 473 103 699 303 957 703 583 639 913
447 283 463  29  23 487 463 993 119 883 327 493 423 159 743
217 623   3 399 853 407 103 983  89 463 290 516 212 462 350
960 376 682 962 300 780 486 502 912 800 250 346 172 812 350
870 456 192 162 593 473 915  45 989 873 823 965 425 329 803
973 965 905 919 133 673 665 235 509 613 673 815 165 992 326
322 148 972 962 286 255 941 541 265 323 925 281 601  95 973
445 721  11 525 473  65 511 164 138 672  18 428 154 448 848
414 456 310 312 798 104 566 520 302 248 694 976 430 392 198
184 829 373 181 631 101 969 613 840 740 778 458 284 760 390
821 461 843 513  17 901 711 993 293 157 274  94 192 156 574
 34 124   4 878 450 476 712 914 838 669 875 299 823 329 699
815 559 813 459 522 788 168 586 966 232 308 833 251 631 107
813 883 451 509 615  77 281 613 459 205 380 274 302  35 805

请

先

不

要

直

接

看

答

案

，

最

好

自

己

先

尝

试

一

下

解题过程：

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

第一步：先解决5x5矩阵的简单问题

首先能够想到的是递归的办法，5x5的矩阵和可以转换成求4x4的矩阵和，递归下去就可以找到答案。

依次划掉第0列的5个元素，可以得到相应的剩余矩阵，可以看到剩余矩阵的维数变成4x4，一直递推到1x1。

为了代码写起来方便，还有将来的性能考虑，用一维数组。

先求剩余矩阵。

fn residual_matrix(mat: &Vec<usize>, row:usize, col:usize) -> Vec<usize> {
    let dim = (mat.len() as f64).sqrt() as usize;
    let mut m2 = vec![];
    for index in 0..mat.len() {
        let r = index / dim;
        let c = index % dim;
        if row != r && col != c {
            m2.push(mat[index]);
        }
    }
    m2
}

递归程序很容易写出来，运行后得到3315，正确答案。

fn main() {
    let mat = vec![
        7, 53, 183, 439, 863, // row 0
        497, 383, 563, 79, 973, // row 1
        287, 63, 343, 169, 583, // row 2
        627, 343, 773, 959, 943, // row 3
        767, 473, 103, 699, 303, // row 4
    ];
    println!("{}", matrix_sum(&mat));
}

fn matrix_sum(mat: &Vec<usize>) -> usize {
    if mat.len() == 1 {
        return mat[0];
    }

    let dim = (mat.len() as f64).sqrt() as usize;
    let mut max_sum = 0;
    for row in 0..dim {
        let m2 = residual_matrix(&mat, row, 0);
        let temp = mat[row*dim] + matrix_sum(&m2);
        if temp > max_sum {
            max_sum = temp;
        }
    }
    max_sum
}

把15x15的矩阵换上，程序运行非常慢，打印了几个中间变量，发现几个小时也算不出来。

第二步：优化性能

15层递归，中间有大量重复的计算，例如下面的两种情况，剩余矩阵是相同的，但会重复计算，类似的情况非常多，所以性能非常差。得对中间结果进行缓存处理。

首先想到一种缓存办法，生成一个与原始矩阵一样大小的矩阵，里面充满很大的随机数。每个元素也相应位置上的随机数进行异或运算，生成哈希值。这样会保证相同的矩阵生成的哈希值肯定是相同的，只要哈希表的容量足够大（0xFFFFFFFF），就不会有哈希冲突的情况发生。

fn main() {
    let mat = vec![
        7, 53, 183, 439, 863, 497, 383, 563, 79, 973, 287, 63, 343, 169, 583, 627, 343, 773, 959,
        943, 767, 473, 103, 699, 303, 957, 703, 583, 639, 913, 447, 283, 463, 29, 23, 487, 463,
        993, 119, 883, 327, 493, 423, 159, 743, 217, 623, 3, 399, 853, 407, 103, 983, 89, 463, 290,
        516, 212, 462, 350, 960, 376, 682, 962, 300, 780, 486, 502, 912, 800, 250, 346, 172, 812,
        350, 870, 456, 192, 162, 593, 473, 915, 45, 989, 873, 823, 965, 425, 329, 803, 973, 965,
        905, 919, 133, 673, 665, 235, 509, 613, 673, 815, 165, 992, 326, 322, 148, 972, 962, 286,
        255, 941, 541, 265, 323, 925, 281, 601, 95, 973, 445, 721, 11, 525, 473, 65, 511, 164, 138,
        672, 18, 428, 154, 448, 848, 414, 456, 310, 312, 798, 104, 566, 520, 302, 248, 694, 976,
        430, 392, 198, 184, 829, 373, 181, 631, 101, 969, 613, 840, 740, 778, 458, 284, 760, 390,
        821, 461, 843, 513, 17, 901, 711, 993, 293, 157, 274, 94, 192, 156, 574, 34, 124, 4, 878,
        450, 476, 712, 914, 838, 669, 875, 299, 823, 329, 699, 815, 559, 813, 459, 522, 788, 168,
        586, 966, 232, 308, 833, 251, 631, 107, 813, 883, 451, 509, 615, 77, 281, 613, 459, 205,
        380, 274, 302, 35, 805,
    ];

    let mut hash_mat = vec![];
    for _i in 0..mat.len() {
        hash_mat.push(rand::random::<usize>());
    }
    //println!("{:?}", hash_mat );

    let mut cache = vec![0; 0xFFFFFFFF];
    println!("{}", matrix_sum(&mat, &hash_mat, &mut cache));
}

fn matrix_sum(mat: &Vec<usize>, hash_mat: &Vec<usize>, cache: &mut Vec<usize>) -> usize {
    if mat.len() == 1 {
        return mat[0];
    }

    let hash = hash_code(mat, hash_mat);
    if cache[hash] > 0 {
        return cache[hash]; // 缓存命中
    }

    let dim = (mat.len() as f64).sqrt() as usize;
    let mut max_sum = 0;
    for row in 0..dim {
        let m2 = residual_matrix(&mat, row, 0);
        let hash_mat2 = residual_matrix(&hash_mat, row, 0);
        let temp = mat[row * dim] + matrix_sum(&m2, &hash_mat2, cache);
        if temp > max_sum {
            max_sum = temp;
        }
    }
    cache[hash] = max_sum; // 写入缓存
    max_sum
}

// 剩余矩阵
fn residual_matrix(mat: &Vec<usize>, row: usize, col: usize) -> Vec<usize> {
    let dim = (mat.len() as f64).sqrt() as usize;
    let mut m2 = vec![];
    for index in 0..mat.len() {
        let r = index / dim;
        let c = index % dim;
        if row != r && col != c {
            m2.push(mat[index]);
        }
    }
    m2
}

fn hash_code(mat: &Vec<usize>, hash_mat: &Vec<usize>) -> usize {
    let mut hash = 0;
    for i in 0..mat.len() {
        hash ^= mat[i] ^ hash_mat[i];
    }
    hash & 0xFFFFFFFF
}

改进的程序速度很快，可以在1、2秒内得到结果。

第三步：改进一下哈希算法

这里的哈希算法不太美观，引入了一个哈希表矩阵，为了不引起哈希冲突，需要较大内存。这个算法能够计算出最大的矩阵和，但具体选择了哪些位置上的元素？没有显示的答案。

这里引入一种简便的记号，[4, 1, 2, 3, 0]表示选中了(第0列, 第4行)、(1, 1)、(2, 2)、(3, 3)和(4, 0)这几个位置上的元素，这样用一维向量就可以记录选中的位置。

如果某个位置被选中，可以用二进制1表示，这样哈希编码后很省空间，2^15就够用。

fn hash_code(path: &Vec<usize>) -> usize {
    let mut hash = 0;
    for p in path {
        hash += 1_usize << p;
    }
    hash
}

当然，也可以写成一行

fn hash_code(path: &Vec<usize>) -> usize {
    path.iter().map(|x| 1_usize << x).sum()
}

最后的程序：

fn main() {
    let mat = vec![
        7, 53, 183, 439, 863, 497, 383, 563, 79, 973, 287, 63, 343, 169, 583, 627, 343, 773, 959,
        943, 767, 473, 103, 699, 303, 957, 703, 583, 639, 913, 447, 283, 463, 29, 23, 487, 463,
        993, 119, 883, 327, 493, 423, 159, 743, 217, 623, 3, 399, 853, 407, 103, 983, 89, 463, 290,
        516, 212, 462, 350, 960, 376, 682, 962, 300, 780, 486, 502, 912, 800, 250, 346, 172, 812,
        350, 870, 456, 192, 162, 593, 473, 915, 45, 989, 873, 823, 965, 425, 329, 803, 973, 965,
        905, 919, 133, 673, 665, 235, 509, 613, 673, 815, 165, 992, 326, 322, 148, 972, 962, 286,
        255, 941, 541, 265, 323, 925, 281, 601, 95, 973, 445, 721, 11, 525, 473, 65, 511, 164, 138,
        672, 18, 428, 154, 448, 848, 414, 456, 310, 312, 798, 104, 566, 520, 302, 248, 694, 976,
        430, 392, 198, 184, 829, 373, 181, 631, 101, 969, 613, 840, 740, 778, 458, 284, 760, 390,
        821, 461, 843, 513, 17, 901, 711, 993, 293, 157, 274, 94, 192, 156, 574, 34, 124, 4, 878,
        450, 476, 712, 914, 838, 669, 875, 299, 823, 329, 699, 815, 559, 813, 459, 522, 788, 168,
        586, 966, 232, 308, 833, 251, 631, 107, 813, 883, 451, 509, 615, 77, 281, 613, 459, 205,
        380, 274, 302, 35, 805,
    ];

    let mut cache = vec![0; 2_usize.pow(15)];
    let mut path = vec![];
    println!("{}", matrix_sum(&mut cache, &mat, &mut path));
}

fn matrix_sum(cache: &mut [usize], mat: &Vec<usize>, path: &mut Vec<usize>) -> usize {
    // 总列数
    let dim = (mat.len() as f64).sqrt() as usize;
    // 准备填充的列号
    let cur_col = path.len();
    if cur_col == dim {
        return 0;
    }

    let hash = hash_code(path);
    if cache[hash] > 0 {
        return cache[hash];
    }

    let mut max_sum = 0;
    for row in 0..dim {
        if !path.contains(&row) {
            path.push(row);
            let temp = mat[row * dim + cur_col] + matrix_sum(cache, &mat, path);
            if temp > max_sum {
                max_sum = temp;
            }
            path.pop();
        }
    }
    cache[hash] = max_sum;
    //println!("{:?} {} {}", &path, hash, max_sum);
    max_sum
}

fn hash_code(path: &Vec<usize>) -> usize {
    path.iter().map(|x| 1_usize<<x).sum()
}

这里还是没有得到最后选中了哪些元素，我想到了一个笨办法，由于缓存里有大量的中间结果，可以从第0列到第n列计算一遍，把这条路径找出来。这个输出路径的函数不太美观，但能用。

fn matrix_output(cache: &Vec<usize>, mat: &Vec<usize>) {
    let dim = (mat.len() as f64).sqrt() as usize;
    let mut path = vec![];
    for col in 0..dim {
        let mut selected_row = 999;
        let mut max = 0;
        for row in 0..dim {
            if !path.contains(&row) {
                path.push(row);
                let hash = hash_code(&path);
                if mat[row * dim + col] + cache[hash] > max {
                    max = mat[row * dim + col] + cache[hash];
                    selected_row = row;
                }
                path.pop();
            }
        }
        path.push(selected_row);
    }
    println!("path: {:?}", path);
}

第四步：遗传算法

除了暴力的递归算法外，还有其它算法，比如遗传算法，每次生成一定数量的后代，对矩阵上相应元素求和，和较大的属于优质后代，淘汰掉那些矩阵和较小的，不断迭代之后，就可以得到最优结果。

任意交换两个元素就可以产生一种选择元素的方法。

算法并不复杂，不需严格排序，将矩阵和最大的放在队列的最前面（用insert语句），其它的扔在后面就行（用push语句），是金子早晚会发光的，经过几代淘汰，优质的后代会跑到队列的前面。

fn main() {
    let mat = vec![
        7, 53, 183, 439, 863, 497, 383, 563, 79, 973, 287, 63, 343, 169, 583, 627, 343, 773, 959,
        943, 767, 473, 103, 699, 303, 957, 703, 583, 639, 913, 447, 283, 463, 29, 23, 487, 463,
        993, 119, 883, 327, 493, 423, 159, 743, 217, 623, 3, 399, 853, 407, 103, 983, 89, 463, 290,
        516, 212, 462, 350, 960, 376, 682, 962, 300, 780, 486, 502, 912, 800, 250, 346, 172, 812,
        350, 870, 456, 192, 162, 593, 473, 915, 45, 989, 873, 823, 965, 425, 329, 803, 973, 965,
        905, 919, 133, 673, 665, 235, 509, 613, 673, 815, 165, 992, 326, 322, 148, 972, 962, 286,
        255, 941, 541, 265, 323, 925, 281, 601, 95, 973, 445, 721, 11, 525, 473, 65, 511, 164, 138,
        672, 18, 428, 154, 448, 848, 414, 456, 310, 312, 798, 104, 566, 520, 302, 248, 694, 976,
        430, 392, 198, 184, 829, 373, 181, 631, 101, 969, 613, 840, 740, 778, 458, 284, 760, 390,
        821, 461, 843, 513, 17, 901, 711, 993, 293, 157, 274, 94, 192, 156, 574, 34, 124, 4, 878,
        450, 476, 712, 914, 838, 669, 875, 299, 823, 329, 699, 815, 559, 813, 459, 522, 788, 168,
        586, 966, 232, 308, 833, 251, 631, 107, 813, 883, 451, 509, 615, 77, 281, 613, 459, 205,
        380, 274, 302, 35, 805,
    ];

    let reproduce_num = 7;
    let mut paths = vec![vec![0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14]];
    let mut max = 0;
    for _generation in 0..20000 {
        for _i in 0..reproduce_num * reproduce_num {
            let path = &paths[rand::random::<usize>() % paths.len()];
            let path_new = swap(
                &path,
                rand::random::<usize>() % 15,
                rand::random::<usize>() % 15,
            );
            if !paths.contains(&path_new) {
                let ev = eval(&mat, &path_new);
                if ev > max {
                    max = ev;
                    paths.insert(0, path_new);
                } else {
                    paths.push(path_new);
                }
            }
        }
        // 只保留前面几个优质的
        if paths.len() > reproduce_num {
            for p in (reproduce_num..paths.len()).rev() {
                paths.remove(p);
            }
        }
    }
    println!("{} {:?}", max, paths[0]);
}

fn swap(path: &Vec<usize>, i: usize, j: usize) -> Vec<usize> {
    let mut path2 = path.clone();
    path2[i] = path[j];
    path2[j] = path[i];
    return path2;
}

fn eval(mat: &Vec<usize>, path: &Vec<usize>) -> usize {
    let dim = (mat.len() as f64).sqrt() as usize;
    let mut sum = 0;
    for (col, row) in path.iter().enumerate() {
        let index = row * dim + col;
        sum += mat[index];
    }
    sum
}

算法中繁殖系数用了7，尽量随机生成下一代，繁衍2万次之后，99%的情况下都会得到正确答案，提高繁衍次数会99.99%的得到正确答案。

遗传算法的缺点是算法不稳定，可能无法得到最优解，只能得到局部较优解。但这种算法也有优点，可以进行规模化，当矩阵更大时，算法仍然有效。只需考虑如何让后代有一定概率的变异性，后代更优异。

第五步：更复杂的算法

这类问题还可以用动态规划，匈牙利算法等，算法比较复杂，感兴趣的朋友自行探索吧。

--- END ---

我把解决这些问题的过程记录了下来，写成了一本《用欧拉计划学 Rust 编程》PDF电子书，请随意下载。

链接：https://pan.baidu.com/s/1NRfTwAcUFH-QS8jMwo6pqw

提取码：qfha

本文参与腾讯云自媒体同步曝光计划，分享自微信公众号。

原始发表：2020-02-15，如有侵权请联系 cloudcommunity@tencent.com 删除

编程算法