Проблема распараллеливания F # при вычислении совершенных чисел?

Question

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

Программа работает (почти) идеально, но когда я открываю диспетчер задач, он все еще работает в одном потоке.Это означает, что я, должно быть, делаю что-то не так, но мои знания F # все еще находятся в «начальной» фазе.

Я постараюсь изложить этот вопрос настолько ясно, насколько это возможно, но если мне не удастся это сделатьпоэтому, пожалуйста, дайте мне знать.

Совершенное число - это число, в котором сумма всех его делителей (кроме самого числа) равна самому числу (например, 6 является идеальным, поскольку сумма егоделителями 1, 2 и 3 являются 6).

Я использую простые числа для ускорения вычислений, то есть меня не интересуют (огромные) списки, в которых хранятся все делители.Для этого я использую формулу, которая оказалась правильной для Евклида: (2 * (степень числа - 1)) * (2 * (степень числа - 1)), где последний - простое число Мерсенна.Я использовал очень быстрый алгоритм из stackoverflow (автор @Juliet), чтобы определить, является ли данное число простым числом.

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

   let perfectNumbersTwo (n : int) =  
    seq { for i in 1..n do 
           if (PowShift i) - 1I |> isPrime 
           then yield PowShift (i-1) * ((PowShift i)-1I)
        } 

Вспомогательная функция PowShift реализована следующим образом:

    let inline PowShift (exp:int32) = 1I <<< exp ;;

Я использую битовый сдвигоператор, так как основа всех степенных расчетов от 2, следовательно, это может быть легкий путь.Конечно, я все еще благодарен за вклад в ответ на вопрос, который я задал по этому вопросу: F # Проблемы с питанием, которые принимают оба аргумента как bigints> F # Проблемы с питанием, которые принимают оба аргумента как bigints

Созданная Джульеттой функция ( заимствованная здесь ) выглядит следующим образом:

let isPrime ( n : bigint) = 
    let maxFactor = bigint(sqrt(float n))
    let rec loop testPrime tog =
        if testPrime > maxFactor then true
        elif n % testPrime = 0I then false
        else loop (testPrime + tog) (6I - tog)
    if n = 2I || n = 3I || n = 5I then true
    elif n <= 1I || n % 2I = 0I || n % 3I = 0I || n % 5I = 0I then false
    else loop 7I 4I;;

Используя этот код без параллельных операций, на моем ноутбуке требуется около 9 минут, чтобы найти 9-е совершенное число (который состоит из 37 цифр и может быть найден со значением 31 для показателя степени).Так как мой ноутбук имеет ЦП с двумя ядрами, и только одно из них работает с 50 процентами (полная загрузка для одного ядра), я подумал, что мог бы ускорить вычисления, параллельно вычисляя результаты.

Поэтому я изменил свойФункция perfectnumber выглядит следующим образом:

//Now the function again, but async for parallel computing
let perfectNumbersAsync ( n : int) =
    async {
        try
            for x in 1.. n do
                if PowShift x - 1I |> isPrime then
                    let result = PowShift (x-1) * ((PowShift x)-1I)
                    printfn "Found %A as a perfect number" result
        with
            | ex -> printfn "Error%s" (ex.Message);
    }

Чтобы вызвать эту функцию, я использую небольшую вспомогательную функцию для ее запуска:

 let runPerfects n =
    [n]
        |> Seq.map perfectNumbersAsync
        |> Async.Parallel
        |> Async.RunSynchronously
        |> ignore

Где результат асинхронного вычисления игнорируется, так как яотображение его в функции perfectNumbersAsync.

Приведенный выше код компилируется и запускается, однако он все еще использует только одно ядро ​​(хотя при вычислении 9-го идеального числа он работает на 10 секунд быстрее).Я боюсь, что это связано с вспомогательными функциями PowShift и isPrime, но я не уверен.Должен ли я поместить код этих вспомогательных функций в асинхронный блок perfectNumbersAsync?Это не улучшает читабельность ...

Чем больше я играю с F #, тем больше я учусь ценить этот язык, но, как и в этом случае, иногда мне нужны некоторые эксперты:).

Заранее спасибо за прочтение, я только надеюсь, что прояснил себя немного ...

Роберт.

Answer 1

Один быстрый комментарий по скорости и параллелизуемости,

Ваш isPrime равен O (sqrt (n)), и каждое последующее n примерно в 2 раза больше, чем последнее, поэтому для его вычисления потребуется примерно 1,5 раза, что означает, что вычисление последних чисел будет займет гораздо больше времени

Я сделал несколько взломов с тестированием на простоту, и нашел несколько полезных вещей:

1. Для больших N (вы тестируете числа с 20 цифрами), простая плотность на самом деле довольно низкая, поэтому вы будете делать много делений по составным числам. Лучшим подходом является предварительный расчет таблицы простых чисел (с использованием сита) до некоторого максимального предела (вероятно, определяемого количеством памяти). Обратите внимание, что вы, скорее всего, найдете факторы с небольшими числами. Как только вам не хватит памяти для вашей таблицы, вы можете проверить остальные числа с помощью существующей функции с большей начальной точкой.

2. Другой подход заключается в использовании нескольких потоков в проверке. Например, вы в настоящее время проверяете x,x+4,x+6... как факторы. Будучи немного умнее, вы можете сделать число конгруэнтным 1 мод 3 в 1 потоке, а числа конгруэнтному 2 мод 3 в другом потоке.

Нет. 2 является самым простым, но № 1 является более эффективным и обеспечивает потенциал для выполнения потока управления с исключениями OutOfMemoryException, которые всегда могут быть интересны

