Как (эффективно) найти самый большой простой множитель числа с помощью Haskell?

Я пытаюсь практиковать Haskell, решая некоторые задачи в Project Euler. В задаче 3 нам нужно найти самый большой простой делитель числа 600851475143, что я уже делал в Java несколько лет назад.

Я придумал следующее:

primes :: [Int]
primes = sieve [2..]
    where sieve (p:xs) = p : sieve (filter (\x -> x `rem` p /= 0) xs)

biggestPrimeFactor :: Int -> Int
biggestPrimeFactor 1 = 0
biggestPrimeFactor x =
    if x `elem` takeWhile (< x + 1) primes
        then x
        else last (filter (\y -> x `rem` y == 0) (takeWhile (< x `div` 2) primes))

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

Обратите внимание, что я ищу не «оптимальное решение», а скорее такое, которое, по крайней мере, умеренно эффективно для больших чисел и просто для понимания и реализации, поскольку я все еще новичок в Haskell.

🤔 А знаете ли вы, что...
Haskell является хорошим выбором для параллельного и многозадачного программирования.


1
71
1

Ответ:

Решено

Здесь у вас есть два основных источника медлительности. Легче всего обратиться к граничному условию в biggestPrimeFactor. Проверка до p > x `div` 2 асимптотически хуже, чем проверка до p^2 > x. Но даже это очень неоптимально, когда число имеет много факторов. Наибольший фактор может быть намного меньше, чем sqrt x. Если вы постоянно уменьшаете целевое число по мере нахождения факторов, вы можете учесть это и значительно ускорить обработку случайных входных данных.

Вот пример этого, включая примечание Даниэля Вагнера из комментариев:

-- Naive trial division against a list of primes. Doesn't do anything
-- intelligent when asked to factor a number less than 2.
factorsNaive :: [Integer] -> Integer -> [Integer]
factorsNaive primes@(p : ps) x
    | p * p > x = [x]
    | otherwise = case x `quotRem` p of
        (q, 0) -> p : factorsNaive primes q
        _ -> factorsNaive ps x

Несколько заметок:

  1. Я решил передать список primes. Это относится к следующему разделу, но это также позволило мне написать это без помощника.
  2. Я специализировался на Integer вместо Int, потому что хотел подбрасывать большие числа, не заботясь о том, что такое maxBound :: Int. Это медленнее, но я решил сначала использовать корректность по умолчанию.
  3. Я убрал обход списка ввода. Выполнение этого за один проход немного эффективнее, но в основном это чище.
  4. Строго говоря, это правильно, даже если входной список содержит непростые числа, при условии, что список начинается с 2, монотонно не убывает и в конечном итоге содержит все простые числа.
  5. Обратите внимание, что когда он рекурсивно, он либо отбрасывает простое число, либо создает его. Это никогда не будет делать и то, и другое одновременно. Это простой способ убедиться, что он не пропускает повторяющиеся факторы.

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

Или, по крайней мере, с факторингом все в порядке, если у вас есть удобный список простых чисел. Оказывается, это вторая основная причина замедления вашего кода. Ваш список простых чисел генерируется медленно, поскольку он становится длиннее.

Ваше определение primes, по сути, объединяет кучу фильтров в списке ввода. Каждое произведенное простое число должно пройти проверку фильтра для каждого предыдущего простого числа. Это может показаться знакомым — для генерации первых n простых чисел требуется как минимум O(n^2). (На самом деле это больше, потому что деление становится более затратным по мере увеличения чисел, но давайте пока проигнорируем это.) Это известный (математикам, я должен был проверить это, чтобы быть уверенным) результат, что количество простых чисел меньше или равно n приближается к n/ln n по мере того, как n становится большим. Это приближается к линейному, когда n становится большим, поэтому создание списка от primes до n приближается к O (n ^ 2), когда n становится большим.

(Да, с этим аргументом полный бардак. Формальная его версия представлена ​​в статье Мелиссы О'Нил «Подлинное сито Эратосфена». См. там гораздо более строгую аргументацию результата.)

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

-- trial division. let's work in Integer for predictable correctness
-- on positive numbers
trialPrimes :: [Integer]
trialPrimes = 2 : sieve [3, 5 ..]
  where
    sieve (p : ps) = p : sieve (filter (\x -> x `rem` p /= 0) ps)

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

В ghci (без компиляции или оптимизации, а они действительно могут иметь значение) это было достаточно быстро, чтобы разложить произведение двух пятизначных простых чисел за несколько секунд.

ghci> :set +s
ghci> factorsNaive trialPrimes $ 84761 * 60821
[60821,84761]
(5.98 secs, 4,103,321,840 bytes)

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

ghci> factorsNaive trialPrimes $ 84761 * 60821
[60821,84761]
(0.01 secs, 6,934,688 bytes)

Это также показывает, что во время выполнения абсолютно доминирует создание списка простых чисел. Наивная факторизация почти мгновенна в том масштабе, когда список простых чисел уже находится в памяти.

Но не стоит особо доверять производительности интерпретируемого кода.

main :: IO ()
main = print (factorsNaive trialPrimes $ 84761 * 60821)

дает

carl@DESKTOP:~/hask/2023$ ghc -O2 -rtsopts factor.hs
[1 of 2] Compiling Main             ( factor.hs, factor.o )
[2 of 2] Linking factor
carl@DESKTOP:~/hask/2023$ ./factor +RTS -s
[60821,84761]
   1,884,787,896 bytes allocated in the heap
      32,303,080 bytes copied during GC
          89,072 bytes maximum residency (2 sample(s))
          29,400 bytes maximum slop
               7 MiB total memory in use (0 MB lost due to fragmentation)

                                     Tot time (elapsed)  Avg pause  Max pause
  Gen  0       326 colls,     0 par    0.021s   0.021s     0.0001s    0.0002s
  Gen  1         2 colls,     0 par    0.000s   0.000s     0.0002s    0.0004s

  INIT    time    0.000s  (  0.000s elapsed)
  MUT     time    0.523s  (  0.522s elapsed)
  GC      time    0.021s  (  0.022s elapsed)
  EXIT    time    0.000s  (  0.007s elapsed)
  Total   time    0.545s  (  0.550s elapsed)

  %GC     time       0.0%  (0.0% elapsed)

  Alloc rate    3,603,678,988 bytes per MUT second

  Productivity  96.0% of total user, 94.8% of total elapsed

Это сократило время выполнения с шести секунд до полсекунды. (Да, +RTS -s довольно многословно для этого, но это быстро и просто.) Я думаю, что это разумное место, чтобы остановиться на коде начального уровня.

Если вы хотите изучить более эффективную генерацию простых чисел, пакет primes on hackage содержит реализацию алгоритма из статьи О'Нила и реализацию наивного факторинга, эквивалентную приведенной здесь.