如何选择位运算序列生成指定比特置位概率的64位长整数？

1. 公式简化与问题转化 #

你提出的目标概率表达式可以简化为：
p ≈ 4⁻ᵏ - 4^(-(k+n))
其中k是AND运算次数，n是OR运算次数，二者均为非负整数。令m = k + n，则表达式也可写成p ≈ 4⁻ᵏ - 4⁻ᵐ（要求m ≥ k）。我们的核心任务就是找到整数m ≥ k ≥ 0，让该式与输入p的误差最小。

2. 筛选候选k值 #

由于4⁻ᵏ随k增大指数递减，我们只需在有限范围内筛选k的候选值：

• 计算基准值 k0 = floor( -ln(p) / (2*ln(2)) )，这是解方程4⁻ᵏ = p得到的近似整数解。
• 候选k取k0-1, k0, k0+1（若k0-1为负数则只保留非负候选），偏离这个范围的k会导致4⁻ᵏ与p的量级差距过大，误差必然更高。

3. 对每个候选k匹配最优m #

针对每个候选k，分两种情况处理：

• 情况1：4⁻ᵏ ≤ p：此时4⁻ᵏ -4⁻ᵐ始终小于p，且m越大，4⁻ᵐ越小，结果越接近4⁻ᵏ。直接取足够大的m（比如20，此时4⁻²⁰≈9e-13，对double精度可忽略误差），对应的n = m -k。
• 情况2：4⁻ᵏ > p：解方程4⁻ᵐ = 4⁻ᵏ - p得到m0 = ceil( -ln(4⁻ᵏ - p)/(2*ln(2)) )，取floor(m0)和ceil(m0)作为m的候选值（需保证m ≥k），计算两个候选对应的误差，选更小的那个。

4. 确定全局最优解 #

遍历所有(k,m)候选组合，计算每个组合的误差|p - (4⁻ᵏ -4⁻ᵐ)|，最终选择误差最小的一组，对应的n = m -k就是所需的OR运算次数。

示例演示 #

以输入p=0.6为例：

• 计算k0 = floor( -ln(0.6)/(2*ln2) )≈0，候选k为0、1。
  • k=0时：4⁻⁰=1，target=1-0.6=0.4，m0≈ceil(0.661)=1。候选m=1时，结果为1-0.25=0.75，误差0.15；m=0时结果为0，误差0.6。
  • k=1时：4⁻¹=0.25 <0.6，取m=20，结果≈0.25，误差0.35。
• 最优组合为k=0, m=1，对应n=1，误差0.15。

伪代码实现 #

import math

def find_optimal_kn(p):
    min_error = float('inf')
    best_k = 0
    best_n = 0

    if p <= 1e-16:
        return (0, 0)
    if p >= 1 - 1e-16:
        return (0, 20)

    # 生成k的候选值
    k0 = math.floor(-math.log(p) / (2 * math.log(2)))
    k_candidates = set()
    for d in (-1, 0, 1):
        k = k0 + d
        if k >= 0:
            k_candidates.add(k)

    for k in k_candidates:
        four_neg_k = 4 ** (-k)
        target = four_neg_k - p

        if target <= 0:
            # 取足够大的m，忽略4^-m的影响
            m = 20
            current_val = four_neg_k - (4 ** (-m))
            error = abs(p - current_val)
            if error < min_error:
                min_error = error
                best_k = k
                best_n = m - k
        else:
            m0 = math.log(1 / target) / (2 * math.log(2))
            m_candidates = [math.floor(m0), math.ceil(m0)]
            for m in m_candidates:
                if m < k:
                    m = k
                current_val = four_neg_k - (4 ** (-m))
                error = abs(p - current_val)
                if error < min_error:
                    min_error = error
                    best_k = k
                    best_n = m - k

    # 单独处理p接近1的情况
    if p > 0.9:
        m = math.ceil(-math.log(1 - p) / (2 * math.log(2)))
        current_val = 1 - (4 ** (-m))
        error = abs(p - current_val)
        if error < min_error:
            best_k = 0
            best_n = m

    return (best_k, best_n)

内容的提问来源于stack exchange，提问作者Andrew Kornder