mardi 21 juillet 2015

Le serpent Vietnamien et profilage sous Scilab

voici un  petit problème de mathématiques qui a attiré l'attention de nombreux internautes il y a 2 mois de cela : il s'agit de trouver comment placer  les nombres de 1 à 9 (une et une seule fois chacun)  dans les 9 cases  vides  du "serpent" ci-dessous (notées a(i))  pour que le résultat  de la formule obtenue  soit égal à  66 :

a(1)+13*a(2)/a(3)+a(4)+12*a(5)-a(6)-11+a(7)*a(8)/a(9)-10=66
les 128 solutions du serpent vietnamien
L'origine du buzz autour de ce casse tête  vient du fait que cet exercice aurait été donné  par un instituteur Vietnamien  à ses élèves de 8 ans, élèves qui seraient donc en classe de CE2 en France. Etant donné l'effort considérable qu'ont du fournir  les quelques adultes à avoir essayé de trouver "manuellement" une solution au problème (qui en comporte exactement 128)  ce puzzle a vite été considéré  comme une preuve des faiblesses de notre système éducatif en mathématiques  !?!?

On pourrait discuter de l'intérêt pédagogique de cet exercice, dont le but est probablement plus de faire manipuler les opérations arithmétiques de base à de jeunes élèves plutôt que de trouver un processus de résolution exhaustif, mais cet exercice m'a de suite fait penser à l'anecdote du calcul de la somme des entiers de 1 à n par Carl Friedrich Gauss : à l'âge de neuf ans Gauss aurait stupéfié son maître d'école en calculant très rapidement la somme des entiers de 1 à 100, alors que le maître s'attendait à ce que ce calcul occupe toute la classe un long moment.

Un peu à la manière d'un Gauss des temps modernes ne serait-il pas plus intéressant de contourner le problème ... mais en utilisant un petit programme pour le résoudre? L'idée de base est d'utiliser la rapidité d'un ordinateur pour tester les 9!=362880 combinaisons possibles et trouver ainsi toutes les solutions. C'est ce qu'on appelle dans le jargon informatique une recherche par Force Brute (souvent anglicisée en "brute force"). Cela peut apparaître comme une inadmissible tricherie mais, tout comme le calcul de Gauss, la mise en place de cette stratégie s'avère être une véritable méthode mathématique très utilisée en théorie des graphes !

premiers niveaux de l'arbre  représentant  les différentes combinaisons 

En effet pour construire toutes les combinaisons possibles il faut trouver toutes les manières de ranger les valeurs de 1 à 9. Si chaque combinaison est représentée par une liste [a(1), a(2), a(3), a(4), a(5), a(6), a(7), a(8), a(9)] la construction de cette liste peut se visualiser sur un arbre dont :
  • le premier niveau consiste à choisir une des 9 valeurs possibles pour a(1),
  • au niveau suivant choisir une des 8 valeurs restantes pour a(2) ,
  • ..., 
  • au dernier niveau il ne reste qu'une valeur pour a(9)
Construire et tester toutes les combinaisons possibles du problème revient donc à parcourir toutes les branches de cet arbre, or la théorie des graphes nous fournies de nombreuses méthodes efficaces pour parcourir un graphe. Ici je vais utiliser l'algorithme de parcours en profondeur qui a la particularité de s'écrire de manière très concise sous forme récursive  dans le cas des arbres :

fonction Parcours_profondeur(A,a)
             A est un arbre (orienté)  et a un sommet de A
             Pour chaque b successeur de a dans A faire
                      Parcours_profondeur(A,b)
             fin_faire

