help matrix
//-------------------------------------------------------
5
// en tapant le chiffre 5, Scilab attribue la valeur 5 à la variable ans
// (pour answer)
ans^2
// ans élevé au carré donne 25
abs(-5)
// valeur absolue
m=10^6
sqrt(m)
// racine carrée 
y=%e
log(y)
// %e est la constante e
sin(%pi)
// noter que le résultat n'est pas exactement zéro
1+%eps
// %eps est la précision : 1+%eps est indistinguable de 1
1+%eps==1   // est vrai !
1+2*%eps==1 // est faux.
//-------------------------------------------------------
//-------------------------------------------------------
v=[3.8,-4,%pi/6]
size(v)
// dimensions de v : 1 ligne et 3 colonnes
w=v'
// transposition

x=[1,2]; y=[3,4,5]
z=[x,y];
// construction d'un vecteur par blocs

t=[4:9]
// vecteur des réels entre 4 et 9 par pas de 1 (pas implicite)
t=[4:1:9]
// vecteur des réels entre 4 et 9 par pas de 1 (pas explicite)
t=[0:0.1:1]
// vecteur des réels entre 0 et 1 par pas de 0.1

u=sqrt(2)*[1:2:8]'
size(u)

t=rand(1,5)
// un vecteur à 1 ligne et 5 colonnes 
// et contenant des nombres au hasard dans [0,1] 
//-------------------------------------------------------
[1,2;3,4]
[11,12,13,14,15,16,17,18,19;21,22,23,24,25,...
26,27,28,29;31,32,33,34,35,36,37,38,39]
// dans une instruction trop longue pour tenir dans une ligne,
// mettre ... avant de passer à la ligne
diag([5 4 3])
// matrice diagonale
eye(6,6)
// des 1 sur la diagonale
B=eye(6,7)
A=ones(3,7)
C=[A;(-6)*B]
D=zeros(2,5)
E=rand(D)
rand(2,5)
//-------------------------------------------------------
M=<À COMPLÉTER>
M=<À COMPLÉTER>
//-------------------------------------------------------
u=2*%pi*rand() // un nombre au hasard dans [0,2*pi]
w=[cos(u) sin(u)]
norm(w)   // norme classique |$L^2$|
norm(w,1) // norme |$L^1$|

t=[0:%pi/2:2*%pi]
v=sin(t)
[m,k]=max(v)
// la valeur maximale des éléments du vecteur v est m
// et elle est atteinte pour l'élément d'indice k : m=v(k)
[m,k]=min(v)
sign(v)
// signe 1 (+) ou -1 (-) et sign(0)=0
//-------------------------------------------------------
1==0    // la réponse à l'assertion ``1 égale 0'' est F false
1~=0    // la réponse à l'assertion ``1 différent de 0'' est F false
1==0 & 1~=0   // et : la réponse est F false
1==0 | 1~=0   // ou : la réponse est T true

t=[0:%pi/2:2*%pi]
v=sin(t)
v>0
// renvoie un vecteur de T (true) ou F (false) selon que 
// l'élément correspondant de v est ou non >0
v>=0
// convertit les T et F en 1 et 0
v(v>=0)
// extrait les élément positifs ou nuls de v
//-------------------------------------------------------
w=1:9
sum(w)
// somme de tous les éléments de w
cumsum(w)
// vecteur donnant les sommes cumulées
//-------------------------------------------------------
A=rand(2,3)
B=sin(A)
A+B
//-------------------------------------------------------
G=[ones(1,4); 2*ones(1,4)] 
sum(G,'c')
// somme sur les lignes : le résultat est un vecteur colonne ('c' pour column)
sum(G,'r')
// somme sur les colonnes : le résultat est un vecteur  ligne ('r' pour row)
//-------------------------------------------------------
A'
//-------------------------------------------------------
//-------------------------------------------------------
A'*A
A*A'
C=rand(2,3)
A'*C
//-------------------------------------------------------
//-------------------------------------------------------
//-------------------------------------------------------
w=1:2:9
w(2)
w($) // dernier élément
w($-1) // avant-dernier élément
//-------------------------------------------------------
E=[11:19;21:29;31:39;41:49;51:59;61:69]
E(1,1) // l'élément de la ligne 1 colonne 1
E(3,4) // l'élément de la ligne 3 colonne 4
E(1,:) // la ligne 1
E(:,5) // la colonne 5
E(2:4,:)
// la sous-matrice formée des lignes allant de 2 à 4
E(2:3,7:9)
// la sous-matrice formée des éléments appartenant 
// aux lignes allant de 2 à 3 et aux colonnes de 7 à 9
E([1,3,5],[2,4,6,8])
// la sous-matrice formée des éléments appartenant 
// aux lignes 1 3 5 et aux colonnes 2 4 6 8
E(:,$) // dernière colonne
E(:,$-1) // avant-dernière colonne
E(2:$,:) // les lignes de la deuxième à la dernière
E(2:($-1),:) // les lignes de la deuxième à l'avant-dernière
//-------------------------------------------------------
A=int(20*rand(1,10))
// partie entière
[sa,ia]=gsort(A)
// tri : sa est le résultat du tri, ia les indices correspondants
//-------------------------------------------------------
//-------------------------------------------------------
x=[1 2 3]
x.*x

