Mon aventure sur SPOJ et autres CodeForces

Son programme doit sortir la liste des nombres premiers entre $a$ et $b$. Il s’agit donc d’un crible d’Eratosthène jusque $b$, duquel il ignorera les valeurs inférieures à $a$.

Ce que je dis, c’est que pour générer par exemple le crible des nombres jusque 50, il peut en effet limiter sa recherche des multiples à partir de 7, car il sait qu’aucun nombre de l’intervalle $[8, 50]$ ne sera multiple d’au moins deux nombres premiers de ce même intervalle.

Mes messages portent juste sur le fait que l’amélioration en question peut être directement apportée à prime_sieve, pas besoin d’avoir une fonction englobante pour cela et de répéter deux fois l’opération.

Edit : Bon, vous commenciez quand même à me faire douter, alors j’ai vérifié.

Ce que je propose, donc, c’est de remplacer prime_sieve par une fonction similaire à la suivante :

1
2
3
4
5
6
7
8

def prime_sieve(start, end):
    sieve = [False, False] + [True] * (end - 2)
    for i in range(2, int(end**0.5)):
        if not sieve[i]:
            continue
        for j in range(i*2, end, i):
            sieve[j] = False
    return [i for i, prime in enumerate(sieve) if prime and i >= start]

J’ai peut-être juste usé d’un abus de langage sur « crible d’Eratosthène », je pensais qu’on pouvait qualifier ainsi toute méthode générant une liste de nombres premiers inférieurs à une borne maximum, alors qu’il semble s’agir d’une méthode bien particulière (qui reste plus ou moins celle utilisée dans mon exemple, si ce n’est que la recherche des multiples est stoppée plus tôt).

Merci pour le code en exemple, ça m’a permis de comprendre ce que tu voulais dire. Et c’est effectivement ce à quoi je pensais quand je disais :

@entwanne, tu es sûr que ce que tu proposes n’est pas au moins linéaire en la valeur de la borne supérieure de l’intervalle ?

Je ne sais pas si tu as lu le sujet (c’est vrai que @Ricocotam n’a pas rappelé les contraintes), mais end peut aller jusqu’à $10^9$, donc l’algo que tu proposes est trop lent (et prend trop de mémoire).

On peut s’amuser à calculer la complexité temporelle exacte de ta fonction prime_sieve executée pour a = start et b = end (en notant $P_x$ l’ensemble des premiers $\le x$) :

L’avantage de la méthode de @Ricocotam est qu’elle ne remplit les cases du tableau du crible que jusqu’à $\sqrt{b}$, donc sa complexité doit être quelque chose comme (cf. théorème des nombres premiers) :

Ceci dit, ta méthode n’est pas inintéressante et c’est plutôt sur cette base qu’il faudrait partir pour résoudre ce problème proprement (je crois que des algos analogues à celui de @Ricocotam peuvent passer sur le sujet du SPOJ s’ils sont optimisés à la main, mais on peut faire mieux en terme de complexité).

J’ai peut-être juste usé d’un abus de langage sur « crible d’Eratosthène », je pensais qu’on pouvait qualifier ainsi toute méthode générant une liste de nombres premiers inférieurs à une borne maximum, alors qu’il semble s’agir d’une méthode bien particulière (qui reste plus ou moins celle utilisée dans mon exemple, si ce n’est que la recherche des multiples est stoppée plus tôt).

Pour moi y a pas d’abus de langage, ce que tu as écrit est bien un crible d’Eratosthène.

Non, je ne l’ai pas lu, j’ai demandé à ce qu’il explique le problème et ses contraintes, et pensais qu’il l’avait ensuite fait.

L’algo que je propose est la reprise du sien, si j’avais à résoudre ce problème, je ne procèderai pas de la sorte et utiliserai d’autres structures de données.

Ce qui me gène dans son code initial, c’est que pour les nombres entre $start$ et $\sqrt{end}$, déjà présents dans le tableau primes, il s’amuse à en chercher les diviseurs. Alors qu’on sait ces nombres premiers !

Il a effectivement mis un cas particulier pour les nombres entre start et sqrt(end) , mais ce cas est marginal. Dans l’exercice on nous dit qu’on a jusqu’à 10 intervalles [a,b] à traiter. Chaque intervalle a une amplitude b-a inférieure à 100000, et b peut aller jusqu’à 10^9

