Как проверить, может ли число быть представлено первичной степенью (n-й корень является простым или нет)

Question

Я пытаюсь эту проблему некоторое время, но получаю неправильный ответ снова и снова. номер может быть очень большой < = 2^2014. 22086. Prime Power Test

Объяснение о моем алгоритме:

1. Для заданного числа я Проверяется число может быть представлено как форма основной мощности или нет.
2. Таким образом, максимальный предел для проверки основной мощности является лог п база 2.
3. Наконец задача сводится к нахождению п-й корень из числа, и если она первична у нас есть ответ еще проверить все i до log (n base 2) и exit ,
4. Я использовал всевозможные оптимизации и протестировал огромные тестовые сценарии, и для всего моего алгоритма дается правильный ответ
5. , но судья говорит неправильный ответ.
6. SPOJ есть еще одна аналогичная проблема с небольшими ограничениями п < = 10^18, для которых я уже получил принятый с Python и C++ (Best решателя в C++)

Вот мой питон код Пожалуйста, предложите мне, если я что-то не так, я не очень разбираюсь в 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() 

Answer 1

это не совсем ответ, но у меня нет достаточно места, чтобы записать его в качестве комментария.

Итак, если проблема все еще не решена, вы можете попробовать следующую функцию для 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 

Answer 2

Я не буду читать весь этот код, хотя я подозреваю, что проблема - неточность с плавающей запятой. Вот моя программа, чтобы определить, является ли число n главной силой; она возвращает простое р и мощность к:

# 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 

Программа использует Ферма малая теорема и использует свидетельство к составленности из п, который находится в алгоритме Миллера-Рабина , Он приводится как Алгоритм 1.7.5 в книге Анри Коэна Курс по вычислительной теории алгебраических чисел. Вы можете увидеть программу в действии по адресу http://ideone.com/cNzQYr.

Answer 3

ваш код выглядеть немного усложненной для решения этой задачи, я не буду беспокоить, чтобы проверить это, но то, что вам нужно, являются следующие

1. is_prime, естественно
2. простой генератор, дополнительный
3. вычислить п-й корень из числа в точном пути

для первого я рекомендую детерминированную форму Miller-Rabin test с подходящим Набор witness для задела на точный результат до 1543267864443420616877677640751301 (1,543 х 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() 

И выше код должен быть сам объяснительное

Заполнение черноту

Как вы знакомы с тестом Миллера-Рабина, я остановлюсь только на том, что если вас это интересует, вы можете найти реализацию детерминистской версии here, просто обновите список witness, и вы готовы к работе.

Для sieve, просто изменить тот, который вы используете, чтобы вернуть список с номером простых чисел, как это, например [ p for p,is_p in enumerate(sieve) if is_p ]

С теми из пути, единственное, что левая вычислить п-й корень из номер и сделать это точным образом, нам нужно получить отрывок из этой надоедливой арифметики с плавающей запятой, которая дает только головные боли, а ответ реализует Nth root algorithm, используя только целочисленную арифметику, которая очень похожа на ту, что у вас есть 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 

и все выше, вы можете решить эту проблему без неудобного