Возможна оптимизация или использование параллельных вычислений

Question

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

вход 100 2 даст выход 3 потому 100 = 10^2 = 6^2 + 8^2 = 1^2 + 3^2 + 4^2 + 5^2 + 7^2 и ввод 100 3 даст выход 1, потому что 100 = 1^3 + 2^3 + 3^3 + 4^3

Так моя функция для решения этой проблемы:

findNums :: Int -> Int -> Int 
findNums a b = length [xs | xs <- (drop 1 .) subsequences [pow x b | x <- [1..c]], foldr (+) (head xs) (tail xs) == a] where c = root a b 0 

root :: Int -> Int -> Int -> Int 
root n a i 
    | pow i a <= n && pow (i+1) a > n = i 
    | otherwise = root n a (i+1) 

pow :: Int -> Int -> Int 
pow _ 0 = 1 
pow n a = n * pow n (a - 1) 

Я нахожу все возможные значения, которые могут быть в моем наборе чисел, которые будут содержать нужный номер. Затем я нахожу все возможные подсписки внутри этого списка и вижу, сколько из них складывается до нужного числа. Это работает, но поскольку это исчерпывающий поиск, требуется много времени для ввода, например, 800 2. Можно ли оптимизировать последовательности так, чтобы возвращались только «правдоподобные» подпоследовательности? Или лучше посмотреть на параллельные вычисления для такого рода проблем?

Answer 1

Давайте рассмотрим несколько вещей.

Бенчмаркинг

Первое: давайте удостоверимся, что мы на самом деле делаем улучшение, как мы идем!Для этого нам понадобятся некоторые ориентиры. Пакет criterion идеально подходит для этого. Мы также обязательно скомпилируем с оптимизацией (так -O2 по всем вызовам GHC). Вот как простая настройка вверх тест может быть:

import Criterion.Main 

-- your code goes here 

main = defaultMain 
    [ bench "findNums 100 2" (nf (uncurry findNums) (100, 2)) 
    , bench "findNums 800 2" (nf (uncurry findNums) (800, 2)) 
    ] 

Можно также осуществить тест, как nf (findNums 100) 2, но я выбираю этот путь, так что мы не можем «обмануть», предварительно рассчитав таблицу поиска для 100, тем самым толкая вся работа в контрольной установке, а не в той части, где фактически выполняется эталон. Вот результат для первоначальной реализации:

benchmarking 100 2 
time     762.7 ns (757.4 ns .. 768.5 ns) 
        1.000 R² (1.000 R² .. 1.000 R²) 
mean     762.5 ns (760.4 ns .. 765.3 ns) 
std dev    7.706 ns (6.378 ns .. 10.59 ns) 

benchmarking 800 2 
time     29.17 s (28.28 s .. 29.87 s) 
        1.000 R² (1.000 R² .. 1.000 R²) 
mean     29.26 s (29.08 s .. 29.35 s) 
std dev    159.2 ms (0.0 s .. 165.2 ms) 
variance introduced by outliers: 19% (moderately inflated) 

Используйте библиотеки

Теперь, низко висящие плоды использовать существующие реализации вещей и надеются, что их авторы сделали что-то лучше, чем у нас. С этой целью мы будем использовать стандартную функцию (^) вместо pow и integerRoot от arithmoi вместо root. Кроме того, мы заменим ленивый foldr на строгий номер foldl. Для моего собственного здравомыслия я также переформатировал очень длинную линию на короткие. Полный результат теперь выглядит следующим образом:

import Criterion.Main 
import Data.List 
import Math.NumberTheory.Powers 

sum' :: Num a => [a] -> a 
sum' = foldl' (+) 0 

findNums :: Int -> Int -> Int 
findNums a b = length 
    [ xs 
    | xs <- drop 1 . subsequences $ [x^b | x <- [1..c]] 
    , sum' xs == a 
    ] where c = integerRoot b a 

main = defaultMain 
    [ bench "100 2" (nf (uncurry findNums) (100, 2)) 
    , bench "800 2" (nf (uncurry findNums) (800, 2)) 
    ] 

Результаты тестирования выглядят, как это сейчас:

benchmarking 100 2 
time     722.8 ns (721.3 ns .. 724.3 ns) 
        1.000 R² (1.000 R² .. 1.000 R²) 
mean     722.6 ns (721.4 ns .. 724.1 ns) 
std dev    4.440 ns (3.738 ns .. 5.674 ns) 

