Что является узким местом в этом предикате, связанном с простыми числами?

Question

Итак, вот оно: я пытаюсь вычислить сумму всех простых чисел ниже двух миллионов (для эта проблема ), но моя программа очень медленная.Я знаю, что сам по себе алгоритм ужасно плох и грубоват, но мне кажется, что он намного медленнее, чем следовало бы.
Здесь я ограничиваю поиск до 20000, чтобы результат не был слишком долгим.
Я не думаю, что этот предикат трудно понять, но я все равно объясню: я вычисляю список всех простых чисел ниже 20000 и затем суммирую их.Часть суммы в порядке, часть простых чисел действительно медленная.

problem_010(R) :-
    p010(3, [], Primes),
    sumlist([2|Primes], R).
p010(20001, Primes, Primes) :- !.
p010(Current, Primes, Result) :-
    (
        prime(Current, Primes)
    ->  append([Primes, [Current]], NewPrimes)
    ;   NewPrimes = Primes
    ),
    NewCurrent is Current + 2,
    p010(NewCurrent, NewPrimes, Result).
prime(_, []) :- !.
prime(N, [Prime|_Primes]) :- 0 is N mod Prime, !, fail.
prime(ToTest, [_|Primes]) :- prime(ToTest, Primes).

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

РЕДАКТИРОВАТЬ: я уже что-то нашел, добавив новые простые числа вместо того, чтобы поместить их в начало списка,У меня есть простые числа, которые встречаются в начале чаще, так что это примерно в 3 раза быстрее.Тем не менее, нужно немного понимания:)

Answer 1

Во-первых, Пролог здесь не ошибается.

Есть очень умные способы генерирования простых чисел.Но в качестве дешевого начала просто накапливайте простые числа в обратном порядке!(7,9 с -> 2,6 с) Таким образом, меньшие тестируются быстрее.Затем рассмотрите возможность проверки только с простыми числами до 141. Большие простые числа не могут быть фактором.

Тогда вместо того, чтобы переходить только через числа, не делимые на 2, вы можете добавить 3, 5, 7.

Есть люди, пишущие статьи по этой "проблеме".См., Например, этот документ , хотя это немного сложная дискуссия о том, каким на самом деле был «подлинный» алгоритм, 22 века назад, когда последний выпуск абакуса отмечался как таблетки Саламина .

Answer 2

Рассмотрим, например, метод сита («Сито Эратосфена»): сначала создайте список [2,3,4,5,6, .... N], используя, например, numlist / 3. Первое число в списке - простое число, сохраните его. Удалите его кратные из остальной части списка. Следующее число в оставшемся списке снова является простым. Снова устраните его кратные. И так далее. Список будет сокращаться довольно быстро, и у вас останутся только простые числа.

Answer 3

Прежде всего, добавление в конец списка с помощью append/3 происходит довольно медленно. Если необходимо, используйте вместо этого списки различий. (Лично я стараюсь максимально избегать append/3)

Во-вторых, ваша prime/2 всегда проверяет весь список при проверке простого числа. Это излишне медленно. Вместо этого вы можете просто проверить идентификатор, найти интегральный коэффициент вплоть до квадратного корня числа, которое вы хотите проверить.

problem_010(R) :-
    p010(3, 2, R).
p010(2000001, Primes, Primes) :- !.
p010(Current, In, Result) :-
    ( prime(Current) -> Out is In+Current ; Out=In ),
    NewCurrent is Current + 2,
    p010(NewCurrent, Out, Result).

prime(2).
prime(3).
prime(X) :-
    integer(X),
    X > 3,
    X mod 2 =\= 0,
    \+is_composite(X, 3).   % was: has_factor(X, 3)

is_composite(X, F) :-       % was: has_factor(X, F) 
    X mod F =:= 0, !.
is_composite(X, F) :- 
    F * F < X,
    F2 is F + 2,
    is_composite(X, F2).

Отказ от ответственности: я нашел эту реализацию prime/1 и has_factor/2 путем поиска в Google.

Этот код дает:

?- problem_010(R).
R = 142913828922
Yes (12.87s cpu)

Вот еще более быстрый код:

problem_010(R) :-
    Max = 2000001,
    functor(Bools, [], Max),
    Sqrt is integer(floor(sqrt(Max))),
    remove_multiples(2, Sqrt, Max, Bools),
    compute_sum(2, Max, 0, R, Bools).

% up to square root of Max, remove multiples by setting bool to 0
remove_multiples(I, Sqrt, _, _) :- I > Sqrt, !.
remove_multiples(I, Sqrt, Max, Bools) :-
    arg(I, Bools, B),
    (
        B == 0
    ->
        true % already removed: do nothing
    ;
        J is 2*I, % start at next multiple of I
        remove(J, I, Max, Bools)
    ),
    I1 is I+1,
    remove_multiples(I1, Sqrt, Max, Bools).

remove(I, _, Max, _) :- I > Max, !.
remove(I, Add, Max, Bools) :-
     arg(I, Bools, 0), % remove multiple by setting bool to 0
     J is I+Add,
     remove(J, Add, Max, Bools).

% sum up places that are not zero
compute_sum(Max, Max, R, R, _) :- !.
compute_sum(I, Max, RI, R, Bools) :-
    arg(I, Bools, B),
    (B == 0 -> RO = RI ;  RO is RI + I ),
    I1 is I+1,
    compute_sum(I1, Max, RO, R, Bools).

Это работает на порядок быстрее, чем код, который я дал выше:

?- problem_010(R).
R = 142913828922
Yes (0.82s cpu)

Answer 4

ОК, до редактирования проблема была только в алгоритме (imho).Как вы заметили, эффективнее сначала проверить, делится ли число на меньшие простые числа;в конечном наборе число делится на 3 больше, чем на 32147.

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

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

Теперь, что такое список различий?Вместо создания обычного списка [1,2,3] вы создаете этот [1,2,3 | T].Обратите внимание, что мы оставляем хвост необоснованным.Затем, если мы хотим добавить один элемент (или несколько) в конец списка, мы можем просто сказать T = [4 | NT].awesome?

Следующее решение (накопление простых чисел в обратном порядке, остановка при простом> sqrt (N), списки различий для добавления) занимает 0,063 для простых чисел 20k и 17 секунд для простых чисел 2m, в то время как исходный код занимает 3,7 секунды для20k и append / 3 версии 1.3sec.

problem_010(R) :-
    p010(3, Primes, Primes),
    sumlist([2|Primes], R).
p010(2000001, _Primes,[]) :- !.                                %checking for primes till 2mil
p010(Current, Primes,PrimesTail) :-
    R is sqrt(Current),
    (
        prime(R,Current, Primes)
    ->  PrimesTail = [Current|NewPrimesTail]
    ;   NewPrimesTail = PrimesTail
    ),
    NewCurrent is Current + 2,
    p010(NewCurrent, Primes,NewPrimesTail).
prime(_,_, Tail) :- var(Tail),!.
prime(R,_N, [Prime|_Primes]):-
    Prime>R.
prime(_R,N, [Prime|_Primes]) :-0 is N mod Prime, !, fail.
prime(R,ToTest, [_|Primes]) :- prime(R,ToTest, Primes).

также, рассматривая добавление чисел во время их генерации, чтобы избежать лишних o (n) из-за sumlist / 2
в конце концов, вывсегда можно реализовать алгоритм AKS, который выполняется за полиномиальное время (XD)