PROGRAM Quaternions;

{
Epreuve informatique de l'Ecole Polytechnique           92.013
--------------------------------------------------------------


1/   On cherche à généraliser la notion de nombre complexe en
     définissant le quaternion, quadruplet de nombres réels
     (a,b,c,d), que l'on note :

              q=a+bi+cj+dk

     comme combinaison linéaire de quatre quaternions de base 1,
     i, j  et k

     (1,0,0,0)=1    (0,1,0,0)=i   (0,0,1,0)=j   (0,0,0,1)=k


     L'addition des quaternions correspond exactement à l'ad-
                                                 4
     dition usuelle des vecteurs dans l'espace  R :

     (a,b,c,d) + (a',b',c',d') = (a+a',b+b',c+c',d+d')


     La multiplication des quaternions est une extension de
     celle des nombres complexes, où a est la partie réelle et
     (b,c,d) la partie imaginaire du quaternion. Cette
     multiplication, qui est distributive à gauche et à droite
     sur l'addition, admet 1 comme élément neutre et vérifie :


         i.i = -1     j.j = -1     k.k = -1

         i.j = k    j.i = -k
         j.k = i    k.j = -i
         k.i = j    i.k = -j



     -    Ecrire un programme Pascal capable de calculer la
          somme et le produit de deux quaternions quelconques.




2/   Déterminer les quaternions qui vérifient l'équation :

                     2
                    q  = -1

     Montrer que la notion de module d'un nombre complexe peut
     être étendue aux quaternions et écrire un programme Pascal
     pour calculer le module d'un quaternion.


3/   La multiplication de deux quaternions n'étant pas commu-
     tative, on peut définir le quotient à gauche et le quo-
     tient à droite de deux quaternions.
     Existent-t-ils toujours ? Sont-ils uniques ? dans quels cas
     sont-ils égaux ?
     Ecrire le programme Pascal permetttant de les calculer.



4/   On étend maintenant la notion d'exponentielle aux quater-
     nions avec la définition suivante :

             inf
             ====   n               2     3           n
             \     q               q     q           q
     Exp(q)= /    --- =  1  + q + --- + --- + ... + --- + ...
             ====  n!              2     3!          n!
             n=0

     Cette somme a-t-elle un sens pour tout quaternion q ?
     Ecrire un programme Pascal capable de la calculer.
     Que peut-on dire de son module ?


5/   Tracer sur un graphe la variation en fonction de la va-
     riable réelle t des quatre composantes de la fonction
     y= Exp(qt) pour une valeur quelconque de q.

     A quelle condition cette fonction peut-elle être périodi-
     que, et comment peut-on alors calculer cette période en
     fonction de q ?


+------------------------------------------------------------+
¦              Imprimer tous les résultats                   ¦
¦     en indiquant chaque fois à quoi ils correspondent      ¦
+------------------------------------------------------------+

                            -=-=-=-
 }


uses
  crt, modubase, entrees, matric;

const
     nmax=30;
     eps=1E-5;

Type
    quaternion = vecteur;

VAR
   un,i,j,k    : quaternion ;
   n           :integer;
   question    : char;
   pyt         : array[1..4,1..4] of quaternion;
   d           : shortint ;           {nb décimales}

Function entier(x:real):boolean;
begin
     if abs(x)<maxint then
        entier:=abs(x-round(x))<eps
     else
         entier:=false;
end;


Function strq(q:quaternion):string;
var
   t,tx :string;
   m:shortint;
   x:integer;
begin
     tx:='';
     if q[1]<>0 then
      begin
           if entier(q[1]) then
              begin
                 x:=round(q[1]);
                 str(x:2,tx);
              end
           else
              str(abs(q[1]):d+2:d,tx);
      end;
     for m:=2 to 4 do
         begin
           if abs(q[m])>eps then
              begin
                   if entier(q[m]) then
                      begin
                           x:=round(q[m]);
                           str(abs(x):2,t);
                      end
                   else
                      str(abs(q[m]):d+2:d,t);
                   if q[m]<0 then
                      tx:=tx+'-
                   else
                      tx:=tx+'+';
                   if entier(q[m]) then
                      begin
                           if round(abs(q[m]))<>1 then
                              tx:=tx+t;
                      end
                    else
                       tx:=tx+'('+t+')';
                   tx:=tx+copy('1ijk',m,1)
              end;
        end;
        if tx='' then
         strq :='0'
        else
         if copy(tx,1,1)='+' then
            strq:=copy(tx,2,length(tx)-1)
         else
            strq := tx;
end;


Procedure show(q:quaternion);
begin
     write(strq(q));
end;

Procedure entrer(var q:quaternion);
var
   m:shortint;
