[Axiom-mail] Looking for beginning with axiom.

[fmaltey]

Hello, 

I compile axiom from cvs on my Linux box, because axiom has now an 
open licence.  After some hours I get a running axiom. 

And what I do now ? 
I can't get any useful help file. 
I only find help file in *.dvi format about compiling.

the )hd command doesn't work.

My first idea is to compare axiom with mupad (or maple) for my use 
with my students (first year after Baccalauréat in France).

So I'm looking for a tutorial whith function names, basic types, 
matrix, solver, and so. Do you have any idea ?

If you want to have any idea about my works with my students
I post 2 examples with mupad :
The first one about magic matrix. 
The second about a stupid minpoly algorithm.
May you help me to translate it, if it's possible...

I'm also looking for an interface in emacs because 
texmacs is a little too slow on my PC.

Thanks a lot. 

Francois Maltey

 
// 
// septembre 2003 
// 
// recherche des matrices magiques d'ordre impair 
// à partir de l'épreuve du CAPES externe de 1999 (5-eme partie)
//

export (linalg) ;

//
// la somme des n lignes est 1+2+...+n^2 ; d'où la somme d'une ligne.
//
sommeTheorique := n -> (n*(n^2+1))/2 ; 

//
// Somme d'une ligne, d'une colonne et de la diagonale secondaire
//
sommeLignei := (M, i) -> _plus (M[i,k] $k=1..nrows(M)) ;
sommeColi   := (M, i) -> _plus (M[k,i] $k=1..nrows(M)) ;
antitrace   := M      -> _plus (M[nrows(M)-k+1,k] $k=1..nrows(M)) ;
ensCoeff    := M      -> {M[i,j] $i=1..nrows (M) $j=1..ncols (M)} ; 

//
// Les différents tests à effectuer pour tester si une matrice est magique.
//
estMagique := 
M ->     bool (nrows (M)=ncols (M)) 
    and _and (bool (sommeLignei (M,i) = sommeTheorique (nrows (M))) 
                                                           $i=1..nrows (M))
    and _and (bool (sommeColi (M,i) = sommeTheorique (nrows (M)))
                                                           $i=1..nrows (M))
    and bool (tr (M) = sommeTheorique (nrows (M)))
    and bool (antitrace (M) = sommeTheorique (nrows (M))
    and bool (ensCoeff (M) = {$1..nrows(M)^2}))
        
//
// Peut se simplifier en montrant que l'ensemble des sommes des lignes, des
// colonnes et des diagonales est un singleton.
//
estMagique := 
M ->    bool (nrows (M)=ncols (M)) 
    and bool (card ({sommeLignei (M,i) $i=1..nrows (M),  
                     sommeColi   (M,i) $i=1..nrows (M),
                     tr (M), antitrace (M)}) = 1)
    and bool (ensCoeff (M) = {$1..nrows(M)^2})
//
// Un exemple de matrice magique
//
M3 := matrix ([[4,3,8],[9,5,1],[2,7,6]]) ;
estMagique (M3) ;

//
// Calculs valables pour n premier avec 3 uniquement
// numlib::igcdmult (3,n)[2] est le coefficient de l'équation de Bézout 
// associé à 3 dans 3x + ny = 1 ; c'est l'inverse modulo n de 3.
//
// Le calcul (X - 1 mod n) + 1 calcule la congruence comprise entre 1 et n.
//
// La fonction beta est réservée par le système en analyse complexe.
//
alpha := (i,j,n) -> (numlib::igcdmult(3,n)[2]*(2*i-j) - 1 mod n) + 1 ;
betaa := (i,j,n) -> (numlib::igcdmult(3,n)[2]*(2*j-i) - 1 mod n) + 1 ;

//
// Une autre possibilité consiste à calculer 
// expr (1 / (Dom::IntergerMod (n)(3)) pour calculer m (cf énoncé et corrigé)
// ou solve (3*m=1, m, Domain=Dom::IntergerMod (n))
//
// Une dernière possibilité utilise linsolve en précisant en troisième 
// argument Domain=Dom::IntergerMod (n)
//

//
// Ces sommes sont constantes car la fonction alpha(.,j,n) réordonne
// tous les i entre 1 et n
//
[_plus(alpha (i,j,7) $i=1..7) $j=1..7] ;
          
//
// Construction de la matrice W énumérée dans l'ordre 
//
matriceW := n -> matrix (n,n,(i,j)->n*(i-1)+j) ;

//
// Construction de la matrice magique G
// 
matriceG := n -> matrix (n,n,(i,j)->(alpha(i,j,n)-1)*n+betaa(i,j,n));

//
// Vérification que les sommes sur les sous-diagonales et 
// les anti-sous-diagonales sont constantes. 
//

//
// Construction d'une matrice magique à partir de deux matrices magiques
//
produitMagique := 
(A,B) -> 
  matrix(nrows(A)*nrows(B),nrows(A)*nrows(B),
         (u,v)->            A[(u-1 mod nrows(A))+1,(v-1 mod nrows(A))+1]
              + nrows(A)^2*(B[(u-1 div nrows(A))+1,(v-1 div nrows(A))+1]-1)) ;

//
// Quelques vérifications de matrices magiques
//
estMagique (produitMagique (M3,M3)) ;
estMagique (produitMagique (M3,matriceG(5))) ;
estMagique (produitMagique (matriceG(5),matriceG(7))) ;

////////////////////////////////////////////////////////////////////////

