Les nombres de Stirling de première espèce non signés $\left[n\atop k\right]$ comptent les permutations de $\mathfrak{S}_n$ ayant exactement $k$ cycles.
Ce sont aussi les coefficients des polynômes $$(t)^n = t^{\overline{n}} = t(t+1)(t+2)\cdots (t+n-1) = \sum_{k=1}^n \left[n\atop k\right]t^k.$$
Les nombres de Stirling de première espèce signés $s(n,k)$ sont les coefficients de
$$(t)_n = t^{\underline{n}} = t(t-1)(t-2)\cdots (t-n+1) = \sum_{k=1}^n s(n,k) t^k.$$Coder une fonction stirling1(n,k,signed=True) calculant les deux versions à l'aide de la récurrence du triangle vue en cours.
def stirling1(n,k,signed=True):
if k<=0 or k>n: return 0
elif n==1 and k==1: return 1
else:
e = -1 if signed else 1
return stirling1(n-1,k-1,signed)+e*(n-1)*stirling1(n-1,k,signed)
for n in range(10):
for k in range(n+1): print(stirling1(n,k,signed=False), end=' ')
print()
for k in range(1,11): print (stirling1(10,k), stirling1(10,k,signed=False))
Coder une fonction fallpow(n,variable=t) calculant $(t)_n$.
from sympy import *
var('x t')
from functools import reduce
def fallpow(n,variable=t):
return reduce(lambda a,b:a*b, [variable-i for i in range(n)])
fallpow(5)
Les nombres de Stirling de deuxième espèce $S(n,k)$ comptent les partitions en $k$ blocs d'un ensemble de $n$ éléments.
Ce sont les coefficients des polynômes (dits de Bell ou de Touchard) $B_n(t)$ dont la série génératrice exponentielle est $$\sum_{n\ge 0}B_n(t)\frac{x^n}{n!}=e^{t(e^x-1)}.$$
Calculer les premiers $B_n(t)$ au moyen de cette série, puis coder une fonction stirling2(n,k) calculant directement les coefficients au moyen de la récurrence du triangle vue en cours.
def stirling2(n,k):
if k<=0 or k>n: return 0
elif n==1 and k==1: return 1
else: return stirling2(n-1,k-1)+k*stirling2(n-1,k)
for n in range(10):
for k in range(n+1): print (stirling2(n,k), end=' ')
print()
En utilisant le développement $$t^n = \sum_{k=1}^n S(n,k)(t)_k$$ et la propriété $$\sum_{s=0}^{t-1}(s)_k = \frac{(t)_{k+1}}{k+1},$$ calculer les premiers polynômes $S_p(t)$ donnant les sommes $$0^p+1^p+\cdots+(n-1)^p = S_p(n).$$
def S(p):
return sum([stirling2(p,j)*fallpow(j+1)/(j+1) for j in range(1,p+1)])
for i in range(1,9): print (S(i).factor())
On va maintenant vérifier que le nombre de partitions d'un ensemble à n éléments en k blocs est bien $S(n,k)$. Pour engendrer récursivement les partitions d'une liste ll, on pourra utiliser deux fonctions auxiliaires qui insérent son dernier élément dans les partitions du reste :
p =((1,3,6),(2,),(4,5)) # une partition de [6], représentée par un tuple de tuples
def insert_all(x,p): # p partition ensembliste, on insère x dans un bloc existant de toutes les manières possibles
return [p[:i] + ((x,)+p[i],) + p[i+1:] for i in range(len(p))]
insert_all(7,p)
def add_block(x,p): # ajout d'un nouveau bloc {x}
return ((x,),)+p
add_block(7,p)
from functools import reduce
def set_partitions(ll):
if not ll: return [tuple()]
else:
mm = set_partitions(ll[1:]); x=ll[0]
return [add_block(x,p) for p in mm] + reduce(lambda a,b:a+b, [insert_all(x,p) for p in mm])
set_partitions(range(1,5))# Et voilà le travail ...
On peut maintenant les dénombrer par nombre de blocs :
def block_lengths(n):
pp = set_partitions(range(1,n+1))
ll = list(map(len,pp))
return [ll.count(i) for i in range(1,n+1)]
On retrouve bien les S(n,k) ...
for i in range(1,10): print (block_lengths(i))
Evidemment, notre algorithme n'est pas très efficace. On trouvera ici le meilleur algorithme connu pour engendrer les partitions ensemblistes. Le programmer en Cython, et comparer ses performances avec l'algorithme naif codé précédemment.
# Version pur Python
def d2w(d):
return ''.join([str(d[i+1]) for i in range(len(d))])
def SetPartitions(n):
ll = []
c = {}
def SP(m,p):
if p<n:
for i in range(1,m+1):
c[p] = i
SP(m,p+1)
c[p] = m+1
SP(m+1,p+1)
elif p==n:
for i in range(1,m+2):
c[p] = i
ll.append(d2w(c))
SP(0,1)
return ll
# d'abord, la méthode bête et brutale
from time import time
t = time(); print (len(set_partitions(range(1,11))));print (time()-t)
# ensuite, l'algorithme subtil
t = time(); print (len(SetPartitions(10)));print (time()-t)
# Et maintenant, la version Cython
from sp import SetPartitions as setp
t = time(); print (len(setp(10)));print (time()-t)
Voilà la marche à suivre ...
$ cat sp.pyx
#def d2w(d):
# return ''.join([str(d[i+1]) for i in range(len(d))])
def SetPartitions(int n):
ll = []
c = [0]*n
def SP(int m,int p):
if p<n:
for i in range(1,m+1):
c[p-1] = i
SP(m,p+1)
c[p-1] = m+1
SP(m+1,p+1)
elif p==n:
for i in range(1,m+2):
c[p-1] = i
ll.append(tuple(c))
SP(0,1)
return ll
$ cat setup.py
from distutils.core import setup
from Cython.Build import cythonize
setup( ext_modules = cythonize("sp.pyx"))
$ python setup.py build_ext --inplace
$ ls sp.*
sp.c sp.py sp.pyx sp.pyx~ sp.so*