y=[-6 12 8]
x.*y

A=rand(2,3);
B=ones(2,3);
A.*B
//-------------------------------------------------------
x=[1 2 3]
y=1 ./x
// un blanc suit le 1, car sinon 1. serait interprété en 1.0 et
// l'opération serait / (résolution de système linéaire)
x.*y

A=rand(2,3);
B=rand(A);
A./B
//-------------------------------------------------------
x=[1 2 3]
y=x.^2
z=x.^[5 10 -2]

A=rand(2,3);
A.^3
B=rand(A);
A.^B
//-------------------------------------------------------
//-------------------------------------------------------

//-------------------------------------------------------
function res=f(x)
   res=x*x*x;
endfunction   

f(4) // =  64
//-------------------------------------------------------
function [a,b,c,d]=g(x)
   a=x;b=x*x;c=x*x*x;d=1;
endfunction   

[x,y,z,t]=g(4) // x=4  y=16  z=64  t=1
//-------------------------------------------------------
U=2.5;
i=0;// calcule le plus petit entier > U
while i <= U
    i = i + 1;
end
//-------------------------------------------------------
  U=grand(1,1,'unf',0,1);// tire *une* v.a. de loi uniforme sur [0,1]
  U=grand(1,3,'unf',0,1);// tire un vecteur 1x3 de v.a. de loi uniforme 
                         // sur [0,1] :      U(1), U(2), U(3)
//-------------------------------------------------------




//-------------------------------------------------------
function res=rand_disc_iter(P)
// P est une loi discrète représentée par une vecteur *ligne*
// Version itérative.
  U=grand(1,1,'unf',0,1);// On tire selon une loi uniforme entre 0 et 1
  // "On inverse la fonction de répartition"
  Q=0;i=0;
  while Q <= U
    i = i + 1;
    Q = Q + <À COMPLÉTER>; // on veut |$Q=\P(N\leq i)$|
  end
  res=i;
endfunction;

p=10;
P=ones(1,p)/p;
for i=1:1000 do X(i)=rand_disc_iter(P); end;
mean(X) // la moyenne des tirages
v=variance(X) // la variance des tirages
//-------------------------------------------------------
//-------------------------------------------------------
X=[2,7,8];
max(X)  // = 8
//-------------------------------------------------------
2>7 // = | F |
X > 7 // = | F F T |
cumsum(X) > 7  // = | F T T |
//-------------------------------------------------------
find(cumsum(X) > 7) // = |  2  3 |
X(find(X > 7)) // =  |  8 |
X(find(cumsum(X) > 7))  // =  |  7  8 |
//-------------------------------------------------------
 X=1:10;// 1,2,3, .... , 10
 Z=0;
 for i=1:length(X) do 
    Z = Z + X(i);
 end;
Z  // = |  55 |

 S=0;Z=0;
 for i=1:length(X) do 
    Z = Z + X(i);
    S(i)=Z;
 end;
S  // = |1   3   6  10  15  21  28  36  45  55 |
//-------------------------------------------------------
printf("%d+%d ?=? %d\n",1,1,2);  // Output: 1+1 ?=? 2
printf("%d+%d ?=? %d\n",1,2,2);  // Output: 1+2 ?=? 2
//-------------------------------------------------------
function res=rand_disc(P)
// P est une loi discrète représenté par une vecteur *ligne*
// Version "sans boucle"
  Q=[0,cumsum(P)];// On calcule la fonction de répartition de la loi
                  // |$Q(i)=\P(N\leq i-1)$|, |$i-1$| car |$Q(1)=0$|.
  U=grand(1,1,'unf',0,1);// On tire selon une loi uniforme entre 0 et 1
  res=max(<À COMPLÉTER>);// res=le plus grand des |$i$| tel que |$U \geq Q(i) = \P(N\leq i-1)$|
