Le problème du tas de boulets

On dit que pendant la guerre de 1870, Édouard Lucas a fait le constat qu'il est difficile de mettre en tas pyramidal, un nombre carré de boulets. Ce dessin de David Eppstein illustre le problème:

Ce tas de boulets comprend 16+9+4+1=30 boulets et effectivement 30 n'est pas un carré.

Ce problème est intéressant du point de vue algorithmique: Comment

• Calculer la somme des carrés de entiers consécutifs ?
• Boucler jusqu'à ce que cette somme soit un carré ?

Comment tester si un nombre est un carré ?

L'algorithme utilisé consiste à regarder s'il y a des chiffres après la virgule dans le développement décimal de √(n) :

1. √(1) = 1 : pas de chiffre après la virgule, 1 est un carré ;
2. √(2) ∼ 1,414: après la virgule il y a les chiffres 414 etc, donc 2 n'est pas un carré ;
3. √(3) ∼ 1,732: Là encore il y a des chiffres après la virgule et 3 n'est pas un carré ;
4. √(4) = 2 : Pas de chiffres après la virgule, et 4 est un carré.

Comment tester si un nombre est entier ?

CoffeeScript étant basé sur JavaScript, n'a pas de fonction donnant les éventuels chiffres après la virgule. Donc pour les connaître, on soustrait à x sa partie entière avec x - Math.floor x. Oui mais... L'image de 2,999999999999 par cette fonction est 0,9999999999989999 qui a bien des chiffres après la virgule, alors que 2,999999999999 est entier "aux erreurs de calcul près". Donc,

• Au lieu de tronquer avec Math.floor, il vaut mieux arrondir avec Math.round au cas où x est, de peu, inférieur à un entier;
• Au lieu de dire que x est entier s'il est égal à son arrondi (différence nulle), on dira que x est entier si la différence entre x et l'entier le plus proche est inférieure à 10^-15 :

isInteger = (x) -> (Math.abs(x-Math.round x) < 1e-15) alert (isInteger x for x in [2..3.2] by 0.2)

Lancer cet algorithme

Maintenant que CoffeeScript est enrichi d'une nouvelle fonction isInteger, on va s'en servir pour définir une nouvelle fonction isSquare :

Un entier est un carré si et seulement si sa racine est entière

isInteger = (x) -> (Math.abs(x-Math.round x) < 1e-15) isSquare = (n) -> isInteger Math.sqrt n alert (isSquare n for n in [1..10])

Lancer cet algorithme

Comment additionner des carrés ?

Ce qui précède (création d'algorithmes pour faciliter la suite) ne sert pas à calculer la somme des carrés, mais on s'en servira après. Par contre, l'algorithme pour la somme des carrés est classique:

• On a besoin de deux variables, sum pour la somme et index pour l'indice dans la boucle.
• Ces deux variables sont initialisées à 0.
• Lors de chaque passage dans la boucle,
  • index est incrémenté avec index++
  • Son carré est additionné à la somme courante sum
• Si on boucle 4 fois, on a donc la somme des 4 premiers carrés

sumOfSquares = (n) -> sum = 0 for index in [0..n] sum += index*index sum alert sumOfSquares 4

Lancer cet algorithme

Le fait de boucler avec for plutôt que while permet de raccourcir le code; on peut raccourcir encore plus en inversant les proposition principale et subordonnée :

sumOfSquares = (n) -> sum = 0 sum += index*index for index in [0..n] sum alert (sumOfSquares n for n in [1..10])

Lancer cet algorithme

Il ne reste maintenant plus qu'à assembler les résultats des deux premières parties pour pouvoir chercher l'éventuelle possibilité que la somme des carrés des entiers successifs finisse par être elle-même un carré.

tu ne t'arrêteras que lorsque tu auras fini !

L'algorithme est simple :

On additionne les carrés jusqu'à ce que la somme soit carrée

isInteger = (x) -> (Math.abs(x-Math.round x) < 1e-15) isSquare = (n) -> isInteger Math.sqrt n sumOfSquares = (n) -> sum = 0 sum += index*index for index in [0..n] sum n = 2 n++ until isSquare sumOfSquares n alert n

Lancer cet algorithme

On boucle à partir de n = 2 parce que des tas de 0 ou 1 boulet ne sont pas considérés comme des tas...

Pour aller plus vite

Si on détaille les boucles, on constate qu'on calcule plusieurs fois le carré de 2 :

isInteger = (x) -> (Math.abs(x-Math.round x) < 1e-15) isSquare = (n) -> isInteger Math.sqrt n sumOfSquares = (n) -> sum = 0 index = 0 until index is n index++ sum += index*index sum n = 2 until isSquare sumOfSquares n n++ alert n

Lancer cet algorithme

Ceci parce qu'il y a des boucles imbriquées et que la somme des carrés jusqu'à n+1 est juste (n+1)² de plus que la somme des carrés allant jusqu'à n: Autrement dit, la suite u_n des carrés des nombres inférieurs ou égaux à n, vérifie la relation de récurrence suivante :

u_n+1=u_n+(n+1)²

Ceci permet de simplifier l'algorithme, en calculant moins souvent les carrés, et en cherchant plus vite :

isInteger = (x) -> (Math.abs(x-Math.round x) < 1e-15) isSquare = (n) -> isInteger Math.sqrt n n = 2 sum = 5 until isSquare sum n++ sum += n*n alert n

Lancer cet algorithme

On peut améliorer l'affichage final :

isInteger = (x) -> (Math.abs(x-Math.round x) < 1e-15) isSquare = (n) -> isInteger Math.sqrt n n = 2 sum = 5 until isSquare sum n++ sum += n*n alert "La somme des carrés des #{n} premiers entiers est #{sum} qui est le carré de #{Math.sqrt sum}"

Lancer cet algorithme

Ceci répond à la question de l'existence d'un nombre n (autre que 0 et 1) tel que la somme des carrés des n premiers entiers est elle-même un carré. Mais y en a-t-il d'autres ?

Calcul formel

Pour retrouver algébriquement ces résultats, on peut faire appel à du calcul formel, en effet il existe une forme explicite pour la somme des carrés. Ici on va utiliser Xcas, mais d'autres logiciels de calcul formel font tout autant l'affaire. Tout d'abord, on peut retrouver la somme des 4 premiers carrés avec sum(k^2,k,0,4) qui répond bien 30. Ensuite il suffit de remplacer 4 par n en faisant sum(k^2,k,0,n) pour avoir (2*(n+1)³-3*(n+1)²+n+1)/6. Mieux, avec simplify(sum(k^2,k,0,n)), on apprend que la somme des n premiers carrés est (2*n³+3*n²+n)/6 (formule à montrer par récurrence évidemment). On a donc déjà un moyen de calculer la somme des carrés sans boucler. Par ailleurs, en demandant à Xcas factor((2*n^3+3*n^2+n)/6) on a une forme encore plus simple pour la somme des carrés: n(n+1)(2n+1)/6. Il s'agit donc de résoudre l'équation diophantienne 2x³+3x²+x=6y². Xcas n'y arrive pas, mais Wolfram Alpha y parvient, en lui soumettant le calcul 2*x^3+3*x^2+x=6*y^2.

Wolfram Alpha dessine la courbe elliptique (puisque l'équation qu'on lui a fournie est celle d'une courbe elliptique):

Parmi les résultats affichés, on voit

Integer solutions:

• x=1, y=±1
• x=24, y=±70