用欧拉计划学Rust编程，第323题，随机整数按位或运算

用欧拉计划学Rust编程，第323题，随机整数按位或运算

题目描述：

理解题意，一个数初始值为0，不断与一个随机数进行‘或’运算，直到32个二进制位都为1，所需要步数记为N。进行无数数这样的试验，问N的平均值是多少。

第一步：

先用蒙特卡罗模拟方法试验一下，看看答案的大概范围。

fn trial() -> u32 {
    let mut x = 0_u32;
    for i in 1.. {
        x |= rand::random::<u32>();
        if x == 0xFFFFFFFF {
            return i;
        }
    }
    0
}

运行这个trail()函数10万次，取平均值，可以得到答案在6.35附近。

第二步：

原题的x是32位整数，可以从最简单的情况入手，假设x为1位的整数。
记此时的期望值为E(1)。
进行第一步操作之后，有两种可能：
第一种情况，1/2的概率，直接变成了1，即一步完成
第二种情况，1/2的概率，仍是0，回到了初始情况，即需要1+E(1)步完成。
这时有：E(1) = 1/2 + 1/2 * (1 + E(1))
可以解出：E(1) = 2

第三步：

假设x是2位整数，此时的期望值记为E(2)。
进行第一步操作之后，有三种可能：
第一种情况，1/4的概率，直接变成了(11)binary，即一步完成
第二种情况，2/4的概率，有一位变为1，即(01)binary或(10)binary,这时是E(1)时的情况，期望值为1+E(1)。
第三种情况，1/4的概率，仍为0，这时是1+E(2)。
列出式子，E(2) = 1/4 + 2/4 * (1+ E(1)) + 1/4 * (1 + E(2))
3/4 * E(2) = 8/4
E(2) = 8/3

第四步：

这时，你可能还没有总结出规律，再来一次，假设x是3位整数，此时的期望值记为E(3)。
进行第一步操作之后，有四种可能：
第一种情况，1/8的概率，直接变成了(111)binary，即一步完成
第二种情况，3/8的概率，有两位变为1，即(011)binary或(101)binary或(110)binary,这时是E(1)时的情况，期望值为1+E(1)。
第三种情况，3/8的概率，有一位变为1，即(001)binary或(010)binary或(100)binary,这时是E(2)时的情况，期望值为1+E(2)。
第四种情况，仍为0，这时是1+E(3)。
列出式子，E(3) = 1/8 + 3/8 * (1+ E(1)) + 3/8 * (1 + E(2)) + 1/8 * (1 + E(3))
7/8 * E(3) = 1/8 + 3/8 * (1+ E(1)) + 3/8 * (1 + E(2)) + 1/8
7/8 * E(3) = 1/8 + 3/8 * 3 + 3/8 * 11/3 + 1/8
E(3) = 22/7

第五步，现在可以总结出最一般的规律：

E(N) = 1/2^N + C(N,1)/2^N * (1+E(1)) + C(N,2)/2^N * (1+E(2)) + ... + C(N,N-1)/2^N * (1+E(N-1)) + 1/2^N (1+E(N))
这里的C(N, i)是排列组合数，从N个数里取i个数的组合数。
耐心地把右侧的E(N)挪到左侧：
(2^N-1)/2^N * E(N) = 1/2^N + C(N,1)/2^N * (1+E(1)) + C(N,2)/2^N * (1+E(2)) + ... + C(N,N-1)/2^N * (1+E(N-1)) + 1/2^N
最后可以得到：
E(N) = [1 + C(N,1) * (1+E(1)) + C(N,2) * (1+E(2)) + ... + C(N,N-1) * (1+E(N-1)) + 1] / (2^N-1)
我们知道E(1)=2，不断套用上面的公式，就可以得到最终的答案。
求组合数的函数：

fn comb(n: usize, r: usize) -> f64 {
    let mut c = 1_f64;
    for i in 0..r {
        c *= ((n - i) as f64) / (r - i) as f64;
    }
    c
}

主程序：

let mut e = vec![0.0, 2.0];
for n in 2..=32 {
    let mut f = 1.0 as f64;
    for r in 1..=n {
        f += comb(n, r) * (1.0 + e[n - r]);
    }
    f = f / (2_f64.powf(n as f64) - 1.0);
    e.push(f);
    println!("E({}) = {:.10}", n, f);
}