Есть ли способ оптимизировать этот код для нахождения делителей числа? - PullRequest
Новая
       87

Есть ли способ оптимизировать этот код для нахождения делителей числа?

4 голосов
/ 09 апреля 2020

Я написал программу на Julia для эффективного вычисления делителей числа n. Алгоритм является оригинальным (насколько мне известно), и он основан на Сите Эратосфена . По сути, это работает так:

Для данного простого числа p, пусть p^k || n; каждое число m в списке, удовлетворяющее p^{k+1} | m, удаляется, и этот процесс повторяется для каждого простого p < n.

Простые числа вычисляются на месте с использованием традиционного сита Эратосфена. 

function ν(p, n)     #returns the smallest power of p that does not divide n
    q = 1

    for i = 0:n
        if mod(n, q) != 0
            return (i, q) 
        end

        q *= p
    end
end

function divisors(n)    #returns a list of divisors of n
    dsieve, psieve = BitArray([true for i = 1:n]), BitArray([true for i = 1:n])
    psieve[1] = false

    for i = 1:n
        if psieve[i] && dsieve[i]
            #sieving out the non-primes
            for j = i^2:i:n
                psieve[j] = false
            end

            #sieving out the non-divisors
            v = ν(i, n)[2]
            for j = v:v:n
                dsieve[j] = false
            end
        end
    end
    return dsieve #the code for converting this BitArray to an array of divisors has been omitted for clarity
end

Хотя это прекрасно работает, я считаю неэффективным использование двух сит одновременно. Я думаю, что эту проблему можно решить, разрешив каждому элементу в массиве sieve принимать три разных значения (соответствующих unchecked, divisor и not divisor), но тогда это уже не может быть реализовано как BitArray .

Я также пытался изменить функцию ν, чтобы сделать ее более эффективной: 

function ν₀(p, n)      #the same as ν, but implemented differently
    q = p
    while mod(n, q) == 0
        q = q^2
    end

    q = floor(Int64, √q)
    q < p ? 1 : q * ν₀(p, n÷q)    #change 1 to p to get the smallest power of p that does not divide n
end

Хотя это и сложнее, но немного быстрее, чем в предыдущем алгоритме, особенно когда сила p деления n велика.

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

1 Ответ

3 голосов
/ 09 апреля 2020

Есть пара вещей, на которые я могу указать - 

dsieve, psieve = BitArray([true for i = 1:n]), BitArray([true for i = 1:n])

выделяет дважды для каждого массива (список comp и затем преобразование). Это будет хорошо: (Edit: @DNF указывает здесь на превосходство Vector{Bool}) 

dsieve = fill(true, n)
psieve = fill(true, n)

Далее мы можем убедиться, что используем любой вид более быстрой индексации, используя 

for i in eachindex(psieve)

вместо ручного диапазона. Затем вы можете добавить к -1 oop с помощью 

@inbounds for i in eachindex(psieve)

или go еще больше, если вы используете Julia 1.3 или более позднюю версию, и многопоточность (если вы установили JULIA_NUM_THREADS до запуска) 

@inbounds Threads.@threads for i in eachindex(psieve)
...