endfunction;
//-------------------------------------------------------
function test_0()
// Un exemple de 10 tirages uniforme sur |$\{1,2,3\}$|
  p=3;
  P=ones(1,p)/p;// loi uniforme sur |$\{1,2,3\}$|

  N=1000;
  for i=1:N do // 10 tirages uniformes sur |$\{1,2,3\}$|
    X(i) = rand_disc(P);
    printf("%d.",X(i));
  end
  printf("\n");
endfunction
//-------------------------------------------------------
function [T_A,T_B] = random_times(p,q)
// On simule selon des lois uniformes sur 1:p et 1:q
// On en déduit les tirages de T_A et T_B

  P=ones(1,p)/p;// Loi uniforme sur 1:p
  Q=ones(1,q)/q;// Loi uniforme sur 1:q
  i_1=rand_disc(P);j_1=rand_disc(Q);
  i_2=rand_disc(P);j_2=rand_disc(Q);

  // Calcul de T_A et T_B en fonction de (i_1,j_1) (i_2,j_2)
  t_1_A = <À COMPLÉTER>
  t_2_A = <À COMPLÉTER>
  T_A = min(t_1_A,t_2_A);
  
  t_1_B = <À COMPLÉTER>;
  t_2_B = <À COMPLÉTER>;
  T_B = min(t_1_B,t_2_B);
endfunction;
//-------------------------------------------------------
x=12;
if x==3 then 
   printf("Bravo, vous avez gagné 2 euros.\n");
elseif x==4  then 
   printf("Bravo, vous avez gagné 1 euros.\n");
else 
   printf("Vous avez perdu, sorry.\n");
end;
//-------------------------------------------------------
function qui_gagne_MC(p,q,N)
  statA=0;statB=0;
  for i=[1:N] do
    [T_A,T_B] = random_times(p,q);
    if <À COMPLÉTER> then // B gagne
      <À COMPLÉTER>
    elseif <À COMPLÉTER> then  // A gagne
      <À COMPLÉTER>
    end;
  end; 
  erreur=1/sqrt(N);// erreur maximum probable
  pA=statA/N;pB=statB/N;
  printf("pA=%.3f+-%.3f, pB=%.3f+-%.3f\n", ...
               statA/N, erreur, statB/N, erreur);
  if (pA+erreur < pB-erreur) then
     printf("p=%d,q=%d : ",p,q);
     printf("Il est très probable que B gagne en moyenne.\n");
  elseif (pA-erreur > pB+erreur) then
     printf("p=%d,q=%d : ",p,q);
     printf("Il est très probable que A gagne en moyenne.\n");
  else
     printf("p=%d,q=%d, N=%d : On ne peut pas conclure: ",p,q,N);
     printf("il faut augmenter le nombre de tirages N.\n");
  end
  
endfunction;

function test_1()
  N=100000;
  qui_gagne_MC(2,3,N);
  qui_gagne_MC(3,4,N);
  qui_gagne_MC(4,5,N);
  qui_gagne_MC(4,4,N);
endfunction;
//-------------------------------------------------------
function [H_AB] = qui_gagne_elem(p,q)
  Ascore=0;Bscore=0;nbre_cas_egalite=0;
  valeurs_TAB=zeros(2,1);nAB=0;

  // On génére toutes les positions pour les 2 {\oe}ufs
  for i_1=1:p do
    for j_1=1:q do
      for i_2=1:p do
        for j_2=1:q do
          T_A=min((j_1-1)*p+i_1,(j_2-1)*p+i_2);
          T_B=min((i_1-1)*q+j_1,(i_2-1)*q+j_2);
          if T_A > T_B then Bscore=Bscore+1;
          elseif T_A < T_B then Ascore=Ascore+1;
          else// T_A == T_B
            nbre_cas_egalite=nbre_cas_egalite+1;
          end
          // on garde les valeurs du couple dans un tableau
          nAB=nAB+1;valeurs_TAB(:,nAB)=[T_A;T_B];
        end
      end
    end
  end
      
  // On calcule la loi du couple (T_A,T_B) a partir des valeurs prises
  H_AB=histo(p*q,valeurs_TAB);

  // Verification que les lois marginales sont identiques.
  // Elles doivent l'etre!
  H_A=<À COMPLÉTER>; // comment calculer la loi de T_A ?
  H_B=<À COMPLÉTER>; // comment calculer la loi de T_B ?
  if norm(H_A-H_B) >= 1.e-7 then;
    printf("Warning: les lois marginales doivent être égales.\n");
  end

  // Imprime des informations à l'écran
  printf('p=%d, q=%d, ',p,q);
  report(p*q,Ascore,Bscore,nbre_cas_egalite,H_AB);