Dans notre cas on peut se passer de la variable "A"  (qui est l'arbre fixé plus haut) , le sommet "a " représente  la combinaison en cours de construction. Il faut ajouter à cette trame une variable pour stocker les solutions  trouvée, qui doit être passée (et récupérée) d'appels en appels, et  la vérification de la formule , qui se fait quand on arrive aux feuilles de l'arbre.  Ceci s'écrit très facilement en langage scilab, où l'on peut  simuler une pile à l'aide de tableau à taille variable (grâce à l'allocation dynamique de la mémoire ). En voici une implémentation assez courte :

function sol=brute_force(a,sol)
    // a : liste des 9 valeurs (tableau 1x9)
    // sol : liste des solution (tableau nx9)
    // utilisation : sol=brute_force([],[]);
    if length(a)==9 // la liste a est remplie
        // on teste  si la combinaison donne bien 66
        val=a(1)+13*a(2)/a(3)+a(4)+12*a(5)-a(6)-11+a(7)*a(8)/a(9)-10 
        if val==66 then
            sol=[sol;a]
        end
    else //on continue de remplir la liste 
        valeurs_restantes=1:9
        valeurs_restantes(a)=[]
        for k=valeurs_restantes
            b=[a,k]//  ajout d'une nouvelle valeur  en fin de a
            sol=brute_force(b,sol)
        end
    end
endfunction

avec cette fonction la résolution du problème a pris environ 10 secondes sur mon portable ce qui sans être très long est quand même conséquent.  Je me suis alors demandé s'il était possible d'optimiser un peu ma fonction ?  Pour ce faire scilab dispose de commandes permettant d'analyser le coût des calculs lors de l'exécution d'une fonction , on parle de profilage (profiling en anglais).   Pour profiler  une fonction il faut d'abord indiquer à scilab  qu'on souhaite analyser son exécution avec la commande add_profiling :

add_profiling("brute_force")  //    prépare l'analyse de la  fonction
sol=brute_force([],[]);       //  exécution "normale" de la fonction

La fonction une fois "profilée"  s'exécute apparemment plus lentement (17sec dans mon cas) à cause de la récupération des données. Il existe ensuite plusieurs  manières d'éditer les données de profilage, la commande profile(brute_force)  permet de récupérer les données brutes sous forme d'un tableau à 3 colonnes :donnant respectivement pour chaque ligne de la fonction :
  • le nombre d'appels effectif, 
  • le temps CPU  cumulé , 
  • un indicateur  du coût de calcul pour l'interpréteur  
Le temps CPU  correspond au temps  nécessaire au calcul pour une machine à un seul cœur  exécutant ce seul processus. Sur des machines à plusieurs cœurs  et exécutant plusieurs processus en parallèles ce temps va être notablement différent  du temps réel mais il évalue plus précisément  le travail demandé au processeur. Ces données une fois récoltées peuvent être affichées de manière plus conviviale avec :
  • showprofile(brute_force)   qui permet d'afficher les données de profilage sous forme de chaînes de caractères dans la console
  • plotprofile("brute_force")   qui permet d'afficher les données de profilage  sous forme de figures dans la fenêtre graphiques
voilà ce que donne ces deux   commandes :

-->showprofile(brute_force)
|986410|0.84  |0 |function sol=fun(a, sol)
|362880|0.72  |0 |  if length(a) == 9 then
|362880|2.21  |44|    val = a(1) + 13 * a(2)/a(3) + a(4) + 12 * a(5)
                            - a(6) - 11 + a(7) * a(8)/a(9) - 10
|128   |0     |0 |    if val == 66 then
|128   |0     |4 |      sol = [sol;a]
|362880|0.45  |0 |    end,
|623530|1.22  |0 |  else
|623530|1.2   |4 |    valeurs_restantes = 1:9
|623530|0.93  |3 |    valeurs_restantes(a) = []
|986409|1.34  |0 |    for k = valeurs_restantes,
|986409|1.62  |4 |      b = [a,k]
|986409|135.57|5 |      sol = brute_force(b, sol)
|623530|0.56  |0 |    end,
|986410|0.79  |0 |  end,
|986410|0.73  |0 |endfunction


on constate immédiatement que ce qui prend  le plus de temps c'est :
  • le calcul de la somme  ligne 3  (indice 44) pour la vérification de la condition
    a(1)+13*a(2)/a(3)+a(4)+12*a(5) -a(6)-11+ a(7)*a(8)/a(9)-10=66 
  •  la répétition des 986410 appels récursifs ligne 12 (temps CPU 135s  sur un totale de 148s)   ce chiffre me semble un peu bizarre je m'attendais à voir apparaître  9!  (et étonnamment on a que 986410~ 9! *e  avec e=2.71828...   la base des exponentielles)
On remarque au passage que le stockage des solutions  lignes 4 et 5 n'est exécuté que 128 fois  exactement le nombre de solution (ouf!). Pour accélérer ce programme il faut donc essayer de réduire le nombre de combinaisons testées. Une idée simple pour y arriver est de remarquer que dans la formule  il y a 2 divisions  dans la formule à évaluer  ligne 3, et le  résultat global des deux expressions  comportant des divisions , a(2)/a(3)  et  a(7)*a(8)/a(9), doit forcément être entier si on veut que le résultat final donne 66.  Une optimisation simple consiste donc à choisir les 5 coefficients a(2), a(3), a(7), a(8), a(9)  puis à vérifier que 13*a(2)/a(3)+ a(7)*a(8)/a(9)  est bien entier avant de construire le reste de la combinaison. Voici ce que donne ma nouvelle fonction brute_force1 après profilage  :

-->showprofile(brute_force1)
|96682|0.09 |0 |function sol=fun(a, sol)
|29232|0.06 |0 |  if length(a) == 9 then
|29232|0.08 |20|    a = a([6,1,2,7,8,9,3,4,5])
|29232|0.18 |44|    val = a(1) + 13 * a(2)/a(3) + a(4) + 12 * a(5) 
                          - a(6) - 11 + a(7) * a(8)/a(9) - 10
|128  |0    |0 |    if val == 66 then
|128  |0    |4 |      sol = [sol;a]
|29232|0.04 |0 |    end,
|67450|0.14 |0 |  else
|15120|0.03 |0 |    if length(a) == 5 then
|15120|1.54 |27|     bool = pmodulo(13*a(1)/a(2)+a(3)*a(4)/a(5),1)==0
|52330|0.17 |2 |    else   bool = %t
|67450|0.06 |0 |    end,
|53548|0.05 |0 |    if bool then
|53548|0.1  |4 |      valeurs_restantes = 1:9
|53548|0.08 |3 |      valeurs_restantes(a) = []
|96681|0.14 |0 |      for k = valeurs_restantes,
|96681|0.16 |4 |        b = [a,k]
|96681|23.31|5 |        sol = brute_force1(b, sol)
|53548|0.05 |0 |      end,
|67450|0.06 |0 |    end,
|96682|0.08 |0 |  end,
|96682|0.07 |0 |endfunction


le résultat est impressionnant, même si la ligne ajoutée consomme pas mal de ressources  (1.54s en temps CPU et un indice de 27  pour le travail de l'interpréteur)  elle a permis de réduire   le nombre de combinaisons testées, et donc le  nombre d'appels récursif, par un facteur 10, et on a toujours bien les 128  solutions  du problème! Le temps réel d'exécution  est passé à 2 secondes, c'est bien mais c'est encore très loin du temps mis par un programme écrit en langage C  (quelques dixièmes de secondes sans optimisations). Si on avait envie de tester l'existence de solutions où les a(i) peuvent prendre des valeurs égales (toujours entre 1 et 9) on aurait alors
$$9^9 =387420489 \approx 1067.627\times 9!$$
ce qui  nécessiterait environ 1 heure de calculs ! On aurait pu améliorer encore un peu la rapidité du programme en effectuant des optimisations  comme :
  • déclarer comme variable globale  la ligne valeurs_restantes = 1:9  pour éviter  de la recalculer à chaque fois
  • modifier la formule a(1)+13*a(2)/a(3)+a(4)+12*a(5) -a(6)-11+ a(7)*a(8)/a(9)-10  pour éviter d'avoir à remettre dans l'ordre les coefficients  de la liste a avec la commande a = a([6,1,2,7,8,9,3,4,5])
mais le gain pratique sera extrêmement faible au regard  de la réécriture du programme qui deviendra beaucoup plus obscur ...

Aucun commentaire:

Enregistrer un commentaire

Pour écrire des formules mathématiques vous pouvez utiliser la syntaxe latex en mettant vos formules entre des "dollars" $ \$....\$ $ par exemple :
- $\sum_{n=1}^\infty {1\over n^2}={\pi^2\over 6}$ s'obtient avec \sum_{n=1}^\infty {1\over n^2}={\pi^2\over 6}
- $\mathbb R$ s'obtient avec {\mathbb R} et $\mathcal D$ s'obtient avec {\mathcal D}
- pour les crochets $\langle .,. \rangle$ dans les commentaires utilisez \langle .,. \rangle
vous pouvez écrire du html dans les commentaires :
- italique <i> ... </i> gras <b> ... </b>
- lien <a href="http://adresse "> .... </a>