Le problème
Le but de cette page est d'arriver à exprimer par une formule simple, i.e. sans sommation, l'expression suivante pour k et n donnés :
On connait tous la réponse pour les cas simples k=1 et k=2 :
On peut se demander "Pourquoi ça t'intéresse ?". Ma réponse est très simple : "parce que". C'est quelque chose qui m'attire personnellement et que je n'explique pas. J'ai appris, bien après avoir trouvé comment obtenir la formule pour n'importe quel k, que le travail avait déjà été fait de la même façon par A.W.F. Edwards au vingtième siècle et d'une autre complètement différente par Jacobi en 1631.
�a peut sembler très naïf mais je ne savais pas que cette formule existait depuis un bon moment et s'appelait la formule de Faulhaber, donc il faut considérer le texte ci-dessous comme un voyage d'un aveugle dans un labyrinthe imaginaire...
Pour citer mes sources, c'est Knuth qui m'a appris l'existence de Faulhaber.
La méthode
On veut donc calculer la somme S des N premiers entiers chacun élevé à la puissance k mais ça n'est pas évident à priori. L'idée est de passer par une intégrale, c'est toujours plus simple à calculer qu'une somme discrète, mais ça nous oblige à trouver une mystérieuse fonction f telle que :
Pour trouver f(x) et pour expliquer pourquoi les bornes de l'intégration ne sont pas celles de la somme, on établit l'équation suivante :
On aboutit, en posant =P(k,x))
, un polynôme en x de degré k, à une résolution qui fait intervenir le développement du binôme, par celui de ^j)
, et des inversions de matrices triangulaires LU.
Résolution
L'équation {2} fait intervenir les matrices naturellement : quand on écrit f(x) sous la forme d'un polynôme, l'intégration, longue pour des valeurs de k élevée, mais très facile, donne Q(i+1)-Q(i), où Q(x) est la primitive de f(x), et cette différence contient des termes en (i+1)^j pour j <= k+1.
Cette différence, et donc l'intégrale, est une expression contenant les coefficients du polynôme f recherché. Cette expression devant être égale à 
, on sait tout de suite que le terme de plus haut degré de f est égal à 1. On a donc :
avec w un polynôme de degré k-1
Pour poursuivre la résolution on peut écrire l'équation {2} sous forme matricielle : A c = (1 0 0 ... 0), avec A la matrice contenant les coefficients du triangle de Pascal issus des développements binomiaux précédent, a le vecteur contenant les coefficients du polynôme f recherché et le membre de droite figurant l'égalité à .
Pour résoudre cette équation on va écrire la matrice A sous la forme A = LU, de manière à pouvoir obtenir :
où U et L sont aisément inversables puisques triangulaires. Je ne m'étends d'ailleurs pas ici sur cette inversion dont on peut trouver la description facilement.
Implémentation
On veut donc convertir la méthode ci-dessus en C++ mais il faut bien faire attention aux problèmes de précisions sur les nombres flottants qui peuvent surgir quand on fait des inversions de matrices ou des multiplications de polynômes successives. Ma méthode pour éviter ça ? Utiliser une classe représentant des fractions et manipuler des entiers ! Les calculs deviennent exacts (ce qui est le but) mais effectivement un peu plus lent (ce qui pour une fois ne m'embète pas du tout !).
Pour le code ci-dessous, les classes utilisées (fraction, polynôme, matrice et vecteur) appartiennent à ma librairie mathématiques en cours de construction.
void Faulhaber(int k)
{
//On remplit une matrice avec les coefficients du triangle de pascal
QiMatrix<QiFraction<long long int>> Pascal(k+1);
Pascal.null();
QiPolynomial< QiFraction<long long int> > Basic(k+1);
Basic.setRaw(1, 1.0f, 1.0f); //x+1
QiPolynomial< QiFraction<long long int> > Current(Basic);
int i, j;
for(i=0; i<=k; i++) {
for(j=0; j<=i; j++) Pascal.set(k-j, k-i, Current.get(j));
Current *= Basic;
}
Pascal.inverse();
//On construit la matrice diagonale T avec k, k-1, ..., 2, 1 (qui est l'inverse de la matrice 1/k, 1/(k-1), ..., 1)
QiMatrix< QiFraction<long long int> > Diag(k+1);
Diag.null();
for(i=0; i<=k; i++) Diag.set(i, i, (long long int)k+1-i);
Diag *= Pascal; //On fait L^-1 * T^-1
//On multiplie par (1, 0, ..., 0) pour récupérer la première colonne
QiVector< QiFraction<long long int> > BaseVector(k+1), Result(k+1);
BaseVector.base(0);
Result = Diag * BaseVector;
//On transforme ce vecteur en polynôme
QiPolynomial< QiFraction<long long int> > Formula(k+1);
Formula.null();
for(i=0; i<=k; i++) Formula.set(k-i, Result.get(i));
//On intègre ce polynôme entre 0 et N+1 en calculant sa primitive et en générant un deuxième polynôme dans lequel on remplace x par n+1
Formula.primitive();
QiPolynomial< QiFraction<long long int> > NPlus1(1);
NPlus1.setRaw(1, 1.0f, 1.0f);
Formula.substitute(NPlus1);
//Présentation des résultats
printf("Formule finale :\n"); Formula.dumpObjects();
printf("On divise par 1/%d :\n", k+1);
Formula /= QiFraction<long long int>(1,k+1);
Formula.dumpObjects();
}
La sortie de ce programme pour quelques valeurs de k ressemble à ce qui suit :
Formule finale : 1X+1
On divise par 1/1 : 1X+1
Résultats pour k=1 :
Formule finale : 1/2X^2+1/2X
On divise par 1/2 : 1X^2+1X
Résultats pour k=2 :
Formule finale : 1/3X^3+1/2X^2+1/6X
On divise par 1/3 : 1X^3+3/2X^2+1/2X
Résultats pour k=3 :
Formule finale : 1/4X^4+1/2X^3+1/4X^2
On divise par 1/4 : 1X^4+2X^3+1X^2
Résultats pour k=4 :
Formule finale : 1/5X^5+1/2X^4+1/3X^3-1/30X
On divise par 1/5 : 1X^5+5/2X^4+5/3X^3-1/6X
Résultats pour k=5 :
Formule finale : 1/6X^6+1/2X^5+5/12X^4-1/12X^2
On divise par 1/6 : 1X^6+3X^5+5/2X^4-1/2X^2
Résultats pour k=6 :
Formule finale : 1/7X^7+1/2X^6+1/2X^5-1/6X^3+1/42X
On divise par 1/7 : 1X^7+7/2X^6+7/2X^5-7/6X^3+1/6X
Résultats pour k=7 :
Formule finale : 1/8X^8+1/2X^7+7/12X^6-7/24X^4+1/12X^2
On divise par 1/8 : 1X^8+4X^7+14/3X^6-7/3X^4+2/3X^2
Résultats pour k=8 :
Formule finale : 1/9X^9+1/2X^8+2/3X^7-7/15X^5+2/9X^3-1/30X
On divise par 1/9 : 1X^9+9/2X^8+6X^7-21/5X^5+2X^3-3/10X
Résultats
Voici les formules des sommes calculées et factorisées à la fois à la main et par ordinateur pour les 8 premières valeurs de k :
=n)
 = x - \frac{1}{2})
