Philippe.Dumas@inria.fr


http://algo.inria.fr/dumas/Maple/

Retour en page principale Table des matières Index

Maple V.4 worksheet, bis Maple V.5 worksheet, bis

Réduction à une équation de Fermat | Réduction à une équation auxilaire
Première tentative de résolution | Fractions continuées
Exponentiation binaire

Le dieu du soleil avait un troupeau de bétail composé de taureaux et de vaches ; une partie des bêtes était blanche, une autre était noire, une troisième était tachetée et une quatrième était brune.

Pour ce qui est des taureaux, le nombre de blancs surpassait le nombre de noirs d'un demi plus un tiers du nombre de bruns ; le nombre de noirs était plus grand que le nombre de bruns d'un quart plus un cinquième du nombre de tachetés ; le nombre de tachetés dépassait le nombre de bruns d'un sixième plus un septième du nombre de blancs.

En ce qui concerne les vaches, le nombre de blanches était égal à un tiers plus un quart du nombre total d'animaux noirs ; le nombre de noires égalait un quart plus un cinquième du nombre de bêtes tachetées ; le nombre de tachetées valait un sixième plus un septième du nombre d'animaux blancs ; le nombre vaches brunes se montait à un sixième plus un septième du nombre total d'animaux blancs.

Quelle était la composition du troupeau ?

C'est ainsi que le problème des boeufs (ou taureaux) d'Archimède est posé dans un manuscrit en grec trouvé en 1773 dans la bibliothèque de Wolfenbüttel par G. E. Lessing (cf. [ Doerrie65 ] pour une version anglaise du texte). Ce problème est usuellement attribué à Archimède. On en trouve une discussion assez approfondie dans [ Dickson66, chap. XII, pp. 342-345].

Réduction à une équation de Fermat. Appelons respectivement W, X, Y, Z le nombre de taureaux blancs, noirs, tachetés et bruns, puis w, x, y, z le nombre de vaches des couleurs correspondantes. Le problème précédent est linéaire et se résout sans peine, même si l'on ne cherche que des solutions entières. Cependant le problème comporte aussi deux contraintes non linéaires :

Nous allons voir que la résolution de ce problème amène des calculs assez importants et nous donnons donc quelques temps de calcul à titre indicatif.

> global_start:=time():

Nous définissons les inconnues et les équations.

> inc1:={W,w,X,x,Y,y,Z,z}:

  inc2:={n,m}:

> equW:=W=(1/2+1/3)*X+Z:

  equX:=X=(1/4+1/5)*Y+Z:

  equY:=Y=(1/6+1/7)*W+Z:

  equw:=w=(1/3+1/4)*(X+x):

  equx:=x=(1/4+1/5)*(Y+y):

  equy:=y=(1/5+1/6)*(Z+z):

  equz:=z=(1/6+1/7)*(W+w):

  equm:=m^2=W+X:

  equn:=Y+Z=n*(n+1)/2:

Nous résolvons d'abord la partie linéaire du système.

> S1:=isolve({seq(equ.i,i=inc1 minus {Z})},k);

Nous substituons le résultat obtenu dans l'équation qui exprime que la somme W+X est un carré.

> S2:=[isolve(subs(S1,equm),p)];

Le résultat obtenu n'est pas celui qui est annoncé page 108. Cependant il s'agit d'une illusion.

> for i to nops(S2) do

    kk.i:=factor(subs(S2[i],k));

    mm.i:=factor(subs(S2[i],m));

  od;

> for i to nops(S2) do

    k.i:=factor(subs(p=(i mod 2+pp)/2,kk.i));

    m.i:=factor(subs(p=(i mod 2+pp)/2,mm.i));

  od;

Nous écrivons donc explicitement le résultat sous la forme annoncée.

> S2:={k=4456749*p^2,m=-8913498*p};

Ensuite nous substituons dans l'équation qui exprime que Y+Z est un nombre triangulaire. Il ne reste plus que deux inconnues, p et n.

> subs(S1,S2,equn);

