python - 如何检查数字是否可以表示素数（第n个根是否为素数）

Question

我正在尝试这个问题一段时间，但一次又一次地得到错误的答案。数字可能非常大 <=2^2014。 22086. 主功率测试

关于我的算法的解释：

对于给定的数字，我正在检查该数字是否可以表示为素数的形式。
因此，检查主功率的最大限制是 log n base 2。
最后问题简化为找到一个数字的第 n 个根，如果它是素数，我们有我们的答案，否则检查所有i直到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 的《计算代数数论课程》一书中作为算法 1.7.5 给出。您可以在以下位置查看正在运行的程序 http://ideone.com/cNzQYr。

score 1 · Accepted Answer

这不是一个真正的答案，但我没有足够的空间来写它作为评论。

所以，如果问题还没有解决，你可以试试下面的 nth_rootp 函数，虽然它有点难看（它只是一个二进制搜索来找到函数的精确值）：

def nth_rootp(n,k):
    r = int(round(math.log(n,2)/k))
    left = 2**(r-1)
    right = 2**(r+1)
    if left**k == n:
        return left
    if right**k == n:
        return right
    while left**k < n and right**k > n:
        tmp = (left + right)/2
        if tmp**k == n:
            return tmp
        if tmp == left or tmp == right:
            return tmp
        if tmp**k < n:
            left = tmp
        else:
            if tmp**k > n:
                right = tmp

score 1 · Accepted Answer

您的代码对于此任务看起来有点过于复杂，我不会费心检查它，但您需要的是以下内容

is_prime, 自然
主发生器，可选
以精确的方式计算数字的 n 次方根

对于第一个，我建议使用Miller-Rabin 测试的确定性形式，并使用一组适当的见证来保证准确的结果，直到 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 测试，我只会提一下，如果您有兴趣，可以在此处找到确定性版本的实现，只需更新见证人列表即可。

对于sieve，只需更改您正在使用的返回一个带有质数的列表，例如[ p for p,is_p in enumerate(sieve) if is_p ]

有了这些，剩下的唯一一件事就是计算数字的第 n 个根，并以精确的方式做到这一点，我们需要摆脱只会产生头痛的讨厌的浮点运算，答案是实现第N 个仅使用整数算术的根算法，这与您已经使用的算法非常相似，我使用Mark Dickinson为立方根isqrt制作的算法来指导自己并对其进行概括，我得到了

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

有了以上所有，您就可以解决问题而不会带来不便

score 0 · Accepted Answer

from sympy.ntheory import factorint

q=int(input("Give me the number q=")) 

fact=factorint(q)        #We factor the number q=p_1^{n_1}*p_2^{n_2}*...

p_1=list(fact.keys())    #We create a list from keys to be the the numbers p_1,p_2,...

n_1=list(fact.values())  #We create a list from values to be the the numbers n_1,n_2,...

p=int(p_1[0])            

n=int(n_1[0])            

if q!=p**n:              #Check if the number q=p_{1}[0]**n_{1}[0]=p**n.

    print("The number "+str(q)+" is not a prime power")

else:

    print("The number "+str(q)+" is a prime power")

    print("The prime number p="+str(p))

    print("The natural number n="+str(n))

python - 如何检查数字是否可以表示素数（第n个根是否为素数）

4 回答 4

Related

Reference