endfunction;
//-------------------------------------------------------
function H_AB=histo(v_max,samples)
// Calcule la loi du couple (T_A,T_B) a partir des valeurs prises
// Les valeurs de samples sont supposées entières
// comprises entre 1 et v_max
  H_AB=0
  nbre=size(samples);nbre=nbre(2);
  for k=1:v_max
    for l=1:v_max
      // Calcul du nbre de tirages valant (k,l) / Taille
      H_AB(k,l)= ...
          length( find((samples(1,:)==k) & (samples(2,:)==l))) ./ nbre; 
    end;
  end; 
endfunction;
//-------------------------------------------------------
function report(Taille,Ascore,Bscore,nbre_cas_egalite,H_AB)
  out_fmt="Taille=%d, diff=%d, egalite=%d : ";
  if(Bscore>Ascore)      then out_fmt=out_fmt+"C''est B qui gagne\n";
  elseif (Ascore>Bscore) then out_fmt=out_fmt+"C''est A qui gagne\n";
                         else out_fmt=out_fmt+"Egalité\n"; end
  printf(out_fmt,Taille,Bscore-Ascore,nbre_cas_egalite);    
endfunction;
//-------------------------------------------------------
function test_2()
  p=2;q=p+1;
  loi_TAB = qui_gagne_elem(p,q);

  for p=[2:9] do qui_gagne_elem(p,p+1);end;
  for p=[2:8] do qui_gagne_elem(p,p+2);end;
  for p=[2:7] do qui_gagne_elem(p,p+3);end;
endfunction;
//-------------------------------------------------------
A=[1,2;2,4];
// |  1  2 |
// |  2  4 |
A .*. A 
// |   1   2   2   4 |
// |   2   4   4   8 |
// |   2   4   4   8 |
// |   4   8   8  16 |
//-------------------------------------------------------
function qui_gagne_optim(p,q)
// Copyright JPC
    pq=p*q;
    // Les temps d'arrivée dans les cases
    A=1:pq;A=matrix(A,p,-1); //  pour le joueur A 
    B=A;B=matrix(B,q,-1);B=B';  // ... et pour le joueur B
    // On genere toutes les positions du couple d'oeufs. 
    // Une position c'est un couple de couple (O1(k),O2(k))
    O1=ones(1,pq).*.(1:pq); // .*. : produit de Kronecker entre matrices.
    O2=(1:pq).*.ones(1,pq); 
    // pour A et B on calcule le temps min pour atteindre le premier oeuf 
    // A(O1) [A(O2)]: temps mis par A pour atteindre le 1er [2eme] oeuf
    mA=min(A(O1),A(O2));
    // B(O1) [B(O2)]: temps mis par B pour atteindre le 1er [2eme] oeuf
    mB=min(B(O1),B(O2));
    // les cas ou A et B ont mis le meme temps à atteindre un oeuf
    I=find(mA==mB);
    nbre_cas_egalite=size(I,'*');
    mA(I)=[];mB(I)=[];
    // On regarde qui est le vainqueur
    [m,I]=min(mA,mB);
    // On compte les cas gagnants pour les 2 joueurs
    Ascore=length(find(I == 1)); 
    Bscore=length(find(I == 2));
  
    printf('p=%d, q=%d,',p,q);
    report(p*q,Ascore,Bscore,nbre_cas_egalite,0);
endfunction;
//-------------------------------------------------------
function test_6()
  // Ca va beaucoup plus vite ...
  for p=[2:20] do qui_gagne_optim(p,p+1);end;
  for p=[2:20] do qui_gagne_optim(p,p+2);end;
  for p=[2:20] do qui_gagne_optim(p,p+3);end;
endfunction;
//-------------------------------------------------------
function [H_AB] = qui_gagne_permutation(taille,sigma_A_1,sigma_B_1)
  Ascore=0;
  Bscore=0;
  nbre_cas_egalite=0;
  valeurs_TAB=zeros(2,taille);
  nAB=0;

  // On génère toutes les positions
  for i_1=1:taille do
    for i_2=1:taille do
      T_A=min(sigma_A_1(i_1),sigma_A_1(i_2));
      T_B=min(sigma_B_1(i_1),sigma_B_1(i_2));
      if T_A > T_B then Bscore=Bscore+1;
      elseif T_A < T_B then Ascore=Ascore+1;
      else// T_A == T_B
        nbre_cas_egalite=nbre_cas_egalite+1;
      end
      nAB=nAB+1;valeurs_TAB(:,nAB)=[T_A;T_B];
    end
  end

  // Génère la loi du couple T_A, T_B
  H_AB=histo(taille,valeurs_TAB);
  report(taille,Ascore,Bscore,nbre_cas_egalite,H_AB);