begin
     for m:=1 to 4 do
      begin
          write('terme en ',copy('1ijk',m,1),'=');
          Entre_reel(q[m]);
       end;
     writeLN('On vient d''entrer q=',strq(q));
end;


Procedure genq(VAR q:quaternion;s,a,b,c:real);
begin
     q[1]:=s; q[2]:=a; q[3]:=b; q[4]:=c;
end;


Procedure scal(var q:quaternion;s:real;q1:quaternion);
begin
     genq(q,s*q1[1],s*q1[2],s*q1[3],s*q1[4])
end;

procedure add(var q:quaternion;q1,q2:quaternion);
var
   m:shortint;
begin
     for m:=1 to 4 do
         q[m]:=q1[m]+q2[m];
end;

procedure mul(var q:quaternion;q1,q2:quaternion);
var
   x:quaternion;
   m1,m2:shortint;
begin
     genq(q,0,0,0,0);
     for m1:=1 to 4 do
         for m2:=1 to 4 do
          begin
             scal(x,q1[m1]*q2[m2],pyt[m1,m2]);
             add(q,q,x);
          end;
end;


Procedure showmat(x:matrice);
var
   m1,m2: shortint;
begin
     for m1:=1 to 4 do
      begi
       for m2:=1 to 4 do
        begin
           write(x[m1,m2]:9:5,' : ');
        end;
        writeLN('');
      end
end;

Procedure mat(var x:matrice;q:quaternion);
begin
     x[1,1]:= q[1]; x[1,2]:=-q[2]; x[1,3]:=-q[3]; x[1,4]:=-q[4];
     x[2,1]:= q[2]; x[2,2]:= q[1]; x[2,3]:= q[4]; x[2,4]:=-q[3];
     x[3,1]:= q[3]; x[3,2]:=-q[4]; x[3,3]:= q[1]; x[3,4]:= q[2];
     x[4,1]:= q[4]; x[4,2]:= q[3]; x[4,3]:=-q[2]; x[4,4]:= q[1];
end;

Procedure Conjugue(Var q:quaternion);
var
   m:shortint;
begin
     for m:=2 to 4 do
         q[m]:=-q[m];
end;

Function Module(Var q:quaternion):real;
var
   m:shortint;
   sum:real;
begin
     sum:=0;
     for m:=1 to 4 do
         sum:=sum+sqr(q[m]);
     Module:=sqrt(sum);
end;


Procedure inverse(var q1:quaternion;q:quaternion);
begin
     scal(q1,1/sqr(module(q)),q);
     Conjugue(q1);
end;


Procedure quotientg(var q:quaternion;dividende,diviseur:quaternion);
begin
     inverse(q,diviseur);
     mul(q,q,dividende);
end;

Procedure quotientd(var q:quaternion;dividende,diviseur:quaternion);
begin
     inverse(q,diviseur);
     mul(q,dividende,q);
end;

Procedure Expq(var q1:quaternion;q:quaternion);
var
   x:quaternion;
   n:shortint;
begin
     scal(x,1,un)
     scal(q1,1,un);
     for n:=1 to 50 do
         begin
              mul(x,x,q);
              scal(x,1/n,x);
              add(q1,q1,x);
          end;
end;


Procedure Courbe(q:quaternion;tmin,tmax,ymax:real);
var
   m : shortint;
   t,t1,pas       : real;
   a,b,c,d,h,u,r  : real;
   y,y1,z,z1,qt :quaternion;
begin
      ModeGraphique;
      Couleur(Rouge);
      Fenetre(tmin,tmax,-ymax,ymax);
      X_Axe(0,0,1) ;
      Y_Axe(0,0,1);
      pas := (tmax-tmin)/200;
      t := tmin;
      couleur(-Brillant);
      Deplace(tmin,-ymax/2);
      Ecris('exp(qt) avec q='+strq(q));
      t1:=t;
      a:=q[1];
      b:=q[2];
      c:=q[3];
      d:=q[4];
      r:=sqrt(sqr(b)+sqr(c)+sqr(d));
      While (t<tmax) and not KeyPressed do
      begin
           scal(qt,t,q);
           Deplace(tmin,ymax/2);
           EcrisReel(t);
           Ecris('->');
           u:=exp(t*a);
           h:=u*sin(t*r)/r;
           genq(z,u*cos(t*r),b*h,c*h,d*h);
           Expq(y,qt);
           for m:=1 to 4 do
            begin
                 case m of
                 1:  couleur(Brillant);
                 2:  couleur(Vert);
                 3:  couleur(Jaune);
                 4:  couleur(Rouge);
                 end;
                 Deplace(t1,y1[m]);
                 Trace(t,y[m]);
                 Deplace(t1,z1[m]);
                 Trace(t,z[m]);
            end;
           scal(y1,1,y);
           scal(z1,1,z);
       {   couleur(-Jaune);
           Deplace(t1,z1);
           Trace(t,z); }
           z1:=z;
           t1:=t;
           t:=t+pas;
      end ;