Donc la probabilité d’avoir a >= sqrt(end) est très faible. Idem, la probabilité que 2 intervalles se superposent est très faible.

Pour moi, il faut : - 1. Calculer end = le max des 10 plafonds. - 2. Prendre sqrt(end) - 3. Recenser tous les nombres premiers de 2 jusqu’à ce nombre. - 4. Dans chaque intervalle, rechercher les nombres premiers, en s’appuyant sur le tableau construit en 3.

Dans cette dernière étape, on a intérêt à faire une boucle du genre for i in range ( start, end, 2).

Ca permet d’économiser beaucoup de divisions par 2.

Ou encore :

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12

st30  = math.ceil(i/30)*30
for i in range( start, st30) 
   teste (i)
for i in range( st30, end, 30) 
   teste (i+1)
   teste (i+7)
   teste (i+11)
   teste (i+13)
   teste (i+17)
   teste (i+19)
   teste (i+23)
   teste (i+29)

Et dans la fonction teste(), on teste si i est divisible par un nombre premier entre 7 et …

On économise 3 Millions de divisions.

J’ai voulu mettre en œuvre les différentes solutions, pour voir un peu ce qu’elles avaient dans le ventre. La solution telle que je la propose est en moyenne 3 à 4 fois plus rapide que celle donnée dans le post initial du sujet.

Je n’ai pas testé au-delà de $10^6$, et aucune des deux n’est viable pour $10^9$.

Le quiproquo vient justement du fait que c’est un algo très différent que tu proposes, puisque faire le crible jusqu’à $\sqrt{b}$ c’est pas pareil que faire l’optimisation de la boucle extérieure du crible jusqu’à $\sqrt{b}$. Je dirais sans trop m’avancer que @Ricocotam a déjà pensé à ta solution, qui est la plus naturelle pour ce problème.

D’ailleurs si ta dernière remarque n’est pas juste une généralité, je serais intéressé de savoir quel point de vue/structures de données tu aurais choisi pour ce problème. C’est pas un défi hein, juste de la curiosité. (J’avoue que je vois pas trop ce qu’on peut faire de vraiment différent.)

Je rejoins @elegance sur ce point : c’est une optimisation qui me paraît négligeable. Calculer les diviseurs d’un nombre $\le 10^4$, c’est rapide par rapport au reste de l’algo.

J’ai voulu mettre en œuvre les différentes solutions, pour voir un peu ce qu’elles avaient dans le ventre. La solution telle que je la propose est en moyenne 3 à 4 fois plus rapide que celle donnée dans le post initial du sujet.

Comment tu choisis tes intervalles pour faire ta comparaison ? Parce que si tu vas jusqu’à $10^6$ tu respectes pas vraiment l’esprit du problème initial.

Pas d’accord avec ça. Je ne sais pas si, cette fois, tu as lu le sujet que j’ai posté dans mon dernier message, mais il y a la contrainte supplémentaire $b-a\le 10^6$ (en plus de $a\le b\le 10^9$). En ce sens, l’algo de @Ricocotam est certes assez limite en temps, mais il est testable pour $b\approx 10^9$, au contraire d’un algo qui alloue un tableau de taille $10^9$…

Ta méthode @elegance est bonne (c’est plus ou moins l’algo attendu), mais je trouve toujours un peu dommage de donner direct une solution sans laisser l’OP chercher.

J’ai pas le temps de développer mais j’avais effectivement mal compris @Entwanne, mea culpa. J’ai implémenté initialement ce que tu proposes mais c’est beaucoup trop lent. Je dois faire passer l’algo en moins de 6s avec comme pire cas cette entrée :

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11

10
999000000 1000000000
999000000 1000000000
999000000 1000000000
999000000 1000000000
999000000 1000000000
999000000 1000000000
999000000 1000000000
999000000 1000000000
999000000 1000000000
999000000 1000000000

Sinon pour la superposition d’intervalle, c’est probablement la raison de la mauvaise sortie de l’algo, merci, ça va m’aider

