1

这个问题源于我几乎写在这个问题下面的评论，其中 Zack 正在计算大量模数的阶乘（为了这个问题，我们将假设它是素数）。Zack 使用传统的阶乘计算，在每次乘法时取余数。

我几乎评论说要考虑的替代方案是Montgomery multiplication，但仔细想想，我只看到这种技术用于加速同一个被乘数的多次乘法（特别是加快n mod p 的计算）。

我的问题是：蒙哥马利乘法可以用来加速n的计算！大 n 和 p 的 mod p？

4

1 回答 1

3

天真地，不；您需要将乘积的 n 项中的每一项转换为“蒙哥马利空间”，因此您有 n 个完全约简 mod m，与“通常”算法相同。

然而，阶乘不仅仅是 n 项的任意乘积；它更有条理。特别是，如果你已经有了“蒙哥马化” kr mod m，那么你可以用一个非常便宜的减价来获得(k+1)r mod m。

所以这是完全可行的，虽然我以前没见过。我继续写了一个快速而肮脏的实现（非常未经测试，我根本不会相信它）：

// returns m^-1 mod 2**64 via clever 2-adic arithmetic (http://arxiv.org/pdf/1209.6626.pdf)
uint64_t inverse(uint64_t m) {
    assert(m % 2 == 1);
    uint64_t minv = 2 - m;
    uint64_t m_1 = m - 1;
    for (int i=1; i<6; i+=1) { m_1 *= m_1; minv *= (1 + m_1); }
    return minv;
}

uint64_t montgomery_reduce(__uint128_t x, uint64_t minv, uint64_t m) {
    return x + (__uint128_t)((uint64_t)x*-minv)*m >> 64;
}

uint64_t montgomery_multiply(uint64_t x, uint64_t y, uint64_t minv, uint64_t m) {
    return montgomery_reduce(full_product(x, y), minv, m);
}

uint64_t montgomery_factorial(uint64_t x, uint64_t m) {
    assert(x < m && m % 2 == 1);
    uint64_t minv = inverse(m); // m^-1 mod 2**64
    uint64_t r_mod_m = -m % m;  // 2**64 mod m
    uint64_t mont_term = r_mod_m;
    uint64_t mont_result = r_mod_m;
    for (uint64_t k=2; k<=x; k++) {
        // Compute the montgomerized product term: kr mod m = (k-1)r + r mod m.
        mont_term += r_mod_m;
        if (mont_term >= m) mont_term -= m;
        // Update the result by multiplying in the new term.
        mont_result = montgomery_multiply(mont_result, mont_term, minv, m);
    }
    // Final reduction
    return montgomery_reduce(mont_result, minv, m);
}

并将其与通常的实现进行基准测试：

__uint128_t full_product(uint64_t x, uint64_t y) {
    return (__uint128_t)x*y;
}

uint64_t naive_factorial(uint64_t x, uint64_t m) {
    assert(x < m);
    uint64_t result = x ? x : 1;
    while (x --> 2) result = full_product(result,x) % m;
    return result;
}

并针对使用一些内联汇编的通常实现来解决轻微的低效率问题：

uint64_t x86_asm_factorial(uint64_t x, uint64_t m) {
    assert(x < m);
    uint64_t result = x ? x : 1;
    while (x --> 2) {
        __asm__("mov %[result], %%rax; mul %[x]; div %[m]"
                : [result] "+d" (result) : [x] "r" (x), [m] "r" (m) : "%rax", "flags");
    }
    return result;
}

对于相当大的 x，在我的 Haswell 笔记本电脑上的结果如下：

implementation   speedup
---------------------------
naive            1.00x
x86_asm          1.76x
montgomery       5.68x

所以这看起来确实是一场不错的胜利。蒙哥马利实现的代码生成相当不错，但也可以通过手写汇编进一步改进。

对于“适度”的 x 和 m，这是一种有趣的方法。一旦 x 变大，在 x 中具有次线性复杂度的各种方法必然会胜出；阶乘具有如此多的结构，这种方法没有利用。

于 2014-08-23T14:04:56.070 回答