//
// Calcul du polynôme minimal de la matrice A 
//
// 1/ Étude du rang des familles de matrices (A^k) pour k variant de 1 à p.
// 2/ Construction du polynôme minimal unitaire
//
// 1/ Recherche du rang = max (p tq (A^k)k=0..p est libre) 
//
// Correspond à la bibliothèque de calcul de l'algèbre linéaire
//
export (linalg) ;

//
// Un exemple de matrice
//
M := matrix ([[1, 2, 3], [3, 3, 3], [3, 2, 1]]) ; 

//
// Transforme une matrice en vecteur et réciproquement pour appliquer
// les fonctions usuelles de l'algèbre linéaire : calcul du rang et co.
//
matrice2vecteur := M -> matrix ([M[i,j] $j=1..ncols(M) $i=1..nrows(M)]) ;
vecteur2matricebr := (V, n) -> [[V[i+k*n] $i=1..n] $k=0..n-1] ; 
vecteur2matrice := V -> matrix (vecteur2matricebr (V, sqrt (vecdim (V)))) ;
//
// Test de vérification
//
bool (vecteur2matrice (matrice2vecteur (M)) = M) ; 
//
// Cette fonction calcule successivement le rang des famille (A^k)k=0..p
// et renvoie la dernière (c'est-à-dire la plus grande) des valeurs de p
// pour laquelle la famille est libre : c'est la valeur précédente p pour
// lequel la famille est liée.
//
rangM := proc (M) local p ;
  begin
    p := 0 ;
    while (nops (basis ([matrice2vecteur (M^k) $k=0..p])) = p+1) 
       do p := p+1 end_while ; 
    p-1 
  end ; 
//
// Tests de vérification
//
rangM (M) ;                                      // réponse 2
rangM (matrix ([[i+j+1 $i=1..5] $j=1..5])) ;     // réponse 2
rangM (matrix ([[i+j+1 $i=1..5] $j=1..5])) ;     // réponse 2
rangM (matrix ([[(i+j)^3 $i=1..5] $j=1..5])) ;   // réponse 4
//
// Autre méthode de calcul de rangM
// La liste L contient la liste des vecteurs associés aux matrices M^k
//
rangMa  := M -> rangMbr (M, [matrice2vecteur (M^0)]) ;
rangMbr := (M, L) -> 
  if nops (L) = nops (basis (L)) 
    then rangMbr (M, [matrice2vecteur (M^(nops (L)))].L)
    else nops (L)-2 
  end_if ; 
//
// Tests de vérification
//
rangMa (M) ;                                      // réponse 2
rangMa (matrix ([[i+j+1 $i=1..5] $j=1..5])) ;     // réponse 2
rangMa (matrix ([[i+j+1 $i=1..5] $j=1..5])) ;     // réponse 2
rangMa (matrix ([[(i+j)^3 $i=1..5] $j=1..5])) ;   // réponse 4
//
//
// 2/ Construction du polynôme minimal unitaire
//
//
// But : rechercher de proche en proche le plus petit des p 
// tel que la matrice M^p est combinaison linéaire des M^k précédents
// Les calculs n'ont pas lieu sur les matrices mais sur les vecteurs associés
//
PolMiniM := proc (M) local p, S, b, a_continuer ;
  begin
    p := 1 ;
    a_continuer := TRUE ; 
    S := matrice2vecteur (M^0) ; 
    RES := [] ; 
    while RES = [] 
      do 
        b := matrice2vecteur (M^p) ; 
        RES  := matlinsolve (S, b, Unique) ;   
        S := S.b ;
        p := p+1 ;
     end_while ;
    X^(p-1) - _plus (RES[i]*X^(i-1) $i=1..vecdim (RES)) ; 
  end_proc ;
//
// Le paramètre Unique à matlinsolve renvoie soit la seule 
// solution du système linéaire soit [].
// L'opérateur . concatène matrices et vecteurs 
// et augmente le nombre de colonnes.
// Les coordonnées du vecteur solution sont les coefficients de X^(k-1)
// dans le polynôme minimal.
//
Res := PolMiniM (M) ; 
_plus (coeff (Res, X, k) * M^k $k=0..ncols (M)) ; 

M1 := matrix ([[i+j+1 $i=1..5] $j=1..5]) ;
M2 := matrix ([[i+j+1 $i=1..5] $j=1..5]) ;
M3 := matrix ([[(i+2*j)^3 $i=1..5] $j=1..5]) ;  
//
// Vérification : le résultat est la matrice nulle.
// Le polynôme calculé est bien le polynôme minimal.
//
Res := PolMiniM (M1) ; 
_plus (coeff (Res, X, k) * M1^k $k=0..ncols (M1)) ; 

Res := PolMiniM (M2) ; 
_plus (coeff (Res, X, k) * M2^k $k=0..ncols (M2)) ; 

Res := PolMiniM (M3) ; 
_plus (coeff (Res, X, k) * M3^k $k=0..ncols (M3)) ; 

PolMiniMb := M -> PolMiniMbr (M, [], matrice2vecteur (M^0), 1) ;

PolMiniMbr := (M, R, S, p) 
  -> if R=[] 
        then 
          PolMiniMbr (M, matlinsolve (S, matrice2vecteur (M^p), Unique),
                      S.matrice2vecteur(M^p), p+1)
        else R
      end_if ; 

PolMiniMb (M1) ;

Res := PolMiniM (M1) ; 
_plus (coeff (Res, X, k) * M1^k $k=0..ncols (M1)) ;