Je vous répondrais un peu plus en détail plus tard. J’ai essayé pas mal de chose avant d’arriver à ce code (il peut encore être amélioré, j’en proposerais une version finale dans un billet), notamment j’ai commencé par de la récursivité mais ça s’est mal passé

D’après l’énoncé, $n - m <= 100000$. L’écart dans ton exemple est 10 fois plus grand.

Grâce à ton sujet j’ai découvert SPOJ, j’ai donc fait le problème PRIME1 (en Go). Je n’ai pas suivi ton sujet par contre donc je ne sais pas si on a utilisé la même méthode. Je peux partager mon code si tu veux (MP ou ici ?).

Oh oui vas y ici pas de problème !

J’avais pas vu que l’écart était plus petit que ça ^^, mais bon mon algo fonctionne quand même en un temps correct avec cet écart

Je me permet de double poste puisque j’ai enfin réussit !

J’ai réglé divers problème et fait de belles améliorations grâce à vos remarques. Je vous poste le code ici en brut pour les curieux. Je vais essayer, maintenant que je comprend un peu mieux le problème et les pièges, de le refaire de manière récursive n’étant pas très à l’aise avec. Une fois ceci fait je ferais un billet (article ou tribune ?) en expliquant où on été les difficultés. Une fois que le billet sera publié je changerais le premier post.

Le code :

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56

# Prime1 on SPOJ - SpheOnlineJudge.com

import sys
import math

def prime_sieve(end):
    nb = [False, False] + [True for i in range(2, end + 1)]  # 0 and 1 to False
    i = 2
    max = math.ceil(math.sqrt(len(nb)))
    while i < max:
        if nb[i]:
            for j in range(2*i, len(nb), i):
                nb[j] = False
        i += 1
    return [n for n in range(len(nb)) if nb[n]]


def main():
    nb_test = int(sys.stdin.readline())
    cases = []
    max_end = 0
    for e in range(nb_test):
        start, end = map(int, sys.stdin.readline().split(" "))
        cases.append({"start": start, "end": end})
        if max_end < end:
            max_end = end

    primes = prime_sieve(math.ceil(math.sqrt(max_end)))

    for case in cases:
        start, end = case["start"], case["end"]
        if start == 1 or start == 2:
            start = 3
            if end >= 2:
                print(2)
            else:
                print()

        for i in range(start, end + 1):
            sqrt_i = math.ceil(math.sqrt(i))
            is_prime = True
            for p in primes:
                if p > sqrt_i:
                    break
                if i % p == 0:
                    is_prime = False
                    break
            if is_prime:
                print(str(i))
        print()  # New line and each test case


if __name__ == '__main__':
    main()