L'équation se ramène à une équation de Fermat (dite aussi de Pell-Fermat, mais cette appelation est non fondée et due à une erreur d'Euler [ Dickson66, chap. XII, p. 341]).

> fermat:=op(1,%)-op(2,%);

> subs(n=N/2-1/2,fermat);

> fermat1:=expand(8*%);

Réduction à une équation auxilaire. Il ne reste plus qu'à résoudre cette équation. On pourrait évidemment tenter de la résoudre par isolve comme à la page 105, mais les nombres en jeu sont trop grands que pour qu'une telle instruction réussisse. Regardons mieux le nombre qui détermine les propriétés des solutions.

> d:=coeff(fermat1,p,2);

> ifactor(d);

Le nombre d n'est pas libre de carré.

> df:=readlib(ifactors)(d)[2];

La taille du nombre d est excessive et nous allons commencer par résoudre l'équation de Fermat associée à la partie libre de carré de d, que nous notons dd.

> sqrtd:=mul(i[1]^iquo(i[2],2),i=df)

                         *sqrt(mul(i[1]^irem(i[2],2),i=df));

> dd:=mul(i[1]^irem(i[2],2),i=df);

> sqrtdd:=sqrt(dd);

Première tentative de résolution. À ce point, Paul Zimmermann me fait remarquer que l'on peut poursuivre la résolution comme suit. On applique d'abord la procédure isolve à l'équation de Fermat disons réduite.

> isolve(N^2-dd*P^2=1);

La réponse est assez illisible. En fait les quatre solutions fournies ne différent que par des changements de signe (ceci est plus clair en V.3, où le prettyprinter écrit le résultat avec les abréviations %1, %2, ...). Celle qui nous importe est celle qui donne des valeurs positives.

> 1/9458988*4729494^(1/2)*

 ((109931986732829734979866232821433543901088049

   +50549485234315033074477819735540408986340*4729494^(1/2))^_N1

 -(109931986732829734979866232821433543901088049

   -50549485234315033074477819735540408986340*4729494^(1/2))^_N1);

Il faut maintenant trouver les valeurs divisibles par 9314. On remarque que l'on a

P(n)=sqrt(dd)/9458988*((a+b*sqrt(dd))^n-(a-b*sqrt(dd))^n)

où

a=109931986732829734979866232821433543901088049,
b=50549485234315033074477819735540408986340

et toujours dd=4729494. Autrement dit, on a

P(n)=sqrt(dd)/9458988*(x1^n-x2^n)

où x1, x2 sont les deux racines de l'équation x^2-2*a*x+1=0, donc les P(n) vérifient la récurrence P(n) = 2*a*P(n-1)-P(n-2), et si on note Pmod(n)=P(n) mod 9314, les Pmod(n) vérifient la même récurrence. On peut donc calculer comme suit le premier exposant pour lequel l'entier P est divisible par 9314.

> a:=109931986732829734979866232821433543901088049:

  b:=50549485234315033074477819735540408986340:

> Pproc := proc(n)

    expand(sqrtdd/9458988*((a+b*sqrtdd)^n-(a-b*sqrtdd)^n))

  end:

> (2*a) mod 9314;

> Pmodproc := proc(n)

     option remember;

    (8812*Pmodproc(n-1)-Pmodproc(n-2)) mod 9314

  end:

> Pmodproc(0):=Pproc(0) mod 9314:

> Pmodproc(1):=Pproc(1) mod 9314:

> start:=time():

  for nu while Pmodproc(nu)<>0 do od:

  time()-start;

  nu-1; 

Il suffit donc d'élever à la puissance 2328 la solution de l'équation de Fermat associée à dd pour obtenir la plus petite solution de l'équation de Fermat associée à d. Cependant il n'est pas clair que le calcul se termine vu la taille des nombres mis en jeu. D'autre part la procédure isolve apparaît jusqu'ici comme une boîte noire, ce qui assez désagréable. Nous allons donc reprendre le calcul par une autre voie. Cela nous permettra de comprendre la théorie sous-jacente et d'obtenir effectivement la plus petite solution du problème.

