Найдите простые числа, используя теорему Уилсона

Следующий код является воплощением теоремы Вильсона. Он использует простой факториал и теорему о делении для проверки простоты. Моя проблема связана с производительностью кода, и есть две наиболее важные проблемы. Во-первых, время, за которое берется остаток, намного превышает остальную часть кода; Во-вторых, при добавлении 1 к числителю возникает заметный дефицит времени из-за того, что я считаю проблемой при обновлении ссылки на другой объект для очень больших целых чисел.

Есть ли способ исправить эти проблемы, кроме, конечно, ускорения этой программы путем кодирования в cpp?

n = 100000


def absolute(x: int):
    if (x == 2):
        a = Diction[x - 1] * x * (x + 1)
        Diction.pop(x - 1)
        Diction.update({(x + 1): a})
        return 2
    try:
      a = Diction[x-1] * x
      return a
    finally:
        a *= (x+1)
        Diction.pop(x-1)
        Diction.update({(x+1):a})

def main():
    for i in range(3, n+1,2):
        if (((absolute(i - 1) + 1 ) % i) == 0):
            prime.append(i)

if __name__ == '__main__':
    import time
    import cProfile
    start = time.perf_counter()

    Diction = {1: 1}
    prime = [2]
    cProfile.run('main()')
    end = time.perf_counter()
    print(len(prime), end-start)

Я ожидаю, что какой-то другой возможный держатель данных, который использует более оптимальный формат, будет работать лучше в этом случае — возможно, numpy или реализация numba.

Общие комментарии к обзору

Помимо производительности, есть несколько общих моментов обзора:

• Не ставьте код после if __name__ == '__main__': кроме как просто позвонить main().
• Поместите импорт в начало кода.
• Не полагайтесь на глобальные переменные для обмена данными между функциями (Diction, prime). Используйте аргументы функции и возвращаемые значения для обмена данными.
• Используйте осмысленные, описательные имена переменных. prime содержит primes (множественное число), так что лучше называть его primes. Diction относится к типу данных, который не является описательным именем переменной.
• Используйте прописные буквы для общемодульных констант (или, в этом случае: используйте параметры функции вместо константы n).
• cProfile уже дает вам время выполнения. Нет необходимости вычислять это снова, используя time.

Основная тема: производительность

Если я правильно понимаю ваш код, вы используете Diction в комбинации с absolute следить за \$(я — 1)!\$. Мне это кажется слишком сложным. Вы можете отслеживать факториал внутри цикла. Это экономит max_prime / 2 - 1 призывает dict операции pop а также update.

Реализуя все приведенные выше советы, это дает следующий код (подсказка типа и задокументированная строка документации):

import cProfile


def generate_primes(max_prime: int) -> list:
    """Returns a list of primes until max_prime (included)

    This function uses Wilkinson's theorem.

    Args:
        max_prime (int): upper limit for the primes to be calculated

    Returns:
        list: list of primes from 2 up until max_primes (included)
    """
    primes = [2]
    factorial = 2  # factorial of i-1

    for i in range(3, max_prime + 1, 2):
        if (factorial + 1) % i == 0:
            primes.append(i)
        factorial *= i * (i + 1)

    return primes


def main():
    max_prime = 100000

    profiler = cProfile.Profile()
    primes = profiler.runcall(generate_primes, max_prime)
    profiler.print_stats()

    print(len(primes))


if __name__ == '__main__':
    main()

Выполнение этого дает примерно ту же производительность, что и исходный код. 🙁

Как вы уже подозревали: есть два вычисления, которые, вероятно, составляют большую часть времени вычислений: модуль и произведение факториала (которое быстро становится очень и очень большим). Чтобы измерить это, я создал две функции, которые выполняют эти вычисления, поэтому cProfile может их профилировать. Это дает следующий результат:

109591 function calls in 10.326 seconds

Ordered by: standard name

ncalls  tottime  percall  cumtime  percall filename:lineno(function)
     1    0.872    0.872   10.326   10.326 slow_prime_test.py:12(generate_primes)
 49999    5.592    0.000    5.592    0.000 slow_prime_test.py:4(modulo_is_0)
 49999    3.862    0.000    3.862    0.000 slow_prime_test.py:8(multiply)
  9591    0.000    0.000    0.000    0.000 {method 'append' of 'list' objects}
     1    0.000    0.000    0.000    0.000 {method 'disable' of '_lsprof.Profiler' objects}

