Модуль быстрой простой факторизации

Я ищу алгоритм реализации или очистки для получения простой факторизации N в любом питоне, псевдокоде или прочем хорошо читаемом. Есть несколько требований/фактов:

  • N находится между 1 и ~ 20 цифрами
  • Нет предварительно вычисленной таблицы поиска, memoization в порядке.
  • Не нужно быть математически доказанным (например, при необходимости может полагаться на гипотезу Голдбаха)
  • Не обязательно быть точным, разрешено быть вероятностным/детерминированным, если необходимо

Мне нужен алгоритм быстрой простой факторизации не только для себя, но и для использования во многих других алгоритмах, таких как вычисление Euler phi (n).

Я пробовал другие алгоритмы из Википедии и такие, но я не мог их понять (ECM), или я не мог создать рабочую реализацию из алгоритма (Pollard-Brent).

Мне действительно интересен алгоритм Поллард-Брент, поэтому любая информация/реализации на нем была бы действительно приятной.

Спасибо!

ИЗМЕНИТЬ

После небольшого общения я создал довольно быстрый модуль prime/factorization. Он сочетает в себе оптимизированный алгоритм пробного деления, алгоритм Полларда-Брент, тест примитивности мельника-рабина и самый быстрый в прямом эфире, который я нашел в Интернете. gcd - регулярная реализация GCD Евклида (двоичный Euclid GCD значительно медленнее, чем обычный).

Bounty

О, радость, можно получить щедрость! Но как я могу его выиграть?

  • Найдите в моем модуле оптимизацию или ошибку.
  • Предоставить альтернативные/лучшие алгоритмы/реализации.

Ответ, который является наиболее полным/конструктивным, получает награду.

И, наконец, сам модуль:

import random

