python - 数が素数べき乗で表せるか調べる方法（n乗根が素数か否か）

Question

私はしばらくこの問題を試していますが、何度も間違った答えを得ています。数が非常に大きくなる可能性があります <=2^2014. 22086. プライムパワーテスト

私のアルゴリズムについての説明:

与えられた数値について、その数値を素数べき乗の形式で表現できるかどうかを確認しています。
したがって、素数べき乗をチェックする最大制限は log n base 2 です。
i最後に、問題は数値の n乗根log (n base 2)を見つけることになり、それが素数の場合は答えが得られますexit。
私はあらゆる種類の最適化を使用し、膨大なテストケースをテストしましたが、すべてのアルゴリズムで正しい答えが得られました
しかし、ジャッジは間違った答えを言います。
Spoj には、小さな制約 n<=10^18 に関する別の同様の問題があり、Python と C++ (C++ で最高のソルバー) で既に受け入れられています。

これが私のpythonコードです。何か間違ったことをしている場合は、私に提案してください。私はPythonにあまり精通していないので、私のアルゴリズムは少し長いです。前もって感謝します。

私のアルゴリズム:

import math
import sys
import fractions
import random
import decimal
write = sys.stdout.write
def sieve(n):
    sqrtn = int(n**0.5)
    sieve = [True] * (n+1)
    sieve[0] = False
    sieve[1] = False
    for i in range(2, sqrtn+1):
        if sieve[i]:
            m = n//i - i
            sieve[i*i:n+1:i] = [False] * (m+1)
    return sieve
def gcd(a, b):
    while b:
        a, b = b, a%b
    return a
def mr_pass(a, s, d, n):
    a_to_power = pow(a, d, n)
    if a_to_power == 1:
        return True
    for i in range(s-1):
        if a_to_power == n - 1:
            return True
        a_to_power = (a_to_power * a_to_power) % n
    return a_to_power == n - 1
isprime=sieve(1000000)
sprime= [2,3,5,7,11,13,17,19,23,29,31,37,41,43,47,53,59,61,67,71,73,79,83,89,97,101,103,107,109,113,127,131,137,139,149,151,157,163,167,173,179,181,191,193,197,199,211,223,227,229,233,239,241,251,257,263,269,271,277,281,283,293,307,311,313,317,331,337,347,349,353,359,367,373,379,383,389,397,401,409,419,421,431,433,439,443,449,457,461,463,467,479,487,491,499,503,509,521,523,541,547,557,563,569,571,577,587,593,599,601,607,613,617,619,631,641,643,647,653,659,661,673,677,683,691,701,709,719,727,733,739,743,751,757,761,769,773,787,797,809,811,821,823,827,829,839,853,857,859,863,877,881,883,887,907,911,919,929,937,941,947,953,967,971,977,983,991,997]
def smooth_num(n):
    c=0
    for a in sprime:
        if(n%a==0):
            c+=1
        if(c>=2):
            return True;
    return False
def is_prime(n):
    if(n<1000000):
        return isprime[n]
    if any((n % p) == 0 for p in sprime):
        return False
    if n==2:
        return True
    d = n - 1
    s = 0
    while d % 2 == 0:
        d >>= 1
        s += 1
    for repeat in range(10):
        a=random.randint(1,n-1)
        if not mr_pass(a, s, d, n):
            return False
    return True
def iroot(n,k):
    hi = 1
    while pow(hi, k) < n:
        hi *= 2
    lo = hi // 2
    while hi - lo > 1:
        mid = (lo + hi) // 2
        midToK = (mid**k)
        if midToK < n:
            lo = mid
        elif n < midToK:
            hi = mid
        else:
            return mid
    if (hi**k) == n:
        return hi
    else:
        return lo