EDIT: Таким образом, я реализовал обе эти идеи, он обнаружил 2305843008139952128 почти мгновенно, а 2658455991569831744654692615953842176 занимает 7 минут на моем компьютере (четырехъядерный процессор AMD 3200). Большую часть времени тратится на проверку простых чисел 2 ^ 61, поэтому для проверки простых чисел лучше использовать лучший алгоритм: код здесь

let swatch = new System.Diagnostics.Stopwatch()
swatch.Start()
let inline PowShift (exp:int32) = 1I <<< exp ;;
let limit = 10000000 //go to a limit, makes table gen slow, but should pay off
printfn "making table"
//returns an array of all the primes up to limit
let table =
    let table = Array.create limit true //use bools in the table to save on memory
    let tlimit = int (sqrt (float limit)) //max test no for table, ints should be fine
    table.[1] <- false //special case
    [2..tlimit] 
    |> List.iter (fun t -> 
        if table.[t]  then //simple optimisation
            let mutable v = t*2
            while v < limit do
                table.[v] <- false
                v <- v + t)
    let out = Array.create (50847534) 0I //wolfram alpha provides pi(1 billion) - want to minimize memory
    let mutable idx = 0
    for x in [1..(limit-1)] do
        if table.[x] then
            out.[idx] <- bigint x
            idx <- idx + 1
    out |> Array.filter (fun t -> t <> 0I) //wolfram no is for 1 billion as limit, we use a smaller number
printfn "table made"

let rec isploop testprime incr max n=
    if testprime > max then true
    else if n % testprime = 0I then false
    else isploop (testprime + incr) incr max n

let isPrime ( n : bigint) = 
    //first test the table
    let maxFactor = bigint(sqrt(float n))
    match table |> Array.tryFind (fun t -> n % t = 0I && t <= maxFactor) with
    |Some(t) -> 
        false
    |None -> //now slow test
        //I have 4 cores so
        let bases = [|limit;limit+1;limit+3;limit+4|] //uses the fact that 10^x congruent to 1 mod 3
        //for 2 cores, drop last 2 terms above and change 6I to 3I
        match bases |> Array.map (fun t -> async {return isploop (bigint t) 6I maxFactor n}) |> Async.Parallel |> Async.RunSynchronously |> Array.tryFind (fun t -> t = false) with
        |Some(t) -> false
        |None -> true


let pcount = ref 0
let perfectNumbersTwo (n : int) =  
    seq { for i in 2..n do 
           if (isPrime (bigint i)) then
                if (PowShift i) - 1I |> isPrime then
                    pcount := !pcount + 1
                    if !pcount = 9 then
                        swatch.Stop()
                        printfn "total time %f seconds, %i:%i m:s"  (swatch.Elapsed.TotalSeconds) (swatch.Elapsed.Minutes) (swatch.Elapsed.Seconds)
                    yield PowShift (i-1) * ((PowShift i)-1I)
        } 


perfectNumbersTwo 62 |> Seq.iter (printfn "PERFECT: %A") //62 gives 9th number

printfn "done"
System.Console.Read() |> ignore

Answer 2

@ Комментарий Джеффри Сакса определенно интересен, поэтому я потратил некоторое время на небольшой эксперимент. Тест Лукаса-Лемера записывается следующим образом:

let lucasLehmer p =
    let m = (PowShift p) - 1I
    let rec loop i acc =
        if i = p-2 then acc
        else loop (i+1) ((acc*acc - 2I)%m)
    (loop 0 4I) = 0I

С помощью теста Лукаса-Лемера я могу очень быстро получить первые несколько идеальных чисел:

let mersenne (i: int) =     
    if i = 2 || (isPrime (bigint i) && lucasLehmer i) then
        let p = PowShift i
        Some ((p/2I) * (p-1I))
    else None

let runPerfects n =
    seq [1..n]
        |> Seq.choose mersenne
        |> Seq.toArray

let m1 = runPerfects 2048;; // Real: 00:00:07.839, CPU: 00:00:07.878, GC gen0: 112, gen1: 2, gen2: 1

Тест Лукаса-Лемера помогает сократить время проверки простых чисел. Вместо проверки делимости 2 ^ p-1, которая занимает O(sqrt(2^p-1)), мы используем тест на простоту, который не более O(p^3). С n = 2048 я могу найти первые 15 чисел Мерсенна за 7,83 секунды. 15-й номер Мерсенна - это номер i = 1279, состоящий из 770 цифр.

Я попытался распараллелить runPerfects, используя модуль PSeq в F # Powerpack . PSeq не сохраняет порядок исходной последовательности, поэтому, если честно, я отсортировал выходную последовательность. Поскольку тест на простоту довольно сбалансирован между индексами, результат весьма обнадеживающий:

#r "FSharp.Powerpack.Parallel.Seq.dll"    
open Microsoft.FSharp.Collections

let runPerfectsPar n =
    seq [1..n]
        |> PSeq.choose mersenne
        |> PSeq.sort (* align with sequential version *)
        |> PSeq.toArray 

let m2 = runPerfectsPar 2048;; // Real: 00:00:02.288, CPU: 00:00:07.987, GC gen0: 115, gen1: 1, gen2: 0

При том же входе параллельная версия заняла 2,28 секунды, что эквивалентно ускорению в 3,4 раза на моей четырехъядерной машине. Я считаю, что результат может быть улучшен, если вы используете Parallel.For конструировать и разумно разделить входной диапазон.