Fractions continuées. On sait que les solutions de l'équation de Fermat sont liées au développement en fraction continuée de sqrt(d). Nous revenons donc à la méthode classique qui passe par les fractions continuées. On pourra consulter à ce sujet [Perron54a], [RoSz92] ou [Faisant91].

Nous avons besoin du développement en fraction continuée de sqrt(4729494), car 4729494 est congru à 2 modulo 4. Le logiciel propose pour cela la procédure convert/confrac. Elle prend en entrée une approximation décimale du nombre à développer et renvoie une troncature du développement. Cependant son emploi n'est pas satisfaisant car nous sommes dans un cas où un calcul exact peut être mené. De plus nous ne voulons pas seulement déterminer quelques quotients partiels de la fraction. Nous savons que la fraction est périodique, car nous avons affaire à un algébrique de degré 2, et nous voulons déterminer la période, car sa connaissance est fondamentale pour obtenir les solutions de l'équation de Fermat.

Nous écrivons donc une procédure algtoconfrac qui prend en entrée un nombre algébrique exprimé par radicaux et une borne sur le nombre d'itérations de l'algorithme des fractions continuées. Dans le cas d'un algébrique de degré 2 la procédure renvoie une séquence de deux listes ou bien une liste. Dans le premier cas, la première liste, qui peut être vide, fournit le motif initial et la seconde liste fournit le motif périodique. Dans le second cas, la liste fournit une troncature de la suite des quotients partiels (et la période n'a pas été trouvée car la première occurrence du motif périodique se termine après la borne qui a été passée en deuxième argument). On notera comment est utilisé assigned pour reconnaître qu'un quotient complet a déjà été rencontré. Ceci permet d'avoir un coût de calcul linéaire en le nombre de quotients. D'autre part l'algorithme demande de reconnaître zéro et la procédure de normalisation adéquate est radnormal. Enfin l'emploi de assigned n'est correct que si l'on dispose d'une forme canonique, c'est-à-dire d'une écriture unique, pour les quotients complets. Ceci est fournit dans le cas qui nous intéresse, celui d'une extension quadratique Q(sqrt(d)) avec d libre de carré, par radnormal/rationalized.

> algtoconfrac:=proc(x,N::nonnegint)

    local n,z,T,history,i,delta;

    z:=radnormal(x,rationalized);

    for n to N do

      T[n]:=floor(z);

      if assigned(history[z]) then

        i:=history[z];

        RETURN([seq(T[j],j=1..i-1)],[seq(T[j],j=i..n-1)])

      else

        history[z]:=n;

      fi;

      delta:=radnormal(z-T[n]);

      if delta=0 then RETURN([seq(T[i],i=1..n)]) 

      else z:=radnormal(1/delta,rationalized)

      fi;

    od;

    [seq(T[n],n=1..N)]

  end: # algtoconfrac

Nous appliquons la procédure à sqrt(4729494).

> start:=time():

  CFe:=algtoconfrac(sqrtdd,1000):

  time()-start;

> CFe;

La procédure a fourni la période et les motifs. On note la forme très particulière (et attendue) du développement.

> m1:=nops(cFe[1]);

> m2:=nops(CFe[2]);

La période vaut 92 (et non 91 comme indiqué dans [ Dickson66, chap. XII, p. 344]).

Les réduites de la fraction vont fournir les solutions de l'équation de l'équation de Fermat. Plus précisément comme la période m2 est paire les réduites A_nu/B_nu d'indice nu=k*m2-1 fournissent les solutions de x^2-4729494*y^2=1. Si cette période avait été impaire, nous aurions eu au second membre (-1)^k. Nous employons donc le nombre m3 ci-dessous.

> if type(m2,even) then m3:=m2 else m3:=2*m2 fi;

Vérifions que la réduite d'indice m3-1 fournit l'unité fondamentale de l'anneau Z[sqrt(dd)].

> L:=[op(CFe[1]),op(CFe[2])]:

