Changeset 1d5418 in git


Ignore:
Timestamp:
Dec 15, 2004, 8:02:33 PM (20 years ago)
Author:
Viktor Levandovskyy <levandov@…>
Branches:
(u'fieker-DuVal', '117eb8c30fc9e991c4decca4832b1d19036c4c65')(u'spielwiese', 'd25190065115c859833252500a64cfb7b11e3a50')
Children:
501f7c553e083fd3e3bfad340f57a827228a7911
Parents:
55e3087801d49cf32ed31572abf9343bd66a4d0a
Message:
*minor changes:
 *usage of inner procedures for computing next Ad images has been changed
 +instead the procedure ApplyAd is introduced.


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

Legend:

Unmodified
Added
Removed
  • Singular/LIB/center.lib

    r55e308 r1d5418  
    11///////////////////////////////////////////////////////////////////////////////
    2 version="$Id: center.lib,v 1.11 2004-12-07 16:02:02 plural Exp $";
     2version="$Id: center.lib,v 1.12 2004-12-15 19:02:33 plural Exp $";
    33category="Noncommutative";
    44info="
    55LIBRARY:  center.lib      computation of central elements of G-algebras and Q-algebras.
    6 AUTHORS:  Oleksandr Motsak,        motsak@mathematik.uni-kl.de.
     6AUTHOR:  Oleksandr Motsak,        motsak@mathematik.uni-kl.de.
    77OVERVIEW:
    88 This is a library for computing the central elements and centralisators of elements in various noncommutative algebras.
    99 Implementation is mainly based on algorithms, written in the frame of the diploma thesis by Oleksandr Motsak (advisor: Prof. S.A. Ovsienko, using the definitions from diploma thesis by V.Levandovskyy), at Kyiv Taras Shevchenko University:  'An algorithm for the computation of the center of noncommutative polynomial algebra'.
    10  @* The last version of this library can be found via internet.
     10 @* The last version of this library can be found via internet at author's home page.
    1111
    1212SUPPORT: Forschungsschwerpunkt 'Mathematik und Praxis'
     
    2323
    2424
    25 ///////////////////////////////////////////////////////////////////////////////
     25/******************************************************/
    2626// stuff
    2727
    28 ///////////////////////////////////////////////////////////////////////////////
     28/******************************************************/
    2929static proc myValue ( def s, list # )
    30 "
     30  "
    3131        return s or (typeof(s))(#)
    3232"
    3333{
    34         def @p = s;
    35         if ( size(#) > 0 )
    36         {
    37                 if ( typeof( #[1] ) == typeof(s) )
    38                 {
    39                         @p = #[1];
    40                 };
    41         };
    42         return (@p);
    43 };
    44 
    45 ///////////////////////////////////////////////////////////////////////////////
     34  def @p = s;
     35  if ( size(#) > 0 )
     36    {
     37      if ( typeof( #[1] ) == typeof(s) )
     38        {
     39          @p = #[1];
     40        };
     41    };
     42  return (@p);
     43};
     44
     45/******************************************************/
    4646static proc myInt ( list # )
    47 "
     47  "
    4848        return 0 or int(#)
    4949"
    5050{
    51         int @p = 0;
    52         return (myValue ( @p, # ));
    53 };
    54 
    55 ///////////////////////////////////////////////////////////////////////////////
     51  int @p = 0;
     52  return (myValue ( @p, # ));
     53};
     54
     55/******************************************************/
    5656static proc myPoly ( list # )
    57 "
     57  "
    5858        return 0 or poly(#)
    5959"
    6060{
    61         poly @p = 0;
    62         return (myValue( @p, # ));
    63 };
    64 
    65 ///////////////////////////////////////////////////////////////////////////////
     61  poly @p = 0;
     62  return (myValue( @p, # ));
     63};
     64
     65/******************************************************/
    6666static proc myRing ( list # )
    67 "
     67  "
    6868        return basring or ring(#)
    6969"
    7070{
    71         def @p = basering;
    72         return (myValue ( @p, # ));
    73 };
    74 
    75 ///////////////////////////////////////////////////////////////////////////////
     71  def @p = basering;
     72  return (myValue ( @p, # ));
     73};
     74
     75/******************************************************/
    7676static proc myString ( list # )
    77 "
     77  "
    7878        return basring or ring(#)
    7979"
    8080{
    81         string @p = "";
    82         return (myValue ( @p, # ));
    83 };
    84 
    85 ///////////////////////////////////////////////////////////////////////////////
     81  string @p = "";
     82  return (myValue ( @p, # ));
     83};
     84
     85/******************************************************/
    8686static proc myIdeal ( list # )
    87 "
     87  "
    8888        return 0 or ideal(#)
    8989"
    9090{
    91         ideal @p;
    92         return (myValue ( @p, # ));
    93 };
    94 
    95 ///////////////////////////////////////////////////////////////////////////////
    96 /// for debug
    97 
    98 ///////////////////////////////////////////////////////////////////////////////
     91  ideal @p;
     92  return (myValue ( @p, # ));
     93};
     94
     95/******************************************************/
     96// for debug
     97
     98/******************************************************/
    9999static proc toprint()
    100100{
    101     return (0);
    102 };
    103 
    104 ///////////////////////////////////////////////////////////////////////////////
    105 static proc Print( def a )
    106 {
    107         if ( toprint() )
    108         {
    109                 print (a);
    110         };
     101  return (0);
     102};
     103
     104/******************************************************/
     105static proc Print( list # )
     106{
     107  if ( toprint() )
     108    {
     109      print (#);
     110    };
    111111};   
    112112
    113 ///////////////////////////////////////////////////////////////////////////////
     113/******************************************************/
    114114static proc BCall(string Name, list #)
    115115{
    116         if( toprint() )
    117         {
    118                 "CALL: ", Name, "( ", string(#), " )";
    119         };
    120 };
    121 
    122 ///////////////////////////////////////////////////////////////////////////////
     116  if( toprint() )
     117    {
     118      "CALL: ", Name, "( ", string(#), " )";
     119    };
     120};
     121
     122/******************************************************/
    123123static proc ECall(string Name, list #)
    124124{
    125    if( toprint() )
    126    {
    127                 "RET: ", Name, "(...)={ ", string(#), " }";
    128    };
    129 };
    130 
    131 
    132 ////////////////////////////////////////////////////////////////////
    133 ////////////////////////////////////////////////////////////////////
     125  if( toprint() )
     126    {
     127      "RET : ", Name, "(...)={ ", string(#), " }";
     128    };
     129};
     130
     131
     132/******************************************************/
    134133// "my" degres functions and lie procedures...
    135134
    136 ////////////////////////////////////////////////////////////////////
    137 static proc maxDeg( proc givendeg, def z )
    138 "
     135/******************************************************/
     136static proc maxDeg( def z )
     137  "
    139138        returns: int : max of givendeg( z_i ), among all z_i \in z
    140139                returns -1 if z is empty or contain only 0 polynomial
    141140"
    142141{
    143         int max = -1; int d;
    144 
    145         for(int i=size(z); i>0; i-- )
    146         {
    147                 d = givendeg(z[i]);
    148                 if( d > max )
    149                 {
    150                         max = d;
    151                 };
    152         };
    153 
    154         return(max);
    155 };
    156 
    157 ////////////////////////////////////////////////////////////////////
    158 static proc minDeg( proc givendeg, def z )
    159 "
    160         returns: int : min of givendeg( z_i ), among all z_i \in z
    161                 returns -1 if z is empty or contain only 0 polynomial
    162 "
    163 {
    164         if( size(z) == 0 )
    165         {
    166             return(-1);
    167         };
    168 
    169         int min = -2; int d;
    170 
    171         for(int i = size(z); i>0; i-- )
    172         {
    173                 d = givendeg(z[i]);
    174                 if( (min == -2) or (d < min) )
    175                 {
    176                         min = d;
    177                 };
    178 
    179         };
    180        
    181         return(min);
    182 };
    183 
    184 ////////////////////////////////////////////////////////////////////
    185 static proc myMonomDeg( poly p)
    186 "
    187         input: poly p
    188         return: just number of veriables in lead monomial of p!
    189         Note: for use by maxList and minList
    190 "
    191 {
    192         if( p == 0 )
    193         {
    194             return (-1);
    195         };
    196        
    197         intvec e = leadexp(p);
    198         int    d = 0;
    199        
    200         for ( int i = size(e); i>0; i-- )
    201         {
    202                 d = d + e[i];
    203         };
    204 
    205         return (d);
    206 };
    207 
    208 ////////////////////////////////////////////////////////////////////
    209 static proc stdDeg ( poly p)
    210 "
    211         input: poly p
    212         return: degree of p
    213         Note: just in order to have it as 'proc' (for use by maxList)
    214 "
    215 {
    216         return (deg(p));
    217 };
    218 
    219 ////////////////////////////////////////////////////////////////////
    220 static proc myDeg ( poly p)
    221 "
    222         input: poly p
    223         return: 'my degree' of p
    224         no real use... just for testing...
    225 "
    226 {
    227         return( maxDeg(myMonomDeg,p) );
    228 };
    229 
    230 ///////////////////////////////////////////////////////////////////////////////
    231 ///////////////////////////////////////////////////////////////////////////////
    232 // some bisiness for my_var and lie_rank
    233 
    234 ////////////////////////////////////////////////////////////////////
     142  int max = -1;
     143  int d;
     144
     145  for(int i=size(z); i>0; i-- )
     146    {
     147      d = deg(z[i]);
     148      if( d > max )
     149        {
     150          max = d;
     151        };
     152    };
     153
     154  return(max);
     155};
     156
     157/******************************************************/
     158// some bisiness for my_var
     159
     160/******************************************************/
    235161static proc myCoeff ( poly f, poly m )
    236 "
     162  "
    237163        input: poly f,
    238164        return: coeff at m
    239165"
    240166{
    241         m = leadmonom(m);
    242 
    243         int i = size(f);
    244        
    245         while ( (i>0) and (leadmonom(f[i])<m) )
    246         {
    247             i--;
    248         };
    249         if( i == 0 )
    250         {
    251             return ( 0 );           
    252         };
    253 
    254         if(leadmonom(f[i]) == m)
    255         {
    256             return ( leadcoef(f[i]) );
    257         };
    258        
    259         return ( 0 );
    260 };
    261 
    262 
    263 ///////////////////////////////////////////////////////////////////////////////
     167  m = leadmonom(m);
     168
     169  int i = size(f);
     170       
     171  while ( (i>0) and (leadmonom(f[i])<m) )
     172    {
     173      i--;
     174    };
     175  if( i == 0 )
     176    {
     177      return ( 0 );         
     178    };
     179
     180  if(leadmonom(f[i]) == m)
     181    {
     182      return ( leadcoef(f[i]) );
     183    };
     184       
     185  return ( 0 );
     186};
     187
     188/******************************************************/
    264189static proc my_vars ()
    265190{
    266         int  N = nvars( basering );
    267         list V = list();
    268 
    269         for ( int k = 1; k <= N; k++ )
    270         {
    271                 V[k] = var(k);
    272         };
    273        
    274         return (V);
    275 };
    276 
    277 ///////////////////////////////////////////////////////////////////////////////
     191  int  N = nvars( basering );
     192  list V = list();
     193
     194  for ( int k = 1; k <= N; k++ )
     195    {
     196      V[k] = var(k);
     197    };
     198       
     199  return (V);
     200};
     201
     202/******************************************************/
    278203static proc my_commutators ()
    279204{
    280         list V = my_vars();
    281         int  N = size( V );
    282 
    283         matrix M[N][N] = 0;
    284 
    285         poly v, d; int i;
    286        
    287         for ( int k = 2; k <= N; k++ )
    288         {
    289                 v = V[k];
    290                 for ( i = 1; i < k; i++ )
    291                 {
    292                         d = V[i];
    293                         M[k,i] =  v*d - d*v;    // [var(k),var(i)]
    294                         M[i,k] = -M[k,i];       // [var(i),var(k)] ==  -[var(k),var(i)]
    295                 };
    296         };
    297 
    298         return (M);
    299 };
    300 
    301 ///////////////////////////////////////////////////////////////////////////////
     205  list V = my_vars();
     206  int  N = size( V );
     207
     208  matrix M[N][N] = 0;
     209
     210  poly v, d; int i;
     211       
     212  for ( int k = 2; k <= N; k++ )
     213    {
     214      v = V[k];
     215      for ( i = 1; i < k; i++ )
     216        {
     217          d = V[i];
     218          M[k,i] =  v*d - d*v;  // [var(k),var(i)]
     219          M[i,k] = -M[k,i];     // [var(i),var(k)] ==  -[var(k),var(i)]
     220        };
     221    };
     222
     223  return (M);
     224};
     225
     226/******************************************************/
    302227static proc my_associated ()
    303228"
     
    308233"
    309234{
    310         int i, j; poly a; poly d;
    311 
    312         int N    = nvars( basering );
    313         list V   = my_vars();
    314         matrix M = my_commutators();
    315 
    316         matrix A[N][N];
    317 
    318         int cartan;
    319        
    320         list RESULT = list();
    321         int  r_begin = 1;
    322         int  r_end   = N;
    323        
    324 
    325         for ( int k = 1; k <= N; k++ )
    326         {
    327                 A = 0;
    328                 cartan = 1;
    329                 for ( i = 1; i <= N; i++ )
     235  int i, j; poly a; poly d;
     236
     237  int N          = nvars( basering );
     238  list V         = my_vars();
     239  matrix M = my_commutators();
     240
     241  matrix A[N][N];
     242
     243  int cartan;
     244       
     245  list RESULT = list();
     246  int  r_begin = 1;
     247  int  r_end   = N;
     248       
     249
     250  for ( int k = 1; k <= N; k++ )
     251    {
     252      A = 0;
     253      cartan = 1;
     254      for ( i = 1; i <= N; i++ )
     255        {
     256          d = M[k,i];
     257          for ( j = 1; j <= N; j++ )
     258            {
     259              a = myCoeff( d, V[j] );   
     260              A[i,j] =  a;
     261
     262              if( (i!=j) and (a!=0) )
    330263                {
    331                         d = M[k,i];
    332                         for ( j = 1; j <= N; j++ )
    333                         {
    334                                 a = myCoeff( d, V[j] ); A[i,j] =  a;
    335 
    336                                 if( (i!=j) and (a!=0) )
    337                                 {
    338                                     cartan = 0;
    339                                 };                                     
    340                         };
    341                 };
    342                
    343                 if ( cartan )
    344                 {       
    345                     RESULT[r_begin] = list( V[k], cartan, A );
    346                     r_begin++;
    347                 } else
    348                 {
    349                     RESULT[r_end] = list( V[k], cartan, A );
    350                     r_end--;           
    351                 };
     264                  cartan = 0;
     265                };                                     
     266            };
     267        };
     268
     269      if ( cartan )
     270        {       
     271          RESULT[r_begin] = list( V[k], cartan, A );
     272          r_begin++;
     273        } else
     274          {
     275            RESULT[r_end] = list( V[k], cartan, A );
     276            r_end--;           
     277          };
    352278                               
    353         };
    354        
    355         return (RESULT);
    356 
    357 };
    358 
    359 /******************************************************************************/
     279    };
     280       
     281  return (RESULT);
     282
     283};
     284
     285/******************************************************/
    360286static proc my_var_init()
    361287"
     
    364290{
    365291  if( defined( @@@SORTEDVARARRAY ) )
    366   {
    367         kill ( @@@SORTEDVARARRAY );
    368   };
     292    {
     293      kill ( @@@SORTEDVARARRAY );
     294    };
    369295 
    370296  list @@@SORTEDVARARRAY;
     
    374300 
    375301  for ( int i = 1; i<= size(V); i++ )
    376   {
     302    {
    377303      @@@SORTEDVARARRAY[i] = V[i][1];
    378   };
     304    };
    379305 
    380306  Print( "@@@SORTEDVARARRAY: " + string(@@@SORTEDVARARRAY) );
     
    383309};
    384310
    385 /******************************************************************************/
     311/******************************************************/
    386312static proc my_var(int @number )
    387 "
     313  "
    388314  'my' var (with my order of variables),
    389315  before use call my_var_init()
     
    391317{
    392318  if( defined( @@@SORTEDVARARRAY ) )
    393   {
    394         return( @@@SORTEDVARARRAY[@number] );
    395   };
     319    {
     320      return( @@@SORTEDVARARRAY[@number] );
     321    };
    396322
    397323  Print( "Error: my_var_init() was not called before this..." );
     
    400326};
    401327
    402 /******************************************************************************/
     328/******************************************************/
    403329static proc my_var_done()
    404 "
     330  "
    405331  this procedure should be called after finishing using my_var in order to clean up.
    406332"
    407333{
    408334  if( defined( @@@SORTEDVARARRAY ) )
    409   {
    410         kill ( @@@SORTEDVARARRAY );
    411   };
    412 };
    413 
    414 
    415 
    416 ///////////////////////////////////////////////////////////////////////////////
    417 static proc index( int i )
    418 "
    419   i need indexing of arrays from 0 => ARRAY[ index(i) ] - everywhere...
    420 "
    421 {
    422         return (i + 1);
    423 };
    424 
    425 ///////////////////////////////////////////////////////////////////////////////
     335    {
     336      kill ( @@@SORTEDVARARRAY );
     337    };
     338};
     339
     340
     341/******************************************************/
     342// QNF stuff
     343
     344/******************************************************/
    426345static proc myBitParam( int param, int n )
    427 "
     346  "
    428347    return: string: param in bin within no less then n digits (adding zeroes)
    429348"
    430349{
    431         string s = ""; 
    432        
    433         while ( param != 0 )
    434         {
    435             s = string ( param % 2 ) + s;
    436             param = param / 2;
    437             n --;       
    438         };
    439         while ( n > 0 )
    440         {
    441             s = "0" + s;
    442             n --;
    443         };
    444         return ( s );
    445 };
    446 
    447 ///////////////////////////////////////////////////////////////////////////////
     350  string s = "";       
     351       
     352  while ( param != 0 )
     353    {
     354      s = string ( param % 2 ) + s;
     355      param = param / 2;
     356      n --;     
     357    };
     358  while ( n > 0 )
     359    {
     360      s = "0" + s;
     361      n --;
     362    };
     363  return ( s );
     364};
     365
     366/******************************************************/
    448367static proc QNF_init(list #)
    449368{
    450     if( defined(@@@MY_QIDEAL) )
    451     {
    452                 kill (@@@MY_QIDEAL);
    453     };
    454 
    455     if( typeof(basering) == "qring" )
    456     {
    457                 ideal @@@MY_QIDEAL = twostd(myIdeal(#));
    458                 /// setup redSB, redTail options
    459                 export(@@@MY_QIDEAL);
    460        
    461                 option(redSB);
    462                 option(redTail);
    463 
    464     };
    465 };
    466 
    467 ///////////////////////////////////////////////////////////////////////////////
     369  if( defined(@@@MY_QIDEAL) )
     370    {
     371      kill (@@@MY_QIDEAL);
     372    };
     373
     374  if( typeof(basering) == "qring" )
     375    {
     376      ideal @@@MY_QIDEAL = twostd(myIdeal(#));
     377      // setup redSB, redTail options
     378      export(@@@MY_QIDEAL);
     379     
     380      option(redSB);
     381      option(redTail);
     382     
     383    };
     384};
     385
     386/******************************************************/
    468387static proc QNF_done()
    469388{
    470     if( defined(@@@MY_QIDEAL) )
    471     {
    472                 kill (@@@MY_QIDEAL);
    473     };
    474 };
    475 
    476 ///////////////////////////////////////////////////////////////////////////////
    477 static proc QNF_poly ( poly p, list # ) )
    478 {
    479     if( defined(@@@MY_QIDEAL) )
    480     {
    481 /*
     389  if( defined(@@@MY_QIDEAL) )
     390    {
     391      kill (@@@MY_QIDEAL);
     392    };
     393};
     394
     395/******************************************************/
     396static proc QNF_poly ( poly p, list # )
     397{
     398  if( defined(@@@MY_QIDEAL) )
     399    {
     400      /*
    482401        int param = 1 ; // by def: reduce both 1st and 2nd entries
    483402        param = myValue(param, #);
    484403        string bits = myBitParam( param, 8 );
    485 */
    486         BCall( "QNF_poly", p );
    487 
    488         p = NF( p, @@@MY_QIDEAL );
    489 
    490         ECall( "QNF_poly", p );
    491     }; /// QRing
    492 
    493     return(p);
    494 };
    495 
    496 
    497 ///////////////////////////////////////////////////////////////////////////////
     404      */
     405      BCall( "QNF_poly", p );
     406     
     407      p = NF( p, @@@MY_QIDEAL );
     408     
     409      ECall( "QNF_poly", p );
     410    }; // QRing
     411
     412  return(p);
     413};
     414
     415
     416/******************************************************/
    498417static proc QNF_list ( list l, list # )
    499 "
     418  "
    500419    expects list of records with two polynomial entries
    501420    and in case of Qring allpies NF( * , std(0) ) to entries (depends on parameter)
     
    516435{
    517436
    518     if( defined(@@@MY_QIDEAL) )
    519     {
    520 
    521         int param = 1 + 2; // by def: reduce both 1st and 2nd entries
    522         param = myValue(param, #);
    523        
    524         string bits = myBitParam( param, 8 );
    525 
    526         BCall( "QNF_list", "list, " + bits );
    527 
    528         poly temp;
    529        
    530         for ( int i = size(l); i>0 ; i -- )
    531         {
    532 
    533             if ( typeof( l[i] ) == "poly" )
    534             {
    535 
    536                     if ( (bits[8] == "1") or (bits[6] == "1") )
     437  if( defined(@@@MY_QIDEAL) )
     438    {
     439
     440      int param = 1 + 2; // by def: reduce both 1st and 2nd entries
     441      param = myValue(param, #);
     442       
     443      string bits = myBitParam( param, 8 );
     444
     445      BCall( "QNF_list", "list, " + bits );
     446
     447      poly temp;
     448       
     449      for ( int i = size(l); i>0 ; i -- )
     450        {
     451
     452          if ( typeof( l[i] ) == "poly" )
     453            {
     454
     455              if ( (bits[8] == "1") or (bits[6] == "1") )
     456                {
     457                  temp = NF( l[i], @@@MY_QIDEAL );
     458                           
     459                  if ( bits[6] == "1" )
     460                    {// for PBW in qrings: kill out pbw entries where pbw monom was reduced
     461                           
     462                      if( temp != l[i] )
     463                        {
     464                          l = delete( l, i );
     465                          i --;
     466                          continue;
     467                        };
     468                    };
     469                           
     470                  if ( bits[8] == "1" )
    537471                    {
    538                             temp = NF( l[i], @@@MY_QIDEAL );
     472                      l[i] = temp;
     473                    };
     474                };
     475            };
     476
     477          if ( typeof( l[i] ) == "list" )
     478            {   
     479              // 1st
     480             
     481              if ( size(l[i])>0 )
     482                {               
     483                  if ( (bits[8] == "1") or (bits[6] == "1") )
     484                    {
     485                      if( typeof( l[i][1] ) == "poly" )
     486                        {
    539487                           
    540                             if ( bits[6] == "1" )
    541                             {/// for PBW in qrings: kill out pbw entries where pbw monom was reduced
     488                          temp = NF( l[i][1], @@@MY_QIDEAL );
    542489                           
    543                                 if( temp != l[i] )
     490                          if ( bits[6] == "1" )
     491                            {// for PBW in qrings: kill out pbw entries where pbw monom was reduced
     492                               
     493                              if( temp != l[i][1] )
    544494                                {
    545                                     l = delete( l, i );
    546                                     i --;
    547                                     continue;
     495                                  l = delete( l, i );
     496                                  i --;
     497                                  continue;
    548498                                };
    549499                            };
    550500                           
    551                             if ( bits[8] == "1" )
     501                          if ( bits[8] == "1" )
    552502                            {
    553                                 l[i] = temp;
    554                             };
    555                     };
    556             };
    557 
    558             if ( typeof( l[i] ) == "list" )
    559             {   
    560                 /// 1st
    561 
    562                 if ( size(l[i])>0 )
    563                 {               
    564                     if ( (bits[8] == "1") or (bits[6] == "1") )
    565                     {
    566                         if( typeof( l[i][1] ) == "poly" )
    567                         {
    568 
    569                             temp = NF( l[i][1], @@@MY_QIDEAL );
    570                            
    571                             if ( bits[6] == "1" )
    572                             {/// for PBW in qrings: kill out pbw entries where pbw monom was reduced
    573                            
    574                                 if( temp != l[i][1] )
    575                                 {
    576                                     l = delete( l, i );
    577                                     i --;
    578                                     continue;
    579                                 };
     503                              l[i][1] = temp;
    580504                            };
    581505                           
    582                             if ( bits[8] == "1" )
    583                             {
    584                                 l[i][1] = temp;
    585                             };
    586 
    587506                        };
    588507                    };
    589508                };
    590 
    591                 // 2nd
     509             
     510              // 2nd
    592511               
    593                 if ( size(l[i])>1 )
     512              if ( size(l[i])>1 )
    594513                {               
    595                     if ( bits[7] == "1" )
     514                  if ( bits[7] == "1" )
    596515                    {
    597                         if( typeof( l[i][2] ) == "poly" )
     516                      if( typeof( l[i][2] ) == "poly" )
    598517                        {
    599                             temp = NF( l[i][2], @@@MY_QIDEAL );
    600                            
    601                             l[i][2] = temp;
     518                          temp = NF( l[i][2], @@@MY_QIDEAL );
     519                         
     520                          l[i][2] = temp;
    602521                        };
    603522                    };
     
    605524            };
    606525        };
    607 
    608         ECall( "QNF_list", "list" );
    609 
     526       
     527      ECall( "QNF_list", "list" );
     528       
    610529    }; // Qring
    611530   
    612     return ( l );
    613 };
    614 
    615 ///////////////////////////////////////////////////////////////////////////////
    616 /// things for zCommon
    617 
    618 /******************************************************************************/
     531  return ( l );
     532};
     533
     534/******************************************************/
     535// things for zCommon
     536
     537/******************************************************/
     538static proc index( int i )
     539  "
     540  i need indexing of arrays from 0 => ARRAY[ index(i) ] - everywhere...
     541"
     542{
     543  return (i + 1);
     544};
     545
     546/******************************************************/
    619547static proc uni_poly( poly p )
    620 "
     548  "
    621549        returns polynomial with the same monomials but without coefficients.
    622550"
    623551{
    624         poly @tt = poly(0);   
    625         for ( int @k = size(p); @k > 0; @k -- )
    626         {
    627                 @tt = @tt + leadmonom(p[@k]);
    628         };   
    629         return (@tt);   
    630 };
    631 
    632 /******************************************************************************/
     552  poly @tt = poly(0);   
     553  for ( int @k = size(p); @k > 0; @k -- )
     554    {
     555      @tt = @tt + leadmonom(p[@k]);
     556    };   
     557  return (@tt);
     558};
     559
     560/******************************************************/
    633561static proc get_max ( poly @t, number @def )
    634562{
    635     int @n = size( @t );
     563  int @n = size( @t );
     564 
     565  if ( @n == 0 )
     566    {
     567      return (@def);
     568    };
    636569   
    637     if ( @n == 0 )
    638     {
    639         return (@def);
    640     };
    641    
    642     number @max = leadcoef(@t[1]);
    643     number @mm;
    644    
    645     if ( @n > 1)
    646     {
    647         for ( int @i = 2; @i <= @n ;@i ++ )
    648         {
    649             @mm = leadcoef ( @t[@i] );
    650             if ( @mm < 0 )
    651             {
    652                 @mm = -@mm;
    653             };
    654            
    655             if( @mm > @max )
    656             {
    657                 @max = @mm;
    658             };
    659         };
    660     };
    661 
    662     @max = @max + 1;
    663     if ( @max == 0 )
    664     {
    665         @max = @max + 1;
    666     };
    667     return ( 1/@max );
    668 };
    669 
    670 
    671 /******************************************************************************/
     570  number @max = leadcoef(@t[1]);
     571  number @mm;
     572 
     573  if ( @n > 1)
     574    {
     575      for ( int @i = 2; @i <= @n ;@i ++ )
     576        {
     577          @mm = leadcoef ( @t[@i] );
     578          if ( @mm < 0 )
     579            {
     580              @mm = -@mm;
     581            };
     582         
     583          if( @mm > @max )
     584            {
     585              @max = @mm;
     586            };
     587        };
     588    };
     589 
     590  @max = @max + 1;
     591  if ( @max == 0 )
     592    {
     593      @max = @max + 1;
     594    };
     595  return ( 1/@max );
     596};
     597
     598
     599/******************************************************/
    672600static proc get_uni_sum(list @given, int @num)
    673601{
    674     int @l, @k;
    675 
    676     poly @t, @tt;
    677 
    678     @l = size (@given);
    679 
    680     @t = poly(0);
    681     for ( @k = @l; @k > 0; @k -- )
    682     {
    683         if (@num == 1)
    684         {
    685             @tt = @given[@k];
     602  int @l, @k;
     603 
     604  poly @t, @tt;
     605
     606  @l = size (@given);
     607 
     608  @t = poly(0);
     609  for ( @k = @l; @k > 0; @k -- )
     610    {
     611      if (@num == 1)
     612        {
     613          @tt = @given[@k];
    686614        } else
    687         {
     615          {
    688616            @tt = @given[@k][2];       
    689         };
    690        
    691         @t = @t + uni_poly( @tt );
    692     };
    693 
    694     return ( uni_poly(@t) );
    695 };
    696 
    697 
    698 
    699 ////////////////////////////////////////////////////////////////////
     617          };
     618     
     619      @t = @t + uni_poly( @tt );
     620    };
     621 
     622  return ( uni_poly(@t) );
     623};
     624
     625
     626
     627/******************************************************/
    700628static proc LM ( intvec exp )
    701 "
     629  "
    702630        input : given exponent
    703631        return: monom with this exponent...
    704632"
    705633{
    706         poly @f = 1; int @deg;
    707         for( int @i = 1; @i <= size(exp); @i ++ )
    708         {
    709                 @deg = exp[@i];
    710                 if ( @deg > 0 )
    711                 {
    712                         @f = @f * (var(@i)^(@deg));
    713                 };
    714         };
    715 
    716         return (@f);
    717 };
    718 
    719 
    720 
    721 
    722 /////////////////////////////////////////////////////////////////////////////
    723 /////////////////////////////////////////////////////////////////////////////
    724 /// reduced centers - computation of "minimal" set of generators
    725 
    726 LIB "general.lib"; /// for sort
    727 
    728 //////////////////////////////////////////////////////////////////////////
     634  poly @f = 1; int @deg;
     635  for( int @i = 1; @i <= size(exp); @i ++ )
     636    {
     637      @deg = exp[@i];
     638      if ( @deg > 0 )
     639        {
     640          @f = @f * (var(@i)^(@deg));
     641        };
     642    };
     643
     644  return (@f);
     645};
     646
     647/******************************************************/
     648// reduced centers - computation of "minimal" set of generators
     649
     650LIB "general.lib"; // for sort
     651
     652/******************************************************/
    729653static proc zSort ( list @z )
    730 "
     654  "
    731655        Sort elements of a list of polynoms,
    732656        and normalize
    733657"
    734658{
    735         int n = size(@z);
    736         if ( n < 1 )
    737         {// empty list
    738                 return( @z );
    739         };
    740 
    741         if ( n == 1 )
    742         {
    743                 if ( @z[1] == 0 ) // if zero => empty list
    744                 {
    745                     return(list());
    746                 };         
    747 
    748                 @z[1] =  @z[1] * ( 1/leadcoef(@z[1]) ) ;
    749 
    750                 return( @z );
    751         };
    752 
    753         int i = 1;
    754        
    755         while ( i<=n )
    756         {
    757             if (size( @z[i] ) != 0)
    758             {
    759                 break;
    760             };
    761             i++;
    762         };
    763        
    764         if ( i > n )
    765         { /// all zeroes
    766             return(list());
    767         };
    768        
    769         ideal id = @z[i];
    770         i++;
    771        
    772         while ( i<= n )
    773         {
    774             if( @z[i] != 0 )
    775             {
    776                 id=@z[i],id;
    777             };
    778             i++;
    779         };
    780 
    781         // can sort only ideal generators:-(
    782         list srt = sort(id);
    783 
    784         poly p;
    785         // convert srt[1]:ideal to result:list.
    786         list result = list();
    787         for ( i = size(srt[1]); i>0; i-- )
    788         {
    789                 p = srt[1][i];
    790                 if ( p == 0 )
    791                 {
    792                     i --;
    793                     continue;
    794                 };
    795                 p = p* (1/leadcoef(p));         /// normalize
    796                 result = list(p) + result;
    797         };
    798 
    799 //      "OUT SORT::";
    800 //      result;
    801        
    802         return ( result );
    803 };
    804 
    805 
    806 
    807 //////////////////////////////////////////////////////////////////////////
     659  int n = size(@z);
     660  if ( n < 1 )
     661    {// empty list
     662      return( @z );
     663    };
     664
     665  if ( n == 1 )
     666    {
     667      if ( @z[1] == 0 ) // if zero => empty list
     668        {
     669          return(list());
     670        };         
     671
     672      @z[1] =  @z[1] * ( 1/leadcoef(@z[1]) ) ;
     673
     674      return( @z );
     675    };
     676
     677  int i = 1;
     678       
     679  while ( i<=n )
     680    {
     681      if (size( @z[i] ) != 0)
     682        {
     683          break;
     684        };
     685      i++;
     686    };
     687       
     688  if ( i > n )
     689    { // all zeroes
     690      return(list());
     691    };
     692       
     693  ideal id = @z[i];
     694  i++;
     695       
     696  while ( i<= n )
     697    {
     698      if( @z[i] != 0 )
     699        {
     700          id=@z[i],id;
     701        };
     702      i++;
     703    };
     704
     705  // can sort only ideal generators:-(
     706  list srt = sort(id);
     707
     708  poly p;
     709  // convert srt[1]:ideal to result:list.
     710  list result = list();
     711  for ( i = size(srt[1]); i>0; i-- )
     712    {
     713      p = srt[1][i];
     714      if ( p == 0 )
     715        {
     716          i --;
     717          continue;
     718        };
     719      p = p* (1/leadcoef(p));           // normalize
     720      result = list(p) + result;
     721    };
     722
     723  //    "OUT SORT::";
     724  //    result;
     725       
     726  return ( result );
     727};
     728
     729
     730
     731/******************************************************/
    808732static proc zRefine ( list @z )
    809 "
     733  "
    810734        kill terms by leading monomials...
    811735        Note: based on myCoeff 
    812736"
    813737{
    814         BCall("zRefine", "list");
    815 
    816         int @i, @j;
    817 
    818         poly @f  = 0; // leadmonom
    819 
    820         poly @ff = 0; // polyes
    821         poly @gg = 0;
    822         number @nf, @ng;
    823 
    824         int flag = 1;
    825 
    826         while ( flag == 1 )
    827         {
    828                 flag = 0;
    829 
    830                 @z = zSort(@z); // sort out, < ...
    831 
    832                 if( size(@z) < 2 )
     738  BCall("zRefine", "list");
     739
     740  int @i, @j;
     741
     742  poly @f  = 0; // leadmonom
     743
     744  poly @ff = 0; // polyes
     745  poly @gg = 0;
     746  number @nf, @ng;
     747
     748  int flag = 1;
     749
     750  while ( flag == 1 )
     751    {
     752      flag = 0;
     753
     754      @z = zSort(@z); // sort out, < ...
     755
     756      if( size(@z) < 2 )
     757        {
     758          return (@z);
     759        };
     760
     761      for ( @i = size(@z); @i > 0; @i -- ) // at 1st the biggest ... then smallest...
     762        {
     763
     764          @ff    = @z[@i];
     765
     766          if( size(@ff) == 0 )
     767            {
     768              @z = delete ( @z , @i );
     769              @i --;
     770              continue;
     771            };
     772
     773          @ff    = @ff*(1/leadcoef(@ff));
     774          @z[@i] = @ff;
     775          @f    = leadmonom(@ff);
     776
     777          for ( @j = (@i-1); (@j>0); @j --  )
     778            {
     779              @gg = @z[@j];
     780              @ng = myCoeff(@gg, @f); // leads?
     781              if( @ng!=0 )
    833782                {
    834                         return (@z);
    835                 };
    836 
    837                 for ( @i = size(@z); @i > 0; @i -- ) // at 1st the biggest ... then smallest...
     783                  @z[@j] = @gg - @ng * @ff;
     784                  flag = 1;
     785                };             
     786            };
     787          for ( @j = (@i+1); (@j<=size(@z)); @j ++ )
     788            {
     789              @gg = @z[@j];
     790              @ng = myCoeff(@gg, @f);
     791              if( @ng!=0 )
    838792                {
    839 
    840                         @ff    = @z[@i];
    841 
    842                         if( size(@ff) == 0 )
    843                         {
    844                                 @z = delete ( @z , @i );
    845                                 @i --;
    846                                 continue;
    847                         };
    848 
    849                         @ff    = @ff*(1/leadcoef(@ff));
    850                         @z[@i] = @ff;
    851                         @f      = leadmonom(@ff);
    852 
    853                         for ( @j = (@i-1); (@j>0); @j --  )
    854                         {
    855                                 @gg = @z[@j];
    856                                 @ng = myCoeff(@gg, @f); /// leads?
    857                                 if( @ng!=0 )
    858                                 {
    859                                         @z[@j] = @gg - @ng * @ff;
    860                                         flag = 1;
    861                                 };             
    862                         };
    863                         for ( @j = (@i+1); (@j<=size(@z)); @j ++ )
    864                         {
    865                                 @gg = @z[@j];
    866                                 @ng = myCoeff(@gg, @f);
    867                                 if( @ng!=0 )
    868                                 {
    869                                         @z[@j] = @gg - @ng * @ff;
    870                                         flag = 1;
    871                                 };             
    872                         };
    873 
    874                 };
    875         };
    876 
    877         ECall("zRefine", "list");
    878 
    879         return  ( @z );
    880 
    881 };
    882 
    883 
    884 ////////////////////////////////////////////////////////////////////////////////
    885 /// procedures for building "bad leadmonomials" set
    886 
    887 
    888 ////////////////////////////////////////////////////////////////////////////////
     793                  @z[@j] = @gg - @ng * @ff;
     794                  flag = 1;
     795                };             
     796            };
     797
     798        };
     799    };
     800
     801  ECall("zRefine", "list");
     802
     803  return        ( @z );
     804
     805};
     806
     807
     808/******************************************************/
     809// procedures for building "bad leadmonomials" set
     810
     811
     812/******************************************************/
    889813static proc checkPolyUniq( list l, poly p )
    890 "
     814  "
    891815        check wheather p sits already in l, assume l to be size-sorted <
    892816        return: -1 if present
     
    894818"
    895819{
    896         BCall( "checkPolyUniq", "{ " + string(l) + " }, " + string(p) );
    897         int n = size(l);
    898         int i = 1;
    899 
    900         int s = size(p);
    901        
    902         while( i<= n )
    903         {
    904             if ( size(l[i]) >= s )
    905             {
    906                 break;
    907             };
     820  BCall( "checkPolyUniq", "{ " + string(l) + " }, " + string(p) );
     821  int n = size(l);
     822  int i = 1;
     823
     824  int s = size(p);
     825       
     826  while( i<= n )
     827    {
     828      if ( size(l[i]) >= s )
     829        {
     830          break;
     831        };
    908832           
    909             i ++;
    910         };
    911 
    912         /// now: size(l[i]) >= s
    913         while( i<= n )
    914         {
    915             if ( size(l[i]) == s )
    916             {
    917                 break;
     833      i ++;
     834    };
     835
     836  // now: size(l[i]) >= s
     837  while( i<= n )
     838    {
     839      if ( size(l[i]) == s )
     840        {
     841          break;
    918842               
    919             };
    920             if ( l[i] == p )
    921             {
    922                 ECall( "checkPolyUniq", -1 );
    923                 return (-1);
    924             };
    925             i ++;
    926         };
    927 
    928         ECall( "checkPolyUniq", 1 );
    929         return ( 1 );
    930 
    931 };
    932 
    933 ////////////////////////////////////////////////////////////////////////////////
     843        };
     844      if ( l[i] == p )
     845        {
     846          ECall( "checkPolyUniq", -1 );
     847          return (-1);
     848        };
     849      i ++;
     850    };
     851
     852  ECall( "checkPolyUniq", 1 );
     853  return ( 1 );
     854
     855};
     856
     857/******************************************************/
    934858static proc addPolyUniq( list l, poly p )
    935 "
     859  "
    936860        add p into l uniquely, and keep l size-sorted <
    937861"
    938862{
    939         BCall( "addPolyUniq", " { " + string(l) + " }, " + string(p) );
    940        
    941         int n = size(l);
    942 
    943         if ( n == 0 )
    944         {
    945                 l = list(p);
     863  BCall( "addPolyUniq", " { " + string(l) + " }, " + string(p) );
     864       
     865  int n = size(l);
     866
     867  if ( n == 0 )
     868    {
     869      l = list(p);
    946870               
    947                 ECall( "addPolyUniq", l );
     871      ECall( "addPolyUniq", l );
    948872               
    949                 return (l);
    950         };
    951 
    952         int s = size(p);
    953        
    954         int i = 1;
    955         while( i<= n )
    956         {
    957                 if( size(l[i]) > s )
    958                 {
    959                         l = insert( l, p, i-1 ) ;
    960                         break;
    961                 };
     873      return (l);
     874    };
     875
     876  int s = size(p);
     877       
     878  int i = 1;
     879  while( i<= n )
     880    {
     881      if( size(l[i]) > s )
     882        {
     883          l = insert( l, p, i-1 ) ;
     884          break;
     885        };
    962886               
    963                 if( size(l[i]) == s )
    964                 {
    965                         if ( l[i] == p )
    966                         {
    967                                 break;
    968                         };
    969                 };
     887      if( size(l[i]) == s )
     888        {
     889          if ( l[i] == p )
     890            {
     891              break;
     892            };
     893        };
    970894               
    971                 i++;
    972         };
    973        
    974         if( i > n )
    975         {
    976             l = l + list(p);
    977         };
    978        
    979         ECall( "addPolyUniq", l );
    980         return(l);
    981 };
    982 
    983 
    984 ////////////////////////////////////////////////////////////////////////////////
     895      i++;
     896    };
     897       
     898  if( i > n )
     899    {
     900      l = l + list(p);
     901    };
     902       
     903  ECall( "addPolyUniq", l );
     904  return(l);
     905};
     906
     907
     908/******************************************************/
    985909static proc mergePolysUniq( list a, list b )
    986 "
     910  "
    987911        merge lists uniq
    988912"
    989913{
    990         BCall( "mergePolysUniq", "{ " + string(a) + " }, { " + string(b) + "} " );
    991        
    992         for( int i = 1; i <= size(b); i ++ )
    993         {
    994                 a = addPolyUniq(a, b[i]);
    995         };
    996        
    997         ECall( "mergePolysUniq", a );
    998        
    999         return (a);
    1000 };
    1001 
    1002 
    1003 ////////////////////////////////////////////////////////////////////////////////
     914  BCall( "mergePolysUniq", "{ " + string(a) + " }, { " + string(b) + "} " );
     915       
     916  for( int i = 1; i <= size(b); i ++ )
     917    {
     918      a = addPolyUniq(a, b[i]);
     919    };
     920       
     921  ECall( "mergePolysUniq", a );
     922       
     923  return (a);
     924};
     925
     926
     927/******************************************************/
    1004928static proc sortPolysUniq( list a )
    1005 "
     929  "
    1006930        sort list uniq
    1007931"
    1008932{
    1009         BCall( "sortPolysUniq", a );
    1010        
    1011         if( size(a) < 2 )
    1012         {
    1013                 ECall( "sortPolysUniq", a );
    1014                 return(a);
    1015         };
    1016        
    1017         list b = list(a[1]);
    1018        
    1019         for( int i = 2; i <= size(a); i ++ )
    1020         {
    1021                 b = addPolyUniq(b, a[i]);
    1022         };
    1023        
    1024         ECall( "sortPolysUniq", b );
    1025        
    1026         return (b);
    1027 };
    1028 
    1029 ////////////////////////////////////////////////////////////////////////////////
     933  BCall( "sortPolysUniq", a );
     934       
     935  if( size(a) < 2 )
     936    {
     937      ECall( "sortPolysUniq", a );
     938      return(a);
     939    };
     940       
     941  list b = list(a[1]);
     942       
     943  for( int i = 2; i <= size(a); i ++ )
     944    {
     945      b = addPolyUniq(b, a[i]);
     946    };
     947       
     948  ECall( "sortPolysUniq", b );
     949       
     950  return (b);
     951};
     952
     953/******************************************************/
    1030954static proc addRecordUniq ( list leadD, list newD, intvec texp, poly tm, int kind, list prs )
    1031 "
     955  "
    1032956        if kind = 0 => for PBW - no products
    1033957        if kind = 1 => with products
    1034958"
    1035959{
    1036         BCall ("addRecordUniq", "{old leads}, {new lead}, [" + string(texp) + "], " + string(tm) + ", " + string(kind) + ", {" + string(prs) + "}");
    1037 
    1038         int i;
    1039         int f_add = 0;
    1040 
    1041         prs = sortPolysUniq(prs);
    1042        
    1043         /// trick:
    1044         ///  check for presens of a monomial @tm in current index poly of @leads (=> in list @leads)
    1045         ///  if char = 0 then new_size > (if not present) or =  (if already present)
    1046         ///  if char > 0 then new_size > (if not present) or =< (if already present)
    1047         /// !!!!!
    1048         if( size(tm + leadD[2]) > size(leadD[2]) )
    1049         {
    1050                 f_add = 1;
     960  BCall ("addRecordUniq", "{old leads}, {new lead}, [" + string(texp) + "], " + string(tm) + ", " + string(kind) + ", {" + string(prs) + "}");
     961
     962  int i;
     963  int f_add = 0;
     964
     965  prs = sortPolysUniq(prs);
     966       
     967  // trick:
     968  //  check for presens of a monomial @tm in current index poly of @leads (=> in list @leads)
     969  //  if char = 0 then new_size > (if not present) or =  (if already present)
     970  //  if char > 0 then new_size > (if not present) or =< (if already present)
     971  // !!!!!
     972  if( size(tm + leadD[2]) > size(leadD[2]) )
     973    {
     974      f_add = 1;
     975    } else
     976      {
     977        if ( kind != 0 )
     978          {
     979            for ( i = 1; i<= size(leadD[1]); i++ )
     980              {
     981                if ( leadD[1][i][2] == tm )
     982                  {
     983                    for ( i = size(prs); i>0; i-- )
     984                      {
     985                        f_add = checkPolyUniq( leadD[1][i][3], prs[i] );
     986                        if( f_add == -1 )
     987                          {
     988                            prs = delete(prs, i);
     989                          };
     990                      };
     991
     992                    break;
     993                  };
     994              };                       
     995          };
     996      };
     997
     998  if ( f_add == 1 )
     999    {
     1000
     1001      list newlist ;
     1002      if ( kind != 0 )
     1003        {
     1004          newlist =  list ( list ( texp, tm, prs ) );
    10511005        } else
    1052         {
    1053                 if ( kind != 0 )
     1006          {
     1007            newlist =  list ( list ( texp, tm ) );
     1008          };
     1009
     1010
     1011      if ( size(newD[1]) == 0 )
     1012        {
     1013          newD[1] =  newlist;
     1014          newD[2] =  tm;
     1015        } else
     1016          {
     1017
     1018            if( size(tm + newD[2]) > size(newD[2]) )
     1019              {
     1020                newD[1] = newD[1] + newlist;
     1021                newD[2] = newD[2] + tm;
     1022              } else
    10541023                {
    1055                         for ( i = 1; i<= size(leadD[1]); i++ )
     1024                  if ( kind != 0 )
     1025                    {
     1026                      for ( i = 1; i<= size(newD[1]); i++ )
    10561027                        {
    1057                                 if ( leadD[1][i][2] == tm )
     1028                          if ( newD[1][i][2] == tm )
     1029                            {
     1030                              newD[1][i][3] = mergePolysUniq( newD[1][i][3], prs );
     1031                              break;
     1032                            };
     1033                        };
     1034                    };
     1035                };
     1036          };
     1037    };
     1038
     1039  ECall("addRecordUniq", "{new leads}");
     1040  return (newD);
     1041};
     1042
     1043
     1044/******************************************************/
     1045static proc mergeRecordsUniq ( list old, list new, int kind )
     1046{
     1047  BCall ("mergeRecordsUniq", "{old leads}, {new leads}, " + string(kind) );
     1048       
     1049  if(  size(new[1]) > 0 )
     1050    {
     1051      if ( size (old[1]) == 0 )
     1052        {
     1053          old = new;
     1054        } else
     1055          {
     1056            if ( kind != 0 )
     1057              {
     1058                int io;
     1059                for ( int in = 1; in <= size( new[1] ); in ++ )
     1060                  {
     1061                    if( size( new[1][in][2] + old[2] ) > size( old[2] ) )
     1062                      {
     1063                        old[1] = old[1] + list(new[1][in]);
     1064                        old[2] = old[2] + new[1][in][2];
     1065                      } else
     1066                        {
     1067                          for( io = 1; io <= size( old[1] ); io ++ )
     1068                            {
     1069                              if( old[1][io][2] == new[1][in][2] )
    10581070                                {
    1059                                         for ( i = size(prs); i>0; i-- )
    1060                                         {
    1061                                                 f_add = checkPolyUniq( leadD[1][i][3], prs[i] );
    1062                                                 if( f_add == -1 )
    1063                                                 {
    1064                                                         prs = delete(prs, i);
    1065                                                 };
    1066                                         };
    1067 
    1068                                         break;
     1071                                  old[1][io][3] = mergePolysUniq( old[1][io][3], new[1][in][3] );
     1072                                  break;
    10691073                                };
    1070                         };                     
     1074                            };                                         
     1075                        };
     1076                  };
     1077              } else
     1078                {
     1079                  old[1] = old[1] + new[1];
     1080                  old[2] = old[2] + new[2];
    10711081                };
    1072         };
    1073 
    1074         if ( f_add == 1 )
    1075         {
    1076 
    1077                 list newlist ;
    1078                 if ( kind != 0 )
    1079                 {
    1080                         newlist =  list ( list ( texp, tm, prs ) );
    1081                 } else
    1082                 {
    1083                         newlist =  list ( list ( texp, tm ) );
    1084                 };
    1085 
    1086 
    1087                 if ( size(newD[1]) == 0 )
    1088                 {
    1089                         newD[1] =  newlist;
    1090                         newD[2] =  tm;
    1091                 } else
    1092                 {
    1093 
    1094                         if( size(tm + newD[2]) > size(newD[2]) )
    1095                         {
    1096                                 newD[1] = newD[1] + newlist;
    1097                                 newD[2] = newD[2] + tm;
    1098                         } else
    1099                         {
    1100                                 if ( kind != 0 )
    1101                                 {
    1102                                         for ( i = 1; i<= size(newD[1]); i++ )
    1103                                         {
    1104                                                 if ( newD[1][i][2] == tm )
    1105                                                 {
    1106                                                         newD[1][i][3] = mergePolysUniq( newD[1][i][3], prs );
    1107                                                         break;
    1108                                                 };
    1109                                         };
    1110                                 };
    1111                         };
    1112                 };
    1113         };
    1114 
    1115         ECall("addRecordUniq", "{new leads}");
    1116         return (newD);
    1117 };
    1118 
    1119 
    1120 ////////////////////////////////////////////////////////////////////////////////
    1121 static proc mergeRecordsUniq ( list old, list new, int kind )
    1122 {
    1123         BCall ("mergeRecordsUniq", "{old leads}, {new leads}, " + string(kind) );
    1124        
    1125         if(  size(new[1]) > 0 )
    1126         {
    1127                 if ( size (old[1]) == 0 )
    1128                 {
    1129                         old = new;
    1130                 } else
    1131                 {
    1132                         if ( kind != 0 )
    1133                         {       
    1134                                 int io;
    1135                                 for ( int in = 1; in <= size( new[1] ); in ++ )
    1136                                 {
    1137                                         if( size( new[1][in][2] + old[2] ) > size( old[2] ) )
    1138                                         {
    1139                                                 old[1] = old[1] + list(new[1][in]);
    1140                                                 old[2] = old[2] + new[1][in][2];
    1141                                         } else
    1142                                         {
    1143                                                 for( io = 1; io <= size( old[1] ); io ++ )
    1144                                                 {
    1145                                                         if( old[1][io][2] == new[1][in][2] )
    1146                                                         {
    1147                                                                 old[1][io][3] = mergePolysUniq( old[1][io][3], new[1][in][3] );
    1148                                                                 break;
    1149                                                         };
    1150                                                 };                                             
    1151                                         };
    1152                                 };
    1153                         } else
    1154                         {
    1155                                 old[1] = old[1] + new[1];
    1156                                 old[2] = old[2] + new[2];
    1157                         };
    1158                 };
    1159         };
    1160 
    1161         ECall ("mergeRecordsUniq", "{new leads}");
    1162         return (old);
    1163 };
    1164 
    1165 
    1166 
    1167 //////////////////////////////////////////////////////////////////////////
     1082          };
     1083    };
     1084
     1085  ECall ("mergeRecordsUniq", "{new leads}");
     1086  return (old);
     1087};
     1088
     1089
     1090
     1091/******************************************************/
    11681092static proc init_bads(int @deg)
    11691093"
     
    11711095"
    11721096{
    1173         list @l = list();
    1174         for( int @i = 0; @i <= @deg; @i ++ )
    1175         {
    1176                 @l[index(@i)]   = list( list() , 0);
    1177                         /// new items:
    1178                         /// {
    1179                         ///     list of leads,    - empty list
    1180                         ///     sum poly "index", - poly(0)
    1181                         /// }
    1182         };
    1183         return (@l);
    1184 };
    1185 
    1186 //////////////////////////////////////////////////////////////////////////
     1097  list @l = list();
     1098  for( int @i = 0; @i <= @deg; @i ++ )
     1099    {
     1100      @l[index(@i)]     = list( list() , 0);
     1101      // new items:
     1102      // {
     1103      //        list of leads,    - empty list
     1104      //        sum poly "index", - poly(0)
     1105      // }
     1106    };
     1107  return (@l);
     1108};
     1109
     1110/******************************************************/
    11871111static proc zAddBad ( list @leads, list @newz, int @maxDeg, int kind )
    1188 "
     1112  "
    11891113        input:
    11901114                @leads: graded by deg list of
     
    12111135"
    12121136{
    1213         BCall ("zAddBad", "{old leads}, { " + string(@newz) + " }, " + string(@maxDeg) + ", " + string(kind) + "");
    1214 
    1215         int @newSize = size(@newz);
    1216         if ( @newSize < 1 )
    1217         {
    1218                 return (@leads);
    1219         };
    1220 
    1221         int @i, @j, @k, @dd, @deg;
    1222 
    1223 
    1224         poly    @m, @tm, @ttm;
    1225         intvec @exp, @texp, @ttexp;
    1226 
    1227         list @newzz;
    1228         int @size, @newdeg, @mydeg;
    1229 
    1230         poly @sum_old, @sum_new;
    1231        
    1232         poly a, b, @temp;
    1233 
    1234 /*
    1235         if kind = 0 => for PBW - no products
    1236         if kind = 1 => for zReduce - with products
    1237 */
    1238 
    1239         poly pr = 0; int pri; list prs;
    1240 
    1241         for ( @i = @newSize; @i > 0; @i -- )
    1242         {// for every new poly (@newz[@i]) do
    1243 
    1244                 @m   = leadmonom( @newz[@i] );
    1245                 @deg = deg(@m);
    1246                 @exp = leadexp(@m);
    1247 
    1248                 //// newzz void new list
    1249                 @newzz = init_bads(@maxDeg );
    1250 
    1251                 for( @mydeg = @deg; @mydeg <= @maxDeg;  @mydeg = @mydeg + @deg )
    1252                 {/// adding all possiblities for @newz[@i]^@j;
    1253 
    1254                         if ( @mydeg > @deg )
     1137  BCall ("zAddBad", "{old leads}, { " + string(@newz) + " }, " + string(@maxDeg) + ", " + string(kind) + "");
     1138
     1139  int @newSize = size(@newz);
     1140  if ( @newSize < 1 )
     1141    {
     1142      return (@leads);
     1143    };
     1144
     1145  int @i, @j, @k, @dd, @deg;
     1146
     1147
     1148  poly  @m, @tm, @ttm;
     1149  intvec @exp, @texp, @ttexp;
     1150
     1151  list @newzz;
     1152  int @size, @newdeg, @mydeg;
     1153
     1154  poly @sum_old, @sum_new;
     1155       
     1156  poly a, b, @temp;
     1157
     1158  /*
     1159    if kind = 0 => for PBW - no products
     1160    if kind = 1 => for zReduce - with products
     1161  */
     1162
     1163  poly pr = 0; int pri; list prs;
     1164
     1165  for ( @i = @newSize; @i > 0; @i -- )
     1166    {// for every new poly (@newz[@i]) do
     1167
     1168      @m   = leadmonom( @newz[@i] );
     1169      @deg = deg(@m);
     1170      @exp = leadexp(@m);
     1171
     1172      // newzz void new list
     1173      @newzz = init_bads(@maxDeg );
     1174
     1175      for( @mydeg = @deg; @mydeg <= @maxDeg;  @mydeg = @mydeg + @deg )
     1176        {// adding all possiblities for @newz[@i]^@j;
     1177
     1178          if ( @mydeg > @deg )
     1179            {
     1180              @texp     = @texp + @exp;
     1181              @tm       = LM ( @texp );
     1182              if ( kind != 0)
     1183                {
     1184                  pr = QNF_poly( pr * @newz[@i] ); // degrees must be there!!!
     1185                };
     1186            } else
     1187              {
     1188                @texp   = @exp;
     1189                @tm     = @m;
     1190                if ( kind != 0)
     1191                  {
     1192                    pr = @newz[@i];
     1193                  };
     1194              };
     1195
     1196          @temp =  QNF_poly(@tm);
     1197          if( @temp != @tm )
     1198            {
     1199              break;
     1200            };
     1201
     1202          /*!!*/                        @newzz[index(@mydeg)] =
     1203                                          /*!!*/        addRecordUniq( @leads[index(@mydeg)], @newzz[index(@mydeg)], @texp, @tm, kind, list(pr) );
     1204
     1205          for ( @dd = 1; (@dd <= @maxDeg) and ((@dd + @mydeg) <= @maxDeg ); @dd ++ )
     1206            { // for every good "deg"
     1207                               
     1208              @newdeg = @dd + @mydeg; // any deg should be additive!!!
     1209                               
     1210              for ( @k = size(@leads[index(@dd)][1]); @k > 0; @k -- )
     1211                {
     1212
     1213                  @ttexp        = (@leads[index(@dd)][1][@k][1]) + @texp;
     1214                  @ttm  = LM (@ttexp);
     1215                                       
     1216                  if ( kind != 0 )
     1217                    {
     1218                      prs = list();
     1219
     1220                      for( pri = 1; pri <= size(@leads[index(@dd)][1][@k][3]); pri++)
    12551221                        {
    1256                                 @texp   = @texp + @exp;
    1257                                 @tm     = LM ( @texp );
    1258                                 if ( kind != 0)
    1259                                 {
    1260                                         pr = QNF_poly( pr * @newz[@i] ); /// degrees must be there!!!
    1261                                 };
    1262                         } else
    1263                         {
    1264                                 @texp   = @exp;
    1265                                 @tm     = @m;
    1266                                 if ( kind != 0)
    1267                                 {
    1268                                         pr = @newz[@i];
    1269                                 };
     1222                          // to do products into list and add one list !!!
     1223                          a = QNF_poly( pr*@leads[index(@dd)][1][@k][3][pri]);
     1224                          b = QNF_poly( @leads[index(@dd)][1][@k][3][pri]*pr);
     1225
     1226                          prs= prs + list(a);
     1227
     1228                          if ( a!=b )
     1229                            {
     1230                              prs= prs + list(b);
     1231                            };
    12701232                        };
    12711233
    1272                         @temp =  QNF_poly(@tm);
    1273                         if( @temp != @tm )
    1274                         {
    1275                                 break;
    1276                         };
    1277 
    1278 /*!!*/                  @newzz[index(@mydeg)] =
    1279 /*!!*/  addRecordUniq( @leads[index(@mydeg)], @newzz[index(@mydeg)], @texp, @tm, kind, list(pr) );
    1280 
    1281                         for ( @dd = 1; (@dd <= @maxDeg) and ((@dd + @mydeg) <= @maxDeg ); @dd ++ )
    1282                         { //// for every good "deg"
    1283                                
    1284                                 @newdeg = @dd + @mydeg; /// any deg should be additive!!!
    1285                                
    1286                                 for ( @k = size(@leads[index(@dd)][1]); @k > 0; @k -- )
    1287                                 {
    1288 
    1289                                         @ttexp  = (@leads[index(@dd)][1][@k][1]) + @texp;
    1290                                         @ttm    = LM (@ttexp);
    1291                                        
    1292                                         if ( kind != 0 )
    1293                                         {
    1294                                                 prs = list();
    1295 
    1296                                                 for( pri = 1; pri <= size(@leads[index(@dd)][1][@k][3]); pri++)
    1297                                                 {
    1298                                                         /// to do products into list and add one list !!!
    1299                                                         a = QNF_poly( pr*@leads[index(@dd)][1][@k][3][pri]);
    1300                                                         b = QNF_poly( @leads[index(@dd)][1][@k][3][pri]*pr);
    1301 
    1302                                                         prs= prs + list(a);
    1303 
    1304                                                         if ( a!=b )
    1305                                                         {
    1306                                                             prs= prs + list(b);
    1307                                                         };
    1308                                                 };
    1309 
    1310                                         } else
    1311                                         {
    1312                                                 prs = list(pr);
    1313                                         };
    1314 
    1315                                         @temp =  QNF_poly(@ttm);
    1316                                         if( @temp == @ttm )
    1317                                         {
    1318 
    1319 /*!!*/                                          @newzz[index(@newdeg)] =
    1320 /*!!*/  addRecordUniq( @leads[index(@newdeg)], @newzz[index(@newdeg)], @ttexp, @ttm, kind, prs );
    1321                                         };
    1322 
    1323 
    1324                                 };                             
    1325                         }; /// for
    1326 
    1327                         if ( @deg == 0 )
    1328                         {
    1329                                 break;
    1330                         };
    1331                 };
    1332 
    1333                 for ( @mydeg = 0; @mydeg <= @maxDeg; @mydeg ++ )
    1334                 { /// adding currently generated to result
    1335                         @leads[index(@mydeg)] = mergeRecordsUniq ( @leads[index(@mydeg)], @newzz[index(@mydeg)], kind );
    1336                 };
    1337 
    1338         };
    1339 
    1340         ECall ("zAddBad", "{new leads}");
    1341 
    1342         return (@leads);
    1343 };
    1344 
    1345 
    1346 //////////////////////////////////////////////////////////////////////////
    1347 /// procedure for reducing a given poly wrt already calculated "badleadings"
    1348 //////////////////////////////////////////////////////////////////////////
    1349 
    1350 
    1351 //////////////////////////////////////////////////////////////////////////
     1234                    } else
     1235                      {
     1236                        prs = list(pr);
     1237                      };
     1238
     1239                  @temp =  QNF_poly(@ttm);
     1240                  if( @temp == @ttm )
     1241                    {
     1242
     1243                      /*!!*/                                            @newzz[index(@newdeg)] =
     1244                                                                          /*!!*/        addRecordUniq( @leads[index(@newdeg)], @newzz[index(@newdeg)], @ttexp, @ttm, kind, prs );
     1245                    };
     1246
     1247
     1248                };                             
     1249            }; // for
     1250
     1251          if ( @deg == 0 )
     1252            {
     1253              break;
     1254            };
     1255        };
     1256
     1257      for ( @mydeg = 0; @mydeg <= @maxDeg; @mydeg ++ )
     1258        { // adding currently generated to result
     1259          @leads[index(@mydeg)] = mergeRecordsUniq ( @leads[index(@mydeg)], @newzz[index(@mydeg)], kind );
     1260        };
     1261
     1262    };
     1263
     1264  ECall ("zAddBad", "{new leads}");
     1265
     1266  return (@leads);
     1267};
     1268
     1269
     1270/******************************************************/
     1271// procedure for reducing a given poly wrt already calculated "badleadings"
     1272
     1273/******************************************************/
    13521274static proc zReducePoly ( list leads, poly new, poly anfang )
    1353 "
     1275  "
    13541276        reduce poly new wrt found leads,
    13551277        return: list of all possible reductions...
    13561278"
    13571279{
    1358         poly temp = new;
    1359         poly rest = anfang;
    1360 
    1361         poly   lm;
    1362         number lc;
    1363 
    1364         list LEADS;
    1365 
    1366         int i, n;
    1367         list prs;
    1368 
    1369         while ( size(temp) > 0 )
    1370         {
    1371                 // do for every term... beginning with leading... 'till smallest
    1372                 lm = leadmonom( temp );
    1373 
    1374                 LEADS = leads[index(deg( lm ))]; // currently users bad leading products...
    1375 
    1376                 if( size(LEADS[2] + lm ) <=  size(LEADS[2]) )
    1377                 { // need to reduce, since leacmonom already in LEADS
     1280  poly temp = new;
     1281  poly rest = anfang;
     1282
     1283  poly   lm;
     1284  number lc;
     1285
     1286  list LEADS;
     1287
     1288  int i, n;
     1289  list prs;
     1290
     1291  while ( size(temp) > 0 )
     1292    {
     1293      // do for every term... beginning with leading... 'till smallest
     1294      lm = leadmonom( temp );
     1295
     1296      LEADS = leads[index(deg( lm ))]; // currently users bad leading products...
     1297
     1298      if( size(LEADS[2] + lm ) <=  size(LEADS[2]) )
     1299        { // need to reduce, since leacmonom already in LEADS
    13781300                   
    1379                         for ( i = 1; i <= size(LEADS[1]); i++ )
    1380                         {
    1381                                 if( LEADS[1][i][2] == lm )
     1301          for ( i = 1; i <= size(LEADS[1]); i++ )
     1302            {
     1303              if( LEADS[1][i][2] == lm )
     1304                {
     1305                               
     1306                  lc = leadcoef( temp ); // no need be the unit
     1307
     1308                  prs = LEADS[1][i][3]; // shouldbe generated by zAddBad with kind == 1
     1309                  n = size(prs) ;
     1310
     1311                  if ( n == 1 )
     1312                    { // no recursion
     1313
     1314                      temp = leadcoef(prs[1]) * temp - lc * prs[1]; // leadmonom goes down
     1315
     1316                    } else
     1317                      { // do recursion
     1318
     1319                        list result = list();
     1320                        poly newnew;
     1321                        int f_addrest = 0;
     1322                                               
     1323                        for( int pri = 1; pri <= n ; pri ++ )
     1324                          {     
     1325                            newnew = leadcoef(prs[pri]) * temp - lc * prs[pri]; // leadmonom goes down
     1326                                                       
     1327                            if( size( newnew ) > 0 )
     1328                              {
     1329                                result = result + zReducePoly(leads, newnew, rest);
     1330                              } else
    13821331                                {
    1383                                
    1384                                         lc = leadcoef( temp ); // no need be the unit
    1385 
    1386                                         prs = LEADS[1][i][3]; // shouldbe generated by zAddBad with kind == 1
    1387                                         n = size(prs) ;
    1388 
    1389                                         if ( n == 1 )
    1390                                         { // no recursion
    1391 
    1392                                                 temp = leadcoef(prs[1]) * temp - lc * prs[1]; // leadmonom goes down
    1393 
    1394                                         } else
    1395                                         { // do recursion
    1396 
    1397                                                 list result = list();
    1398                                                 poly newnew;
    1399                                                 int f_addrest = 0;
    1400                                                
    1401                                                 for( int pri = 1; pri <= n ; pri ++ )
    1402                                                 {       
    1403                                                         newnew = leadcoef(prs[pri]) * temp - lc * prs[pri]; // leadmonom goes down
    1404                                                        
    1405                                                         if( size( newnew ) > 0 )
    1406                                                         {
    1407                                                             result = result + zReducePoly(leads, newnew, rest);
    1408                                                         } else
    1409                                                         {
    1410                                                             f_addrest = 1;
    1411                                                         };
    1412                                                 };
    1413 
    1414                                                 if ( f_addrest == 1 )
    1415                                                 {
    1416                                                     result = result + list(rest);
    1417                                                 };
    1418                                                 return ( result );
    1419 
    1420                                         };
    1421                                         break;
     1332                                  f_addrest = 1;
    14221333                                };
    1423                         };
    1424 
    1425                 } else
    1426                 { // no such leadmonom in leads
    1427        
    1428                         rest = rest + lead ( temp );
    1429                         temp = temp - lead ( temp ); // leadcoeff goes down
     1334                          };
     1335
     1336                        if ( f_addrest == 1 )
     1337                          {
     1338                            result = result + list(rest);
     1339                          };
     1340                        return ( result );
     1341
     1342                      };
     1343                  break;
    14301344                };
    1431         };
    1432 
    1433         return (list(rest));
    1434 
    1435 };
    1436 
    1437 
    1438 //////////////////////////////////////////////////////////////////////////
     1345            };
     1346
     1347        } else
     1348          { // no such leadmonom in leads
     1349       
     1350            rest = rest + lead ( temp );
     1351            temp = temp - lead ( temp ); // leadcoeff goes down
     1352          };
     1353    };
     1354
     1355  return (list(rest));
     1356
     1357};
     1358
     1359
     1360/******************************************************/
    14391361static proc zCancel ( list @tt, list @leads )
    1440 "
     1362  "
    14411363        just kill entries of plane PBW base with leading monomials in @leads...
    14421364"
    14431365{
    1444         int @i;
    1445 
    1446         if ( (size(@tt) == 0) )
    1447         {
    1448                 return (@tt);
    1449         };
    1450 
    1451         list result = list();
    1452         poly g, f;
    1453 
    1454         for ( @i = size(@tt); @i > 0 ; @i -- )
    1455         /// for all PBW entries:
    1456         {
    1457                 g = leadmonom(@tt[@i][1]);    /// pbw entry
    1458                 f = @leads[index(deg(g))][2]; /// 'index' poly from @leads
    1459 
    1460                 if ( size(f + g) > size(f) ) /// if g not in @leads (works only for monomials)
    1461                 {
    1462                         result = list( @tt[@i] ) + result;
    1463                 };
    1464         };
    1465 
    1466         return (result);
    1467 
    1468 };
    1469 
    1470 //////////////////////////////////////////////////////////////////////////
    1471 /// procedure for computing "minimal" set of generators...
    1472 //////////////////////////////////////////////////////////////////////////
    1473 
    1474 //////////////////////////////////////////////////////////////////////////
     1366  int @i;
     1367
     1368  if ( (size(@tt) == 0) )
     1369    {
     1370      return (@tt);
     1371    };
     1372
     1373  list result = list();
     1374  poly g, f;
     1375
     1376  for ( @i = size(@tt); @i > 0 ; @i -- )
     1377    // for all PBW entries:
     1378    {
     1379      g = leadmonom(@tt[@i][1]);    // pbw entry
     1380      f = @leads[index(deg(g))][2]; // 'index' poly from @leads
     1381
     1382      if ( size(f + g) > size(f) ) // if g not in @leads (works only for monomials)
     1383        {
     1384          result = list( @tt[@i] ) + result;
     1385        };
     1386    };
     1387
     1388  return (result);
     1389
     1390};
     1391
     1392/******************************************************/
     1393// procedure for computing "minimal" set of generators...
     1394
     1395/******************************************************/
    14751396static proc zReduce( list @z )
    1476 "
     1397  "
    14771398        reduce a set @z - base of Center as V.S.
    14781399        into a minimal set!!
    14791400"
    14801401{
    1481         BCall( "zReduce", @z );
    1482        
    1483         @z = zRefine ( @z );
    1484         int n = size( @z );
    1485 
    1486         if ( n < 2 )
    1487         {
    1488                 return (@z);
    1489         };
    1490 
    1491         int d = maxDeg( stdDeg, @z );
    1492 
    1493         list leads = init_bads( d );
    1494 
    1495         list @red;
    1496 
    1497         int f_add; int j;
    1498 
    1499         list result = list();
    1500 
    1501         poly p;
    1502        
    1503         for ( int i = 1; i <= n ; i++ )
    1504         {// in this order... from least to maximal...
    1505 
    1506                 p = @z[i];
    1507 
    1508                 @red = zReducePoly ( leads, p, 0 );
    1509 
    1510                 f_add = 1;
     1402  BCall( "zReduce", @z );
     1403       
     1404  @z = zRefine ( @z );
     1405  int n = size( @z );
     1406
     1407  if ( n < 2 )
     1408    {
     1409      return (@z);
     1410    };
     1411
     1412  int d = maxDeg( @z );
     1413
     1414  if( d== -1 )
     1415    {
     1416      return (@z);
     1417    };
     1418
     1419  list leads = init_bads( d );
     1420
     1421  list @red;
     1422
     1423  int f_add; int j;
     1424
     1425  list result = list();
     1426
     1427  poly p;
     1428       
     1429  for ( int i = 1; i <= n ; i++ )
     1430    {// in this order... from least to maximal...
     1431
     1432      p = @z[i];
     1433
     1434      @red = zReducePoly ( leads, p, 0 );
     1435
     1436      f_add = 1;
    15111437               
    1512                 for( j = 1; j <= size(@red); j++ )
    1513                 {
    1514                         if ( @red[j] == 0 )
    1515                         {
    1516                                 f_add = 0;
    1517                                 break; // nothing new....
    1518                         };
    1519                 };
    1520 
    1521                 @red = zRefine( @red ); /// there will be no zeroes after ???
    1522 
    1523                 if( size(@red) > 0 )
    1524                 {                   
    1525                     leads = zAddBad( leads, @red, d, 1);
    1526                 };
    1527 
    1528                 if ( f_add == 1 )
    1529                 {// we got something new....
    1530                         result = result + list(@red); // ??? which one to add? => trying to add all
    1531                 };
    1532         };
    1533        
    1534         ECall( "zReduce", result );
    1535        
    1536         return (result);
    1537 };
    1538 
    1539 
    1540 ////////////////////////////////////////////////////////////////////////////
    1541 ////////////////////////////////////////////////////////////////////////////
     1438      for( j = 1; j <= size(@red); j++ )
     1439        {
     1440          if ( @red[j] == 0 )
     1441            {
     1442              f_add = 0;
     1443              break; // nothing new....
     1444            };
     1445        };
     1446     
     1447      @red = zRefine( @red ); // there will be no zeroes after ???
     1448
     1449      if( size(@red) > 0 )
     1450        {                   
     1451          leads = zAddBad( leads, @red, d, 1);
     1452        };
     1453     
     1454      if ( f_add == 1 )
     1455        {
     1456          // we got something new....
     1457          result = result + list(@red); // ??? which one to add? => trying to add all
     1458        };
     1459    };
     1460       
     1461  ECall( "zReduce", result );
     1462       
     1463  return (result);
     1464};
     1465
     1466
     1467/******************************************************/
    15421468// PBW procedures ....
    15431469
     
    15471473*/
    15481474
    1549 
    1550 /////////////////////////////////////////////////////////////////////////////
     1475/******************************************************/
    15511476static proc PBW_init()
    1552 "
     1477  "
    15531478PBW[ index(0) ] = PBW part of degree 0:
    15541479record:
     
    15581483"
    15591484{
    1560         BCall( "PBW_init" );
    1561         return (list(list(1, 0, 1)));
    1562 };
    1563 
    1564 
    1565 /////////////////////////////////////////////////////////////////////////////
     1485  BCall( "PBW_init" );
     1486  return (list(list(1, 0, 1)));
     1487};
     1488
     1489
     1490/******************************************************/
    15661491static proc PBW_sys(list VARS, poly p)
    1567 "
     1492  "
    15681493        calculate the array [1..NVARS] of records:
    15691494record[i] =
     
    15761501{
    15771502
    1578         BCall( "PBW_sys" );
    1579 
    1580         poly t;
    1581         for (int v = size(VARS); v>0; v -- )
    1582         {
    1583                 t       = VARS[v];
    1584                 VARS[v] = list( t, QNF_poly( p*t - t*p ), deg(t) ) ;
    1585         };
     1503  BCall( "PBW_sys" );
     1504
     1505  poly t;
     1506  for (int v = size(VARS); v>0; v -- )
     1507    {
     1508      t         = VARS[v];
     1509      VARS[v] = list( t, QNF_poly( p*t - t*p ), deg(t) ) ;
     1510    };
    15861511               
    1587         return (VARS);
    1588 };
    1589 
    1590 /////////////////////////////////////////////////////////////////////////////
     1512  return (VARS);
     1513};
     1514
     1515/******************************************************/
    15911516static proc PBW_done( list PBW )
    1592 "
     1517  "
    15931518        collects all together, from graded lists into plane list.       
    15941519       
     
    15961521"
    15971522{
    1598         BCall( "PBW_done" );
    1599 
    1600         list result = list();
    1601 
    1602         int n = size(PBW);
    1603         for ( int i = 1; i <= n; i++)
    1604         {
    1605                 result = result + PBW[i];
    1606         };
    1607 
    1608         return (result);
    1609 };
    1610 
    1611 /////////////////////////////////////////////////////////////////////////////
     1523  BCall( "PBW_done" );
     1524
     1525  list result = list();
     1526
     1527  int n = size(PBW);
     1528  for ( int i = 1; i <= n; i++)
     1529    {
     1530      result = result + PBW[i];
     1531    };
     1532
     1533  return (result);
     1534};
     1535
     1536/******************************************************/
    16121537static proc PBW_next ( list PBW, int k, list sys)
    16131538{
    1614         BCall( "PBW_next", k );
    1615 
    1616         list temp;
    1617         int kk, nn, ii;
    1618         list result = list(); // PBW[ index(k) ] ought to be empty ??
    1619         int N = size(sys);
    1620        
    1621         for ( int i = 1; i <= N; i ++ ) // for all vars
    1622         {
    1623                 kk = (k - sys[i][3]);   // any deg should be additive function wrt multiplication
    1624                 if ( kk >= 0 )
     1539  BCall( "PBW_next", k );
     1540
     1541  list temp;
     1542  int kk, nn, ii;
     1543  list result = list(); // PBW[ index(k) ] ought to be empty ??
     1544  int N = size(sys);
     1545       
     1546  for ( int i = 1; i <= N; i ++ ) // for all vars
     1547    {
     1548      kk = (k - sys[i][3]);   // any deg should be additive function wrt multiplication
     1549      if ( kk >= 0 )
     1550        {
     1551          nn = size( PBW[ index(kk) ] );
     1552          for ( ii = 1; ii <= nn; ii ++ )
     1553            {
     1554              temp = PBW[ index(kk) ][ii];
     1555
     1556              if ( temp[3] <= i )
    16251557                {
    1626                         nn = size( PBW[ index(kk) ] );
    1627                         for ( ii = 1; ii <= nn; ii ++ )
    1628                         {
    1629                                 temp = PBW[ index(kk) ][ii];
    1630 
    1631                                 if ( temp[3] <= i )
    1632                                 {
    1633                                         temp[2] = temp[2]*sys[i][1] + temp[1]*sys[i][2]; // recursive [,]
    1634                                         temp[1] = temp[1]*sys[i][1];
    1635                                         temp[3] = i;
     1558                  temp[2] = temp[2]*sys[i][1] + temp[1]*sys[i][2]; // recursive [,]
     1559                  temp[1] = temp[1]*sys[i][1];
     1560                  temp[3] = i;
    16361561                                       
    1637                                         result = result + list ( temp );
    1638                                 };
    1639                         };
     1562                  result = result + list ( temp );
    16401563                };
    1641         };
    1642 
    1643         result = QNF_list( result, 2 + 4);
    1644         ECall( "PBW_next", "list result" );
    1645         return (result);
    1646        
    1647 };
    1648 
    1649 /////////////////////////////////////////////////////////////////////////////
     1564            };
     1565        };
     1566    };
     1567
     1568  result = QNF_list( result, 2 + 4);
     1569  ECall( "PBW_next", "list result" );
     1570  return (result);
     1571       
     1572};
     1573
     1574/******************************************************/
    16501575static proc PBW_base( int MaxDeg, poly p )
    16511576{
    1652         BCall( "PBW_base", MaxDeg, p );
    1653 
    1654         list PBW = list();
    1655         list SYS = PBW_sys( my_vars(), p );
    1656 
    1657         PBW[index(0)] = PBW_init();
    1658        
    1659         for (int k = 1; k <= MaxDeg; k ++ )
    1660         {               
    1661                 PBW[index(k)] = PBW_next( PBW, k, SYS );
    1662         };
    1663 
    1664         return (PBW_done( PBW ));
     1577  BCall( "PBW_base", MaxDeg, p );
     1578
     1579  list PBW = list();
     1580  list SYS = PBW_sys( my_vars(), p );
     1581
     1582  PBW[index(0)] = PBW_init();
     1583       
     1584  for (int k = 1; k <= MaxDeg; k ++ )
     1585    {           
     1586      PBW[index(k)] = PBW_next( PBW, k, SYS );
     1587    };
     1588
     1589  return (PBW_done( PBW ));
    16651590};
    16661591
    16671592
    1668 /////////////////////////////////////////////////////////////////////////////
    1669 /////////////////////////////////////////////////////////////////////////////
    1670 // standard center procedures... "zCommon"
    1671 
    1672 
    1673 /******************************************************************************/
    1674 static proc my_calc (list @given, poly @BASE_POLY, list #)
     1593/******************************************************/
     1594// standard center procedures...
     1595
     1596/******************************************************/
     1597static proc FSS (list @given, poly @BASE_POLY)
    16751598"
    16761599    Gauss with computation of kernel v.s basis
    1677 "
    1678 {
    1679     list @nums = list ();
    1680     intvec @ones;
    1681     int @j, @k,
    1682         @n, @v,
    1683         @a, @nn;
     1600    to be optimizes
     1601"
     1602{
     1603  BCall( "FSS", "list", "basepoly: ", @BASE_POLY);
     1604
     1605  list @nums = list ();
     1606  intvec @ones;
     1607  int @j, @k,
     1608    @n, @v,
     1609    @a, @nn;
    16841610   
    1685     @n = size( @BASE_POLY );
    1686     @v = size( @given );
     1611  @n = size( @BASE_POLY );
     1612  @v = size( @given );
    16871613   
    1688     if ( @v == 0 )
     1614  if ( (@v == 0) or (@n == 0) )
    16891615    {   
    1690         return ( list() );
    1691     };
    1692 
    1693     if ( size(#)>0 )
    1694     {
    1695         poly @Next_Ad_var;
    1696         if ( typeof(#[1]) == "int")
    1697         {
    1698             int @z_next = #[1];
    1699             @Next_Ad_var = var (@z_next);
    1700         };
    1701        
    1702         if ( typeof(#[1]) == "poly")
    1703         {
    1704             @Next_Ad_var = #[1];
    1705         };
    1706     };
     1616      return (@given);
     1617    };
     1618
     1619  matrix MD[@n][@v];
     1620
     1621  number @max;
     1622  poly @test;
     1623  poly @t;
     1624
     1625  list LM = list(); // rec: { exp, poly }
     1626
     1627  for( @k = @v; @k > 0; @k -- )
     1628    {
     1629      LM[@k] = list();
     1630      @t =  @given[@k][2];
     1631      LM[@k][2]= @t;
     1632      LM[@k][1]= leadexp( @t );
     1633    };
     1634
     1635
     1636  intvec @w;
     1637  for ( @a = 1 ; @a <= @n ; @a ++ )
     1638    {
     1639      @w = leadexp(@BASE_POLY[@a]);
     1640      for( @k = @v; @k > 0; @k -- )
     1641        {
     1642          if( @w == LM[@k][1] )
     1643            {
     1644              @t = LM[@k][2];
     1645              MD[@a, @k] = leadcoef( @t );
     1646              @t = @t - lead ( @t );
     1647              LM[@k][2]  = @t;
     1648
     1649              LM[@k][1]  = leadexp(@t);
     1650            };
     1651        };
     1652
     1653    };
     1654
     1655  int @x, @y;
     1656  number @min;
     1657
     1658  number @div;
     1659  //    @nums = list();
    17071660   
    1708     poly @t;
     1661  // Gauss
     1662//  print("Before Gauss, Matrix: ");
     1663//  print( MD );
     1664
    17091665   
    1710     if ( @n == 0 )
    1711     {   
    1712         if (defined(@Next_Ad_var))
    1713         {
    1714             for ( @k = @v; @k > 0; @k --)
    1715             {
    1716                 @t = @given[@k][1];
    1717                 @given[@k][2] = @Next_Ad_var * @t - @t * @Next_Ad_var;
    1718             };
    1719         };
    1720         return (@given);
    1721     };
    1722    
    1723 
    1724     matrix MD[@n][@v];
    1725 
    1726     number @max;
    1727     poly @test;
    1728 
    1729     list LM = list(); /// rec: { exp, poly }
    1730 
    1731 
    1732 
    1733 
    1734     for( @k = @v; @k > 0; @k -- )
    1735     {
    1736           LM[@k] = list();
    1737           @t =  @given[@k][2];
    1738           LM[@k][2]= @t;
    1739           LM[@k][1]= leadexp( @t );
    1740     };
    1741 
    1742 
    1743     intvec @w;
    1744     for ( @a = 1 ; @a <= @n ; @a ++ )
    1745     {
    1746           @w = leadexp(@BASE_POLY[@a]);
    1747           for( @k = @v; @k > 0; @k -- )
    1748         {
    1749                   if( @w == LM[@k][1] )
    1750                   {
    1751                     @t = LM[@k][2];
    1752                         MD[@a, @k] = leadcoef( @t );
    1753                         @t = @t - lead ( @t );
    1754                         LM[@k][2]  = @t;
    1755 
    1756                         LM[@k][1]  = leadexp(@t);
    1757                   };
    1758         };
    1759 
    1760     };
    1761 
    1762     int @x, @y;
    1763     number @min;
    1764 
    1765     number @div;
    1766 ///    @nums = list();
    1767 
    1768     /// Gauss
    1769 
    1770     @x = 1;
    1771     @y = 1;
    1772     while ( (@y <= @n) && (@x <= @v))
     1666  @x = 1;
     1667  @y = 1;
     1668  while ( (@y <= @n) && (@x <= @v))
    17731669    // main Gauss loop...
    17741670    {
    1775         @min =  leadcoef( MD[@y, @x] ); /// curr diag.
    1776 
    1777         if ( @min == 0 )
     1671      @min =  leadcoef( MD[@y, @x] ); // curr diag.
     1672       
     1673      if ( @min == 0 )
    17781674        // if zero on diag...
    17791675        {
    1780             @j = 0;
    1781 
    1782             /// let's find the minimal
    1783             for ( @k = @y+1; @k <= @n; @k ++ )
    1784             {
    1785                 @max = leadcoef( MD[ @k, @x ] );
    1786                 if ( @max != 0 )
     1676          @j = 0;
     1677               
     1678          // let's find the minimal
     1679          for ( @k = @y+1; @k <= @n; @k ++ )
     1680            {
     1681              @max = leadcoef( MD[ @k, @x ] );
     1682              if ( @max != 0 )
    17871683                {
    1788                     @j = @k;                   
    1789                     @min = @max;
    1790 //////              @k ++;
    1791                     break; //// this pure for
    1792                     /// continue; // for find min
     1684                  @j = @k;                     
     1685                  @min = @max;
     1686                  //                @k ++;
     1687                  break; // this pure for
     1688                  // continue;
     1689                  // for find min
    17931690                };
    1794             }; /// for (@k) // found minimal
    1795        
    1796             if ( @j == 0 )
     1691            }; // for (@k) // found minimal
     1692               
     1693          if ( @j == 0 )
    17971694            {       
    1798                 /// das ist gut! ///           
    1799                 @nums = @nums + list ( list ( @x, @y ) );
    1800                 @x ++;
    1801                 continue; // while
     1695              // das ist gut! //               
     1696              @nums = @nums + list ( list ( @x, @y ) );
     1697              @x ++;
     1698              continue; // while
    18021699            } ;
    18031700           
    1804             for ( @k = @x ; @k <= @v ; @k ++ )
    1805             {
    1806                 @t =  MD[ @j, @k ];
    1807                 MD[ @j, @k ] = MD[ @y, @k ];
    1808                 MD[ @y, @k ] = @t;
    1809             };
     1701          for ( @k = @x ; @k <= @v ; @k ++ )
     1702            {
     1703              @t =  MD[ @j, @k ];
     1704              MD[ @j, @k ] = MD[ @y, @k ];
     1705              MD[ @y, @k ] = @t;
     1706            };
     1707               
     1708        }; // if diag === zero.
     1709      @ones[@y] = @x;
    18101710           
    1811         }; /// if diag === zero.
    1812         @ones[@y] = @x;
    1813        
    1814         if ( @min == 0 )
    1815         {
    1816 //          write_latex_str ( " ******************* ERROR ******************* " );
    1817             quit;
    1818         };
    1819 
     1711      if ( @min == 0 )
     1712        {
     1713          //        write_latex_str ( " ******************* ERROR ******************* " );
     1714          quit;
     1715        };
     1716           
     1717           
     1718      if ( @min != 1) // let's norm the row... to make the '1' !
     1719        {
     1720          @min = 1 / @min;
     1721          for ( @k = @x; @k <= @v; @k++ )
     1722            {
     1723              @max = leadcoef( MD[@y, @k] );
     1724                   
     1725              if ( @max == 0)
     1726                {
     1727                  @k ++ ;
     1728                  continue; // for : norming the row...
     1729                };
     1730                   
     1731              MD[@y, @k] = @max * @min;  // here must be Field ...
     1732            };
     1733               
     1734          //        @min = 1;
     1735               
     1736        };
     1737           
     1738      // here must be @min != 0;
     1739      for ( @k = 1; @k <= @n; @k++ )
     1740        {
     1741          @max = leadcoef( MD[@k, @x] );
     1742               
     1743          if ( ( @k == @y) || (@max == 0) )
     1744            {
     1745              @k ++ ;
     1746              continue;
     1747            };
     1748               
     1749          for ( @j = @x; @j <= @v ; @j ++ )
     1750            {
     1751              MD[ @k, @j ] = MD[ @k, @j ] - @max * MD[ @y, @j ];
     1752            };             
     1753         
     1754        }; //killing
     1755           
     1756      @x ++;
     1757      @y ++;
     1758    }; //main while.
     1759  /******************************************************/
    18201760   
    1821         if ( @min != 1) /// let's norm the row... to make the '1' !
    1822         {
    1823             @min = 1 / @min;
    1824             for ( @k = @x; @k <= @v; @k++ )
    1825             {
    1826                 @max = leadcoef( MD[@y, @k] );
    1827 
    1828                 if ( @max == 0)
    1829                 {
    1830                     @k ++ ;
    1831                     continue; /// for : norming the row...
    1832                 };
    1833            
    1834                 MD[@y, @k] = @max * @min;  ////here must be Field ...
    1835             };
    1836        
    1837 ////        @min = 1;
    1838 
    1839         };
    1840        
    1841         //// here must be @min != 0;
    1842         for ( @k = 1; @k <= @n; @k++ )
    1843         {
    1844             @max = leadcoef( MD[@k, @x] );
    1845 
    1846             if ( ( @k == @y) || (@max == 0) )
    1847             {
    1848                 @k ++ ;
    1849                 continue;
    1850             };
    1851 
    1852             for ( @j = @x; @j <= @v ; @j ++ )
    1853             {
    1854                 MD[ @k, @j ] = MD[ @k, @j ] - @max * MD[ @y, @j ];
    1855             };             
    1856            
    1857         }; ///killing
    1858        
    1859         @x ++;
    1860         @y ++;
    1861     }; //////main while.
    1862 //////////////////////////////////////////////////////
    1863 
    1864 
    1865     /// computation of kernel's basis
    1866 
    1867    
    1868     if ( @x <= @v )
    1869     {
    1870         for ( @k = @x; @k <= @v ; @k ++ )
    1871         {
    1872             @nums = @nums + list ( list ( @k, @n+1 ) );
     1761//  print("Gaussed Matrix: ");
     1762//  print( MD );
     1763
     1764  // computation of kernel's basis
     1765
     1766  if ( @x <= @v )
     1767    {
     1768      for ( @k = @x; @k <= @v ; @k ++ )
     1769        {
     1770          @nums = @nums + list ( list ( @k, @n+1 ) );   
    18731771        }
    18741772    }
    18751773   
    1876     if ( @y <= @n )
    1877     {
    1878         @n = @y - 1;
    1879     };
     1774//  print("Nums: " );
     1775//  print (@nums);
    18801776   
    1881 
    1882     list result = list();
     1777  list result = list();
    18831778   
    1884 /// real calculations of the Base of a Ker as V.S.
    1885 
    1886     for ( @k = 1; @k <= size(@nums) ; @k ++ )
    1887     {
    1888         @x = @nums[@k][1];
    1889         @j  = @nums[@k][2];
    1890        
    1891         @t = @given[@x][1];
    1892        
    1893         for ( @y = 1; @y < @j ; @y ++ )
    1894         /// for every "@x" column
    1895         {
    1896             @max = leadcoef( MD[@y, @x] );
    1897             if ( (@max != 0) )
    1898             {
    1899                 @a = @ones[@y];
    1900                 @t = @t - @max * @given[@a][1];
    1901             };
    1902         };
    1903 
    1904         if ( defined(@Next_Ad_var) )
    1905         {
    1906             result[@k] = list ( @t, @Next_Ad_var * @t - @t * @Next_Ad_var ) ;
     1779  // real calculations of the Base of a Ker as V.S.
     1780
     1781  for ( @k = 1; @k <= size(@nums) ; @k ++ )
     1782    {
     1783      @x = @nums[@k][1];
     1784      @j  = @nums[@k][2];
     1785           
     1786      @t = @given[@x][1];
     1787           
     1788      for ( @y = 1; @y < @j ; @y ++ )
     1789        // for every "@x" column
     1790        {
     1791          @max = leadcoef( MD[@y, @x] );
     1792          if ( (@max != 0) )
     1793            {
     1794              @a = @ones[@y];
     1795              @t = @t - @max * @given[@a][1];
     1796            };
     1797        };
     1798     
     1799      result[@k] = @t;
     1800    };
     1801
     1802  ECall( "FSS", "list" );
     1803  return ( result );
     1804};
     1805
     1806/******************************************************/
     1807static proc reduce_one_per_row(list @given)
     1808{
     1809  BCall( "reduce_one_per_row" );
     1810       
     1811  int @k;
     1812  int @l = size (@given);
     1813       
     1814  if( @l == 0 )
     1815    {
     1816      return (@given);
     1817    };
     1818       
     1819  int @char = char(basering);
     1820       
     1821  intvec @marks;
     1822
     1823  if ( @char == 0 )
     1824    {
     1825      poly @t = poly(0);
     1826    };
     1827
     1828  list @unis = list ();
     1829  poly p;
     1830       
     1831  for ( @k = @l; @k > 0; @k -- )
     1832    {
     1833      p = uni_poly( @given[@k][2] );
     1834      @unis[@k] = p;
     1835
     1836      if( p == 0 )
     1837        {
     1838          @marks[@k] = 2;
    19071839        } else
    1908         {
    1909             result[@k] = @t;
    1910         }
    1911     };
    1912 
    1913     return ( result );
    1914 };
    1915 
    1916 
    1917 /******************************************************************************/
    1918 /*
    1919 static proc add2list( poly l, poly q )
    1920 "
    1921         l is a sorted list of monomials (polynom with coefs: 0 or 1)
    1922         adding monomials from q, checking also for presence
    1923        
    1924         ?? NEED DEBUG ??
    1925 "
    1926 {
    1927         BCall( "add2list", "{" + string(l) + "} , " + string(q) );
    1928 
    1929         poly t;
    1930         for( int i = size(q); i>0; i-- )
    1931         {
    1932 
    1933         /// trick:
    1934         ///  check for presens of a monomial q[i] in l
    1935         ///  if char = 0 then new_size > (if not present) or =  (if already present)
    1936         ///  if char > 0 then new_size > (if not present) or =< (if already present)
    1937                 t = l + leadmonom(q[i]);
    1938                 if ( size(t) > size(l) )
     1840          {
     1841            @marks[@k] = 1;
     1842            if ( @char == 0 )
     1843              {
     1844                @t = @t + p;
     1845              };
     1846          };
     1847    };
     1848       
     1849  if ( @char != 0 )
     1850    {
     1851      def save = basering;
     1852      execute( "ring NewRingWithGoodField = (0), (" + varstr(save) + "), (" + ordstr(save) + "); ");
     1853      poly @t = 0;
     1854               
     1855      if(! defined (@unis) )
     1856        {
     1857          list @unis = imap( save, @unis );
     1858        };
     1859               
     1860      for ( @k = @l; @k > 0; @k -- )
     1861        {
     1862          if( @marks[@k] == 1 )
     1863            {
     1864              @t = @t + @unis[@k];
     1865            };
     1866        };
     1867    };
     1868               
     1869  int @loop_size = size(@t);
     1870  poly @for_delete, @tt;
     1871  int @ll;
     1872       
     1873  while( @loop_size > 0 )
     1874    {
     1875      @for_delete = poly(0);
     1876      @ll = size(@t);   
     1877               
     1878      for ( @k = @ll; @k > 0; @k -- )
     1879        {
     1880          if ( leadcoef(@t[@k]) == 1 )
     1881            {
     1882              @for_delete = @for_delete + @t[@k];
     1883            };
     1884        };
     1885
     1886      @loop_size = size( @for_delete );
     1887
     1888      if ( @loop_size>0 )
     1889        {
     1890          for( @k = @l ; @k > 0 ; @k -- )
     1891            {       
     1892              if ( @marks[@k] == 1)
    19391893                {
    1940                         l = t;
    1941                 };     
    1942         };
    1943        
    1944         ECall( "add2list", "{" + string(l) + "}" );
    1945         return (l);
    1946 };
    1947 */
    1948 
    1949 /******************************************************************************/
    1950 /*
    1951 static proc rm4list( poly l, poly q )
    1952 "
    1953         l is a sorted list of monomials (polynom with coefs: 0 or 1)
    1954         removing from l monomials occuring in q
    1955        
    1956         ?? NEED DEBUG ??
    1957 "
    1958 {
    1959         BCall( "rm4list", "{" + string(l) + "} , " + string(q) );
    1960         poly t;
    1961         for( int i = size(q); i>0; i-- )
    1962         {
    1963 
    1964         /// trick:
    1965         ///  check for presens of a monomial q[i] in l
    1966         ///  if char = 0 then new_size > (if not present) or =  (if already present)
    1967         ///  if char > 0 then new_size > (if not present) or =< (if already present)
    1968                 t = leadmonom(q[i]);
    1969                 if ( size(l+t) <= size(l) )
    1970                 {
    1971                         l = t - t;
    1972                         if( size(l) == 0 )
    1973                         {
    1974                                 break;
    1975                         };
    1976                 };     
    1977         };
    1978        
    1979         ECall( "rm4list", "{" + string(l) + "}" );
    1980         return (l);
    1981 };
    1982 */
    1983 
    1984 /******************************************************************************/
    1985 static proc reduce_one_per_row(list @given)
    1986 {
    1987         BCall( "reduce_one_per_row" );
    1988        
    1989         int @k;
    1990         int @l = size (@given);
    1991        
    1992         if( @l == 0 )
    1993         {
    1994                 return (@given);
    1995         };
    1996        
    1997         int @char = char(basering);
    1998        
    1999         intvec @marks;
    2000 
    2001         if ( @char == 0 )
    2002         {
    2003                 poly @t = poly(0);
    2004         };
    2005 
    2006         list @unis = list ();
    2007         poly p;
    2008        
    2009         for ( @k = @l; @k > 0; @k -- )
    2010         {
    2011                 p = uni_poly( @given[@k][2] );
    2012                 @unis[@k] = p;
    2013 
    2014                 if( p == 0 )
    2015                 {
    2016                         @marks[@k] = 2;
    2017                 } else
    2018                 {
    2019                         @marks[@k] = 1;
    2020                         if ( @char == 0 )
    2021                         {
    2022                                 @t = @t + p;
    2023                         };
     1894                  @tt = @unis[@k];
     1895
     1896                  if( size( @for_delete + @tt ) != ( size( @for_delete )  + size( @tt ) ) )
     1897                    {
     1898                      @t = @t - @tt;
     1899                      @marks[@k] = 0;
     1900                    };
    20241901                };
    2025         };
    2026        
    2027         if ( @char != 0 )
    2028         {
    2029                 def save = basering;
    2030                 execute( "ring NewRingWithGoodField = (0), (" + varstr(save) + "), (" + ordstr(save) + "); ");
    2031                 poly @t = 0;
    2032                
    2033                 if(! defined (@unis) )
    2034                 {
    2035                         list @unis = imap( save, @unis );
    2036                 };
    2037                
    2038                 for ( @k = @l; @k > 0; @k -- )
    2039                 {
    2040                         if( @marks[@k] == 1 )
    2041                         {
    2042                                 @t = @t + @unis[@k];
    2043                         };
    2044                 };
    2045         };
    2046                
    2047         int @loop_size = size(@t);
    2048         poly @for_delete, @tt;
    2049         int @ll;
    2050        
    2051         while( @loop_size > 0 )
    2052         {
    2053                 @for_delete = poly(0);
    2054                 @ll = size(@t);   
    2055                
    2056                 for ( @k = @ll; @k > 0; @k -- )
    2057                 {
    2058                         if ( leadcoef(@t[@k]) == 1 )
    2059                         {
    2060                                 @for_delete = @for_delete + @t[@k];
    2061                         };
    2062                 };
    2063 
    2064                 @loop_size = size( @for_delete );
    2065 
    2066                 if ( @loop_size>0 )
    2067                 {
    2068                         for( @k = @l ; @k > 0 ; @k -- )
    2069                         {           
    2070                                 if ( @marks[@k] == 1)
    2071                                 {
    2072                                         @tt = @unis[@k];
    2073 
    2074                                         if( size( @for_delete + @tt ) != ( size( @for_delete )  + size( @tt ) ) )
    2075                                         {
    2076                                                 @t = @t - @tt;
    2077                                                 @marks[@k] = 0;
    2078                                         };
    2079                                 };
    2080                         };             
    2081                 };             
    2082         };
    2083        
    2084         if ( @char != 0 )
    2085         {
    2086                 setring(save);
    2087                 kill(NewRingWithGoodField);
    2088         };
    2089 
    2090         list @reduced = list();
    2091        
    2092         for ( @k = @l ; @k>0 ; @k --)
    2093         {
    2094                 if (@marks[@k]==2)
    2095                 {
    2096                         @reduced = list ( @given[@k] ) + @reduced ;
    2097                 } else
    2098                 {
    2099                         if (@marks[@k]==1)
    2100                         {
    2101                                 @reduced = @reduced + list ( @given[@k] );
    2102                         };
    2103                 };
    2104         };
     1902            };         
     1903        };             
     1904    };
     1905       
     1906  if ( @char != 0 )
     1907    {
     1908      setring(save);
     1909      kill(NewRingWithGoodField);
     1910    };
     1911
     1912  list @reduced = list();
     1913       
     1914  for ( @k = @l ; @k>0 ; @k --)
     1915    {
     1916      if (@marks[@k]==2)
     1917        {
     1918          @reduced = list ( @given[@k] ) + @reduced ;
     1919        } else
     1920          {
     1921            if (@marks[@k]==1)
     1922              {
     1923                @reduced = @reduced + list ( @given[@k] );
     1924              };
     1925          };
     1926    };
    21051927   
    2106         ECall( "reduce_one_per_row", "{...l...}" );
    2107         return (@reduced);
    2108 };
    2109 
    2110 
    2111 
    2112 /******************************************************************************/
    2113 static proc calc_base (list @AD_GIVEN, list #)
    2114 "
     1928  ECall( "reduce_one_per_row", "structured list" );
     1929  return (@reduced);
     1930};
     1931
     1932
     1933
     1934/******************************************************/
     1935static proc calc_base (list @AD_GIVEN)
     1936  "
    21151937    sort out given 'pbw' into groups due to common monoms in images = ([2])s
    21161938"
    21171939{
    2118         BCall( "calc_base" );
    2119 
    2120     if ( (size(#)>0) )
    2121     {
    2122         int @z_next;
    2123         poly @z_nextp;
    2124        
    2125         if ( typeof(#[1]) == "int")
    2126         {
    2127             @z_next = #[1];
    2128             @z_nextp = var(@z_next);   
    2129         };
    2130        
    2131         if ( typeof(#[1]) == "poly")
    2132         {
    2133             @z_nextp = #[1];
    2134         };
    2135    
    2136     };
    2137    
    2138 
    2139     list @AD_CALC = list();
    2140     int @ll, @k, @j, @loop_size;
    2141 
    2142     poly @t = poly(0);
    2143     poly @tt, @for_delete;
    2144 
    2145     int @t_size;
    2146 
    2147 
    2148     list @GR, @GR_TEMP;
    2149 
    2150     int @number = size(@AD_GIVEN);
    2151 
    2152     while ( @number > 0 )
    2153     {
    2154         @tt = @AD_GIVEN[@number][2];
    2155         if ( size (@tt) == 0)
    2156         {
    2157             @t = @AD_GIVEN[@number][1];
    2158             if (defined(@z_nextp))
    2159             {
    2160                 @tt = @z_nextp * @t - @t * @z_nextp; /// right? (-?)
    2161                 @AD_CALC = @AD_CALC + list ( list ( @t, @tt ) );
    2162             } else
    2163             {
    2164                 @AD_CALC = @AD_CALC + list ( @t );         
    2165             };
    2166            
     1940  BCall( "calc_base" );
     1941
     1942  list @AD_CALC = list();
     1943  int @ll, @k, @j, @loop_size;
     1944
     1945  poly @t = poly(0);
     1946  poly @tt, @for_delete;
     1947
     1948  int @t_size;
     1949
     1950  list @GR, @GR_TEMP;
     1951
     1952  int @number = size(@AD_GIVEN);
     1953
     1954  while ( @number > 0 )
     1955    {
     1956      @tt = @AD_GIVEN[@number][2];
     1957      if ( size (@tt) == 0)
     1958        {
     1959          @AD_CALC = @AD_CALC + list ( @AD_GIVEN[@number][1] );
    21671960        } else
    2168         {
     1961          {
    21691962            @t = uni_poly( @tt );
    21701963            @t_size = size(@t);
     
    21761969            @loop_size = size(@GR);
    21771970            if ( @loop_size == 0 )
    2178             {
     1971              {
    21791972                @GR = list(@GR_TEMP);
    2180             } else
    2181             {
    2182                 for ( @k = @loop_size; @k > 0 ; @k -- )
     1973              } else
    21831974                {
    2184                     @tt = @GR[@k][1];
    2185                     if ( size( @t + @tt ) != ( @t_size + size(@tt) ) )
    2186                     // whether @tt and @i intersencts? ( will not work in char == 2 !!!)
     1975                  for ( @k = @loop_size; @k > 0 ; @k -- )
    21871976                    {
    2188 
    2189                         if ( char(basering) == 0 )
     1977                      @tt = @GR[@k][1];
     1978                      if ( size( @t + @tt ) != ( @t_size + size(@tt) ) )
     1979                        // whether @tt and @i intersencts? ( will not work in char == 2 !!!)
    21901980                        {
    2191                             @GR_TEMP[1] = @GR_TEMP[1] + @tt;
    2192                         } else
    2193                         {
    2194                             @GR_TEMP[1] = uni_poly( @GR_TEMP[1] + @tt );
     1981
     1982                          if ( char(basering) == 0 )
     1983                            {
     1984                              @GR_TEMP[1] = @GR_TEMP[1] + @tt;
     1985                            } else
     1986                              {
     1987                                @GR_TEMP[1] = uni_poly( @GR_TEMP[1] + @tt );
     1988                              };
     1989                       
     1990                          @GR_TEMP[2] = @GR_TEMP[2] + @GR[@k][2];
     1991                          @GR = delete ( @GR, @k );
    21951992                        };
    2196                        
    2197                         @GR_TEMP[2] = @GR_TEMP[2] + @GR[@k][2];
    2198                         @GR = delete ( @GR, @k );
    21991993                    };
     1994                  @GR = @GR + list(@GR_TEMP);
    22001995                };
    2201                 @GR = @GR + list(@GR_TEMP);
    2202             };
    2203         };
    2204         @number --;
    2205     }; ///  main while
     1996          };
     1997      @number --;
     1998    }; //  main while
    22061999   
    2207     list @res;
    2208 
    2209     for ( @k = size(@GR); @k > 0 ; @k -- )
    2210     {
    2211         if ( size (@GR[@k][2]) > 1 ) /// ! zeroes in AD_CALC so here must be non zero
    2212         {
    2213             if(defined(@z_nextp))
    2214             {
    2215                   @res = my_calc ( @GR[@k][2], uni_poly(@GR[@k][1]), @z_nextp);
    2216             } else
    2217             {
    2218                   @res = my_calc ( @GR[@k][2], uni_poly(@GR[@k][1]));   
    2219             }
    2220             if ( size (@res) > 0 )
    2221             {
    2222                   @AD_CALC = @AD_CALC + @res;
     2000  list @res;
     2001
     2002  for ( @k = size(@GR); @k > 0 ; @k -- )
     2003    {
     2004      if ( size (@GR[@k][2]) > 1 ) // ! zeroes in AD_CALC so here must be non zero
     2005        {
     2006          @res = FSS ( @GR[@k][2], uni_poly(@GR[@k][1]));
     2007         
     2008          if ( size (@res) > 0 )
     2009            {
     2010              @AD_CALC = @AD_CALC + @res;
    22232011            };
    22242012        };
    22252013    };   
    22262014
    2227         ECall( "calc_base" );
    2228         return (@AD_CALC);
    2229 
    2230 };
    2231 
    2232 /******************************************************************************/
     2015  ECall( "calc_base" );
     2016  return (@AD_CALC);
     2017
     2018};
     2019
     2020
     2021/******************************************************/
     2022static proc ApplyAd( list l, list # )
     2023"
     2024  Apply Ad_{#} to every element of a list of polynomils
     2025"
     2026{
     2027  BCall( "ApplyAd", l, # );
     2028
     2029  if( (size(l) == 0) or (size(#) == 0) )
     2030    {
     2031      return (l);
     2032    };
     2033
     2034  poly f, t;
     2035 
     2036  if ( typeof(#[1]) == "int")
     2037    {
     2038      int next = #[1];
     2039      f = my_var (next);
     2040    } else
     2041      {
     2042        if ( typeof(#[1]) == "poly")
     2043          {
     2044            f = #[1];
     2045          } else
     2046            {
     2047              print("Error: cannot differentiate with '" + string(#)+"'");
     2048              return ();
     2049            };
     2050      };
     2051
     2052  for ( int k = size(l); k > 0; k --)
     2053    {
     2054      t = l[k];
     2055      l[k] = list( t, f * t - t * f );
     2056    };
     2057
     2058  ECall("ApplyAd", l );
     2059  return (l);
     2060};
     2061
     2062
     2063/******************************************************/
    22332064static proc calc_k_base (list l, list #)
    2234 "
    2235   calculation of Ker as Vector Space
    2236 "
    2237 {
    2238         BCall( "calc_k_base", "list, " + string(#) );
    2239 
    2240         list t = reduce_one_per_row( l, 0); /// optimization (a1)
    2241         return( QNF_list ( calc_base(t, #), 2) );          /// calculation of groupps (a2) + gauss.
    2242 
    2243 };
    2244 
    2245 /******************************************************************************/
    2246 
    2247 ///////////////////////////////////////////////////////////////////////////////
     2065  "
     2066     calculation of Kernel of an Ad operator as a K-Vector Space
     2067"
     2068{
     2069  BCall( "calc_k_base", "list, " + string(#) );
     2070
     2071  list t = reduce_one_per_row( l, 0); // optimization (a1)
     2072  return( QNF_list ( ApplyAd( calc_base(t), # ), 2) );          // calculation of groupps (a2) + gauss.
     2073   
     2074};
     2075
     2076LIB "poly.lib"; // for content in proc makeIdeal
     2077
     2078/******************************************************/
    22482079static proc makeIdeal ( list l )
    2249 "
     2080  "
    22502081        return: ideal: where the generators are polynomials from list, without 1 and zeroes
    22512082"
    22522083{
    2253         ideal I; poly p;
    2254 
    2255         for( int i = 1; i <= size(l); i++ )
    2256         {
    2257                 if ( typeof( l[i] ) == "list" )
    2258                 {
    2259                         p = l[i][1]; /// just take the 1st polynom...
    2260                 } else
    2261                 {
    2262                         p = l[i];
    2263                 };
     2084  ideal I; poly p;
     2085
     2086  for( int i = 1; i <= size(l); i++ )
     2087    {
     2088      if ( typeof( l[i] ) == "list" )
     2089        {
     2090          p = l[i][1]; // just take the 1st polynom...
     2091        } else
     2092          {
     2093            p = l[i];
     2094          };
    22642095               
    2265                 if ( (p != 1) and (p != 0) )
    2266                 {
    2267                         I = I, p;
    2268                 };
    2269         };
    2270        
    2271         I = simplify ( I, 2 ); /// no zeroes
    2272        
    2273         return(I);
    2274 };
    2275 
    2276 /******************************************************************************/
     2096      p = cleardenom( p* (1/content(p)) );
     2097      if ( (p != 1) and (p != 0) )
     2098        {
     2099          I = I, p;
     2100        };
     2101    };
     2102       
     2103  I = simplify ( I, 2 ); // no zeroes
     2104       
     2105  return(I);
     2106};
     2107
     2108/******************************************************/
    22772109static proc inCenter_poly( poly p )
    2278 "
     2110  "
    22792111  if p in center => return 1
    22802112      otherwise     return 0
    22812113"
    22822114{
    2283     poly t;
    2284     for (int k = nvars(basering); k>0; k-- )
    2285     {
    2286         t = var(k);
    2287         if ( QNF_poly(t * p - p * t) != 0 )
    2288         {
    2289                   if( toprint() )
    2290                   {
    2291                         "POLY: ", string (p), " is NOT in center";
    2292                   };
    2293                   return (0);
    2294         };
    2295     };
    2296 
    2297     if( toprint() )
    2298     {
    2299           "POLY: ", string (p), " is in center";
    2300     };
    2301     return (1);
    2302 };
    2303 
    2304 /******************************************************************************/
     2115  poly t;
     2116  for (int k = nvars(basering); k>0; k-- )
     2117    {
     2118      t = var(k);
     2119      if ( QNF_poly(t * p - p * t) != 0 )
     2120        {
     2121          if( toprint() )
     2122            {
     2123              "POLY: ", string (p), " is NOT in center";
     2124            };
     2125          return (0);
     2126        };
     2127    };
     2128
     2129  if( toprint() )
     2130    {
     2131      "POLY: ", string (p), " is in center";
     2132    };
     2133  return (1);
     2134};
     2135
     2136/******************************************************/
    23052137static proc inCenter_list(def l)
    23062138{
    2307     for ( int @i = 1; @i <= size(l); @i++ )
    2308     {
    2309         if ( typeof(l[@i])=="poly" )
    2310         {
    2311             if (! inCenter_poly(l[@i]) )
    2312             {
    2313                 return(0);
     2139  for ( int @i = 1; @i <= size(l); @i++ )
     2140    {
     2141      if ( typeof(l[@i])=="poly" )
     2142        {
     2143          if (! inCenter_poly(l[@i]) )
     2144            {
     2145              return(0);
    23142146            };
    23152147               
    23162148        } else
    2317         {
     2149          {
    23182150            if ( (typeof(l[@i])=="list") or (typeof(l[@i])=="ideal") )
    2319             {
    2320                   if (! inCenter_list(l[@i]) )
     2151              {
     2152                if (! inCenter_list(l[@i]) )
    23212153                  {
    2322                         return(0);
     2154                    return(0);
    23232155                  };
    2324             };
    2325         };
    2326     };
    2327     return(1);
    2328 };
    2329 
    2330 /******************************************************************************/
     2156              };
     2157          };
     2158    };
     2159  return(1);
     2160};
     2161
     2162/******************************************************/
    23312163static proc inCentralizer_poly( poly p, poly f )
    2332 "
     2164  "
    23332165  if p in centralizer(f) => return 1
    23342166      otherwise     return 0
    23352167"
    23362168{
    2337     if ( QNF_poly(f * p - p * f) != 0 )
     2169  if ( QNF_poly(f * p - p * f) != 0 )
    23382170    {     
    2339           if( toprint() )
    2340           {
    2341                 "POLY: ", string (p), " is NOT in centralizer(f)";
    2342           };
    2343           return (0);
    2344     };
    2345        
    2346     if( toprint() )
    2347     {
    2348         "POLY: ", string (p), " is in centralizer(f)";
    2349     };
    2350         return (1);
    2351 };
    2352 
    2353 /******************************************************************************/
     2171      if( toprint() )
     2172        {
     2173          "POLY: ", string (p), " is NOT in centralizer(f)";
     2174        };
     2175      return (0);
     2176    };
     2177       
     2178  if( toprint() )
     2179    {
     2180      "POLY: ", string (p), " is in centralizer(f)";
     2181    };
     2182  return (1);
     2183};
     2184
     2185/******************************************************/
    23542186static proc inCentralizer_list( def l, poly f )
    23552187{
    23562188  for ( int @i = 1; @i <= size(l); @i++ )
    2357   {
    2358         if ( typeof(l[@i])=="poly" )
    2359         {
    2360             if (! inCentralizer_poly(l[@i], f) )
    2361                 {
    2362                   return(0);
    2363                 };
     2189    {
     2190      if ( typeof(l[@i])=="poly" )
     2191        {
     2192          if (! inCentralizer_poly(l[@i], f) )
     2193            {
     2194              return(0);
     2195            };
    23642196               
    23652197        } else
    2366         {
     2198          {
    23672199            if ( (typeof(l[@i])=="list") or (typeof(l[@i])=="ideal") )
    2368             {
    2369                   if (! inCentralizer_list(l[@i], f) )
     2200              {
     2201                if (! inCentralizer_list(l[@i], f) )
    23702202                  {
    2371                         return(0);
     2203                    return(0);
    23722204                  };
    2373             };
    2374         };
    2375   };
     2205              };
     2206          };
     2207    };
    23762208  return(1);
    23772209};
    23782210
    23792211
    2380 ///////////////////////////////////////////////////////////////////////////////
    2381 ///////////////////////////////////////////////////////////////////////////////
    2382 /// main algorithms
    2383 
    2384 
    2385 ///////////////////////////////////////////////////////////////////////////////
    2386 ///////////////////////////////////////////////////////////////////////////////
    2387 /// center's computation
    2388 
    2389 ///////////////////////////////////////////////////////////////////////////////
     2212/******************************************************/
     2213/******************************************************/
     2214// main algorithms
     2215
     2216
     2217/******************************************************/
     2218/******************************************************/
     2219// center's computation
     2220
     2221/******************************************************/
    23902222/* static */ proc center_min_iterative( int MaxDeg, list # )
    2391 "
     2223  "
    23922224        computes the 'minimal' set of central elements (of deg <= MaxDeg) in iterative way
    23932225        Note: based on calc_k_base, zAddBad, zRefine, zCancel, PBW_*
    23942226"
    23952227{
    2396         BCall("center_min_iterative", MaxDeg, #);
    2397        
    2398         int n = myInt(#);
    2399         int m = ( MaxDeg < 0 ); /// no bound on Degree
    2400        
    2401         int MinDeg = 6; /// starting guess for MaxDeg
    2402         int Delta  = 4; /// increment of MaxDeg
    2403        
    2404         if( m )
    2405         {       
    2406                 // minimal guess
    2407                 MaxDeg = MinDeg;               
    2408         };
    2409        
    2410         list @q; int @i; int N = nvars(basering);
    2411 
    2412         my_var_init(); // setup for my_var(n)
    2413         QNF_init();
    2414 
    2415         list PBW = list();
    2416         list PBW_PLAIN = list();
    2417 
    2418         PBW[ index(0) ] = PBW_init();
    2419 
    2420         list SYS = PBW_sys( my_vars(), my_var(1) );
    2421 
    2422        
    2423         list @z = list ();                                      /// center list
    2424         list @l = init_bads( MaxDeg );  /// verbotten loeadexps...
    2425        
    2426         @q = PBW[ index(0) ];
    2427 
    2428         int k = 1;
    2429         while( k <= ( MaxDeg+1 ) )
    2430         {
    2431                 for ( @i = 2; @i <= N; @i ++ )
    2432                 {
    2433                         @q = calc_k_base (@q, my_var(@i));
    2434                 };
    2435 
    2436                 @q = zRefine (calc_k_base(@q)); /// new center!
    2437 
    2438 
    2439                 if ( size(@q) > 0 )
    2440                 {
    2441                         @z = @z + @q; /// computed central elements
     2228  BCall("center_min_iterative", MaxDeg, #);
     2229       
     2230  int n = myInt(#);
     2231  int m = ( MaxDeg < 0 ); // no bound on Degree
     2232       
     2233  int MinDeg = 6; // starting guess for MaxDeg
     2234  int Delta  = 4; // increment of MaxDeg
     2235       
     2236  if( m )
     2237    {   
     2238      // minimal guess
     2239      MaxDeg = MinDeg;         
     2240    };
     2241       
     2242  list @q; int @i; int N = nvars(basering);
     2243
     2244  my_var_init(); // setup for my_var(n)
     2245  QNF_init();
     2246
     2247  list PBW = list();
     2248  list PBW_PLAIN = list();
     2249
     2250  PBW[ index(0) ] = PBW_init();
     2251
     2252  list SYS = PBW_sys( my_vars(), my_var(1) );
     2253
     2254       
     2255  list @z = list ();                                    // center list
     2256  list @l = init_bads( MaxDeg );        // verbotten loeadexps...
     2257       
     2258  @q = PBW[ index(0) ];
     2259
     2260  int k = 1;
     2261  while( k <= ( MaxDeg+1 ) )
     2262    {
     2263      for ( @i = 2; @i <= N; @i ++ )
     2264        {
     2265          @q = calc_k_base (@q, my_var(@i));
     2266        };
     2267
     2268      @q = zRefine (calc_k_base(@q)); // new center!
     2269
     2270
     2271      if ( size(@q) > 0 )
     2272        {
     2273          @z = @z + @q; // computed central elements
    24422274                       
    2443                         if( (n > 0) and (size(@z) > n) )
    2444                         {
    2445                                 break; // found all central elements                   
    2446                         };
    2447                 };
    2448 
    2449                 if( k == ( MaxDeg+1 ) )
    2450                 {
    2451                         if( (n == 0) and ( !m ) )
    2452                         {
    2453                                 break; /// that's all
    2454                         };
     2275          if( (n > 0) and (size(@z) > n) )
     2276            {
     2277              break; // found all central elements                     
     2278            };
     2279        };
     2280
     2281      if( k == ( MaxDeg+1 ) )
     2282        {
     2283          if( (n == 0) and ( !m ) )
     2284            {
     2285              break; // that's all
     2286            };
    24552287                       
    2456                         MaxDeg = MaxDeg + Delta;
     2288          MaxDeg = MaxDeg + Delta;
    24572289                       
    2458                         // renew bad list
    2459                         @l = init_bads( MaxDeg );                       
    2460                         @l = zAddBad( @l, @z, MaxDeg, 0 );
     2290          // renew bad list
     2291          @l = init_bads( MaxDeg );                     
     2292          @l = zAddBad( @l, @z, MaxDeg, 0 );
    24612293                       
    2462                 } else
    2463                 {
     2294        } else
     2295          {
    24642296               
    2465                         if ( size(@q) > 0 )
    2466                         {
    2467                                 @l = zAddBad( @l, @q, MaxDeg, 0 ); /// add all possible 'leadexps' !
    2468                         };
     2297            if ( size(@q) > 0 )
     2298              {
     2299                @l = zAddBad( @l, @q, MaxDeg, 0 ); // add all possible 'leadexps' !
     2300              };
    24692301                       
    2470                 };
    2471 
    2472                 PBW[index(k)] = PBW_next( PBW, k, SYS );
    2473 
    2474                 PBW_PLAIN = PBW_PLAIN + zCancel( PBW[index(k-1)] , @l );
    2475 
    2476                 @q = PBW_PLAIN + zCancel( PBW[index(k)] , @l ); /// kill from @tt all entries with leadexp in @l[@d]
    2477 
    2478                 k ++;
    2479         };
    2480 
    2481         QNF_done();
    2482         my_var_done();
    2483 
    2484         return (makeIdeal(@z));
    2485 };
    2486 
    2487 ///////////////////////////////////////////////////////////////////////////////
     2302          };
     2303
     2304      PBW[index(k)] = PBW_next( PBW, k, SYS );
     2305
     2306      PBW_PLAIN = PBW_PLAIN + zCancel( PBW[index(k-1)] , @l );
     2307
     2308      @q = PBW_PLAIN + zCancel( PBW[index(k)] , @l ); // kill from @tt all entries with leadexp in @l[@d]
     2309
     2310      k ++;
     2311    };
     2312
     2313  QNF_done();
     2314  my_var_done();
     2315
     2316  return (makeIdeal(@z));
     2317};
     2318
     2319/******************************************************/
    24882320static proc center_vectorspace( int MaxDeg )
    2489 "
     2321  "
    24902322        pure calculation of center as a finitely dimmensional Vector Space (deg <= MaxDeg )
    24912323"
    24922324{
    2493         my_var_init();
    2494        
    2495         int  N = nvars( basering );
    2496         list P = PBW_base( MaxDeg, my_var(1) );
    2497        
    2498         for( int v = 2; v <= N; v++ )
    2499         {
    2500                 P = calc_k_base( P, my_var(v) );
    2501         };
    2502 
    2503         my_var_done();
    2504        
    2505         return( calc_k_base ( P ) );
    2506 };
    2507 
    2508 ///////////////////////////////////////////////////////////////////////////////
     2325  my_var_init();
     2326       
     2327  int  N = nvars( basering );
     2328  list P = PBW_base( MaxDeg, my_var(1) );
     2329       
     2330  for( int v = 2; v <= N; v++ )
     2331    {
     2332      P = calc_k_base( P, my_var(v) );
     2333    };
     2334
     2335  my_var_done();
     2336       
     2337  return( calc_k_base ( P ) );
     2338};
     2339
     2340/******************************************************/
    25092341/* static */ proc center_min_vectorspace( int MaxDeg )
    2510 "
     2342  "
    25112343    computes the 'minimal' set of central elements (of deg <= MaxDeg)
    25122344    by reducing the set of it's generators as vector space
     
    25172349{
    25182350
    2519         QNF_init();
    2520         ideal res = makeIdeal( zReduce( center_vectorspace( MaxDeg)));
    2521         QNF_done();
    2522 
    2523         return( res );
    2524 };
    2525 
    2526 
    2527 ///////////////////////////////////////////////////////////////////////////////
    2528 ///////////////////////////////////////////////////////////////////////////////
    2529 /// centralizer's computations
    2530 
    2531 ///////////////////////////////////////////////////////////////////////////////
     2351  QNF_init();
     2352  ideal res = makeIdeal( zReduce( center_vectorspace( MaxDeg)));
     2353  QNF_done();
     2354
     2355  return( res );
     2356};
     2357
     2358
     2359/******************************************************/
     2360/******************************************************/
     2361// centralizer's computations
     2362
     2363/******************************************************/
    25322364static proc centralizer_vectorspace( poly p, int MaxDeg )
    25332365{
    2534         return( calc_k_base( PBW_base( MaxDeg, p )));
    2535 };
    2536 
    2537 ///////////////////////////////////////////////////////////////////////////////
     2366  return( calc_k_base( PBW_base( MaxDeg, p )));
     2367};
     2368
     2369/******************************************************/
    25382370/* static */ proc centralizer_min_vectorspace( poly p, int MaxDeg )
    25392371{
    25402372
    2541         QNF_init();
    2542         ideal res = makeIdeal( zReduce( centralizer_vectorspace( p, MaxDeg)));
    2543         QNF_done();
    2544 
    2545         return( res );
    2546 };
    2547 
    2548 ///////////////////////////////////////////////////////////////////////////////
     2373  QNF_init();
     2374  ideal res = makeIdeal( zReduce( centralizer_vectorspace( p, MaxDeg)));
     2375  QNF_done();
     2376
     2377  return( res );
     2378};
     2379
     2380/******************************************************/
    25492381/* static */ proc centralizer_min_iterative( poly p, int MaxDeg, list # )
    2550 "
     2382  "
    25512383  computes the 'minimal' set of elements (of deg <= MaxDeg) generating centralizer of p in iterative way
    25522384  Note: based on calc_k_base
     
    25562388"
    25572389{
    2558         QNF_init();
    2559 
    2560         int n = myInt(#);
    2561         int m = (MaxDeg < 0);
    2562        
    2563         int MinDeg = 6; /// starting guess for MaxDeg
    2564         int Delta  = 4; /// increment of MaxDeg
    2565        
    2566         if( m )
    2567         {
    2568                 // minimal guess
    2569                 MaxDeg = MinDeg;
    2570         };
    2571 
    2572         list @q;
    2573 
    2574         list PBW = list();
    2575         list PBW_PLAIN = list();
    2576 
    2577         PBW[ index(0) ] = PBW_init();
    2578 
    2579         list SYS = PBW_sys( my_vars(), p );
    2580 
    2581         list @z  = list ();             // result list
    2582         list @l  = init_bads( MaxDeg ); // verbotten loeadexps...
    2583        
    2584         @q = PBW[ index(0) ];
    2585 
    2586         for (int k = 1; k <= ( MaxDeg+1 ); k ++ )
    2587         {
    2588                 @q = zRefine( calc_k_base(@q), 1 );
    2589 
    2590                 if ( size(@q) > 0 )
    2591                 {
    2592                         @z = @z + @q; /// computed desired elements
    2593                        
    2594                         if( (n > 0) and (size(@z) > n) )
    2595                         {
    2596                                 break; // found needed elements
    2597                         };
    2598                 };
    2599 
    2600                 if( k == ( MaxDeg+1 ) )
    2601                 {
    2602                         if( (n == 0) or ( !m ) )
    2603                         {
    2604                                 break; /// that's all
    2605                         };                     
    2606                        
    2607                         MaxDeg = MaxDeg + Delta;
    2608                        
    2609                         // renew bad list
    2610                         @l = init_bads( MaxDeg );                       
    2611                         @l = zAddBad( @l, @z, MaxDeg, 0 );
    2612                        
    2613                 } else
    2614                 {
     2390  QNF_init();
     2391
     2392  int n = myInt(#);
     2393  int m = (MaxDeg < 0);
     2394       
     2395  int MinDeg = 6; // starting guess for MaxDeg
     2396  int Delta  = 4; // increment of MaxDeg
     2397       
     2398  if( m )
     2399    {
     2400      // minimal guess
     2401      MaxDeg = MinDeg;
     2402    };
     2403
     2404  list @q;
     2405
     2406  list PBW = list();
     2407  list PBW_PLAIN = list();
     2408
     2409  PBW[ index(0) ] = PBW_init();
     2410
     2411  list SYS = PBW_sys( my_vars(), p );
     2412
     2413  list @z  = list ();             // result list
     2414  list @l  = init_bads( MaxDeg ); // verbotten loeadexps...
     2415       
     2416  @q = PBW[ index(0) ];
     2417
     2418  for (int k = 1; k <= ( MaxDeg+1 ); k ++ )
     2419    {
     2420      @q = zRefine( calc_k_base(@q), 1 );
     2421
     2422      if ( size(@q) > 0 )
     2423        {
     2424          @z = @z + @q; // computed desired elements
     2425                 
     2426          if( (n > 0) and (size(@z) > n) )
     2427            {
     2428              break; // found needed elements
     2429            };
     2430        };
     2431
     2432      if( k == ( MaxDeg+1 ) )
     2433        {
     2434          if( (n == 0) or ( !m ) )
     2435            {
     2436              break; // that's all
     2437            };                 
     2438                 
     2439          MaxDeg = MaxDeg + Delta;
     2440                 
     2441          // renew bad list
     2442          @l = init_bads( MaxDeg );                     
     2443          @l = zAddBad( @l, @z, MaxDeg, 0 );
     2444                 
     2445        } else
     2446          {
     2447                   
     2448            if ( size(@q) > 0 )
     2449              {
     2450                @l = zAddBad( @l, @q, MaxDeg, 0 ); // add all possible 'leadexps' !
     2451              };
     2452                   
     2453          };
    26152454               
    2616                         if ( size(@q) > 0 )
    2617                         {
    2618                                 @l = zAddBad( @l, @q, MaxDeg, 0 ); /// add all possible 'leadexps' !
    2619                         };
    2620                        
    2621                 };
    2622 
    2623                 PBW[index(k)] = PBW_next( PBW, k, SYS );
    2624                 PBW_PLAIN = PBW_PLAIN + zCancel( PBW[index(k-1)] , @l );
    2625                 @q = PBW_PLAIN + zCancel( PBW[index(k)] , @l );  // kill from @tt all entries with leadexp in @l[@d]
    2626 
    2627         };
    2628 
    2629         QNF_done();
    2630         return (makeIdeal(@z));
    2631 };
    2632 
    2633 
    2634 
    2635 ///////////////////////////////////////////////////////////////////////////////
    2636 ///////////////////////////////////////////////////////////////////////////////
     2455      PBW[index(k)] = PBW_next( PBW, k, SYS );
     2456      PBW_PLAIN = PBW_PLAIN + zCancel( PBW[index(k-1)] , @l );
     2457      @q = PBW_PLAIN + zCancel( PBW[index(k)] , @l );  // kill from @tt all entries with leadexp in @l[@d]
     2458
     2459    };
     2460
     2461  QNF_done();
     2462  return (makeIdeal(@z));
     2463};
     2464
     2465
     2466
     2467/******************************************************/
     2468/******************************************************/
    26372469// main aliases
    26382470
    2639 /******************************************************************************/
     2471/******************************************************/
    26402472proc inCenter( def a )
    2641 "USAGE:   inCenter(a); a poly/list/ideal
     2473  "USAGE:   inCenter(a); a poly/list/ideal
    26422474RETURN:  integer (1 if a in center, 0 otherwise)
    26432475EXAMPLE: example inCenter; shows examples"
    26442476{
    2645         QNF_init();
    2646         def res;
    2647        
    2648         if ( typeof(a) == "poly" )
    2649         {
    2650                 res = inCenter_poly(a);
    2651         } else
    2652         {
    2653                 if ( (typeof(a)=="list") or (typeof(a)=="ideal") )
    2654                 {
    2655                         res = inCenter_list(a);
    2656                 } else
    2657                 {
    2658                         res = a;
    2659                 };
    2660         };
    2661        
    2662         QNF_done();
    2663         return (res);
     2477  QNF_init();
     2478  def res;
     2479       
     2480  if ( typeof(a) == "poly" )
     2481    {
     2482      res = inCenter_poly(a);
     2483    } else
     2484      {
     2485        if ( (typeof(a)=="list") or (typeof(a)=="ideal") )
     2486          {
     2487            res = inCenter_list(a);
     2488          } else
     2489            {
     2490              res = a;
     2491            };
     2492      };
     2493       
     2494  QNF_done();
     2495  return (res);
    26642496}
    26652497example
     
    26852517/******************************************************************************/
    26862518proc inCentralizer( def a, poly f )
    2687 "USAGE:   inCentralizer(a, f); a poly/list/ideal, f poly
     2519  "USAGE:   inCentralizer(a, f); a poly/list/ideal, f poly
    26882520RETURN:  integer (1 if a in centralizer(f), 0 otherwise)
    26892521EXAMPLE: example inCentralizer; shows examples"
    26902522{
    2691         QNF_init();
    2692         def res;
    2693 
    2694         if ( typeof(a) == "poly" )
    2695         {
    2696                 res = inCentralizer_poly(a, f);
    2697         } else
    2698         {
    2699                 if ( (typeof(a)=="list") or (typeof(a)=="ideal") )
    2700                 {
    2701                         res = inCentralizer_list(a, f);
    2702                 } else
    2703                 {
    2704                         res = a;
    2705                 };
    2706         };
    2707 
    2708         QNF_done();
    2709         return (res);
     2523  QNF_init();
     2524  def res;
     2525
     2526  if ( typeof(a) == "poly" )
     2527    {
     2528      res = inCentralizer_poly(a, f);
     2529    } else
     2530      {
     2531        if ( (typeof(a)=="list") or (typeof(a)=="ideal") )
     2532          {
     2533            res = inCentralizer_list(a, f);
     2534          } else
     2535            {
     2536              res = a;
     2537            };
     2538      };
     2539
     2540  QNF_done();
     2541  return (res);
    27102542}
    27112543example
     
    27192551  ncalgebra(1,D); // this is U(sl_2)
    27202552  poly f = x^2;
    2721   poly a = z; /// lies in center
     2553  poly a = z; // lies in center
    27222554  poly b = y^2;
    27232555  inCentralizer(a, f);
     
    27312563
    27322564
    2733 ///////////////////////////////////////////////////////////////////////////////
     2565/******************************************************/
    27342566proc center(int MaxDeg, list # )
    2735 "USAGE:      center(MaxDeg[, N]); int MaxDeg, int N
     2567  "USAGE:      center(MaxDeg[, N]); int MaxDeg, int N
    27362568RETURN:     ideal generated by found elements
    27372569NOTE:       computes the 'minimal' set of central elements.
     
    27442576{
    27452577
    2746         if ( myInt(#) > 0 ) /// given a number of central elements to compute
    2747         {
    2748                 return ( center_min_iterative( MaxDeg, # ) );
    2749         };
    2750 
    2751         if( MaxDeg >= 0 ) /// standard computation - the fastest one (often)
    2752         {
    2753 //              return ( center_min_iterative( MaxDeg, # ) );
    2754                 return ( center_min_vectorspace ( MaxDeg ) );
    2755         };
    2756 
    2757         print( "Error: wrong arguments." );
    2758         return();
     2578  if ( myInt(#) > 0 ) // given a number of central elements to compute
     2579    {
     2580      return ( center_min_iterative( MaxDeg, # ) );
     2581    };
     2582
     2583  if( MaxDeg >= 0 ) // standard computation - the fastest one (often)
     2584    {
     2585      // return ( center_min_iterative( MaxDeg, # ) );
     2586      return ( center_min_vectorspace ( MaxDeg ) );
     2587    };
     2588
     2589  print( "Error: wrong arguments." );
     2590  return();
    27592591       
    27602592}
    27612593example
    27622594{ "EXAMPLE:"; echo = 2;
    2763    ring a=0,(x,y,z),dp;
    2764    matrix D[3][3]=0;
    2765    D[1,2]=-z;
    2766    D[1,3]=2*x;
    2767    D[2,3]=-2*y;
    2768    ncalgebra(1,D); // this is U(sl_2)
    2769    ideal Z = center(2); /// find all central elements of degree <= 2
    2770    Z;
    2771    inCenter(Z);
    2772    ideal ZZ = center(-1, 1 ); /// find the first non trivial central element
    2773    ZZ;
    2774    inCenter(ZZ);
    2775 };
    2776 
    2777 
    2778 ///////////////////////////////////////////////////////////////////////////////
     2595 ring a=0,(x,y,z),dp;
     2596 matrix D[3][3]=0;
     2597 D[1,2]=-z;
     2598 D[1,3]=2*x;
     2599 D[2,3]=-2*y;
     2600 ncalgebra(1,D); // this is U(sl_2)
     2601 ideal Z = center(2); // find all central elements of degree <= 2
     2602 Z;
     2603 inCenter(Z);
     2604 ideal ZZ = center(-1, 1 ); // find the first non trivial central element
     2605 ZZ; "";
     2606 inCenter(ZZ);
     2607};
     2608
     2609
     2610/******************************************************/
    27792611proc centralizer( poly p, int MaxDeg, list # )
    2780 "USAGE:      centralizer(F, MaxDeg[, N]); poly F, int MaxDeg, int N
     2612  "USAGE:      centralizer(F, MaxDeg[, N]); poly F, int MaxDeg, int N
    27812613RETURN:     ideal generated by found elements
    27822614NOTE:       computes the 'minimal' set of elements of centralizer(F).
     
    27892621{
    27902622       
    2791         if( myInt(#) > 0 )
    2792         {
    2793                 return( centralizer_min_iterative( p, MaxDeg, # ) );
    2794         };
    2795        
    2796         if( MaxDeg >= 0 )
    2797         {
    2798 //              return( centralizer_min_iterative( p, MaxDeg, # ) );
    2799                 return( centralizer_min_vectorspace( p, MaxDeg ) );
    2800         };
    2801        
    2802         print( "Error: wrong arguments." );
    2803         return();
     2623  if( myInt(#) > 0 )
     2624    {
     2625      return( centralizer_min_iterative( p, MaxDeg, # ) );
     2626    };
     2627       
     2628  if( MaxDeg >= 0 )
     2629    {
     2630      // return( centralizer_min_iterative( p, MaxDeg, # ) );
     2631      return( centralizer_min_vectorspace( p, MaxDeg ) );
     2632    };
     2633       
     2634  print( "Error: wrong arguments." );
     2635  return();
    28042636
    28052637}
    28062638example
    28072639{ "EXAMPLE:"; echo = 2;
    2808    ring a=0,(x,y,z),dp;
    2809    matrix D[3][3]=0;
    2810    D[1,2]=-z;
    2811    D[1,3]=2*x;
    2812    D[2,3]=-2*y;
    2813    ncalgebra(1,D); // this is U(sl_2)
    2814    poly f = 4*x*y+z^2-2*z; /// central polynomial
    2815    f;
    2816    ideal c = centralizer(f, 2); /// find all elements of degree <= 2 which lies in centralizer of f
    2817    c;
    2818    inCentralizer(c, f);
    2819    ideal cc = centralizer(f, -1, 2 ); /// find at least first two non trivial elements of the centralizer of f
    2820    cc;
    2821    inCentralizer(cc, f);
    2822    poly g = z^2-2*z; /// any polynomial
    2823    g;
    2824    c = centralizer(g, 2); /// find all elements of degree <= 2 which lies in centralizer of g
    2825    c;
    2826    inCentralizer(c, g);
    2827    cc = centralizer(g, -1, 2 ); /// find at least first two non trivial elements of the centralizer of g
    2828    cc;
    2829    inCentralizer(cc, g);
    2830 };
    2831 
    2832 
    2833 ///////////////////////////////////////////////////////////////////////////////
    2834 ///////////////////////////////////////////////////////////////////////////////
    2835 ///////////////////////////////////////////////////////////////////////////////
     2640 ring a=0,(x,y,z),dp;
     2641 matrix D[3][3]=0;
     2642 D[1,2]=-z;
     2643 D[1,3]=2*x;
     2644 D[2,3]=-2*y;
     2645 ncalgebra(1,D); // this is U(sl_2)
     2646 poly f = 4*x*y+z^2-2*z; // central polynomial
     2647 f;
     2648 ideal c = centralizer(f, 2); // find all elements of degree <= 2 which lies in centralizer of f
     2649 c;
     2650 inCentralizer(c, f);
     2651 ideal cc = centralizer(f, -1, 2 ); // find at least first two non trivial elements of the centralizer of f
     2652 cc;
     2653 inCentralizer(cc, f);
     2654 poly g = z^2-2*z; // any polynomial
     2655 g; "";
     2656 c = centralizer(g, 2); // find all elements of degree <= 2 which lies in centralizer of g
     2657 c; "";
     2658 inCentralizer(c, g);
     2659 cc = centralizer(g, -1, 2 ); // find at least first two non trivial elements of the centralizer of g
     2660 cc; "";
     2661 inCentralizer(cc, g);
     2662};
Note: See TracChangeset for help on using the changeset viewer.