source: git/Singular/LIB/center.lib @ 91c978

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