Remarque.

Les $C_{k+1}^t$ représentent les coefficients du binôme de Newton, c'est-à-dire le nombre de combinaisons de $t$ éléments d'un ensemble à $k+1$ éléments.

La définition des nombres de Bernoulli n'est pas tout à fait standardisée: il y traîne chez certains auteurs des facteurs $(-1)^t$; chez d'autres les indices $t$ deviennent $t/2$, et j'en passe. Il convient donc de toujours bien regarder la définition adoptée.

Avec celle-ci, on a:

$$\begin{array}{c@=l} B_0 & 1\ B_1 & -\frac{1}{2} \ B_2 & ... ... & \frac{7}{6} \ B_3 & B_5 = B_7 = \cdots = 0 \ \end{array}$$

On connaît bien sûr des techniques de calcul rapide des nombres de Bernoulli, la plupart récurrentes.

À propos de formules explicites pour calculer rapidement ces nombres, on dispose tout de même du théorème de von Staudt - Clausen qui dit que $B_{2k} + \sum 1/p$ étendue aux nombre premiers $p$ tels que $(p-1)$ divise $(2k)$ est un entier.

Sachant par ailleurs que, pour $k > 0$, on a

$$B_{2k} = (-1)^{k-1} 2 (2k)! \frac{\zeta(2k)}{(2\pi)^{2k}}$$

$\zeta$ étant bien sûr la fonction de Riemann, on peut programmer le problème par un procédé hybride sans utiliser de récurrence (d'abord approcher, en multi-précision, puis obtenir la valeur rationnelle exacte grâce à von Staudt).

On a alors comme valeurs pour $k$ variant de $2$ à $9$ (inclus):

$$\begin{array}{c@=l} \sum_{i=1}^n i^2 & \frac{1}{6} n (n + 1)... ...- 1) (2 n^4 + 4 n^3 - n^2 - 3 (n+3) (n + 1)^2 \ \end{array}$$

On peut également dévélopper une approche algébrique assez détaillée. On montre dans un premier temps que cette somme $S(r,n)$ peut s'écrire $\sum_{i=0}^r H(r,i)n^{i+1}$.

Le fait qu'un tel polynôme existe découle de l'observation suivante: Dans le sous espace vectoriel des polynômes de degrés inférieurs ou égaux à $r$, les $u_i$ ($i$ variant de $0$ à $r$) définis par: $u_i = (1+X)^{i+1} - X^{i+1}$ constituent une base.

Si l'on nomme $(H(r,i) ; i=0,\ldots,r)$ le système des coordonnées de $(1+X)^r$ on retrouve que $H(r,i) n^{i+1} = \sum_{p=1}^n p^r$

L'utilisation de cette base permet de déduire que les $H(r,i)$ sont solution du système $Mv=w$.

où:

• $C_y^x$ désignant le coefficient binômial bien connu...
• $M$ est la matrice triangulaire supérieure dont le terme $m(l,c)$ vaut $C_{l-1}^c$ pour $c$ entre $1$ et $r+1$ et $l$ entre $1$ et $c$ par exemple dans le cas $r=2$, la matrice est:
  
  $$\left( \begin{array}{ccc} C_0^1 & C_0^2 & C_0^3 \ 0 & C_1^2 & C_1^3 \ 0 & 0 & C_2^3 \ \end{array} \right)$$
  
  
• $v$ est le vecteur dont les composantes sont $H(r,0)$ à $H(r,r)$
• $w$ est le vecteur dont les composantes sont $C_0^r$ à $C_r^r$.

Le pivot de Gauss donne rapidement les premières valeurs:

$$\begin{array}{l@=cl} H(r,r) & \frac{1}{r+1} \ H(r,r-1) & \... ...H(r,r-6) & \frac{ C_5^3 }{252} \ H(r,r-7) & 0 \ \end{array}$$

Il y a là de quoi traiter rapidement jusqu'au cas $r=7$ mais...En dérivant l'égalité polynômiale somme des $u_i = (1+X)^r$ on trouve $H(r,i) = \frac{r}{i+1} H(r-1,i-1)$

On en déduit $H(r,i) = \frac{ C_i^r }{i+1} H(r-i,0)$.

Cette relation limite les recherches aux $H(k,0)$ mais, dans le système initial, c'est celui que le pivot de Gauss donne en dernier...

La relation donnant $H(r,i)$ en fonction de $H(r,i+1),\ldots,H(r,r)$ devient

$$H(r,0) = 1 - \sum_{j=0}^{r-1} \frac{ C_j^{r+1} }{r+1} H(j,0)$$

et l'on a $H(0,0) =1$.

Cette dernière relation permet de programmer rapidement le calcul des coefficients dans les cas où $r$ est assez grand...

Par ailleurs, on a remarqué qu'il y avait pas mal de $0$...En fait, on a: pour tout $k$: $H(2k+3,0)=0$. Cela peut se montrer de plusieurs façons...'égalité polynômiale avec $X=1$ puis $X=-1$ donne:

$$\sum_{\stackrel{j=0}{j\text{ pair}}}^k H(k,k-j) = \sum_{\stackrel{j=0}{j\text{ impair}}}^k H(k,k-j) = 1/2$$

en bricolant un peu avec les termes connus on trouve que

$$\sum_{\stackrel{j=0}{j\text{ impair}}} C_j^{k+1} H(j,0)$$

est nulle. cette égalité est valable pour tout $k$ (au moins $3$...) et l'on sait que $H(3,0)$ est nul...donc, de proche en proche...