Таким образом, мы тратим 54% времени на расчет по модулю и 37% времени на расчет (крупных) произведений.

В идеале существует алгоритм, который может вычислить \$((n — 1)! + 1) % n\$ без фактического вычисления факториала. + 1 делает это очень сложным, и я не знаю никакого алгоритма, который может это сделать. Без + 1 мы могли бы использовать формула Лежандра вычислять \$н! % п\$.

Вероятно, вы можете немного ускорить свой код, используя C++ или любой другой скомпилированный высокопроизводительный язык, но факт остается фактом: теорема Уилсона не очень хорошо работает для генерации простых чисел. Об этом также говорится на странице теоремы в Википедии:

Теорема Уилсона использовалась для построения формул для простых чисел, но они слишком медленные, чтобы иметь практическое значение.

Используйте gmpy2 для ускорения вычислений больших целых чисел

Я нашел модуль Python gmpy2, что увеличивает производительность больших целочисленных вычислений. Просто добавление from gmpy2 import mpz и замена factorial = 2 с factorial = mpz('2') повышает производительность моего компьютера в 4,5 раза.

Using gmpy2: True
9592
         109591 function calls in 4.035 seconds

   Ordered by: standard name

   ncalls  tottime  percall  cumtime  percall filename:lineno(function)
        1    0.533    0.533    4.035    4.035 slow_prime_test.py:13(generate_primes)
    49999    3.077    0.000    3.077    0.000 slow_prime_test.py:5(modulo_is_0)
    49999    0.424    0.000    0.424    0.000 slow_prime_test.py:9(multiply)
     9591    0.001    0.000    0.001    0.000 {method 'append' of 'list' objects}
        1    0.000    0.000    0.000    0.000 {method 'disable' of '_lsprof.Profiler' objects}


Using gmpy2: False
         109591 function calls in 18.458 seconds

   Ordered by: standard name

   ncalls  tottime  percall  cumtime  percall filename:lineno(function)
        1    1.438    1.438   18.458   18.458 slow_prime_test.py:13(generate_primes)
    49999   14.895    0.000   14.895    0.000 slow_prime_test.py:5(modulo_is_0)
    49999    2.123    0.000    2.123    0.000 slow_prime_test.py:9(multiply)
     9591    0.001    0.000    0.001    0.000 {method 'append' of 'list' objects}
        1    0.000    0.000    0.000    0.000 {method 'disable' of '_lsprof.Profiler' objects}

На tio.run я даже получаю улучшение производительности примерно в 10 раз. Я думаю, это зависит от компьютерной платформы. Эта ссылка TIO включает окончательный код, который я использовал:

Попробуйте онлайн!

Редактировать (2): используя формулу Лежандра

Поскольку вы вычисляете все простые числа, мы можем использовать тот факт, что при тестировании, если \$n\$ является простым, у нас уже есть все простые числа, меньшие \$n\$. Мы можем использовать это, используя формулу Лежандра. Вместо расчета \$((n — 1)! + 1) % n\$мы можем использовать простую факторизацию \$n — 1\$ и powmod функция gmpy2. Таким образом, мы избегаем использования очень больших целых чисел и (возможно) повышаем скорость выполнения. Я реализовал это, но, к сожалению, он работает намного медленнее, чем оптимизированный для gmpy2 код, упомянутый выше.

def wilsons_prime_test(n: Union[int, gmpy2.mpz],
                       primes: list) -> bool:
    """Checks if ((n - 1)! + 1 ) % n == 0 using the given set of primes < n

    Args:
        n (int or mpz): number of which to check if it is prime
        primes (list): prime list which must give all primes < n

    Returns:
        bool: _description_
    """

    result = gmpy2.mpz('1')

    # Loop over all primes
    for p in primes:
        # Find largest power of p that divides n - 1
        tmp_n = n - 1
        largest_power = 0
        while tmp_n:
            tmp_n //= p
            largest_power += tmp_n

        # Multiply the result with prime ^ largest_power modulo n
        result *= gmpy2.powmod(p, largest_power, n)

    # Return if (result + 1) % n == 0
    return (result + 1) % n == 0