benchmarking 800 2 
time     17.16 s (16.93 s .. 17.64 s) 
        1.000 R² (1.000 R² .. 1.000 R²) 
mean     17.05 s (16.99 s .. 17.15 s) 
std dev    88.10 ms (0.0 s .. 94.58 ms) 

Чуть менее в два раза быстрее с очень небольшим усилием. Ницца!

Лучше алгоритм

Существенная проблема с subsequences является то, что, даже если мы вычисляем, что sum' [x,y,z] > a, мы по-прежнему смотреть на все более длинные подпоследовательности, которые начинаются с [x,y,z]. Учитывая структуру возвращаемого типа 44843689053679157146564308888, мы мало что можем с этим поделать; поэтому давайте разработаем реализацию, которая дает нам немного больше структуры. Мы построим дерево, где пути от корня до любого внутреннего узла дают нам подпоследовательность.

import Data.Tree 

subsequences :: [a] -> Forest a 
subsequences [] = [] 
subsequences (x:xs) = Node x rest : rest where 
    rest = subsequences xs 

(Просто для удовольствия, это приводит к экспоненциально большие семантические деревья с очень малым использования пространства - примерно столько же пространства, как в первоначальном списке -. Из-за агрессивного обмена поддерева) Что здорово об этом представлении, если мы прекратите поиск, мы отсекаем огромные полосы неинтересных результатов. Это может быть реализовано путем реализации что-то вроде takeWhile для списков:

