1

Lucas 疑似素数検定を使用して、数値nが素数か合成数かを判断する関数を作成しようとしています。現時点では、標準テストで作業していますが、それが機能するようになったら、強力なテストを作成します。私はBaillie と Wagstaff による論文を読んでおり、 trn.cファイルで Thomas Nicely による実装に従っています。

完全なテストにはいくつかのステップが含まれることを理解しています: 小さな素数による割り算の試行、nが平方でないことの確認、基数 2 に対する強力な疑似素数テストの実行、そして最後にルーカスの疑似素数テスト。他のすべての部分は処理できますが、Lucas 疑似素数検定に問題があります。Pythonでの私の実装は次のとおりです。

def gcd(a, b):
    while b != 0:
        a, b = b, a % b
    return a

def jacobi(a, m):
    a = a % m; t = 1
    while a != 0:
        while a % 2 == 0:
            a = a / 2
            if m % 8 == 3 or m % 8 == 5:
                t = -1 * t
        a, m = m, a # swap a and m
        if a % 4 == 3 and m % 4 == 3:
            t = -1 * t
        a = a % m
    if m == 1:
        return t
    return 0

def isLucasPrime(n):
    dAbs, sign, d = 5, 1, 5
    while 1:
        if 1 < gcd(d, n) > n:
            return False
        if jacobi(d, n) == -1:
            break
        dAbs, sign = dAbs + 2, sign * -1
        d = dAbs * sign
    p, q = 1, (1 - d) / 4
    print "p, q, d =", p, q, d
    u, v, u2, v2, q, q2 = 0, 2, 1, p, q, 2 * q
    bits = []
    t = (n + 1) / 2
    while t > 0:
        bits.append(t % 2)
        t = t // 2
    h = -1
    while -1 * len(bits) <= h:
        print "u, u2, v, v2, q, q2, bits, bits[h] = ",\
               u, u2, v, v2, q, q2, bits, bits[h]
        u2 = (u2 * v2) % n
        v2 = (v2 * v2 - q2) % n
        if bits[h] == 1:
            u = u2 * v + u * v2
            u = u if u % 2 == 0 else u + n
            u = (u / 2) % n
            v = (v2 * v) + (u2 * u * d)
            v = v if v % 2 == 0 else v + n
            v = (v / 2) % n
        if -1 * len(bits) < h:
            q = (q * q) % n
            q2 = q + q
        h = h - 1
    return u == 0

これを実行すると、83 や 89 などの素数isLucasPrimeが返されますが、これは正しくありません。Falseまた、コンポジット 111 に対しても返されますFalseが、これは正しいです。そして、複合 323 を返します。これは、を返す必要Falseがあるルーカス疑似素数であることがわかっています。実際、私がテストしたすべてのnを返します。isLucasPrimeTrueisLucasPseudoprimeFalse

いくつか質問があります:

1) 私は C/GMP の専門家ではありませんが(n+1)/2、他の著者がビットを左から右に実行するのに対し、Nicely は右から左 (最下位から最上位) のビットを実行しているように思えます。上記のコードは左から右にビットを実行しますが、右から左にビットを実行しようとしても同じ結果になりました。正しい順番は？

2) Nicelyが 1 ビットのu変数とv変数のみを更新するのは奇妙に思えます。これは正しいです？チェーンのインデックスは各ステップで増加するため、ループのたびに 4 つのルーカス チェーン変数すべてを更新することを期待していました。

3) 私は何を間違えましたか?

4

1 に答える 1

3

1) 私は C/GMP の専門家ではありませんが(n+1)/2、他の著者がビットを左から右に実行するのに対し、Nicely は右から左 (最下位から最上位) のビットを実行しているように思えます。上記のコードは左から右にビットを実行しますが、右から左にビットを実行しようとしても同じ結果になりました。正しい順番は？

確かに、ナイスリーは最下位ビットから最上位ビットに移動します。U(2^k)彼はand V(2^k)(もちろんQ^(2^k)、すべてモジュロ) をand変数で計算し、 respをandに格納します。1 ビットが検出されると、残りが変化し、それに応じて更新されます。NmpzU2mmpzV2mU((N+1) % 2^k)V((N+1) % 2^k)mpzUmpzV(N+1) % 2^kmpzUmpzV

もう 1 つの方法はU(p)、 、U(p+1)、V(p)および (オプションで)V(p+1)の接頭辞pを計算しN+1、それらを組み合わせて計算U(2*p+1)し、接頭辞の後の次のビットが 0 か 1U(2*p)かに応じてU(2*p+2)[ の場合は同上] を計算することです。Vp

どちらの方法も正しいです。たとえば、x^N左から右へ、x^pおよびx^(p+1)状態として、または右から左へ、x^(2^k)状態x^(N % 2^k)として、および [および、計算U(n)およびU(n+1)基本的ζ^nにどこを計算するζ = (1 + sqrt(D))/2] を使用して電力を計算できます。

私も他の人も、どうやら左から右への順序の方が簡単だと思います。私は分析を行ったり読んだりしていません。右から左の方が平均して計算コストが低く、そのために右から左を選択した可能性があります。

2) Nicelyが 1 ビットの変数uと変数のみを更新するのは奇妙に思えます。vこれは正しいです？チェーンのインデックスは各ステップで増加するため、ループのたびに 4 つのルーカス チェーン変数すべてを更新することを期待していました。

はい、正解です。ビットが 0(N+1) % 2^k == (N+1) % 2^(k-1)の場合は余りです。2^k

3) 私は何を間違えましたか?

最初の小さなタイプミス:

if 1 < gcd(d, n) > n:

する必要があります

if 1 < gcd(d, n) < n:

もちろん。

もっと実質的には、Nicely の走査順序 (右から左) の更新を使用しますが、反対方向に走査します。もちろん、それは間違った結果をもたらします。

また、更新時にはv

    if bits[h] == 1:
        u = u2 * v + u * v2
        u = u if u % 2 == 0 else u + n
        u = (u / 2) % n
        v = (v2 * v) + (u2 * u * d)
        v = v if v % 2 == 0 else v + n
        v = (v / 2) % n

の新しい値をu使用しますが、古い値を使用する必要があります。

def isLucasPrime(n):
    dAbs, sign, d = 5, 1, 5
    while 1:
        if 1 < gcd(d, n) < n:
            return False
        if jacobi(d, n) == -1:
            break
        dAbs, sign = dAbs + 2, sign * -1
        d = dAbs * sign
    p, q = 1, (1 - d) // 4
    u, v, u2, v2, q, q2 = 0, 2, 1, p, q, 2 * q
    bits = []
    t = (n + 1) // 2
    while t > 0:
        bits.append(t % 2)
        t = t // 2
    h = 0
    while h < len(bits):
        u2 = (u2 * v2) % n
        v2 = (v2 * v2 - q2) % n
        if bits[h] == 1:
            uold = u
            u = u2 * v + u * v2
            u = u if u % 2 == 0 else u + n
            u = (u // 2) % n
            v = (v2 * v) + (u2 * uold * d)
            v = v if v % 2 == 0 else v + n
            v = (v // 2) % n
        if h < len(bits) - 1:
            q = (q * q) % n
            q2 = q + q
        h = h + 1
    return u == 0

動作します（保証はありませんが、正しいと思います。いくつかのテストを行ったところ、すべて合格しました）。

于 2013-02-24T17:08:00.690 に答える