`

L’image de la réussite (3.9s)

Si vous voyez encore des améliorations possibles (je ne les détails pas puisque je le ferais d’ici un jour ou deux) n’hésitez pas.

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131

// http://www.spoj.com/problems/PRIME1/
// 0.49 8.6M
package main

import (
  "bufio"
  "fmt"
  "math"
  "os"
  "strconv"
  "strings"
)

type Interval struct {
  Min int
  Max int
}

func NewInterval(min, max int) Interval {
  return Interval{min, max}
}

func NewIntervalFromString(stringInterval string) Interval {
  interval := strings.Fields(stringInterval)

  min, _ := strconv.Atoi(interval[0])
  max, _ := strconv.Atoi(interval[1])

  return NewInterval(min, max)
}

type Primes struct {
  Map  map[int]struct{}
  Keys []int
}

func NewPrimes() Primes {
  return Primes{make(map[int]struct{}), []int{}}
}

func (p *Primes) Add(n int) {
  p.Map[n] = struct{}{}
  p.Keys = append(p.Keys, n)
}

func (p Primes) Contains(n int) bool {
  _, ok := p.Map[n]
  return ok
}

func (p Primes) IsUnknownPrime(n int) bool {
  if n <= 1 {
      return false
  }

  max := intSqrt(n)

  for _, prime := range p.Keys {
      if prime > max {
          break
      }

      // it isn't prime
      if floatIsInt(float64(n) / float64(prime)) {
          return false
      }
  }

  return true
}

func (p *Primes) Find(i Interval, print bool) {
  for j := i.Min; j <= i.Max; j++ {
      if p.Contains(j) {
          if print {
              fmt.Println(strconv.Itoa(j))
          }
      } else if p.IsUnknownPrime(j) {
          p.Add(j)

          if print {
              fmt.Println(strconv.Itoa(j))
          }
      }
  }
}

func main() {
  s := bufio.NewScanner(os.Stdin)

  // ignore the number of test cases, we don't need it
  s.Scan()

  var intervals []Interval
  for s.Scan() {
      intervals = append(intervals, NewIntervalFromString(s.Text()))
  }

  run(intervals)
}

func run(intervals []Interval) {
  primes := NewPrimes()
  primes.Find(NewInterval(2, intSqrt(biggestMax(intervals))), false)

  for k, interval := range intervals {
      if k != 0 {
          fmt.Println()
      }

      primes.Find(interval, true)
  }
}

func biggestMax(intervals []Interval) (max int) {
  for _, interval := range intervals {
      if interval.Max > max {
          max = interval.Max
      }
  }

  return max
}

func floatIsInt(n float64) bool {
  return n == float64(int(n))
}

func intSqrt(n int) int {
  return int(math.Sqrt(float64(n)))
}

Le fonctionnement est plutôt simple :

• primes est la liste contenant tout les nombres premiers que je connais.
• On remplit primes de tout les nombres premiers entre 2 et int(sqrt(max)) inclusif avec max le plus grand nombre de tout les intervalles.
• Ensuite, on boucle sur tout les intervalles et on affiche les nombres qui sont premiers.
• Pour savoir si n est premier, je regarde d’abord si primes contient n, sinon je boucle sur primes jusqu’à int(sqrt(n)) pour voir si ce n’est pas un nombre premier que je ne connais pas déjà.
• Le type Primes possède une map[int]struct{} pour un lookup rapide (cf Primes.Contains()) et un slice []int pour boucler sur les nombres premiers dans l’ordre (cf ligne 58).
• J’arrive à un temps je pense correct de 0.49s et une utilisation mémoire de 8.6M (avec un slice et sans map, j’étais dans les 0.60 et 6.1M).

Je suis sûr quand gardant le même algorithme il y a moyen de faire des micro-optimisations et qu’il existe d’autres algos plus rapides. Je ne sais même pas comment s’appelle celui que j’ai implémenté.

@Ricocotam : ton prog python prend combien de temps pour PRIME1 ? Tu as essayé de te frotter à PRINT  ?

@tleb : En gros t’as le même algo que @Ricocotam, sauf que tu génères les premiers jusqu’à sqrt(maxend) de manière plus naïve qu’un crible (l’optimisation consistant à rajouter les premiers que tu trouves dans la boucle for k, interval à la liste p est inutile dans le pire cas où en entrée tu as des intervalles disjoints et dont les bornes sont toutes très grandes).

Pour b = maxend, sans compter les opérations sur les map, ton premier appel à find est en $\Theta\left(\frac{b^{3/4}}{\ln b}\right)$, ce qui est moins bien que le crible en $\Theta(\sqrt{b}\ln\ln b)$, mais un peu mieux que l’algo naïf consistant à tester les divisibilités avec tous les nombres (pas forcément premiers) entre 0 et la racine, qui est en $\Theta(b^{3/4})$. Comme $b$ est pas très grand, en fait toutes ces méthodes se valent (le gros du temps est passé dans la phase "online", pas dans le précalcul). En clair, je pense que tu obtiendrais quasiment les mêmes perfs avec un algo comme :

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17

maxend = 10^9
primes = []
i = 2
while i * i <= maxend:
    j = 2
    isprime = true
    while j * j <= i:
        if j % i == 0:
            isprime = false
    if isprime: primes += i
for a, b in intervals:
    for i in {a, ..., b}:
        isprime = true
        for p in primes:
            if i % p == 0:
                isprime = false
        if isprime: print(i)

Dans les deux cas, la phase de recherche dans les intervalles est en $\Theta\left((b-a)\frac{\sqrt{b}}{\ln b}\right)$, ce qui est censé être un peu limite mais qui en fait passe sur SPOJ. L’algo que @elegance a plus ou moins décrit, qui consiste à adapter le crible en le faisant tourner exactement sur chaque intervalle, est en $\Theta((b-a)\ln\ln b)$, ce qui est bien mieux !

Si vous voulez vraiment benchmarker, pour vos 2 algos le pire cas est de la forme :

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11

10
990000000 991000000
991000000 992000000
992000000 993000000
993000000 994000000
994000000 995000000
995000000 996000000
996000000 997000000
997000000 998000000
998000000 999000000
999000000 1000000000

Comment vous avez fait pour être si rapide ? xD. Je suis en 3.9s et 27M de mémoire. Bon, il va falloir que je trouve comment faire mieux. Déjà, je vais oublier tout les nombres pairs que je test, c’est inutile. En revanche pour la partie mémoire je ne sais pas comment mieux faire :/

@Lucas-84 Mon temps sur ton benchmark : "Runtime error #stdin #stdout 2.3s 28384KB"

J’ai modifié mon algo, et l’ai réécrit en cpp. Visiblement, ce n’est pas mal du tout.

Je n’ai pas lu a description des algorithmes proposés ci-dessus, mais le mien est relativement simple :

• on pré-calcule a l’aide d’un crible tous les nombres premiers jusqu’à $\sqrt{MAX}$
• ensuite on emploie des cribles sur les intervalles en se servant de la liste des premiers connus (je ne me soucie pas de recalculer 2 fois le même intervalle, car vu l’énoncé c’est relativement négligeable)

Le reste consiste en diverses astuces d’implémentation pour aller plus vite.

J’avoue que j’allais essayé d’implémenter le crible sur l’ensemble à calculer, visiblement c’est plutôt pas mal alors. J’essaierai ça demain !

@Ricocotam

Dans ton code, pour chaque valeur de i entre start et end, tu calcules sqrt(i) puis pour chaque valeur de p premier, tu testes si p est inférieur à sqrt(i).

Donc en gros, tu calcules environ 1 million de fois une racine-carrée, et tu testes p<sqrt(i) environ 10 ou 20 millions de fois.

Tu peux remplacer tout ça par :

• Pour chacun des 10 intervalles à traiter, calculer sqrt(end)
• Faire un tableau des nombres premiers inférieurs à sqrt(end) extrait du tableau primes qui a été initialisé une fois pour toutes.
• Et du coup plus besoin des 2 opérations en question (sqrt(i) et tests d’infériorité)

Toutes les opérations que j’ai recensées au début seront donc remplacées par 10 calculs de racines carrées, et 10 opérations de copie de range, et quelques divisions en plus ( 100000, 200000 ? pas beaucoup plus)

J’ai besoin absolument de la comparaison en fait. C’est ce qui m’évite de me retrouver avec des intervalles qui sont confondus (pour les petits nombres c’est important). Il faudrait que je fasse différemment du coup..

Edit : J’ai répondu un peu vite à cause de la fatigue. Du coup je sais pas si ce que j’ai dit au dessus est valable ^^. En fait le problème c’est que si le nombre i est premier mais inférieur à sqrt(end), mon nombre va être détecté comme non premier. Je peux peut être cependant faire un test d’égalité.. Je test ça et je te dis. S’il n’y apas de réponse j’édit le message

Edit 2: Alors j’ai apporté la modif que tu décris, j’ai enlevé tout les nombres paires de ce que je test et je me prend du TLE. WTF

Mais ça ne coûtait rien d’ajouter aussi l’optimisation à prime_sieve, qui n’a pas besoin d’itérer jusque end.

Je ne l’ai pas fait, mais intuitivement je n’aurais pas fait un tableau de booléens mais plutôt un set contenant les non-premiers. Même si maintenant, après réflexion, c’est loin d’être la meilleure idée.

Je n’ai toujours pas pris/eu le temps de lire le sujet original de l’exercice, donc je ne connais pas l’esprit initial. J’ai pris $[1, 10^6]$ comme intervalle.

Quand je dis que ce n’est pas viable pour $10^9$, c’est que ça prend trop de temps à mon goût (je n’ai pas eu le courage de le laisser tourner).

@Entwanne J’ai fait ce que tu proposais dans le crible et cela fonctionne beaucoup plus rapidement, me suffisait de comprendre ce que tu disais

À ce propos, j’ai vu que tu avais repris le range avec i*2, mais il serait préférable de commencer à i**2.