#! /usr/bin/env python3
# -*- coding: utf-8; mode: python -*-
""" TP4 au Lycée Lakanal (sujet rédigé par Arnaud Basson).
- *Date :* mercredi 18-11 et jeudi 19-11, puis jeudi 03-11 et 17-11 (2015),
- *Auteur :* Lilian Besson, pour le cours d'informatique pour tous en prépa MP (http://perso.crans.org/besson/infoMP/),
- *Licence :* MIT Licence (http://lbesson.mit-license.org).
"""
from __future__ import print_function, division # Python 2 compatibility
# %% Exercice 1 : Suite définie par récurrence
print("\n\nExercice 1 : Suite definie par recurrence")
[docs]def suite_u(n):
""" Suite :math:`{(u_n)}_n` définie par récurrence :math:`u_0 = 1.0, u_n+1 = u_n + 1 / (2^{u_n})`."""
assert n >= 0
if n == 0:
return 1.0 # u_0 = 1.0
else:
x = suite_u(n - 1) # On ne le calcule qu'une fois
return x + 1.0 / (2 ** x)
print("Premier terme u_0 =", suite_u(0))
print("10 eme terme u_10 =", suite_u(10))
print("100 eme terme u_100 =", suite_u(100), "(= 6.18199... OK)")
print("101 eme terme u_101 =", suite_u(101))
# %% Exercice 2 : Exponentiation rapide
print("\n\nExercice 2 : Exponentiation rapide")
# 2.1.
[docs]def exponentiation_rapide(x, n):
r""" Calcule :math:`x^n` en ne faisant que :math:`O(\log n)` multiplications (flottantes ou entières, selon x)."""
if n < 0:
return exponentiation_rapide(1.0 / x, -n)
elif n == 0: # x**0 = 1 même pour x = 0
return 1
elif n == 1:
return x
elif n % 2 == 0: # n est pair
# x ** (2k) = (x**2)**k
return exponentiation_rapide((x * x), (n / 2))
elif n % 2 == 1: # n est impair
# x ** (2k+1) = x * (x**2)**k
return x * exponentiation_rapide((x * x), ((n - 1) / 2))
else:
raise ValueError("exponentiation_rapide(x, n) valeur invalide pour n.")
print("Exemples avec l'exponentiation rapide :")
for x in [0, 1, 2, 4, 100, 3.1415]:
for n in [0, 1, 2, 3, 10]:
v = exponentiation_rapide(x, n)
assert abs(v - (x**n)) < 1e-9
print(x, "**", n, "=", v)
# 2.2.a.
[docs]def naive_fibonacci(n):
""" Fibonacci(n) :math:`F_n` en :math:`O(2^n)` opérations (et une pile d'appel de taille :math:`O(2^n)`).
.. warning:: C'est très inefficace.
``naive_fibonacci(40)`` prend de l'ordre de :math:`O(2^{40})` operations élémentaires :
- Soit au moins 15 minutes, si on dit qu'un CPU à 3.1 GHz fait :math:`10^9` opérations par seconde,
- ou plutôt de l'ordre d'une journée de calcul avec Python... Beaucoup trop quoi !
"""
assert n >= 0
if n == 0:
return 0
elif n == 1:
return 1
else:
return naive_fibonacci(n - 1) + naive_fibonacci(n - 2)
print("Exemples avec Fibonacci naif :")
for n in [0, 1, 2, 3, 4, 5, 10, 20, 25]:
print("Fibo", n, "=", naive_fibonacci(n))
assert naive_fibonacci(10) == 55
# 2.2.b.
import numpy as np
A = np.array([[0, 1], [1, 1]], dtype=np.int64)
X0 = np.array([0, 1], dtype=np.int64)
# 2.2.c.
[docs]def puissMat(A, n):
r""" Calcule la matrice :math:`A^n` en ne faisant que :math:`O(\log n)` multiplications de matrices
Note : chaque multiplication de matrice est en :math:`O(k^3)` si :math:`A` est de taille :math:`k \times k`.
"""
if n == 0: # A**0 = 1 même pour A = 0
return 1
elif n == 1:
return A
elif n % 2 == 0: # n est pair
# A ** (2k) = (A**2)**k
return puissMat(np.dot(A, A), (n / 2))
elif n % 2 == 1: # n est impair
# A ** (2k+1) = A * (A**2)**k
return np.dot(A, puissMat(np.dot(A, A), ((n - 1) / 2)))
else:
raise ValueError("puissMat(A, n) valeur invalide pour n.")
[docs]def fibo(n):
r""" Fibonacci(n) :math:`F_n` en :math:`O(\log n)` opérations (additions et multiplications de matrices d'entiers 64 bits, de tailles :math:`2 \times 2`).
.. note::
- Utiliser des matrices numpy demande d'utiliser un type de données accepté par numpy, comme ``np.int32`` (entiers 32 bits) ou ``np.int64``, et qui sont limités en taille (jusqu'à :math:`2^{31}` ou :math:`2^{63}` respectivement).
- Donc cette implémentation est exacte tant que :math:`F_n` ne dépasse pas cette capacité maximale.
- On peut corriger soit en réimplémentant nous même le produit matriciel (et ne plus utiliser :py:func:`np.dot`).
"""
An = puissMat(A, n)
Xn = np.dot(An, X0)
return Xn[0]
print("Exemples avec Fibonacci matriciel (fibo) :")
for n in [0, 1, 2, 3, 4, 5, 10, 20, 25, 30, 35, 45]:
print("Fibo", n, "=", fibo(n))
assert fibo(10) == 55
print("Pour de grande valeurs, ça foire :")
for n in [100, 200, 1000]:
print("Fibo", n, "=", fibo(n))
# 2.2.d.
from math import sqrt
phi = (1 + sqrt(5)) / 2
psi = (1 - sqrt(5)) / 2
[docs]def fibo2(n):
r""" Fibonacci(n) :math:`F_n` en :math:`O(1)` opérations (additions et exponentiations **flottantes**).
.. math:: F_n = \frac{1}{\sqrt{5}} \Big( {\left( \frac{1 + \sqrt{2}}{2} \right)}^n - {\left( \frac{1 - \sqrt{2}}{2} \right)}^n \Big)
.. note::
- utiliser des flottants entraîne des erreurs de calcul dues aux arrondis (``fibo2(n)`` n'est presque jamais exacte)
- on pourrait utiliser le module :py:mod:`decimal` pour pouvoir fixer la précision (à 100 chiffres après la virgule au lieu de 14), ou :py:mod:`sympy.mpmath`.
"""
fn = ((phi ** n) - (psi ** n)) / sqrt(5)
return int(fn) # On force à renvoyer un entier
print("Exemples avec Fibonacci flottant via sa formule exacte (fibo2) :")
for n in [0, 1, 2, 3, 4, 5, 10, 20, 25, 30, 35, 45]:
print("Fibo", n, "=", fibo2(n))
assert fibo2(10) == 55
# 2.2.e
print("Comparons ces deux fibo et fibo2 pour n = 78 et n = 95 :")
for n in [78, 95]:
print("Pour n =", n, "fibo(n) =", fibo(n), "et fibo2(n) =", fibo2(n))
# %% Exercice 3 : Version récursive de la recherche dichotomique
print("\n\nExercice 3 : Version recursive de la recherche dichotomique")
print("Exemple de t et x :")
t = [-100005, -451, -78, -54, -6, 0, 7, 9, 10, 25, 67, 31, 78, 911, 5000, 10000, 746154879]
x = 911
print("t =", t, "\nx =", x)
[docs]def rechercheDichoRecopiages(t, x):
r""" (Mauvaise) recherche dichotomique de la valeur x dans le tableau (liste) trié t.
- *entrée* : un tableau d'entiers t trié en ordre croissant et un entier x,
- *sortie* : un indice i tel que t[i]=x ou -1 si x n'est pas dans t.
Cette première implémentation était fausse (il faut ajouter ``k+1`` si on recherche à droite),
et inefficace car le recopiage d'un sous-tableau (``t[:k]`` ou ``t[k+1:]``) prend un temps non-négligeable (:math:`O(n/2)` en l'occurrence),
donc la recherche n'était pas logarithmique en n (la longueur du tableau) mais linéaire en n (:math:`n/2 + n/4 + \dots + 1`)...
"""
if len(t) == 0: # Cas de base
return -1
k = len(t) // 2 # Indice du milieu du tableau
if t[k] == x: # Autre cas de base
return k
elif t[k] < x:
# Recherche dans la moitié droite
return k + 1 + rechercheDichoRecopiages(t[k + 1:], x)
else:
# Recherche dans la moitié gauche
return rechercheDichoRecopiages(t[:k], x)
print("rechercheDichoRecopiages(t, x) =", rechercheDichoRecopiages(t, x))
print("Cette implementation etait fausse, et mauvaise car utilise des recopiages inutiles.")
# Vérification
i = rechercheDichoRecopiages(t, x)
assert t[i] == x
j = t.index(x)
assert i == j
[docs]def rechercheDichoRec(t, x, g, d):
""" Recherche dichotomique récursive de la valeur x dans le tableau (liste) trié t.
- *entrée* : un tableau d'entiers t trié en ordre croissant et un entier x, ainsi que g et d deux entiers (on recherche dans t[g:d+1] = t[g], .., t[d]),
- *sortie* : un indice i tel que t[i]=x ou -1 si x n'est pas dans t.
"""
if (len(t) == 0) or (g > d): # Cas de base
return -1
k = (d + g) // 2 # Indice du milieu du tableau t[g:d+1]
if t[k] == x: # Autre cas de base
return k
elif t[k] < x:
# Recherche dans la moitié droite
return rechercheDichoRec(t, x, k + 1, d)
else:
# Recherche dans la moitié gauche
return rechercheDichoRec(t, x, g, k - 1)
[docs]def rechercheDicho(t, x):
""" Bonne recherche dichotomique de la valeur x dans le tableau (liste) trié t.
- *entrée* : un tableau d'entiers t trié en ordre croissant et un entier x,
- *sortie* : un indice i tel que t[i]=x ou -1 si x n'est pas dans t.
"""
return rechercheDichoRec(t, x, 0, len(t) - 1)
print("Meilleure implementation, effectivement en O(log len(t)) :")
print("rechercheDicho(t, x) =", rechercheDicho(t, x))
# Vérification
i = rechercheDicho(t, x)
assert t[i] == x
j = t.index(x)
assert i == j
# %% Exercice 4 : Courbe du dragon
print("\n\nExercice 4 - Courbe du dragon")
# 4.a.
[docs]def inverse(t):
""" Renvoie le tableau *"inverse"* de t (renverse puis change les signes)."""
# return [-x for x in reversed(t)] # Un peu plus rapide !
n = len(t)
rev_t = []
for k in range(n - 1, -1, -1):
x = t[k] # On part de la fin pour renverser
rev_t.append(-x) # On inverse le signe !
return rev_t
print("Inverse d'un tableau :")
print("Pour t = [1, 1, -1], inverse(t) =", inverse([1, 1, -1]), "(comme voulu).")
# 4.b.
[docs]def listeRec(n):
""" Implémentation récursive de cette liste :math:`L_n` (de taille :math:`2^{(n+1)} - 1`)."""
assert n >= 0
if n == 0:
return [1]
else:
l = listeRec(n - 1)
# On peut concaténer des listes facilement comme ça :
return l + [1] + inverse(l)
for n in [0, 1, 2, 3]:
print("listeRec", n, "= ", listeRec(n))
# 4.c.
# import numpy as np
import matplotlib.pyplot as plt
[docs]def dragon(n, sauvegarde=False):
""" Trace la courbe (récursive) du dragon à l'ordre n (< 22 sinon illisible)."""
angles = listeRec(n)
r = len(angles)
# Création de la liste des affixes complexes M
listeM = [(0 + 0j), (1 + 0j)] # listeM = [0, 1] marche aussi
for k in range(0, r):
# Nouveau point complexe z_{k+2} = z_{k+1} + angles_k i (z_{k+1} - z_k)
zpred = listeM[k]
z = listeM[k + 1]
# print("Pour k =", k, "zpred =", zpred, "et z =", z, "et angles[k] =", angles[k])
nextz = z + angles[k] * 1j * (z - zpred)
listeM.append(nextz)
# Conversion en deux listes de nombres réels
listeX = [u.real for u in listeM]
listeY = [u.imag for u in listeM]
# Affichage de la courbe Dragon(n)
plt.figure() # Création d'une nouvelle figure
plt.plot(listeX, listeY, 'r-') # Tracé de la courbe
plt.xticks([]) # Aucune graduations sur l'axe des x
plt.yticks([]) # Aucune graduations sur l'axe des y
plt.axis('equal') # Pour imposer un repère orthonormal
plt.title('Le dragon fractal Mushu (pour n = {})'.format(n))
plt.show() # On affiche la figure, pour être sûr.
if sauvegarde: # Sauvegarde vers une image PNG
fichier = "TP4_courbe_du_dragon_n={}.png".format(n)
print("Imprime la figure dans le fichier", fichier)
plt.savefig(fichier)
return listeM
if __name__ == '__main__':
print("Demonstration de la courbe du dragon pour n = 3, 4, 5, 15 : (voir graphiques plus bas)")
# dragon(3)
# dragon(4)
# dragon(5)
# dragon(15)
# %% Exercice 5 : Les tours de Hanoï
print("\n\nExercice 5 : Les tours de Hanoï")
# 5.a.
[docs]def hanoi(n, i, j):
""" Affiche les instructions nécessaires pour effectuer le déplacement d’une pyramide de n rondelles du piquet i vers le piquet j.
- Affiche :math:`O(2^n)` lignes (c'est grand !)
- On peut en fait montrer que ce nombre d'opérations est minimal.
Une légende...
.. note::
Dans le grand temple de Bénarès, au-dessous du dôme qui marque le centre du monde, trois aiguilles de diamant, plantées dans une dalle d'airain, hautes d'une coudée et grosses comme le corps d'une abeille.
Sur une de ces aiguilles, Dieu enfila au commencement des siècles, 64 disques d'or pur, le plus large reposant sur l'airain, et les autres, de plus en plus étroits, superposés jusqu'au sommet.
C'est la tour sacrée du Brahmâ.
Nuit et jour, les prêtres se succèdent sur les marches de l'autel, occupés à transporter la tour de la première aiguille sur la troisième, sans s'écarter des règles fixes que nous venons d'indiquer, et qui ont été imposées par Brahma.
Quand tout sera fini, la tour et les brahmes tomberont, et ce sera la fin des mondes !
(`Ce texte vient de cette page Wikipédia <https://fr.wikipedia.org/wiki/Tours_de_Hano%C3%AF#Origine_du_probl.C3.A8me>`_)
.. note:: Une borne sur la durée nécessaire pour déplacer tous les disques ?
Un jeu à 64 disques requiert un *minimum* de :math:`2^{64}-1` déplacements.
En admettant qu'il faille 1 seconde pour déplacer un disque, ce qui fait 86 400 déplacements par jour, la fin du jeu aurait lieu au bout d'environ 213 000 milliards de jours, ce qui équivaut à peu près à 584,5 milliards d'années, soit 43 fois l'âge estimé de l'univers (13,7 milliards d'années selon certaines sources).
"""
if n <= 0:
return # On s'arrête directement !
# Façon rapide de récupérer le troisième k
k = 6 - i - j
# On peut aussi faire k = {1,2,3} - {i,j} avec des ensembles
# Résolution récursive
hanoi(n - 1, i, k) # Étape 1
print(i, "->", j) # Étape 2
hanoi(n - 1, k, j) # Étape 3
print("Exemple, hanoi(2, 1, 2) :")
hanoi(2, 1, 2)
print("Exemple, hanoi(4, 1, 3) :")
hanoi(4, 1, 3)
# %% Exercice 6 : Programmation récursive de l'algorithme d'Euclide
print("\n\nExercice 6 : Programmation recursive de l'algorithme d'Euclide")
# 6.1.a
# On importe une fonction gcd déjà existante, pour comparer
from fractions import gcd
[docs]def pgcd_rec(a, b):
""" Calcule le PGCD de a et b récursivement.
*Note :* elle a un temps logarithmique mais une pile d'appels logarithmique aussi (la version itérative est meilleure, car sans pile d'appel).
"""
if b == 0:
return a
if b > a:
return pgcd_rec(b, a)
r = a % b
return pgcd_rec(b, r)
print("\nQuelques essais de pgcd_rec(a, b) :")
for a in [10, 20, 30, 200, 400]:
for b in [0, 1, 2, 5, 10]:
d = pgcd_rec(a, b)
print("pgcd_rec de", a, "et", b, "=", d)
if d != gcd(a, b): # Check
print("Oups, il y a une erreur pour ce PGCD (gcd(a,b) donne", gcd(a, b), "mais d =", d)
[docs]def pgcd_rec_verbeux(a, b):
""" Calcule le PGCD de a et b récursivement, et affiche les quotiens et restes successifs."""
if b == 0:
print(" b = 0, on finit avec a =", a)
return a
if b > a:
print(" b > a, on les inverse (a =", a, "et b =", b, ")")
return pgcd_rec_verbeux(b, a)
q = a // b
r = a % b
print(" a =", a, "et b =", b, "donne a = bq + r : q =", q, "et r =", r)
print(" --> pgcd_rec(", b, ",", r, ")")
return pgcd_rec_verbeux(b, r)
print("\nQuelques essais de pgcd_rec_verbeux(a, b) :")
for a in [10, 20]:
for b in [2, 3, 7]:
d = pgcd_rec_verbeux(a, b)
print(" => pgcd_rec_verbeux de", a, "et", b, "=", d)
if d != gcd(a, b): # Check
print("Oups, il y a une erreur pour ce PGCD (gcd(a,b) donne", gcd(a, b), "mais d =", d)
# 6.2.a
print("u = v1, et v = u1 - q1 v1 est la relation demandee.")
# 6.2.b
[docs]def euclideEtendu(a, b):
r""" Utilse l'algorithme d'Euclide étendu pour renvoyer un triplet (d, u, v) solution du problème :
- :math:`d = \mathrm{PGCD}(a, b) \,\text{et}\, d = u a + v b`
- Effectue :math:`O( \ln(\min(a, b)) )` divisions entières, et le même nombre de multiplications et soustractions.
"""
if b == 0:
return (a, 1, 0)
if b > a:
d, v, u = euclideEtendu(b, a)
return (d, u, v)
q1, r1 = a // b, a % b
d, u1, v1 = euclideEtendu(b, r1)
return d, v1, (u1 - q1 * v1)
print("\nQuelques essais de euclideEtendu(a, b) :")
for a in [10, 20]:
for b in [2, 3, 7]:
d, u, v = euclideEtendu(a, b)
print("euclideEtendu de", a, "et", b, "donne d =", d, "u =", u, "et v =", v)
if d != u * a + v * b: # Check
print("Oups, il y a une erreur pour ces u et v : d != u * a + v *b !")
if d != gcd(a, b): # Check
print("Oups, il y a une erreur pour ce PGCD (gcd(a,b) donne", gcd(a, b), "mais d =", d)
# %% Exercice 7 : Énumération récursive de permutations
print("\n\nExercice 7 : Énumeration recursive de permutations")
# from itertools import permutations
from math import factorial
[docs]def int_factorial(n):
""" Factorielle de n, obtenue par :py:func:`math.factorial`."""
return int(factorial(n))
# 7.1.a
[docs]def inserer(li, x):
""" Insère l'élément x dans la liste li à toutes les positions possibles (renvoie une liste de listes)."""
resultats = []
for position in range(len(li) + 1):
t = li[:] # Copie fraiche
t.insert(position, x)
resultats.append(t)
return resultats
print("Exemple d'inserer : inserer([1, 2], 3) =", inserer([1, 2], 3))
# 7.1.b
[docs]def permutations(n):
r""" Renvoie la liste des :math:`n!` permutations de l'ensemble :math:`\{1, 2, \dots, n\}`."""
if n == 0:
return [[]] # Liste contenant une liste vide
else:
sous_perm = permutations(n - 1)
resultats = []
for li in sous_perm:
for p in inserer(li, n):
resultats.append(p)
return resultats
print("Exemple de permutations :")
for n in [0, 1, 2, 3, 5]:
l = permutations(n)
assert len(l) == int_factorial(n)
print("\nPour n =", n, "permutations(n) =\n\t", l)
# 7.2.2
[docs]def estSolution(v):
""" Vérifie si la liste des neuf valeurs v[0], .., v[8] (permutation de [1..9]) est solution du *"problème Vietnamien"*."""
# Première solution, avec une division entière //, non satisfaisant
# return v[0] + (13 * v[1] // v[2]) + v[3] + (12 * v[4]) - v[5] - 11 + (v[6] * v[7] // v[8]) - 10 == 66
# Deuxième solution, avec une division /, meilleure mais on rate les calculs qui ont des erreurs d'arrondis
# return v[0] + (13 * v[1] / v[2]) + v[3] + (12 * v[4]) - v[5] - 11 + (v[6] * v[7] / v[8]) - 10 == 66
# Troisième solution : on oublie les arrondis après le 12ème chiffre
resultat = v[0] + (13 * v[1] / v[2]) + v[3] + (12 * v[4]) - v[5] - 11 + (v[6] * v[7] / v[8]) - 10
reponse = abs(resultat - 66) < 1e-12
if reponse:
print("\t Pour v =", v, "resultat =", resultat, "est bonne.")
return reponse
[docs]def compteSolutions():
""" Compte le nombre de solutions de ce "problème Vietnamien" en utilisant notre fonction :py:func:`permutations`.
J'en ai compté exactement **136** : 128 où le calcul numérique donne 66.0, et 8 où ça donne 65.999999999.
(sur 362880 permutations de {1..9} différentes, soit un taux de 0.035% solution).
"""
nb = 0
print("\nDenombrons les solutions de ce \"probleme Vietnamien\" :")
perms = permutations(9)
total = len(perms)
for v in perms:
if estSolution(v):
nb += 1
print("\tLa", nb, "eme solution est", v)
print("==> Finalement, il y a exactement", nb, "solutions (et on a teste", total, "differentes combinaisons).")
tau = nb / total
print("(Soit un taux de solutions de {:.3%} ...)".format(tau))
return nb
compteSolutions()
# %% Exercice 8 : Le problème des huit dames
print("\n\nExercice 8 : Le probleme des huit dames")
# Comptage du nombre total
[docs]def coefficientBinomial(k, n):
r""" Compute n factorial by a direct multiplicative method: (:math:`C^n_k = {k \choose n}`).
- Complexity: O(1) in memory, O(n) in time.
- From https://en.wikipedia.org/wiki/Binomial_coefficient#Binomial_coefficient_in_programming_languages
- From: http://userpages.umbc.edu/~rcampbel/Computers/Python/probstat.html#ProbStat-Combin-Combinations
"""
if k < 0 or k > n:
return 0
if k == 0 or k == n:
return 1
k = min(k, n - k) # take advantage of symmetry
c = 1
for i in range(k):
c *= (n - i) // (i + 1)
return c
[docs]def parmi(k, n):
r""" Calcule k parmi n avec la factorielle (:math:`C^n_k = {k \choose n} = \frac{n!}{k!(n-k)!})`."""
from math import factorial as f
return int(f(n) // (f(k) * f(n - k)))
n = 8
nbTotal = coefficientBinomial(8, n**2)
nbTotal = parmi(8, n**2)
print("Il y a en tout {} parmi {} facon de poser {} reines sur un echiquier, ce qui fait {} (environ {:.2g}).".format(n, n**2, n, nbTotal, float(nbTotal)))
# 8.1 Verifier une liste de positions
[docs]def verifierDames(positions):
""" Vérifie que les dames présentes aux positions (positions[k] = i,j)_k soient valides sur un équiquier de taille n x n.
- Voir https://fr.wikipedia.org/wiki/Problème_des_huit_dames pour une présentation du problème.
On vérifie qu'aucun conflit n'arrive :
- d'abord sur une ligne (nb d'indices de lignes est exactement n),
- puis aucun conflit sur une colonne (nb indices colonnes = n),
- puis sur les diagonales (montantes et descendantes, nb de i+j différents et de i-j différents est exemple n).
.. image:: solution_1_pb_8_reines.png
"""
n = len(positions)
# Aucun autre reine n'est sur la même ligne
if len({i for (i, j) in positions}) < n:
return False
# Aucun autre reine n'est sur la même colonne
if len({j for (i, j) in positions}) < n:
return False
# Aucun autre reine n'est sur la même diagonal descendante
if len({i + j for (i, j) in positions}) < n:
return False
# Aucun autre reine n'est sur la même diagonal montante
if len({i - j for (i, j) in positions}) < n:
return False
# TODO without sets, just count myself
return True
# Un exemple qui n'est pas solution (deux dames colonne 1)
positions = [(0, 1), (1, 3), (2, 5), (3, 7), (4, 2), (5, 0), (6, 6), (7, 1)]
print("Pour les dames aux positions suivantes :", positions)
print("Est-ce une solution valide ?", "Oui" if verifierDames(positions) else "Non.")
# Un exemple qui n'est pas solution (deux dames ligne 0)
positions = [(0, 1), (0, 3), (2, 5), (3, 7), (4, 2), (5, 0), (6, 6), (7, 1)]
print("Pour les dames aux positions suivantes :", positions)
print("Est-ce une solution valide ?", "Oui" if verifierDames(positions) else "Non.")
# Un exemple qui n'est pas solution (toutes les dames sur la même diagonale descendante)
positions = [(0, 0), (7, 7)]
print("Pour les dames aux positions suivantes :", positions)
print("Est-ce une solution valide ?", "Oui" if verifierDames(positions) else "Non.")
# Un exemple qui n'est pas solution (toutes les dames sur la même diagonale montante)
positions = [(0, 7), (5, 2)]
print("Pour les dames aux positions suivantes :", positions)
print("Est-ce une solution valide ?", "Oui" if verifierDames(positions) else "Non.")
[docs]def afficheEchiquier(positions):
""" Affiche une version très simple d'un équiquier."""
n = len(positions)
print("\n" + ("=" * 2 * n))
print("\n".join('_ ' * i + 'Q ' + '_ ' * (n - i - 1) for i in positions))
print(("=" * 2 * n) + "\n")
# Un exemple de solution (depuis la première image sur https://fr.wikipedia.org/wiki/Problème_des_huit_dames)
positions = [(0, 1), (1, 3), (2, 5), (3, 7), (4, 2), (5, 0), (6, 6), (7, 4)]
print("Pour les dames aux positions suivantes :", positions)
print("Est-ce une solution valide ?", "Oui" if verifierDames(positions) else "Non.")
afficheEchiquier([j0 for (i0, j0) in positions])
# 8.1 Brute force par trop brute
[docs]def enumererDames(n=8):
r""" Résolution du problème des n=8 dames, par énumération (presque) force brute.
Pour énumerer la liste des positions potentiellement valide, on se restreint à celles qui n'ont aucune reine sur une même ligne (:math:`i \neq i'`) ou colonne (:math:`j \neq j'`).
On utilise donc notre énumération des permutations (exercice 7) sur :math:`\{1,\dots,n\}`, ce qui est grand (pour :math:`n = 8, 8! = 40320`) mais vraiment petit face à :math:`{n^2}^n` (qui vaut :math:`64^8 \simeq 2 \times 10^{14}` pour 8 reines).
Numériquement, ça prend environ 350 ms pour :math:`n=8` et 3.3 s pour :math:`n=9` :
>>> %timeit enumererDames(n=8)
1 loop, best of 3: 350 ms per loop
>>> %timeit enumererDames(n=9)
1 loop, best of 3: 3.32 s per loop
.. warning::
Numériquement, c'est déjà trop lourd pour :math:`n=10` (environ 40 secondes)...
.. image:: solution_8_pb_8_reines.png
Pour n = 8, il y a 92 solutions, mais seulement 12 *uniques* en tenant compte de transformations telles que des rotations ou des réflexions (via le `lemme de Burnside <https://fr.wikipedia.org/wiki/Lemme_de_Burnside>`_, tient c'était un de `mes développements pour les oraux d'agrégation <http://perso.crans.org/besson/agreg-2014/>`_).
Voici une illustration, un peu dense, de ces 12 solutions :
.. image:: 12_solutions_pb_8_reines.png
"""
perms = permutations(n)
solutions = []
for p in perms:
positions = [(i, p[i] - 1) for i in range(n)] # list(enumerate(p))
if verifierDames(positions):
solutions.append(positions)
return solutions
print("\n\n8.1 : Enumeration complete.")
for n in [1, 2, 3, 4, 5, 8]:
solutions = enumererDames(n)
print("Pour n = {}, il y a {} solution(s) : (par énumération)".format(n, len(solutions)))
for sol in solutions:
print("\t", sol)
# Affiche la dernière :
afficheEchiquier([j0 for (i0, j0) in sol])
# 8.2 Par backtracking
print("\n\n8.2 : par backtracking !")
[docs]def backtrackingDames(n, afficherSolutions=True):
""" Enrobage pour l'initialisation, partant d'une liste ``p`` vide.
- Utilise la function :py:func:`backtrackingDames_rec`.
"""
return backtrackingDames_rec(n, [], afficherSolutions=afficherSolutions)
[docs]def backtrackingDames_rec(n, p, afficherSolutions=True):
r""" Résolution du problème des n dames, par backtracking.
Idée : on construit la permutation p peu à peu.
Explications de l'algorithme :
- Notation : ``p = [a1, ..., a(i)]``,
- on calcule le nb de sol possibles obtenues en complétant ``p``,
- principe : on cherche à placer une reine sur la ligne ``i+1``,
- on doit éviter un certain nombre d'emplacements,
- pour chaque emplacement licite, on poursuit la complétion récursivement.
.. note:: Ça marche bien jusqu'à :math:`n=12`; pour :math:`n=13` ça galère un peu...::
- n = 1 -> 1
- n = 2 -> 0 (réfléchissez-y)
- n = 3 -> 0 (réfléchissez-y)
- n = 4 -> 2
- n = 5 -> 10
- n = 6 -> 4
- n = 7 -> 40
- n = 8 -> 92 (pb d'origine)
- n = 9 -> 352
- n = 10 -> 724
- n = 11 -> 2680
- n = 12 -> 14200
- n = 13 -> 73712
Numériquement, c'est bien plus rapide que l'énumération exhaustive, ça prend environ 9.5 ms pour :math:`n=8` :
>>> %timeit backtrackingDames(n=8, afficherSolutions=False)
100 loops, best of 3: 9.48 ms per loop
>>> %timeit backtrackingDames(n=9, afficherSolutions=False)
100 loops, best of 3: 41.8 ms per loop
Numériquement, pour :math:`n=10` ou :math:`n=12`, c'est encore pas trop lent, même si ça prend déjà 198 ms :
>>> %timeit backtrackingDames(n=10, afficherSolutions=False)
100 loops, best of 3: 198 ms per loop
>>> %timeit backtrackingDames(n=11, afficherSolutions=False)
100 loops, best of 3: 1 s per loop
>>> %timeit backtrackingDames(n=12, afficherSolutions=False)
100 loops, best of 3: 5.52 s ms per loop
.. warning::
Ça commence à galérer pour :math:`n \geq 13` (40 secondes !).
Soyez rassurés : d'autres techniques *encore plus malines* permettent de résoudre ce problème des n dames pour des valeurs de n plus grandes.
Références :
- Voir `sur Wikipédia <https://fr.wikipedia.org/wiki/Problème_des_huit_dames>`_ pour plus de détails (sur la partie backtracking).
- Voir `sur RosettaCode.org <http://rosettacode.org/wiki/N-queens_problem#Python>`_ pour d'autres résolutions du problème des N-reines en Python.
.. note::
*Message de fin* : bonnes vacances et joyeux Noël !!
Meilleurs voeux de réussite et de bonheur pour l'année 2016 !
"""
i = len(p)
# Cas de base
if i == n:
if afficherSolutions:
# FIXME remove afficheEchiquier and put print(p) back!
# print(p)
afficheEchiquier(p)
return 1 # 1 solution trouvée
# Cas récursif : i < n
# emplacements possibles sur la ligne i
emplacements = [True] * n
# il y a une reine ligne k, colonne x
for k, x in enumerate(p):
emplacements[x] = False # pas dans la même colonne
if x >= i - k:
emplacements[x - i + k] = False # pas en diagonale
if x + i - k < n:
emplacements[x + i - k] = False # idem
nbSol = 0
for x in range(n):
if emplacements[x]:
# On ajoute les solutions avec cet emplacement x
nbSol += backtrackingDames_rec(n, p + [x], afficherSolutions)
return nbSol
for n in [1, 2, 3, 4, 5]:
print("Pour n = {}, affichons toutes les solutions : (par backtracking)".format(n))
nbSolutions = backtrackingDames(n, afficherSolutions=True)
print("\n ==> Pour n = {}, il y a {} solution(s) : (par backtracking)".format(n, nbSolutions))
for n in [6, 7, 8, 9, 10, 11, 12]:
nbSolutions = backtrackingDames(n, afficherSolutions=False)
print("\n ==> Pour n = {}, il y a {} solution(s) : (par backtracking)".format(n, nbSolutions))