def primesbelow(N):
    # http://stackoverflow.com/info/2068372/fastest-way-to-list-all-primes-below-n-in-python/3035188#3035188
    #""" Input N>=6, Returns a list of primes, 2 <= p < N """
    correction = N % 6 > 1
    N = {0:N, 1:N-1, 2:N+4, 3:N+3, 4:N+2, 5:N+1}[N%6]
    sieve = [True] * (N // 3)
    sieve[0] = False
    for i in range(int(N ** .5) // 3 + 1):
        if sieve[i]:
            k = (3 * i + 1) | 1
            sieve[k*k // 3::2*k] = [False] * ((N//6 - (k*k)//6 - 1)//k + 1)
            sieve[(k*k + 4*k - 2*k*(i%2)) // 3::2*k] = [False] * ((N // 6 - (k*k + 4*k - 2*k*(i%2))//6 - 1) // k + 1)
    return [2, 3] + [(3 * i + 1) | 1 for i in range(1, N//3 - correction) if sieve[i]]

smallprimeset = set(primesbelow(100000))
_smallprimeset = 100000
def isprime(n, precision=7):
    # http://en.wikipedia.org/wiki/Miller-Rabin_primality_test#Algorithm_and_running_time
    if n < 1:
        raise ValueError("Out of bounds, first argument must be > 0")
    elif n <= 3:
        return n >= 2
    elif n % 2 == 0:
        return False
    elif n < _smallprimeset:
        return n in smallprimeset


    d = n - 1
    s = 0
    while d % 2 == 0:
        d //= 2
        s += 1

    for repeat in range(precision):
        a = random.randrange(2, n - 2)
        x = pow(a, d, n)

        if x == 1 or x == n - 1: continue

        for r in range(s - 1):
            x = pow(x, 2, n)
            if x == 1: return False
            if x == n - 1: break
        else: return False

    return True

# https://comeoncodeon.wordpress.com/2010/09/18/pollard-rho-brent-integer-factorization/
def pollard_brent(n):
    if n % 2 == 0: return 2
    if n % 3 == 0: return 3

    y, c, m = random.randint(1, n-1), random.randint(1, n-1), random.randint(1, n-1)
    g, r, q = 1, 1, 1
    while g == 1:
        x = y
        for i in range(r):
            y = (pow(y, 2, n) + c) % n

        k = 0
        while k < r and g==1:
            ys = y
            for i in range(min(m, r-k)):
                y = (pow(y, 2, n) + c) % n
                q = q * abs(x-y) % n
            g = gcd(q, n)
            k += m
        r *= 2
    if g == n:
        while True:
            ys = (pow(ys, 2, n) + c) % n
            g = gcd(abs(x - ys), n)
            if g > 1:
                break

    return g

smallprimes = primesbelow(1000) # might seem low, but 1000*1000 = 1000000, so this will fully factor every composite < 1000000
def primefactors(n, sort=False):
    factors = []

    for checker in smallprimes:
        while n % checker == 0:
            factors.append(checker)
            n //= checker
        if checker > n: break

    if n < 2: return factors

    while n > 1:
        if isprime(n):
            factors.append(n)
            break
        factor = pollard_brent(n) # trial division did not fully factor, switch to pollard-brent
        factors.extend(primefactors(factor)) # recurse to factor the not necessarily prime factor returned by pollard-brent
        n //= factor

    if sort: factors.sort()

    return factors

def factorization(n):
    factors = {}
    for p1 in primefactors(n):
        try:
            factors[p1] += 1
        except KeyError:
            factors[p1] = 1
    return factors

totients = {}
def totient(n):
    if n == 0: return 1

    try: return totients[n]
    except KeyError: pass

    tot = 1
    for p, exp in factorization(n).items():
        tot *= (p - 1)  *  p ** (exp - 1)

    totients[n] = tot
    return tot

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

def lcm(a, b):
    return abs((a // gcd(a, b)) * b)

Ответы

Ответ 2

Если вы не хотите изобретать колесо, используйте библиотеку sympy

pip install sympy

Используйте функцию sympy.ntheory.factorint

>>> from sympy.ntheory import factorint
>>> factorint(10**20+1)
{73: 1, 5964848081: 1, 1676321: 1, 137: 1}

Вы можете указать несколько очень больших чисел:

>>> factorint(10**100+1)
{401: 1, 5964848081: 1, 1676321: 1, 1601: 1, 1201: 1, 137: 1, 73: 1, 129694419029057750551385771184564274499075700947656757821537291527196801: 1}

Ответ 3

Нет необходимости вычислять smallprimes с помощью primesbelow, используйте smallprimeset для этого.

smallprimes = (2,) + tuple(n for n in xrange(3,1000,2) if n in smallprimeset)

Разделите primefactors на две функции для обработки smallprimes и других для pollard_brent, это может сэкономить пару итераций, поскольку все полномочия smallprimes будут разделены на n.

def primefactors(n, sort=False):
    factors = []

    limit = int(n ** .5) + 1
    for checker in smallprimes:
        print smallprimes[-1]
        if checker > limit: break
        while n % checker == 0:
            factors.append(checker)
            n //= checker


    if n < 2: return factors
    else : 
        factors.extend(bigfactors(n,sort))
        return factors

def bigfactors(n, sort = False):
    factors = []
    while n > 1:
        if isprime(n):
            factors.append(n)
            break
        factor = pollard_brent(n) 
        factors.extend(bigfactors(factor,sort)) # recurse to factor the not necessarily prime factor returned by pollard-brent
        n //= factor

    if sort: factors.sort()    
    return factors

Рассматривая проверенные результаты Померанса, Селфриджа и Вагстаффа и Йешке, вы можете уменьшить повторения в isprime, который использует тест примитивности Миллера-Рабина. Из Wiki.

  • если n < 1,373,653, достаточно проверить a = 2 и 3;
  • если n < 9,080,191, достаточно проверить a = 31 и 73;
  • если n < 4,759,123,141, достаточно проверить a = 2, 7 и 61;
  • если n < 2,152,302,898,747, достаточно проверить a = 2, 3, 5, 7 и 11;
  • если n < 3,474,749,660,383, достаточно проверить a = 2, 3, 5, 7, 11 и 13;
  • если n < 341,550,071,728,321, достаточно проверить a = 2, 3, 5, 7, 11, 13 и 17.

Изменить 1. Исправлен обратный вызов if-else для добавления bigfactors к факторам в primefactors.

Ответ 4

Даже на текущем, есть несколько пятен, которые нужно заметить.

  • Не выполняйте checker*checker каждый цикл, используйте s=ceil(sqrt(num)) и checher < s
  • checher должен плюс 2 каждый раз, игнорировать все четные числа, кроме 2
  • Используйте divmod вместо % и //

Ответ 5

Вероятно, вам нужно сделать какое-то основное обнаружение, которое вы можете посмотреть здесь, Быстрый алгоритм для нахождения простых чисел?

Вы должны прочитать весь блог, хотя есть несколько алгоритмов, которые он перечисляет для проверки правильности.

Ответ 6

Здесь есть библиотека python с набором тестов примитивности (в том числе неправильные для того, что не делать). Он называется pyprimes. Понятно, что стоит упомянуть цель потомства. Я не думаю, что он включает в себя алгоритмы, которые вы упомянули.

Ответ 7

Я просто столкнулся с ошибкой в ​​этом коде при факторизации числа 2**1427 * 31.

  File "buckets.py", line 48, in prettyprime
    factors = primefactors.primefactors(n, sort=True)
  File "/private/tmp/primefactors.py", line 83, in primefactors
    limit = int(n ** .5) + 1
OverflowError: long int too large to convert to float

Этот фрагмент кода:

limit = int(n ** .5) + 1
for checker in smallprimes:
    if checker > limit: break
    while n % checker == 0:
        factors.append(checker)
        n //= checker
        limit = int(n ** .5) + 1
        if checker > limit: break

следует переписать в

for checker in smallprimes:
    while n % checker == 0:
        factors.append(checker)
        n //= checker
    if checker > n: break

который, скорее всего, будет работать быстрее на реальных входах. Квадратный корень медленный - в основном эквивалент многих умножений -, smallprimes имеет только несколько десятков членов, и таким образом мы удаляем вычисление n ** .5 из жесткой внутренней петли, что, безусловно, полезно при факторизации чисел, таких как 2**1427. Просто нет причин для вычисления sqrt(2**1427), sqrt(2**1426), sqrt(2**1425) и т.д. И т.д., Когда все, о чем мы заботимся, это "превышает ли [квадрат] контролера n".

Переписанный код не генерирует исключений при больших цифрах, и примерно в два раза быстрее, чем timeit (на вводе образцов 2 и 2**718 * 31).

Также обратите внимание, что isprime(2) возвращает неверный результат, но это нормально, если мы не полагаемся на него. IMHO вы должны переписать ввод этой функции как

if n <= 3:
    return n >= 2
...

Ответ 8

Вы можете разложить факторизацию до предела, а затем использовать brent, чтобы получить более высокие коэффициенты.

from fractions import gcd
from random import randint

def brent(N):
   if N%2==0: return 2
   y,c,m = randint(1, N-1),randint(1, N-1),randint(1, N-1)
   g,r,q = 1,1,1
   while g==1:             
       x = y
       for i in range(r):
          y = ((y*y)%N+c)%N
       k = 0
       while (k<r and g==1):
          ys = y
          for i in range(min(m,r-k)):
             y = ((y*y)%N+c)%N
             q = q*(abs(x-y))%N
          g = gcd(q,N)
          k = k + m
       r = r*2
   if g==N:
       while True:
          ys = ((ys*ys)%N+c)%N
          g = gcd(abs(x-ys),N)
          if g>1:  break
   return g

def factorize(n1):
    if n1==0: return []
    if n1==1: return [1]
    n=n1
    b=[]
    p=0
    mx=1000000
    while n % 2 ==0 : b.append(2);n//=2
    while n % 3 ==0 : b.append(3);n//=3
    i=5
    inc=2
    while i <=mx:
       while n % i ==0 : b.append(i); n//=i
       i+=inc
       inc=6-inc
    while n>mx:
      p1=n
      while p1!=p:
          p=p1
          p1=brent(p)
      b.append(p1);n//=p1 
    if n!=1:b.append(n)   
    return sorted(b)

from functools import reduce
#n= 2**1427 * 31 #
n= 67898771  * 492574361 * 10000223 *305175781* 722222227*880949 *908909
li=factorize(n)
print (li)
print (n - reduce(lambda x,y :x*y ,li))