=\frac{1}{2}\left(n^2+n\right)=\frac{1}{2}n(n+1))
 = x^2 - x + \frac{1}{6})
=\frac{1}{3}\left(n^3+\frac{3}{2}n^2+\frac{1}{2}n\right)=\frac{1}{3}n(n+1)\left(n+\frac{1}{2}\right))
 = x^3 - \frac{3}{2}x^2 + \frac{1}{2}x)
=\frac{1}{4}\left(n^4+2n^3+n^2\right)=\frac{1}{4}n^2(n+1)^2)
=\frac{n^5}{5}+\frac{n^4}{2}+\frac{n^3}{3}-\frac{n}{30})
=\frac{1}{15}n(n+1)\left(n+\frac{1}{2}\right)\left(3n^2+3n-1\right))
=\frac{1}{15}n(n+1)\left(n+\frac{1}{2}\right)\left(n+\frac{1}{2}-\frac{\sqrt{21}}{6}\right)\left(n+\frac{1}{2}+\frac{\sqrt{21}}{6}\right))
=\frac{1}{5}\left(n^5+\frac{5}{2}n^4+\frac{5}{3}n^3-\frac{1}{6}n\right))
=\frac{n^6}{6}+\frac{n^5}{2}+\frac{5n^4}{12}-\frac{n^2}{12})
=\frac{1}{12}n^2(n+1)^2\left(2n^2+2n-1\right))
=\frac{1}{6}\left(n^6+3n^5+\frac{5}{2}n^4-\frac{1}{2}n^2\right))
=\frac{n^7}{7}+\frac{n^6}{2}+\frac{n^5}{2}-\frac{n^3}{6}+\frac{n}{42})
=\frac{1}{21}n(n+1)\left(n+\frac{1}{2}\right)\left(3n^4+6n^3-3n+1\right))
=\frac{1}{7}\left(n^7+\frac{7}{2}n^6+\frac{7}{2}n^5-\frac{7}{6}n^3+\frac{1}{6}n\right))
=\frac{n^8}{8}+\frac{n^7}{2}+\frac{7n^6}{12}-\frac{7n^4}{24}+\frac{1n^2}{12})
=\frac{1}{24}n^2(n+1)^2\left(3n^4+6n^3-n^2-4n+2\right))
=\frac{1}{24}n^2(n+1)^2\left(3n^2-24n+2\right)\left(n^2+10n+1\right))
=\frac{1}{24}n^2(n+1)^2\left(x-4-\sqrt{\frac{46}{3}}\right)\left(x-4+\sqrt{\frac{46}{3}}\right)\left(x+5-2\sqrt{6}\right)\left(x+5+2\sqrt{6}\right))
=\frac{1}{8}\left(n^8+4n^7+\frac{14}{3}n^6-\frac{7}{3}n^4+\frac{2}{3}n^2\right))
=\frac{n^9}{9}+\frac{n^8}{2}+\frac{2n^7}{3}-\frac{14n^5}{30}+\frac{2n^3}{9}-\frac{n}{30})
=\frac{1}{45}n(n+1)\left(n+\frac{1}{2}\right)\left(5n^6+15n^5+5n^4-15n^3-n^2+9n-3\right))
=\frac{1}{9}\left(n^9+\frac{9}{2}n^8+6n^7-\frac{21}{5}n^5+2n^3-\frac{3}{10}n\right))
On voit des relations entre toutes ces formules à deux indices d'écart (P2 est contenue dans P4, P6 et P8, P3 est contenue dans P5 et P7) mais je ne sais rien en faire pour l'instant. A poursuivre donc. on peut disposer les coefficients des polynômes successifs dans un tableau pour y voir plus clair, en y ajoutant quelques lignes supplémentaires pour donner plus de poids à nos conjectures :
| k | Coeff. | 0 | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | 10 | 11 | 12 | 13 |
 |  |