end;


Function vec(m:shortint):shortint;
begin
     vec:= ((m-2) mod 3)+2;
end;

procedure genpyt;
var
   m,m1,m2,m3:shortint;
begin
     genq(un,1,0,0,0);
     genq(i ,0,1,0,0);
     genq(j ,0,0,1,0);
     genq(k ,0,0,0,1);

     pyt[1,1]:=un;
     pyt[1,2]:=i;
     pyt[1,3]:=j;
     pyt[1,4]:=k;

     for m:=2 to 4 do
       begin
        pyt[m,1]:=pyt[1,m];
        m1:=vec(m)
        m2:=vec(m+1);
        m3:=vec(m+2);
        pyt[m1,m2]:=pyt[1,m3];
        scal(pyt[m2,m1],-1,pyt[1,m3]);
        scal(pyt[m,m],-1,un);
       end;
end;

Procedure Question1;
var
   m1,m2:shortint;
   x,sum : quaternion;
begin
     writeLN('============== Table de Pythagore ========');
     writeLN('');
     write('1=');show(un);writeln('');
     write('i=');show(i );writeln('');
     write('j=');show(j );writeln('');
     write('k=');show(k );writeln('');

     for m1:=1 to 4 do
      begin
       for m2:=1 to 4 do
        begin
           show(Pyt[m1,m2]);
           write(' : ');
        end;
        writeLN('');
      end;


end;

Procedure Question2;
var
   qg,qd,q,produit,dividende,diviseur:quaternion;
begin
     WriteLN('Quaternions à diviser');
     WriteLN('Dividende');
     Entrer(dividende);
     WriteLN('Diviseur');
     Entrer(diviseur);
     Quotientg(qg,dividende,diviseur);
     Quotientd(qd,dividende,diviseur);;
     Mul(produit,qg,diviseur);
     writeLN('A gauche : <',strq(qg),'> * <',strq(diviseur),
                                      '> = <',strq(produit),'>');
     Writeln('modules :',module(qg):8:5,'*',module(diviseur):8:5,
             ' = ',module(qg)*module(diviseur):8:5,' =? ',module(produit):8:5);
     Mul(produit,diviseur,qd);
     writeLN('A droite : <',strq(diviseur),'> * <',strq(qd),
                                      '> = <',strq(produit),'>');
     Writeln('modules :',module(qd):8:5,'*',module(diviseur):8:5,
             ' = ',module(qg)*module(diviseur):8:5,' =? ',module(produit):8:5);
end;



Procedure Question3;
var
   q,q1:quaternion;
begin
     writeLN('---------- Décomposition de Exp(qt) -----------------------');
     writeLN('                                                           ');
     writeLN('                                                           ');
     writeLN('     Exp(qt) = Exp(t*(a+bi+cj+dk))                         ');
     writeLN('                                                           ');
     writeLN('       =      exp(at)*cos(rt)                              ');
     writeLN('       +i * b*exp(at)*sin(rt)/sqr(r)                       ');
     writeLN('       +j * c*exp(at)*sin(rt)/sqr(r)                       ');
     writeLN('       +k * d*exp(at)*sin(rt)/sqr(r)                       ');
     writeLN('                                                           ');
     writeLN('                                                           ');
     writeLN('           en posant r=sqrt(sqr(b)+sqr(c)+sqr(d))          ');
     writeLN('                                                           ');
     writeLN(' --------------------------------------------------------- ');
     d:=8;
     Entrer(q);
     Expq(q1,q);
     writeLN('Exp(',strq(q),') = ',strq(q1));
     Pause;
     Courbe(q,-3,3,2);
end;



Procedure Question4;
begin
end;

Procedure Question5;
begin
end;

Procedure Question6;
begin

end;

Procedure Presentation;
begin
    ModeTexte;
    Efface;
    WriteLN('Les quaternions                        ');
    WriteLN('                                       ');
    WriteLN('                                       ');
    WriteLN(' (1)  Table de multiplication          ');
    WriteLN(' (2)  Quotients                        ');
    WriteLN(' (3)  Exponentielle                    ');
    WriteLN(' (4)                                   ');
    WriteLN(' (5)                                   ');
    WriteLN(' (6)                                   ');
    WriteLN('                                       ');
    GenPyt;
end;



begin

  Initgraphique;
  while true do
  begin
     Presentation;
     write('Question choisie ? ');
     readln(question);
     Efface;
     d :=3;
     writeln('======== Question n°',question,'===============');
     writeln('');

     case question of
          '1':  Question1;
          '2':  Question2;
          '3' : Question3;
          '4' : Question4;
          '5' : Question5;
          '6' : Question6;
     end;
     Pause;
  end;

end.