takeWhileTree :: Monoid m => (m -> Bool) -> Forest m -> Forest m 
takeWhileTree predicate = goForest mempty where 
    goForest m forest = forest >>= goTree m 
    goTree m (Node m' children) = 
     [Node m (goForest (m <> m') children) | predicate m'] 

Давайте дадим ему попробовать. Полный код теперь:

import Criterion.Main 
import Data.Foldable 
import Data.Monoid 
import Data.Tree 
import Math.NumberTheory.Powers 

subsequencesTree :: [a] -> Forest a 
subsequencesTree [] = [] 
subsequencesTree (x:xs) = Node x rest : rest where 
    rest = subsequencesTree xs 

takeWhileTree :: Monoid m => (m -> Bool) -> Forest m -> Forest m 
takeWhileTree predicate = goForest mempty where 
    goForest m forest = forest >>= goTree m 
    goTree m (Node m' children) = let m'' = m <> m' in 
     [Node m' (goForest m'' children) | predicate m''] 

leaves :: Forest a -> [[a]] 
leaves [] = [[]] 
leaves forest = do 
    Node x children <- forest 
    xs <- leaves children 
    return (x:xs) 

findNums :: Int -> Int -> Int 
findNums a b = length 
    [ xs 
    | xs <- leaves 
      . takeWhileTree (<= Sum a) 
      . subsequencesTree 
      $ [Sum (x^b) | x <- [1..c]] 
    , fold xs == Sum a 
    ] where c = integerRoot b a 

main = defaultMain 
    [ bench "100 2" (nf (uncurry findNums) (100, 2)) 
    , bench "800 2" (nf (uncurry findNums) (800, 2)) 
    ] 

Это выглядит как много работы, но с таймингами, это действительно окупается:

benchmarking 100 2 
time     16.67 μs (16.53 μs .. 16.77 μs) 
        0.999 R² (0.999 R² .. 1.000 R²) 
mean     16.60 μs (16.52 μs .. 16.72 μs) 
std dev    325.4 ns (270.5 ns .. 444.1 ns) 
variance introduced by outliers: 17% (moderately inflated) 

benchmarking 800 2 
time     22.59 ms (22.26 ms .. 22.89 ms) 
        0.999 R² (0.999 R² .. 1.000 R²) 
mean     22.44 ms (22.34 ms .. 22.57 ms) 
std dev    260.3 μs (191.6 μs .. 332.2 μs) 

Это фактор ускорения около 1000 на findNums 800 2.

распараллеливания

У меня был пойти на распараллеливание это с помощью concat и parMap в takeWhileTree вместо (>>=), так что отдельные ветви дерева будут изучены параллельно.В каждом случае накладные расходы на распараллеливание намного перевешивают преимущество нескольких потоков. Хорошо, что мы поставили этот бенчмарк в начале!

Answer 2

Как вы полагаете, здесь есть место для оптимизации, не прибегая к распараллеливанию вещей (которые, помните, могут в лучшем случае дать ускорение 4 раза, если вы переходите от одного до четырех параллельных потоков).

Что делает функция subsequences, по существу идет по списку, и для каждого элемента он создает две ветви выполнения: один, где этот элемент включен, и тот, где он отсутствует. Т.е., subsequences [1,2,3] делает:

      start 
        /-------/ \-------\   (take 1 or not) 
      [1,..]     [..] 
      / \    / \  (take 2 or not) 
    [1,2,..]  [1,..]  [2,..] [..] 
    /\   /\  /\ /\ (take 3 or not) 
[1,2,3] [1,2] [1,3] [1] [2,3] [2] [3] [] 

В результате subsequences [1,2,3] затем список, содержащий узлы листа в нижней части.

Теперь, на каждом из промежуточных узлов (то есть [1,2,..]), мы можем проверить результат применения функции значения (т. Е. Суммы квадратов или кубов или т. Д.) К уже принятым номерам. Если мы уже выше цели, нет смысла продолжать эту ветку. Если мы напишем эту логику подпоследовательности поколения самого, мы можем сделать это:

findNums :: Int -> Int -> Int 
findNums a b = findNums' a b 1 0 

findNums' :: Int -> Int -> Int -> Int -> Int 
findNums' a b c s 
    | s + c^b > a = 0 
    | s + c^b == a = 1 
    | otherwise = findNums' a b (c+1) (s + c^b) + 
        findNums' a b (c+1) s 

Здесь c нашего счетчика и s является суммой степеней чисел мы подобрали. Есть три случая в findNums':

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

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

Если ни одно из них не является истинным, мы рекурсируем далее с двумя разными ветвями: один, где мы добавляем c^b к нашей сумме, и тот, где мы воздерживаемся от этого. Мы добавляем результаты вместе, что означает, что результатом здесь будет количество ветвей ниже этой точки, которые нашли правильное решение.

Answer 3

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

Чтобы упростить ситуацию, давайте просто рассмотрим суммы квадратов. Кроме того, мы сначала рассмотрим упорядоченные последовательности (с допустимыми допустимыми значениями ); позже мы рассмотрим, как изменить алгоритм для создания только неупорядоченных последовательностей без каких-либо повторных номеров .

Это наша первая попытка. Алгоритм основан на этой идее:

Идея 1:

Чтобы получить последовательность, сумма квадратов п, первые забрать значение с и последовательность хз, сумма которых квадратов - nc * c и положите эти два вместе.

-- an integer sqrt function 
isqrt n = floor $ (sqrt (fromIntegral n) :: Double) 

pows2a :: Int -> [ [Int] ] 
pows2a n 
    | n < 0  = [] 
    | n == 0 = [ [] ] 
    | otherwise = [ (c:xs) | c <- [start,start-1..1], xs <- pows2a (n-c*c) ] 
    where start = isqrt n 

Это работает, но возвращает перестановки решений, а также растворы с повторяющимися элементами - например, pos2a 6 возвращает [2,1,1], [1,2,1], [1,1,2] и [1,1,1,1,1,1].

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

Идея 2:

Чтобы получить последовательность, сумма квадратов п, первый выбрать значение c и последовательность xs, чья сумма квадратов - nc * c и все элементы которой составляют < c и положите эти два вместе.

Это лишь небольшая модификация нашего первого алгоритма:

pows2b :: Int -> [[Int]] 
pows2b n 
    | n < 0  = [] 
    | n == 0 = [ [] ] 
    | otherwise = [ (c:xs) | c <- [start, start-1..1], xs <- pows2b (n-c*c), all (< c) xs ] 
    where 
    start = isqrt n 

Это работает, но вызов, как pows2b 100 занимает много времени, чтобы закончить, потому что мы делаем вызовы pows2b с несколько раз тот же самый аргумент ,

Мы можем решить эту проблему memoizing результаты, и это то, что pows2c делает:

powslist = map pows2c [0..] 
pows2c n 
    | n == 0 = [ [] ] 
    | otherwise = [ (c:xs) | c <- [s,s-1..1], xs <- powslist !! (n-c*c), all (< c) xs ] 
    where s = isqrt n 

Здесь рекурсивный вызов с аргументом n-c*c заменяется поиск в список, который является способом кэширования ответ.