endfunction;
//-------------------------------------------------------
function test_3()
  taille=6;
  // On tire |$\sigma_A^{-1}$| et |$\sigma_B^{-1}$| directement
  // ce qui ne change pas la loi et évite de calculer l'inverse
  sigma_A_1=grand(1,"prm",[1:taille]');
  sigma_B_1=grand(1,"prm",[1:taille]');
  H_AB=qui_gagne_permutation(taille,sigma_A_1,sigma_B_1)
endfunction;
//-------------------------------------------------------
function sigma=permutation(p,q)
// permutation associée à la transposition d'une matrice |$p \times q$|.
sigma=0;
for i=[1:p]
  for j=[1:q]          
      t_1_A=(j-1)*p+i;// t_1_A = i
      t_1_B=(i-1)*q+j;// t_1_B = sigma_0(i)
      sigma(t_1_A)=t_1_B;// sigma = sigma_0
  end
end
endfunction;
//-------------------------------------------------------
function [H_AB] = qui_gagne_elem_variante(p,q)
  identite=1:p*q;// la permutation identité
  sigma=permutation(p,q);
  printf('p=%d, q=%d, ',p,q);
  H_AB=qui_gagne_permutation(p*q,identite,sigma)
endfunction;
//-------------------------------------------------------
function test_5()
  p=4;q=6;
  H2 = qui_gagne_elem(p,q);// la méthode du début
  H1 = qui_gagne_elem_variante(p,q);// la méthode utilisant |$\sigma$|
  // On vérifie que les 2 méthodes calculent la même chose.
  norm(H1-H2) // on doit trouver 0 (ou qque chose de très petit)
endfunction;
//-------------------------------------------------------
function test_4()
  taille=6;
  sigma_A_1=grand(1,"prm",[1:taille]');

  // la permutation définie par |$\sigma_B^{-1}=\sigma_0(\sigma_A^{-1})$| est toujours meilleure que |$\sigma_A$|
  sigma_0=[taille,1:taille-1]; // [5 1 2 4 6 3]  marcherait aussi
  sigma_B_1=sigma_0(sigma_A_1);

  sigma_A_1'
  // C'est toujours B qui gagne
  H_AB=qui_gagne_permutation(taille,sigma_A_1,sigma_B_1);
endfunction;
//-------------------------------------------------------
function t_A=temps_A_zig_zag(i,j,p,q)
// temps d'atteinte par A de la cas i,j en zig-zag
   if modulo(j,2) == 0 then
       t_A=(j-1)*p+(p-i+1);
   else
       t_A=(j-1)*p+i;
   end
endfunction;

function t_B=temps_B_zig_zag(i,j,p,q)
// temps d'atteinte par B de la cas i,j en zig-zag
    if modulo(i,2) == 0 then
       t_B=(i-1)*q+(q-j+1);
    else
       t_B=(i-1)*q+j;
    end
endfunction;

function sigma=permutation_zig_zag(p,q)
// permutation associée au zig-zag pour A et B.
  sigma=0;
  for i=[1:p]
    for j=[1:q]
      t_1_A=temps_A_zig_zag(i,j,p,q);
      t_1_B=temps_B_zig_zag(i,j,p,q);
      sigma(t_1_A)=t_1_B;// sigma = sigma_0
    end
  end
endfunction;
//-------------------------------------------------------
function [H_AB] = qui_gagne_elem_zz(p,q)
  identite=1:p*q;// la permutation identité
  
  sigma=permutation(p,q);
  printf('\ndirect, p=%d, q=%d, ',p,q);
  H_AB=qui_gagne_permutation(p*q,identite,sigma)
  
  sigma=permutation_zig_zag(p,q);
  printf('zigzag, p=%d, q=%d, ',p,q);
  H_AB=qui_gagne_permutation(p*q,identite,sigma)
endfunction;

function test_7()
  for p=[2:9] do qui_gagne_elem_zz(p,p+1);end;
  for p=[2:8] do qui_gagne_elem_zz(p,p+2);end;
  for p=[2:7] do qui_gagne_elem_zz(p,p+3);end;
endfunction;
//-------------------------------------------------------
function main()
  test_0()
  test_1()
  test_2()
  test_3()
  test_4()
  test_5()
  test_6()
  test_7()
endfunction
//-------------------------------------------------------