Ticket #57: jacobson.tst

  option(prot);

  proc TestJ(matrix m, list #)
  {
  "/////////////////////////////////////////////////";
  "// m: ";
  print(m);
  "// #: [", string(#), "]";
  def J = jacobson(m, #);   
  
  "typeof (result): [", typeof(J), "], size(J): [", size(J), "]!";
  
  "//  a Jacobson Form D for m: ";
  print(J[2]);   
  "// Check: ";
  print(J[1]*m*J[3] - J[2]); // check that U*M*V = D
  
  
  }


  proc Test(matrix m)
  {
  TestJ(m));   

  TestJ(m,-1); 
  TestJ(m, 0); 
  TestJ(m, 1); 
  TestJ(m, 2); 


  TestJ(m, "a"); 

/*
  TestJ(m,-1, "a"); 
  TestJ(m, 0, "a"); 
  TestJ(m, 1, "a"); 
  TestJ(m, 2, "a"); 


  TestJ(m, 0, 2, "a"); 
  TestJ(m,-1, 2, "a"); 
  TestJ(m, 1, 2, "a"); 
  TestJ(m, 2, 2, "a"); 

  TestJ(m,-1, 0); 
  TestJ(m, 0, 0); 
  TestJ(m, 1, 0); 
  TestJ(m, 2, 0); 
*/  



  
  }
  
  LIB "jacobson.lib";
  printlevel = 1;  echo = 2;



  ring r = 0,(x,d),Dp;

  Test(0); 
//  Test(matrix()); 
  Test(matrix(0)); 

  matrix m[2][2] = d,x,0,d; 
  print(m);
  Test(m);

  def R = nc_algebra(1,1);   setring R; // the 1st Weyl algebra

  Test(0); 
//  Test(matrix()); 
  Test(matrix(0)); 

  matrix m[2][2] = d,x,0,d; print(m);
  
  Test(m);
  
  list J = jacobson(m); // returns a list with 3 entries
  print(J[2]); // a Jacobson Form D for m
  print(J[1]*m*J[3] - J[2]); // check that U*M*V = D
  /*   now, let us do the same for the shift algebra  */

  kill R, r;

  ring r2 = 0,(x,s),Dp;

  matrix m[2][2] = s,x,0,s; print(m); // matrix of the same for as above
  Test(m);

  def R2 = nc_algebra(1,s);   setring R2; // the 1st shift algebra
  
  Test(0); 
//  Test(matrix()); 
  Test(matrix(0)); 
  
  matrix m[2][2] = s,x,0,s; print(m); // matrix of the same for as above
  
  Test(m);
  list J = jacobson(m);
  print(J[2]); // a Jacobson Form D, quite different from above
  print(J[1]*m*J[3] - J[2]); // check that U*M*V = D

  kill R2, r2;

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

  ring r = 0,(x,d),Dp;
  def R=nc_algebra(1,1); setring R;
  R; // the Weyl algebra in x and d
  
  Test(matrix(0)); // TODO: BUG HERE!!!

  matrix m[2][2]=d,x,0,d; 
  print(m);
  Test(m);
  list J = jacobson(m); // returns a list with 3 entries
  print(J[2]); // a Jacobson Form D
  // TODO: ARE U AND V REALLY UNIMODULAR?
  J[1]; // U
  J[3]; // V

  print(J[1]*m*J[3] - J[2]); // check that U*M*V = D

  matrix m[3][2]=x, x^4+x^2+21, x^4+x^2+x, x^3+x, 4*x^2+x, x;
  print(m); // M
  Test(m);
  list JJ = jacobson(m); // returns a list with 3 entries

  // TODO: BUG: WHY THE FOLLOWING JACOBSON FORM IS NOT DIAGONAL???!
  print(JJ[2]); // a Jacobson Form D
  print(JJ[1]*m*JJ[3] - JJ[2]); // check that U*M*V = D

  // TODO: ARE U AND V REALLY UNIMODULAR?
  JJ[1]; // U
  JJ[3]; // V

  list S=smith(m,1);

  print(S[2]); // Smith Normal Form S of M
  print(S[1]*m*S[3] - S[2]); // check that U*M*V = S

  // TODO: ARE U AND V REALLY UNIMODULAR?
  S[1]; // U
  S[3]; // V


$$$