118

私はMiller-Rabin primality testを実装しようとしていましたが、なぜ中規模の数値 (~7 桁) にこれほど長い (> 20 秒) かかるのか戸惑いました。最終的に、次のコード行が問題の原因であることがわかりました。

x = a**d % n

(ここでa、 、d、およびnはすべて似ていますが、等しくない中規模の数値です。**は累乗演算子、%はモジュロ演算子です)

次に、次のように置き換えてみました。

x = pow(a, d, n)

それに比べて、それはほとんど瞬時です。

コンテキストについては、元の関数を次に示します。

from random import randint

def primalityTest(n, k):
    if n < 2:
        return False
    if n % 2 == 0:
        return False
    s = 0
    d = n - 1
    while d % 2 == 0:
        s += 1
        d >>= 1
    for i in range(k):
        rand = randint(2, n - 2)
        x = rand**d % n         # offending line
        if x == 1 or x == n - 1:
            continue
        for r in range(s):
            toReturn = True
            x = pow(x, 2, n)
            if x == 1:
                return False
            if x == n - 1:
                toReturn = False
                break
        if toReturn:
            return False
    return True

print(primalityTest(2700643,1))

時限計算の例:

from timeit import timeit

a = 2505626
d = 1520321
n = 2700643

def testA():
    print(a**d % n)

def testB():
    print(pow(a, d, n))

print("time: %(time)fs" % {"time":timeit("testA()", setup="from __main__ import testA", number=1)})
print("time: %(time)fs" % {"time":timeit("testB()", setup="from __main__ import testB", number=1)})

出力 (PyPy 1.9.0 で実行):

2642565
time: 23.785543s
2642565
time: 0.000030s

出力 (Python 3.3.0、2.7.2 で実行すると、非常に類似した時間が返されます):

2642565
time: 14.426975s
2642565
time: 0.000021s

また、関連する質問として、通常は PyPy の方がはるかに高速なのに、Python 2 または 3 で実行すると、PyPy で実行した場合よりもこの計算がほぼ 2 倍速くなるのはなぜですか?

4

4 に答える 4

171

累乗剰余に関するウィキペディアの記事を参照してください。基本的に、 を実行する場合a**d % n、実際には を計算する必要がありa**d、これは非常に大きくなる可能性があります。しかし、それ自体を計算a**d % nせずに計算する方法はあります。オペレーターは、モジュラスをすぐに取得しようとしていることを知るために「未来を見る」ことができないため、これを行うことはできません。a**dpow**

于 2013-01-03T06:03:13.917 に答える
37

BrenBarn が主な質問に答えました。余談:

通常は PyPy の方がはるかに高速なのに、Python 2 または 3 で実行すると PyPy よりもほぼ 2 倍高速なのはなぜですか?

PyPy のパフォーマンス ページを読んだ場合、これはまさに PyPy が得意としない種類のものです。実際、彼らが示した最初の例です。

悪い例には、最適化できないサポート コードによって実行される、大きな long を使用した計算が含まれます。

理論的には、mod が続く巨大なべき乗を (少なくとも最初のパスの後で) 剰余べき乗に変換することは、JIT が行うことができる変換ですが、PyPy の JIT ではできません。

余談ですが、巨大な整数で計算を行う必要がある場合は、 のようなサードパーティ製モジュールを検討することをお勧めします。これらのモジュールはgmpy、場合によっては主流の用途以外で CPython のネイティブ実装よりもはるかに高速であり、多くの機能を備えています。便利ではないという犠牲を払って、そうでなければ自分で書かなければならない追加機能の。

于 2013-01-03T06:17:16.927 に答える
13

a**(2i) mod n剰余累乗を実行するためのショートカットがあります。たとえば、 for every ifrom 1toを見つけて、必要な中間結果をlog(d)一緒に乗算 (mod ) することができます。n3-argument のような専用の累乗剰余関数pow()は、剰余算術を実行していることを認識しているため、このようなトリックを活用できます。Python パーサーは、与えられたそのままの式a**d % nではこれを認識できないため、完全な計算を実行します (これにはさらに時間がかかります)。

于 2013-01-03T06:07:30.083 に答える
3

計算方法x = a**d % nはa、dべき乗し、それを でモジュロすることnです。まず、aが大きい場合、これにより巨大な数が作成され、切り捨てられます。ただし、x = pow(a, d, n)ほとんどの場合、最後の桁のみが追跡されるように最適化されていnます。これは、数を法とする乗算を計算するために必要なすべてです。

于 2013-01-03T06:08:17.120 に答える