Как правильно реализовать факторизацию Ферма в Python? ⇐ Python

[Anonymous]

Я пытаюсь реализовать эффективные алгоритмы факторизации простых чисел в Python. Это не домашнее задание и не работа, это совершенно из любопытства.
Я узнал, что факторизация простых чисел очень сложна:
Я хочу внедрить эффективные алгоритмы для этого как добровольную задачу. Я решил сначала реализовать метод факторизации Ферма, поскольку он кажется достаточно простым.
Код Python, непосредственно транслируемый из псевдокода:

Код: Выделить всё

def Fermat_Factor(n):
a = int(n ** 0.5 + 0.5)
b2 = abs(a**2 - n)
while int(b2**0.5) ** 2 != b2:
a += 1
b2 = a**2 - n

return a - b2**0.5, a + b2**0.5

(Мне нужно использовать abs, иначе b2 легко будет отрицательным, а приведение int завершится с ошибкой TypeError, потому что корень complex)
Как видите, он возвращает два целых числа, произведение которых равно входным данным, но возвращает только два выходных значения и не гарантирует простоту множителей. Я понятия не имею, насколько эффективен этот алгоритм, но факторизация полупростых чисел с использованием этого метода намного эффективнее, чем метод пробного деления, использованный в моем предыдущем вопросе: почему факторизация произведений близких простых чисел происходит намного медленнее, чем произведений разнородных простых чисел.

Код: Выделить всё

In [20]: %timeit FermatFactor(3607*3803)
2.1 μs ± 28.2 ns per loop (mean ± std. dev. of 7 runs, 100,000 loops each)

In [21]: FermatFactor(3607*3803)
Out[21]: [3607, 3803]

In [22]: %timeit FermatFactor(3593 * 3671)
1.69 μs ± 31 ns per loop (mean ± std. dev. of 7 runs, 1,000,000 loops each)

In [23]: FermatFactor(3593 * 3671)
Out[23]: [3593, 3671]

In [24]: %timeit FermatFactor(7187 * 7829)
4.94 μs ± 47.4 ns per loop (mean ± std. dev. of 7 runs, 100,000 loops each)

In [25]: FermatFactor(7187 * 7829)
Out[25]: [7187, 7829]

In [26]: %timeit FermatFactor(8087 * 8089)
1.38 μs ± 12.9 ns per loop (mean ± std. dev. of 7 runs, 1,000,000 loops each)

In [27]: FermatFactor(8087 * 8089)
Out[27]: [8087, 8089]

Поэтому я хочу использовать этот алгоритм для генерации всех простых множителей любого заданного целого числа (конечно, я знаю, что это работает только с нечетными целыми числами, но это не проблема, поскольку степени двойки можно тривиально исключить с помощью бит-хакинга). Самый простой способ, который я могу придумать, — это рекурсивно вызывать Fermat_Factor до тех пор, пока n не станет простым числом. Я не знаю, как проверить, является ли число простым в этом алгоритме, но кое-что заметил:

Код: Выделить всё

In [3]: Fermat_Factor(3)
Out[3]: (1.0, 3.0)

In [4]: Fermat_Factor(5)
Out[4]: (1.0, 3.0)

In [5]: Fermat_Factor(7)
Out[5]: (1.0, 7.0)

In [6]: Fermat_Factor(11)
Out[6]: (1.0, 11.0)

In [7]: Fermat_Factor(13)
Out[7]: (1.0, 13.0)

In [8]: Fermat_Factor(17)
Out[8]: (3.0, 5.0)

In [9]: Fermat_Factor(19)
Out[9]: (1.0, 19.0)

In [10]: Fermat_Factor(23)
Out[10]: (1.0, 23.0)

In [11]: Fermat_Factor(29)
Out[11]: (3.0, 7.0)

In [12]: Fermat_Factor(31)
Out[12]: (1.0, 31.0)

In [13]: Fermat_Factor(37)
Out[13]: (5.0, 7.0)

In [14]: Fermat_Factor(41)
Out[14]: (1.0, 41.0)

Первое число в выводе этого алгоритма для многих простых чисел равно 1, но не для всех, поэтому его нельзя использовать для определения момента, когда рекурсия должна остановиться. Я усвоил это на собственном горьком опыте.
Поэтому я просто решил вместо этого использовать проверку членства предварительно сгенерированного набора простых чисел. Естественно, это приведет к ошибке RecursionError: превышена максимальная глубина рекурсии, когда входное число является простым числом, превышающим максимальное значение набора. Поскольку у меня нет бесконечной памяти, это следует рассматривать как деталь реализации.
Итак, я реализовал рабочую версию (для некоторых входных данных), но для некоторых допустимых входных данных (продуктов простых чисел в пределах лимита) почему-то алгоритм не дает правильного результата:

Код: Выделить всё

import numpy as np
from itertools import cycle

TRIPLE = ((4, 2), (9, 6), (25, 10))
WHEEL = ( 4, 2, 4, 2, 4, 6, 2, 6 )
def prime_sieve(n):
primes = np.ones(n + 1, dtype=bool)
primes[:2] = False
for square, double in TRIPLE:
primes[square::double] = False

wheel = cycle(WHEEL)
k = 7
while (square := k**2)  TEST_LIMIT:
raise ValueError('Number too large')

if n in PRIME_SET:
return [n]

a = int(n ** 0.5 + 0.5)
if a ** 2 == n:
return FermatFactor(a) + FermatFactor(a)

b2 = abs(a**2 - n)
while int(b2**0.5) ** 2 != b2:
a += 1
b2 = a**2 - n

return FermatFactor(factor := int(a - b2**0.5)) + FermatFactor(n // factor)

Это работает для многих входных данных:

Код: Выделить всё

In [18]: FermatFactor(255)
Out[18]: [3, 5, 17]

In [19]: FermatFactor(511)
Out[19]: [7, 73]

In [20]: FermatFactor(441)
Out[20]: [3, 7, 3, 7]

In [21]: FermatFactor(3*5*823)
Out[21]: [3, 5, 823]

In [22]: FermatFactor(37*333667)
Out[22]: [37, 333667]

In [23]: FermatFactor(13 * 37 * 151 * 727 * 3607)
Out[23]: [13, 37, 727, 151, 3607]

Но не все:

Код: Выделить всё

In [25]: FermatFactor(5 * 53 * 163)
Out[25]: [163, 13, 2, 2, 5]

In [26]: FermatFactor(3*5*73*283)
Out[26]: [17, 3, 7, 3, 283]

In [27]: FermatFactor(3 * 11 * 29 * 71 *  137)
Out[27]: [3, 11, 71, 61, 7, 3, 3]

Почему именно так? Как я могу это исправить?

Подробнее здесь: https://stackoverflow.com/questions/793 ... -in-python