Changeset e67578d in git


Ignore:
Timestamp:
Mar 5, 2009, 9:36:11 PM (14 years ago)
Author:
Viktor Levandovskyy <levandov@…>
Branches:
(u'spielwiese', 'd1ba061a762c62d3a25159d8da8b6e17332291fa')
Children:
acff7e71919841ba60543811ceda2b000dbb15d1
Parents:
fe880793d1fe5f78cc5038ae512e18c6cbaab42c
Message:
*levandov: fixes and new proc divideUnits


git-svn-id: file:///usr/local/Singular/svn/trunk@11516 2c84dea3-7e68-4137-9b89-c4e89433aadc
File:
1 edited

Legend:

Unmodified
Added
Removed
  • Singular/LIB/jacobson.lib

    rfe88079 re67578d  
    11//////////////////////////////////////////////////////////////////////////////
    2 version="$Id: jacobson.lib,v 1.9 2009-02-12 20:27:17 levandov Exp $";
     2version="$Id: jacobson.lib,v 1.10 2009-03-05 20:36:11 levandov Exp $";
    33category="System and Control Theory";
    44info="
    55LIBRARY: jacobson.lib     Algorithms for Smith and Jacobson Normal Form
    6 AUTHOR: Kristina Schindelar,     Kristina.Schindelar@math.rwth-aachen.de
    7 
    8 THEORY: We work over a ring R, that is a principal ideal domain.
    9 @*   If R is commutative, we suppose R to be a polynomial ring in one variable.
    10 @* If R is non-commutative, we suppose R to be in two variables, say x and d.
    11 @* We treat then the basering as principal ideal ring with d a polynomial
    12 @* variable and x an invertible one. That is, we work in the Ore localization of R
     6AUTHOR: Kristina Schindelar,     Kristina.Schindelar@math.rwth-aachen.de,
     7@*           Viktor Levandovskyy,      levandov@math.rwth-aachen.de
     8
     9THEORY: We work over a ring R, that is an Euclidean principal ideal domain.
     10@* If R is commutative, we suppose R to be a polynomial ring in one variable.
     11@* If R is non-commutative, we suppose R to have two variables, say x and d.
     12@* We treat then the basering as the Ore localization of R
    1313@* with respect to the mult. closed set S = K[x] without 0.
     14@* Thus, we treat basering as principal ideal ring with d a polynomial
     15@* variable and x an invertible one.
    1416@* Note, that in computations no division by x will actually happen.
    15 @*   Given a rectangular matrix M over R, one can compute unimodular (that is invertible)
    16 @* square matrices U and V, such that U*M*V=D is a diagonal matrix.
    17 @*   Depending on the ring, the diagonal entries of D have certain properties.
     17@*
     18@* Given a rectangular matrix M over R, one can compute unimodular (that is
     19@* invertible) square matrices U and V, such that U*M*V=D is a diagonal matrix.
     20@* Depending on the ring, the diagonal entries of D have certain properties.
    1821
    1922REFERENCES:
    20 @* N. Jacobson, 'The theory of rings', AMS, 1943.
    21 @* Manuel Avelino Insua Hermo, 'Varias perspectives sobre las bases de Groebner: Forma normal de Smith, Algorithme de Berlekamp y algebras de Leibniz'. PhD thesis, Universidad de Santiago de Compostela, 2005.
     23@* (1) N. Jacobson, 'The theory of rings', AMS, 1943.
     24@* (2) Manuel Avelino Insua Hermo, 'Varias perspectives sobre las bases de Groebner :
     25@* Forma normal de Smith, Algorithme de Berlekamp y algebras de Leibniz'.
     26@* PhD thesis, Universidad de Santiago de Compostela, 2005.
    2227
    2328
    2429PROCEDURES:
    2530smith(M[,eng1,eng2]);  compute the Smith Normal Form of M over commutative ring
    26 jacobson(M[,eng]);       compute a weak Jacobson Normal Form of M over non-commutative ring
     31jacobson(M[,eng]); compute a weak Jacobson Normal Form of M over non-commutative ring
     32divideUnits(L);     create ones out of units in the output of smith or jacobson
     33
    2734
    2835SEE ALSO: control_lib
     
    3037
    3138  LIB "poly.lib";
    32   LIB "involut.lib";
     39  LIB "involut.lib"; // involution
     40  LIB "dmodapp.lib"; // engine
     41  LIB "qhmoduli.lib"; // Min
     42
     43proc divideUnits(list L)
     44"USAGE: divideUnits(L); list L
     45RETURN: matrix or list of matrices
     46ASSUME: L is an output of @code{smith} or a @code{jacobson} procedures, that is
     47@* either L contains one rectangular matrix with elements only on the main diagonal
     48@* or L consists of three matrices, where L[1] and L[3] are square invertible matrices
     49@* while L[2] is a rectangular matrix with elements only on the main diagonal
     50PURPOSE: divide out units from the diagonal and reflect this in transformation matrices
     51EXAMPLE: example divideUnits; shows examples
     52"
     53{
     54  // check assumptions
     55  int s = size(L);
     56  if ( (s!=1) && (s!=3) )
     57  {
     58    ERROR("The list has neither 1 nor 3 elements");
     59  }
     60  // check sizes of matrices
     61
     62  if (s==1)
     63  {
     64    list LL; LL[1] = L[1]; LL[2] = L[1];
     65    L = LL;
     66  }
     67 
     68
     69  // divide units out
     70  matrix M = L[2];
     71  int ncM = ncols(M);   int nrM = nrows(M);
     72  //  matrix TM[nrM][nrM]; // m times m matrix
     73  matrix TM = matrix(freemodule(nrM));
     74  int i; int nsize = Min(intvec(nrM,ncM));
     75  poly p; number n; intvec lexp;
     76 
     77  for(i=1; i<=nsize; i++)
     78  {
     79    p = M[i,i];
     80    lexp = leadexp(p);
     81    //    TM[i,i] = 1;
     82    // check: is p a unit?
     83    if (p!=0)
     84    {
     85      if ( lexp == 0)
     86      {
     87        // hence p = n*1
     88        n = leadcoef(p);
     89        TM[i,i] = 1/n;
     90      }
     91    }
     92  }
     93  int ppl = printlevel-voice+2;
     94  dbprint(ppl,"divideUnits: extra transformation matrix is: ");
     95  dbprint(ppl, TM);
     96  L[2] = TM*L[2];
     97  if (s==3)
     98  {
     99    L[1] = TM*L[1];
     100    return(L);
     101  }
     102  else
     103  {
     104    return(L[2]);
     105  }
     106}
     107example
     108{ "EXAMPLE:"; echo = 2;
     109  ring R=(0,m,M,L1,L2,m1,m2,g), D, lp; // two pendula example
     110  matrix P[3][4]=m1*L1*D^2,m2*L2*D^2,(M+m1+m2)*D^2,-1,
     111  m1*L1^2*D^2-m1*L1*g,0,m1*L1*D^2,0,0,
     112  m2*L2^2*D^2-m2*L2*g,m2*L2*D^2,0;
     113  list s=smith(P,1);  // returns a list with 3 entries
     114  print(s[2]); // a diagonal form, close to the Smith form
     115  print(s[1]); // U, left transformation matrix
     116  list t = divideUnits(s);
     117  print(t[2]); // the Smith form of the matrix P
     118  print(t[1]); // U', modified left transformation matrix
     119}
    33120
    34121proc smith(matrix MA, list #)
    35122"USAGE: smith(M[, eng1, eng2]);  M matrix, eng1 and eng2 are optional integers
    36 RETURN: matrix or list
    37 ASSUME: The current ring is assumed to be the commutative polynomial ring in
    38         one variable
    39 PURPOSE: compute the Smith Normal Form of M with transformation matrices (optionally)
    40 NOTE: If the optional integer eng1 is non-zero, returns the list {U,D,V},
    41 where U*M*V = D and the diagonal field entries of D are not normalized.
    42 @*    Otherwise only the matrix D, that is the Smith Normal Form of M, is returned.
    43 @*    Here, U and V are square unimodular (invertible) matrices. The procedure works for rectangular matrix M.
    44 @* The optional integer eng2 determines the engine, that computes the Groebner basis:
    45 @* 0 means 'std' (default), 1 means 'groebner' and 2 means 'slimgb'.
     123RETURN: matrix or list of matrices, depending on arguments
     124ASSUME: Basering is a commutative polynomial ring in one variable
     125PURPOSE: compute the Smith Normal Form of M with (optionally) transformation matrices
     126THEORY: Groebner bases are used for the Smith form like in (2).
     127NOTE: By default, just the Smith normal form of M is returned.
     128@* If the optional integer @code{eng1} is non-zero, the list {U,D,V} is returned
     129@* where U*M*V = D and the diagonal field entries of D are not normalized.
     130@* The normalization of the latter can be done with the 'divideUnits' procedure.
     131@* U and V above are square unimodular (invertible) matrices.
     132@* Note, that the procedure works for a rectangular matrix M.
     133@*
     134@* The optional integer @code{eng2} determines the Groebner basis engine:
     135@* 0 (default) ensures the use of 'slimgb' , otherwise 'std' is used.
    46136DISPLAY: If @code{printlevel}=1, progress debug messages will be printed,
    47137@* if @code{printlevel}>=2, all the debug messages will be printed.
    48138EXAMPLE: example smith; shows examples
     139SEE ALSO: divideUnits, jacobson
    49140"
    50141{
     
    120211  matrix m[3][2]=x, x^4+x^2+21, x^4+x^2+x, x^3+x, 4*x^2+x, x;
    121212  list s=smith(m,1);
    122   print(s[2]);  // Smith form of m
     213  print(s[2]);  // non-normalized Smith form of m
    123214  print(s[1]*m*s[3] - s[2]); // check U*M*V = D
     215  list t = divideUnits(s);
     216  print(t[2]); // the Smith form of m
    124217}
    125218
     
    552645
    553646
    554 
    555 static proc engine(module I, int i)
    556 {
    557   module J;
    558   if (i==0)
    559   {
    560     J = std(I);
    561   }
    562   if (i==1)
    563   {
    564     J = groebner(I);
    565   }
    566   if (i==2)
    567   {
    568     J = slimgb(I);
    569   }
    570   return(J);
    571 }
     647// VL : engine replaced by the one from dmodapp.lib
     648// cases are changed
     649
     650// static proc engine(module I, int i)
     651// {
     652//   module J;
     653//   if (i==0)
     654//   {
     655//     J = std(I);
     656//   }
     657//   if (i==1)
     658//   {
     659//     J = groebner(I);
     660//   }
     661//   if (i==2)
     662//   {
     663//     J = slimgb(I);
     664//   }
     665//   return(J);
     666// }
    572667
    573668proc jacobson(matrix MA, list #)
    574669 "USAGE:  jacobson(M, eng);  M matrix, eng an optional int
    575670RETURN: list
    576 ASSUME: Basering is a noncommutative ring in two variables.
    577 PURPOSE: compute a weak Jacobson Normal Form of M over a noncommutative ring
     671ASSUME: Basering is a (non-commutative) ring in two variables.
     672PURPOSE: compute a weak Jacobson Normal Form of M over the basering
     673THEORY: Groebner bases and involutions are used, generalizing an idea of (2)
    578674NOTE: A list L of matrices {U,D,V} is returned. That is L[1]*M*L[3]=L[2], where
    579675@*      L[2] is a diagonal matrix and L[1], L[3] square invertible (unimodular) matrices.
    580676@*      Note, that M can be rectangular.
    581 @* The optional integer @code{eng} determines the engine, that computes the Groebner basis:
    582 @* 0 means 'std' (default), 1 means 'groebner' and 2 means 'slimgb'.
     677@* The optional integer @code{eng2} determines the Groebner basis engine:
     678@* 0 (default) ensures the use of 'slimgb' , otherwise 'std' is used.
    583679DISPLAY: If @code{printlevel}=1, progress debug messages will be printed,
    584680@* if @code{printlevel}>=2, all the debug messages will be printed.
    585681EXAMPLE: example jacobson; shows examples
     682SEE ALSO: divideUnits, smith
    586683"
    587684{
Note: See TracChangeset for help on using the changeset viewer.