Quel est le goulot d'étranglement dans ce prédicat lié aux primes ?

Question

Alors voilà : J'essaie de calculer la somme de tous les nombres premiers inférieurs à deux millions (pour ce problème ), mais mon programme est très lent. Je sais que l'algorithme en lui-même est terriblement mauvais et qu'il s'agit d'un algorithme de force brute, mais il me semble beaucoup plus lent qu'il ne le devrait.
Ici, je limite la recherche à 20 000 afin que le résultat ne soit pas attendu trop longtemps.
Je ne pense pas que ce prédicat soit difficile à comprendre mais je vais quand même l'expliquer : je calcule la liste de tous les nombres premiers inférieurs à 20 000 et je les additionne. La partie somme est bonne, la partie nombres premiers est très lente.

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).

J'aimerais avoir une idée de la raison pour laquelle il est si lent. S'agit-il d'une bonne implémentation de l'algorithme stupide de force brute, ou y a-t-il une raison quelconque qui fait que Prolog tombe ?

EDIT : J'ai déjà trouvé quelque chose, en ajoutant les nouveaux nombres premiers au lieu de les laisser en tête de liste, j'ai des nombres premiers qui apparaissent plus souvent au début, donc c'est ~3 fois plus rapide. J'ai encore besoin d'un peu d'aide :)

Answer 1

Premièrement, Prolog n'échoue pas ici.

Il existe des méthodes très intelligentes pour générer des nombres premiers. Mais pour commencer, il suffit d'accumuler les nombres premiers dans l'ordre inverse ! (7,9s -> 2,6s) De cette manière, les plus petits sont testés plus tôt. Ensuite, envisagez de ne tester que les nombres premiers jusqu'à 141. Les nombres premiers plus grands ne peuvent pas être un facteur.

Ensuite, au lieu de passer uniquement par les nombres non divisibles par 2, vous pourriez ajouter 3, 5, 7.

Il y a des gens qui écrivent des articles sur ce "problème". Voir, par exemple cet article Bien que la question de savoir ce qu'était l'algorithme "authentique" soit un peu sophistiquée, il y a 22 siècles, lorsque la dernière version du boulier était célébrée comme Comprimés de Salamis .

Answer 2

Tout d'abord, l'ajout à la fin d'une liste à l'aide de la fonction append/3 est assez lent. Si vous devez le faire, utilisez plutôt des listes de différences. (Personnellement, j'essaie d'éviter append/3 autant que possible)

Deuxièmement, votre prime/2 itère toujours sur l'ensemble de la liste lors de la vérification d'une prime. C'est inutilement lent. Vous pouvez plutôt vérifier si vous pouvez trouver un facteur intégral jusqu'à la racine carrée du nombre que vous voulez vérifier.

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).

Avertissement : J'ai trouvé cette mise en œuvre de prime/1 y has_factor/2 en googlant.

Ce code donne :

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

Voici un code encore plus rapide :

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).

Cela fonctionne un ordre de grandeur plus rapide que le code que j'ai donné ci-dessus :

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

Answer 3

Considérons par exemple la méthode du tamis ("tamis d'Eratosthène") : Créez d'abord une liste [2,3,4,5,6,....N], en utilisant par exemple numlist/3. Le premier nombre de la liste est un nombre premier, gardez-le. Éliminez ses multiples du reste de la liste. Le nombre suivant dans la liste restante est à nouveau un nombre premier. Éliminez à nouveau ses multiples. Et ainsi de suite. La liste va se réduire assez rapidement, et vous finirez par ne conserver que les nombres premiers.

Answer 4

OK, avant l'édition, le problème était juste l'algorithme (imho). Comme vous l'avez remarqué, il est plus efficace de vérifier d'abord si le nombre est divisé par les plus petits nombres premiers ; dans un ensemble fini, il y a plus de nombres divisibles par 3 que par 32147.

Une autre amélioration de l'algorithme consiste à arrêter la vérification lorsque les nombres premiers sont supérieurs à la racine carrée du nombre.

Maintenant, après votre changement, il y a effectivement quelques problèmes de prologue : vous utilisez append/3. append/3 est assez lent car vous devez parcourir toute la liste pour placer l'élément à la fin. Au lieu de cela, vous devriez utiliser des listes de différence, ce qui rend le placement de l'élément à la fin vraiment rapide.

Maintenant, qu'est-ce qu'une liste de différence ? Au lieu de créer une liste normale [1,2,3], vous créez celle-ci [1,2,3|T]. Remarquez que nous laissons la queue non renseignée. Ainsi, si nous voulons ajouter un élément (ou plus) à la fin de la liste, nous pouvons simplement dire T=[4|NT]. génial ?

La solution suivante (accumuler les nombres premiers dans l'ordre inverse, s'arrêter lorsque le nombre de nombres premiers>sqrt(N), liste des différences à ajouter) prend 0.063 pour 20k nombres premiers et 17sec pour 2m nombres premiers alors que votre code original prenait 3.7sec pour 20k et la version 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).

également, envisager d'ajouter les nombres pendant que vous les générez pour éviter le o(n) supplémentaire à cause de sumlist/2
au final, vous pouvez toujours implémenter l'algorithme AKS qui fonctionne en temps polynomial (XD)