> oldA:=1:

  A:=CFe[1][1]:

  oldB:=0:

  B:=1:

  for nu to m1 do

    newA:=L[nu+1]*A+oldA;

    newB:=L[nu+1]*B+oldB;

    oldA:=A;

    oldB:=B;

    A:=newA;

    B:=newB;

  od:

  for nu from m1+1 to m3-1 do

    nn:=nu-iquo(nu-m1,m2)*m2;

    newA:=L[nn+1]*A+oldA;

    newB:=L[nn+1]*B+oldB;

    oldA:=A;

    oldB:=B;

    A:=newA;

    B:=newB;

  od:

> A,B;

> A^2-dd*B^2;

Nous obtenons bien 1 comme prévu, ce qui prouve que nous avons affaire à une unité. Le fait que l'unité soit fondamentale est connue par théorème. Cela signifie que toutes les unités s'écrivent (A+B*sqrt(dd))^m. Nous conservons ces valeurs car elles vont resservir.

> XA:=oldA;

  YA:=A;

  XB:=oldB;

  YB:=B;

La difficulté est que nous voulons une réduite A/B où B est divisible par 9314. Pour obtenir un indice adéquat, nous calculons donc les numérateurs et dénominateurs des réduites modulo 9314 jusqu'à voir apparaître 0.

> modulus:=mul(i[1]^iquo(i[2],2),i=df);

> oldA:=1:

  A:=CFe[1][1]:

  oldB:=0:

  B:=1:

  for nu to m1 do

    newA:=L[nu+1]*A+oldA mod modulus;

    newB:=L[nu+1]*B+oldB mod modulus;

    oldA:=A;

    oldB:=B;

    A:=newA;

    B:=newB;

  od:

  for nu from m1+1 to m3-1 do

    nn:=nu-iquo(nu-m1,m2)*m2;

    newA:=L[nn+1]*A+oldA mod modulus;

    newB:=L[nn+1]*B+oldB mod modulus;

    oldA:=A;

    oldB:=B;

    A:=newA;

    B:=newB;

  od:

  for kappa while B<>0 do

    for nu from kappa*m3 to (kappa+1)*m3-1 do

      nn:=nu-iquo(nu-m1,m2)*m2;

      newA:=L[nn+1]*A+oldA mod modulus;

      newB:=L[nn+1]*B+oldB mod modulus;

      oldA:=A;

      oldB:=B;

      A:=newA;

      B:=newB;

    od;

  od:

> right_nu:=nu-1;

Nous avons ainsi le bon indice. Vérifions que le B correspondant est bien multiple de 9314.

> A,B;

Exponentiation binaire. La question qui se pose maintenant est le calcul de cette réduite. En effet l'indice est passablement grand et les deux nombres A et B associés doivent être énormes. Il faut donc un moyen de calcul efficace. Les calculs emploient deux récurrences linéaires d'ordre 2 qui ne diffèrent que par les conditions initiales. Par ailleurs, elles ne sont pas à coefficients constants. À chaque pas de la récurrence, on multiplie un vecteur de Z^2 par une matrice de taille 2x2. Notons T_nu la matrice employée au rang nu. La réduite au rang nu est donc donnée par le vecteur T_nu T_(nu-1) T_(nu-2) ... T_0 V_(-1). Cependant la périodicité permet de grouper les matrices par paquets dans ce produit. Ceci fait apparaître une matrice M et du coup il ne s'agit plus essentiellement que de calculer une puissance de cette matrice. La bonne méthode est alors l'exponentiation binaire.

Regardons la puissance utile de la matrice M.

> qq:=iquo(right_nu,92,'rr');

Il faut calculer M^2328. De plus ce produit devra être appliqué au vecteur V_rr=T_rr T_(rr-1) ... T_0 V_(-1) avec rr le reste dans la division euclidienne.

> rr;

Nous retrouvons la valeur 91, ce qui est normal puisque nous savons que les solutions de l'équation de Fermat sont données par les réduites d'indice k*t-1 si t est la période. C'est pour cela que nous avons conservé les valeurs, qui sont rangées dans XA, YA, XB, YB.