|  | | |
|  | | | |
 |  | | |  | |
 |  | | | | | |
 |  | |  | |  |  | |
 |  | | |  | | | | |
 |  | | | |  | |  | | |
 |  | | |  | |  | |  | | |
 |  | |  | | | |  | | |  | |
 |  | | |  | | | |  | |  | | |
 |  | |  | |  | |  | |  | |  |  | |
|  | | | | |  | |  | | | | | | |
 |  | |  | |  | |  | |  | |  | |  |  | |
Les lignes 9 à 12 contiennent des coefficients calculés par le programme de cette page et n'ont pas été obtenus à la main...
On peut faire les constats suivants :
- le coefficient de degré le plus élevé est toujours égal à 1 :
- celui de plus bas degré est toujours nul, i.e. pas de terme constant ce qui est assez logique car
=0)
:
- la diagonale reliant les coefficients de degré k est une série arithmétique de raison 1/2 :
- la diagonale située en-dessous est également une série avec pour formule :
- quand la colonne 1 contient 0 (ce qui semble être le cas une fois sur deux) la diagonale issue de ce 0 ne contient que des 0 :
- les autres diagonales sont également des séries avec pour formule :
- La somme des coefficients de la ligne k est égale à k+1, ce qui permet de trouver les
par simple soustraction.
Tout ceci semble indiquer qu'il y a une sorte de structure d'une ligne à l'autre. Je dis bien semble car je peux montrer comment calculer une ligne à partir de celle qui la précède mais je ne l'ai pas encore démontré rigoureusement, seulement vérifié jusqu'à un certain indice, c'est loin d'être la même chose.
demonstration needed
� propos de la dernière remarque concernant les coefficients de la première colonne non-nulle je peux effectivement me débrouiller pour les calculer car la somme des autres coefficients est égale à k+1, mais il est important de remarquer que multipliés par le coefficient 1/k+1 ils sont égaux aux nombres de Bernouilli ! On vient donc de retrouver une méthode pour calculer ces nombres analogue à celle du triangle de Pascal pour les coefficients du binôme !formula needed.
Ce qui me fait croire qu'une partie de la démonstration des formules précédentes et de la procédure utilisée pour les obtenir pourrait résider dans ces fameux nombres de Bernouilli mais ça reste un terrain à investiguer.
Usages
Ces formules sont des ingrédients que l'on peut utiliser dans une décomposition polynomiale des termes d'une série (citation neeeded :)).
Questions ouvertes
Qu'en est-il des cas avec k non entier ? A faire ! �a peut sembler intéressant de formuler la sommes des racines carrées des N premiers entiers non ?
Références
Un article complet sur les sommes
La fameuse formule de Faulhaber
L'article de Knuth sur Faulhaber
Un mathématicien/logisticien/généticien qui s'est intéressé à la formule de Faulhaber
Les nombres de Bernouilli