def isqrt(x):
    n = int(x)
    if n == 0:
        return 0
    a, b = divmod(n.bit_length(), 2)
    x = pow(2,(a+b))
    while True:
        y = (x + n//x)>>1
        if y >= x:
            return x
        x = y
maxx=2**1024;minn=2**64
def nth_rootp(n,k):
    return int(round(math.exp(math.log(n)/k),0))
def main():
    for cs in range(int(input())):
        n=int(sys.stdin.readline().strip())
        if(smooth_num(n)):
            write("Invalid order\n")
            continue;
        order = 0;m=0
        power =int(math.log(n,2))
        for i in range(1,power+1):
            if(n<=maxx):
                if i==1:m=n
                elif(i==2):m=isqrt(n)
                elif(i==4):m=isqrt(isqrt(n))
                elif(i==8):m=isqrt(isqrt(isqrt(n)))
                elif(i==16):m=isqrt(isqrt(isqrt(isqrt(n))))
                elif(i==32):m=isqrt(isqrt(isqrt(isqrt(isqrt(n)))))
                elif(i==64):m=isqrt(isqrt(isqrt(isqrt(isqrt(isqrt(n))))))
                elif(i==128):m=isqrt(isqrt(isqrt(isqrt(isqrt(isqrt(isqrt(n)))))))
                elif(i==256):m=isqrt(isqrt(isqrt(isqrt(isqrt(isqrt(isqrt(isqrt(n))))))))
                else:m=int(nth_rootp(n,i))
            else:
                if i==1:m=n
                elif i==2:m=isqrt(n)
                elif(i==4):m=isqrt(isqrt(n))
                elif(i==8):m=isqrt(isqrt(isqrt(n)))
                elif(i==16):m=isqrt(isqrt(isqrt(isqrt(n))))
                elif(i==32):m=isqrt(isqrt(isqrt(isqrt(isqrt(n)))))
                elif(i==64):m=isqrt(isqrt(isqrt(isqrt(isqrt(isqrt(n))))))
                elif(i==128):m=isqrt(isqrt(isqrt(isqrt(isqrt(isqrt(isqrt(n)))))))
                elif(i==256):m=isqrt(isqrt(isqrt(isqrt(isqrt(isqrt(isqrt(isqrt(n))))))))
                else:m=iroot(n,i)
            if m<2:
                order=0
                break
            if(is_prime(m) and n==(m**i)):
                write("%d %d\n"%(m,i))
                order = 1
                break
        if(order==0):
            write("Invalid order\n")
main()

score 4 · Accepted Answer

問題は浮動小数点の不正確さにあると思いますが、すべてのコードを読むつもりはありません。数値nが素数ベキかどうかを判断するプログラムは次のとおりです。素数pと累乗kを返します。

# prime power predicate

from random import randint
from fractions import gcd

def findWitness(n, k=5): # miller-rabin
    s, d = 0, n-1
    while d % 2 == 0:
        s, d = s+1, d/2
    for i in range(k):
        a = randint(2, n-1)
        x = pow(a, d, n)
        if x == 1 or x == n-1: continue
        for r in range(1, s):
            x = (x * x) % n
            if x == 1: return a
            if x == n-1: break
        else: return a
    return 0

# returns p,k such that n=p**k, or 0,0
# assumes n is an integer greater than 1
def primePower(n):
    def checkP(n, p):
        k = 0
        while n > 1 and n % p == 0:
            n, k = n / p, k + 1
        if n == 1: return p, k
        else: return 0, 0
    if n % 2 == 0: return checkP(n, 2)
    q = n
    while True:
        a = findWitness(q)
        if a == 0: return checkP(n, q)
        d = gcd(pow(a,q,n)-a, q)
        if d == 1 or d == q: return 0, 0
        q = d

このプログラムはフェルマーの小定理を使用し、Miller-Rabin アルゴリズムによって検出されたnの複合性に対する証人aを利用します。これは、Henri Cohen の著書A Course in Computational Algebraic Number Theory でAlgorithm 1.7.5 として提供されています。プログラムの動作は次の URL で確認できます。http://ideone.com/cNzQYr .

score 1 · Accepted Answer

あなたのコードは、このタスクに対して少し複雑すぎるように見えます。わざわざチェックするつもりはありませんが、必要なものは次のとおりです。

is_prime、当然
オプションの素数ジェネレータ
正確な方法で数値の n 乗根を計算する

最初のものについては、1543267864443420616877677640751301 (1.543 x 10 ³³ )までの正確な結果を保証するために、適切な証人のセットを使用した決定論的な形式のミラー-ラビンテストをお勧めします。あなたの基準で選ばれた証人

ソリューションのテンプレートは次のとおりです。

import math

def is_prime(n):
    ...

def sieve(n):
    "list of all primes p such that p<n"
    ...

def inthroot(x,n):
    "calculate floor(x**(1/n))"
    ...

def is_a_power(n):
    "return (a,b) if n=a**b otherwise throw ValueError"
    for b in sieve( math.log2(n) +1 ):
        a = inthroot(n,b)
        if a**b == n:
            return a,b
    raise ValueError("is not a power")

def smooth_factorization(n):
    "return (p,e) where p is prime and n = p**e if such value exists, otherwise throw ValueError"
    e=1
    p=n
    while True:
        try:
            p,n = is_a_power(p)
            e = e*n
        except ValueError:
            break
    if is_prime(p): 
        return p,e
    raise ValueError

def main():
    for test in range( int(input()) ):
        try:
            p,e = smooth_factorization( int(input()) )
            print(p,e)
        except ValueError:
            print("Invalid order")

main()

そして、上記のコードは自明であるはずです

黒を埋める

Miller-Rabin テストに精通しているので、興味がある場合は、ここで決定論バージョンの実装を見つけることができ、 witnessのリストを更新するだけで準備完了です。

の場合は、使用しているものを変更して、たとえば次のsieveような素数のリストを返すだけです[ p for p,is_p in enumerate(sieve) if is_p ]

それらが邪魔にならないようにして、残っているのは数のn乗根を計算することだけです。それを正確に行うには、頭痛の種になるだけの厄介な浮動小数点演算を取り除く必要があり、答えはN乗根を実装することです整数演算のみを使用するルートアルゴリズムisqrtは、既に使用しているものと非常によく似ています。Mark Dickinsonが作成したキューブルートを一般化すると、次のようになります。

def inthroot(A, n) :
    "calculate floor( A**(1/n) )"
    #https://en.wikipedia.org/wiki/Nth_root_algorithm
    #https://en.wikipedia.org/wiki/Nth_root#nth_root_algorithm
    #https://stackoverflow.com/questions/35254566/wrong-answer-in-spoj-cubert/35276426#35276426
    #https://stackoverflow.com/questions/39560902/imprecise-results-of-logarithm-and-power-functions-in-python/39561633#39561633
    if A<0:
        if n%2 == 0:
            raise ValueError
        return - inthroot(-A,n)
    if A==0:
        return 0
    n1 = n-1
    if A.bit_length() < 1024: # float(n) safe from overflow
        xk = int( round( pow(A,1.0/n) ) )
        xk = ( n1*xk + A//pow(xk,n1) )//n # Ensure xk >= floor(nthroot(A)).
    else:
        xk = 1 << -(-A.bit_length()//n) # 1 << sum(divmod(A.bit_length(),n))
                                        # power of 2 closer but greater than the nth root of A
    while True:
        sig = A // pow(xk,n1)
        if xk <= sig:
            return xk
        xk = ( n1*xk + sig )//n

上記のすべてで、不便なく問題を解決できます

python - 数が素数べき乗で表せるか調べる方法（n乗根が素数か否か）

4 に答える 4

Related

Reference