Pour déterminer la matrice M, nous appliquons la récurrence à un vecteur générique de façon à effectuer les transformations sur une période, en démarrant au bon indice.

> oldC:=x:

  C:=y:

  for nu from m3 to 2*m3-1 do

    nn:=nu-iquo(nu-m1,m2)*m2;

    newC:=L[nn+1]*C+oldC;

    oldC:=C;

    C:=newC

  od:

Nous récupérons les coefficients et nous les rangeons dans la matrice M.

> M:=array(1..2,1..2,

     [[subs(x=1,y=0,oldC),subs(x=0,y=1,oldC)],

      [subs(x=1,y=0,C),subs(x=0,y=1,C)]]);

Évaluons la taille du résultat avant d'aller plus loin. Nous munissons Z^2 ou R^2 de la norme d'indice infini, la norme du maximum, et les matrices 2x2 de la norme associée. La norme de M a la valeur suivante.

> Mnorm:=max(M[1,1]+M[1,1],M[2,1]+M[2,2]);

Nous en tirons la majoration suivante pour les composantes du résultat.

> evalf(Mnorm^qq*max(XA,XB,YA,YB));

Comme tous les entiers manipulés sont positifs, il y a de forte chance que la borne soit quasiment atteinte.

Passons au calcul. Nous évaluons M^qq par exponentiation binaire. De plus nous faisons apparaître le temps de calcul pour chaque tour de boucle.

> L:=convert(qq,base,2);

> start:=time():

  MM:=M:

  P:=&*():

  for i from 2 to nops(L) do

    MM:=evalm(MM&*MM);

    if L[i]=1 then 

      P:=evalm(P&*MM) 

    fi;

    print(i,time()-start);

  od:

À chaque itération la taille des nombres (comptée par le nombre de chiffres décimaux) double. C'est pourquoi la dernière étape est celle qui domine tout le calcul.

Calculons les deux composantes de la réduite et leur taille.

> A:=evalm(P&*vector([XA,YA]))[2]:

  B:=evalm(P&*vector([XB,YB]))[2]:

> length(A),length(B);

Nous avons bien la taille attendue. Vérifions aussi la cohérence des résultats.

> A mod modulus;

> B mod modulus;

> evalb(A^2-dd*B^2=1);

Revenons maintenant au problème d'Archimède et utilisons le résultat obtenu.

> SOL:=subs(remove(has,S2,m),{N=A,p=iquo(B,modulus)},

                                {n=N/2-1/2} union S1 union S2):

Nous calculons le décuplet solution et nous regardons la taille des nombres en jeu.

> for i in inc1 union inc2 do

    sol.i:=subs(SOL,i);

    print(i,length(sol.i))

  od:

> time()-global_start;

L'information est essentiellement dans le nombre B. Nous regardons sa tête et sa queue.

> evalf(B,50);

> B mod 10^50;

Nous le rangeons dans un fichier archimedes.html ce qui vous permettra de le contempler.

> writeto("archimedes.html"):

  B;

  writeto(terminal):

Par curiosité nous regardons le nombre de bêtes qui constituent le troupeau du dieu du soleil.

> S:=subs({seq(i=sol.i,i=inc1 union inc2)},x+y+z+w+X+Y+Z+W):

> length(S);

> evalf(S);

Le texte qui énonce le problème affirme que ce troupeau paisse dans les prairies de Sicile. L'aire de la Sicile est de 25 608 km^2 !

Pour finir, notons que nous n'avons fourni qu'une solution du problème, en fait la plus petite. Les autres s'obtiennent à l'aide de puissances convenables de la matrice M (les puissances multiples de celle que nous avons employée).

Il apparaît qu' Ilan Vardi a écrit une note [ référence Bibtex, note en postcript] sur ce sujet. Cette note est beaucoup mieux documentée que la présente page. Il est plaisant de constater que les résultats sont similaires.

On peut aussi consulter la page http://www.sciencenews.org/sn_arc98/4_18_98/mathland.htm.

Retour en page principale