1

この質問は、私がほとんどこの質問の下に書いたコメントに由来します。ここで、ザックは大きな数を法とする大きな数の階乗を計算しています (この質問のために素数であると仮定します)。ザックは階乗の従来の計算を使用しており、乗算ごとに剰余を取ります。

検討すべき代替手段はモンゴメリー乗算であるとほとんどコメントしましたが、それについてもっと考えてみると、この手法が同じ被乗数による複数の乗算を高速化するために使用されているのを見ただけです (特に、n mod p の計算を高速化するため)。

私の質問は、モンゴメリー乗算を使用して n の計算を高速化できるかどうかです! 大きな n と p の mod p?

4

1 に答える 1

3

素朴に、いいえ。積の n 項のそれぞれを「モンゴメリー空間」に変換する必要があるため、「通常の」アルゴリズムと同じように、m を法とする n 個の完全削減があります。

ただし、階乗は単なる n 項の任意の積ではありません。それははるかに構造化されています。特に、「Montgomerized」をすでにお持ちの場合は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;
}

マイナーな非効率性を修正するために、いくつかのインライン asm を使用した通常の実装に対して:

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

私の Haswell ラップトップでは、かなり大きな x に対して次のような結果が得られました。

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

ですから、これは本当に素晴らしい勝利のようです。Montgomery 実装の codegen はかなりまともですが、手書きのアセンブリでもさらに改善される可能性があります。

これは、「適度な」x と m の興味深いアプローチです。x が大きくなると、x の線形以下の複雑さを持つさまざまなアプローチが必然的に勝ちます。factorial には非常に多くの構造があるため、この方法では利用できません。

于 2014-08-23T14:04:56.070 に答える