source: git/Singular/LIB/schreyer.lib @ ff12c6

jengelh-datetimespielwiese
Last change on this file since ff12c6 was ff12c6, checked in by Oleksandr Motsak <motsak@…>, 11 years ago
Testing of ClearContent&ClearDenominators add: testing the use of CPolyCoeffsEnumerator by n_ClearDenominators and n_ClearContent add/chg: better generic testing of Clear* functions add: more testing of Clear* functions: cleardenom + ERROR() + exit on wrong test result add: testing of "-leadcoeff" normalization (for TestClearContent) but not for TestClearDenominators... NOTE: ClearContent should make leading coeff. positive! NOTE: testing alg. and trans. extensions of Q aswell as Q itself add: Added testing of Clear* to 'make check'
  • Property mode set to 100644
File size: 29.5 KB
Line 
1//////////////////////////////////////////////////////////////////////////////
2version="$Id$";
3category="General purpose";
4info="
5LIBRARY: schreyer.lib     Helpers for working with the Schreyer induced ordering
6AUTHOR:  Oleksandr Motsak <U@D>, where U={motsak}, D={mathematik.uni-kl.de}
7
8PROCEDURES:
9 Sres(M,l)      Schreyer resolution of module M of maximal length l
10 Ssyz(M)        Schreyer resolution of module M of length 1
11 Scontinue(l)   continue the resolution computation by most l steps
12
13KEYWORDS:  syzygy; Schreyer induced ordering; Schreyer free resolution
14NOTE:  requires the dynamic module: syzextra
15";
16
17static proc prepareSyz( module I, list # )
18{
19  int i;
20  int k = 0;
21  int r = nrows(I);
22  int c = ncols(I);
23
24
25  if( size(#) > 0 )
26  {
27    if( typeof(#[1]) == "int" || typeof(#[1]) == "bigint" )
28    {
29      k = #[1];
30    }
31  }
32
33  if( k < r )
34  {
35    "// *** Wrong k: ", k, " < nrows: ", r, " => setting k = r = ", r;
36    k = r;
37  }
38
39//   "k: ", k;  "c: ", c;   "I: ", I;
40
41  for( i = c; i > 0; i-- )
42  {
43    I[i] = I[i] + gen(k + i);
44  }
45
46//  DetailedPrint(I);
47
48  return(I);
49}
50
51static proc separateSyzGB( module J, int c )
52{
53  module II, G; vector v; int i;
54
55  J = simplify(J, 2);
56
57  for( i = ncols(J); i > 0; i-- )
58  {
59    v = J[i];
60    if( leadcomp(v) > c )
61    {
62      II[i] = v;
63    } else
64    {
65      G[i] = v; // leave only gen(i): i <= c
66    }
67  }
68
69  II = simplify(II, 2);
70  G = simplify(G, 2);
71
72  return (list(G, II));
73}
74
75static proc splitSyzGB( module J, int c )
76{
77  module JJ; vector v, vv; int i;
78
79  for( i = ncols(J); i > 0; i-- )
80  {
81    v = J[i];
82
83    vv = 0;
84   
85    while( leadcomp(v) <= c )
86    {
87      vv = vv + lead(v);
88      v  = v  - lead(v);
89    }
90
91    J[i] = vv;
92    JJ[i] = v;
93  }
94
95  J = simplify(J, 2);
96  JJ = simplify(JJ, 2);
97
98  return (list(J, JJ));
99}
100
101
102static proc Sinit(module M)
103{
104  def @save = basering;
105 
106  int @DEBUG = !system("with", "ndebug");
107  if( @DEBUG )
108  {
109    "Sinit::Input";
110    type(M);
111    DetailedPrint(M);
112    attrib(M);
113  }
114
115  int @RANK = nrows(M); int @SIZE = ncols(M);
116
117  int @IS_A_SB = attrib(M, "isSB"); // ??? only if all weights were zero?!
118
119  if( !@IS_A_SB )
120  {
121    M = std(M); // this should be faster than computing std in S (later on)
122  }
123
124  def S = MakeInducedSchreyerOrdering(1); // 1 puts history terms to the back
125  // TODO: NOTE: +1 causes trouble to Singular interpreter!!!???
126  setring S; // a new ring with a Schreyer ordering
127
128  if( @DEBUG )
129  {
130    "Sinit::StartingISRing";
131    basering;
132//    DetailedPrint(basering);
133  }
134
135  // Setup the leading syzygy^{-1} module to zero:
136  module Z = 0; Z[@RANK] = 0; attrib(Z, "isHomog", intvec(0)); 
137
138  module MRES = Z;
139 
140  list RES; RES[1] = Z;
141
142  module F = freemodule(@RANK);
143  intvec @V = deg(F[1..@RANK]);
144 
145  module M = imap(@save, M);
146 
147  attrib(M, "isHomog", @V);
148  attrib(M, "isSB", 1);
149
150 
151  if( @DEBUG )
152  {
153    "Sinit::SB_Input: ";
154    type(M);
155    attrib(M);
156    attrib(M, "isHomog");
157    DetailedPrint(M);
158  }
159
160  if( @DEBUG )
161  {
162    // 0^th syz. property
163    if( size(module(transpose( transpose(M) * transpose(MRES) ))) > 0 )
164    {
165      transpose( transpose(M) * transpose(MRES) );
166      "transpose( transpose(M) * transpose(MRES) ) != 0!!!";
167      $
168    }
169  }
170
171  RES[size(RES)+1] = M; // list of all syzygy modules
172  MRES = MRES, M;
173
174  attrib(MRES, "isHomog", @V); 
175
176  attrib(S, "InducionLeads", lead(M));
177  attrib(S, "InducionStart", @RANK); 
178 
179  if( @DEBUG )
180  {
181    "Sinit::MRES";
182    DetailedPrint(MRES);
183    attrib(MRES, "isHomog");
184    attrib(S);
185  }
186
187  export RES;
188  export MRES;
189  return (S);
190}
191
192static proc Sstep()
193{
194  int @DEBUG = !system("with", "ndebug");
195
196  if( @DEBUG )
197  {
198    "Sstep::NextInducedRing";
199    DetailedPrint(basering);
200
201    attrib(basering, "InducionLeads");
202    attrib(basering, "InducionStart");
203
204    GetInducedData();
205  }
206
207  // syzygy step:
208
209/*
210  // is initial weights are all zeroes!
211  def L =  lead(M);
212  intvec @V = deg(M[1..ncols(M)]);  @W;  @V;  @W = @V;  attrib(L, "isHomog", @W); 
213  SetInducedReferrence(L, @RANK, 0);
214*/
215
216//  def L =  lead(MRES);
217//  @W = @W, @V;
218//  attrib(L, "isHomog", @W); 
219
220
221  // General setting:
222//  SetInducedReferrence(MRES, 0, 0); // limit: 0!
223  int @l = size(RES);
224
225  module M = RES[@l];
226
227  module L = attrib(basering, "InducionLeads");
228  int limit = attrib(basering, "InducionStart");
229
230//  L;  limit;
231 
232  int @RANK = ncols(MRES) - ncols(M); // nrows(M); // what if M is zero?!
233
234/*
235  if( @RANK !=  nrows(M) )
236  {
237    type(MRES);
238    @RANK;
239    type(M);
240    pause();
241  }
242*/
243 
244  intvec @W = attrib(M, "isHomog");
245  intvec @V = deg(M[1..ncols(M)]);
246  @V = @W, @V;
247   
248  if( @DEBUG )
249  {
250    "Sstep::NextInput: ";
251    M;
252    @V;
253    @RANK;   
254    DetailedPrint(MRES);
255    attrib(MRES, "isHomog");
256  }
257
258 
259     
260  SetInducedReferrence(L, limit, 0);
261 
262  def K = prepareSyz(M, @RANK);
263//  K;
264 
265//   attrib(K, "isHomog", @V);   DetailedPrint(K, 1000);
266
267//  pause();
268 
269  K = idPrepare(K, @RANK); // std(K); // ?
270  K = simplify(K, 2);
271
272//  K;
273
274  module N = separateSyzGB(K, @RANK)[2]; // 1^st syz. module: vectors which start in lower part (comp >= @RANK)
275  attrib(N, "isHomog", @V);
276
277// "N_0: "; N; DetailedPrint(N, 10);
278 
279  N = std(N); // TODO: fix "wrong weights"!!!?
280  attrib(N, "isHomog", @V);
281
282//  N;
283 
284  if( @DEBUG )
285  {
286    if( size(N) > 0 )
287    {
288      // next syz. property
289      if( size(module(transpose( transpose(N) * transpose(MRES) ))) > 0 )
290      {
291        MRES;
292
293        "N: "; N; DetailedPrint(N, 10);
294
295        "K:"; K; DetailedPrint(K, 10);
296
297        "RANKS: ", @RANK;
298
299        "transpose( transpose(N) * transpose(MRES) ) != 0!!!";
300        transpose( transpose(N) * transpose(MRES) );
301
302        "transpose(N) * transpose(MRES): ";
303        transpose(N) * transpose(MRES);
304        DetailedPrint(module(_), 2);
305        $
306      }
307    }
308  }
309 
310  RES[@l + 1] = N; // list of all syzygy modules
311 
312  MRES = MRES, N;
313  attrib(MRES, "isHomog", @V);
314
315
316  L = L, lead(N);
317  attrib(basering, "InducionLeads", L);
318
319  if( @DEBUG )
320  {
321    "Sstep::NextSyzOutput: ";
322    DetailedPrint(N);
323    attrib(N, "isHomog");
324  }
325
326}
327
328proc Scontinue(int l)
329"USAGE:  Scontinue(l)
330RETURN:  nothing, instead it changes RES and MRES variables in the current ring
331PURPOSE: computes further (at most l) syzygies
332NOTE:    must be used within a ring returned by Sres or Ssyz. RES and MRES are
333         explained in Sres
334EXAMPLE: example Scontinue; shows an example
335"
336{
337  def data = GetInducedData();
338           
339  if( (!defined(RES)) || (!defined(MRES)) || (typeof(data) != "list") || (size(data) != 2) )
340  {
341    ERROR("Sorry, but basering does not seem to be returned by Sres or Ssyz");
342  }
343  for (;  (l != 0) && (size(RES[size(RES)]) > 0); l-- )
344  {
345    Sstep();
346  }
347}
348example
349{ "EXAMPLE:"; echo = 2;
350  ring r;
351  module M = maxideal(1); M;
352  def S = Ssyz(M); setring S; S;
353  "Only the first syzygy: ";
354  RES; MRES;
355  "More syzygies: ";
356  Scontinue(10);
357  RES; MRES;
358}
359
360proc Ssyz(module M)
361"USAGE:  Ssyz(M)
362RETURN:  ring, containing a list of modules RES and a module MRES
363PURPOSE: computes the first syzygy module of M (wrt some Schreyer ordering)
364NOTE:    The output is explained in Sres
365EXAMPLE: example Ssyz; shows an example
366"
367{
368  def S = Sinit(M); setring S;
369 
370  Sstep(); // NOTE: what if M is zero?
371
372  return (S);
373}
374example
375{ "EXAMPLE:"; echo = 2;
376  ring r;
377  module M = maxideal(1); M;
378  def S = Ssyz(M); setring S; S;
379  "Only the first syzygy: ";
380  RES;
381  MRES; // Note gen(i)
382  kill S;
383  setring r; kill M;
384
385  module M = 0;
386  def S = Ssyz(M); setring S; S;
387  "Only the first syzygy: ";
388  RES;
389  MRES;
390}
391
392proc Sres(module M, int l)
393"USAGE:  Sres(M, l)
394RETURN:  ring, containing a list of modules RES and a module MRES
395PURPOSE: computes (at most l) syzygy modules of M wrt the classical Schreyer
396         induced ordering with gen(i) > gen(j) if i > j, provided both gens
397         are from the same syzygy level.
398NOTE:    RES contains the images of maps subsituting the beginning of the
399         Schreyer free resolution of baseRing^r/M, while MRES is a sum of
400         these images in a big free sum, containing all the syzygy modules.
401         The syzygy modules are shifted so that gen(i) correspons to MRES[i].
402         The leading zero module RES[0] indicates the fact that coker of the
403         first map is zero. The number of zeroes inducates the rank of input.
404NOTE:    If l == 0 then l is set to be nvars(basering) + 1
405EXAMPLE: example Sres; shows an example
406"
407{
408  def S = Sinit(M); setring S;
409
410  if (l == 0)
411  {
412    l = nvars(basering) + 1; // not really an estimate...?!
413  }
414 
415  Sstep(); l = l - 1;
416 
417  Scontinue(l);
418 
419  return (S);
420}
421example
422{ "EXAMPLE:"; echo = 2;
423  ring r;
424  module M = maxideal(1); M;
425  def S = Sres(M, 0); setring S; S;
426  RES;
427  MRES;
428  kill S;
429  setring r; kill M;
430
431  def A = nc_algebra(-1,0); setring A;
432  ideal Q = var(1)^2, var(2)^2, var(3)^2;
433  qring SCA = twostd(Q);
434  basering;
435
436  module M = maxideal(1);
437  def S = Sres(M, 2); setring S; S;
438  RES;
439  MRES;
440}
441
442
443
444// ================================================================== //
445
446
447LIB "general.lib"; // for sort
448
449/* static proc Tail(def M) // DONE: in C++ (dyn. module: syzextra)!
450{
451  int i = ncols(M); def m;
452  while (i > 0)
453  {
454    m = M[i];
455    m = m - lead(m); // m = tail(m)   
456    M[i] = m;   
457    i--;
458  }
459  return (M);
460}*/
461
462
463/* static */ proc SSinit(def M)
464{
465  if( (typeof(M) != "module") && (typeof(M) != "ideal") )
466  {
467    ERROR("Sorry: need an ideal or a module for input");
468  }
469
470  // TODO! DONE?
471  def @save = basering;
472 
473  int @DEBUG = !system("with", "ndebug");
474  int @SYZCHECK = 1 || @DEBUG; // TODO: only for now!!
475 
476  if( @DEBUG )
477  {
478    "SSinit::Input";
479    type(M);
480//    DetailedPrint(M);
481    attrib(M);
482  }
483
484  int @RANK = nrows(M); int @SIZE = ncols(M);
485
486  int @IS_A_SB = attrib(M, "isSB"); // ??? only if all weights were zero?!
487
488  if( !@IS_A_SB )
489  {
490    def opts = option(get);
491    option(redSB); option(redTail);
492    M = std(M);
493    option(set, opts);
494    kill opts;
495  } else
496  {
497    M = simplify(M, 2 + 4 + 32);
498  }
499
500  def LEAD = lead(M);
501
502  // sort wrt neg.deg.rev.lex!
503  intvec iv_ds = sort(LEAD, "ds", 1)[2]; // ,1 => reversed!
504
505  M = M[iv_ds]; // sort M wrt ds on current leading terms
506  LEAD = LEAD[iv_ds];
507 
508  def TAIL = Tail(M);
509 
510  intvec @DEGS = deg(M[1..@SIZE]); // store actuall degrees of input elements
511 
512  // TODO: what about real modules? weighted ones?
513 
514  list @l = ringlist(@save);
515
516  int @z = 0; ideal @m = maxideal(1); intvec @wdeg = deg(@m[1..ncols(@m)]);
517
518  // NOTE: @wdeg will be ignored anyway :(
519  @l[3] = list(list("C", @z), list("lp", @wdeg));
520
521  kill @z, @wdeg; // since these vars are ring independent!
522
523  def S = ring(@l); // --MakeInducedSchreyerOrdering(1);
524
525  module F = freemodule(@RANK);
526  intvec @V = deg(F[1..@RANK]);
527 
528  setring S; // ring with an easy divisibility test ("C, lex")
529
530  if( @DEBUG )
531  {
532    "SSinit::NewRing(C, lex)";
533    basering;
534//    DetailedPrint(basering);
535  }
536
537  // Setup the leading syzygy^{-1} module to zero:
538  module Z = 0; Z[@RANK] = 0; attrib(Z, "isHomog", intvec(0)); 
539
540  module MRES = Z;
541 
542  list RES;  RES[1] = Z;
543  list LRES; LRES[1] = Z;
544  list TRES; TRES[1] = Z;
545 
546  def M = imap(@save, M);
547
548  attrib(M, "isHomog", @V);
549  attrib(M, "isSB", 1);
550  attrib(M, "degrees", @DEGS); 
551 
552  def LEAD = imap(@save, LEAD);
553 
554  attrib(LEAD, "isHomog", @V);
555  attrib(LEAD, "isSB", 1); 
556 
557  def TAIL = imap(@save, TAIL);
558
559  if( @DEBUG )
560  {
561    "SSinit::(sorted) SB_Input: ";
562    type(M);
563    attrib(M);
564    attrib(M, "isHomog");
565//    DetailedPrint(M);
566  }
567
568  if( @SYZCHECK )
569  {
570    // 0^th syz. property
571    if( size(module(transpose( transpose(M) * transpose(MRES) ))) > 0 )
572    {
573      transpose( transpose(M) * transpose(MRES) );
574      "ERROR: transpose( transpose(M) * transpose(MRES) ) != 0!!!";
575      $
576    }
577  }
578
579  RES[size(RES)+1] = M; // list of all syzygy modules
580  LRES[size(LRES)+1] = LEAD; // list of all syzygy modules
581  TRES[size(TRES)+1] = TAIL; // list of all syzygy modules
582 
583  MRES = MRES, M; //?
584
585  attrib(MRES, "isHomog", @V);
586 
587//  attrib(S, "InducionStart", @RANK);
588  attrib(S, "LEAD2SYZ", 1);
589  attrib(S, "TAILREDSYZ", 0);
590  attrib(S, "DEBUG", @DEBUG);
591  attrib(S, "SYZCHECK", @SYZCHECK);
592 
593  if( @DEBUG )
594  {
595    "SSinit::MRES";
596    MRES;
597//    DetailedPrint(MRES);
598    attrib(MRES, "isHomog");
599    attrib(S);
600  }
601
602  export RES;
603  export MRES;
604  export LRES;
605  export TRES;
606  return (S);
607}
608example
609{ "EXAMPLE:"; echo = 2;
610  ring R = 0, (w, x, y, z), dp;
611
612  def M = maxideal(1);
613  def S = SSinit(M); setring S; S;
614 
615  "Only the first initialization: ";
616  RES; LRES; TRES;
617  MRES;
618
619  kill S; setring R; kill M;
620 
621  ideal M = w^2 - x*z,  w*x - y*z,  x^2 - w*y, x*y - z^2, y^2 - w*z;
622  def S = SSinit(M); setring S; S;
623
624  "Only the first initialization: ";
625  RES; LRES; TRES;
626  MRES;
627
628  kill S; setring R; kill M;
629}
630
631
632LIB "poly.lib"; // for lcm
633
634
635
636/// Compute L(Syz(L))
637proc SSComputeLeadingSyzygyTerms(def L)
638{
639  if( typeof( attrib(basering, "DEBUG") ) == "int" )
640  {
641    int @DEBUG = attrib(basering, "DEBUG");
642  } else
643  {
644    int @DEBUG = !system("with", "ndebug");
645  }
646
647  if( @DEBUG )
648  {
649    "SSComputeLeadingSyzygyTerms::Input: ";
650    L;
651  }
652
653  int i, j, r; intvec iv_ds;
654  int N = ncols(L);
655  def a, b;
656  poly aa, bb;
657
658  bigint c;
659
660  ideal M;
661
662  module S = 0;
663 
664  for(i = 1; i <= N; i++)
665  {
666    a = L[i];
667//    "a: ", a;
668    c = leadcomp(a);
669    r = int(c);
670
671    if(r > 0)
672    {
673      aa = a[r];
674    } else
675    {
676      aa = a;
677    }
678
679    M = 0;
680   
681    for(j = i-1; j > 0; j--)
682    {
683      b = L[j];
684//      "b: ", b;
685
686      if( leadcomp(b) == c )
687      {
688        if(r > 0)
689        {
690          bb = b[r];
691        } else
692        {
693          bb = b;
694        }
695
696        M[j] = (lcm(aa, bb) / aa);
697      }
698    }
699   
700    // TODO: add quotient relations here...
701
702    M = simplify(M, 1 + 2 + 32);
703
704    iv_ds = sort(M, "ds", 1)[2]; // ,1 => reversed!
705
706    M = M[iv_ds];
707   
708    S = S, M * gen(i);
709  }
710
711  S = simplify(S, 2);
712
713  S = sort(S, "ds", 1)[1]; // ,1 => reversed! // TODO: not needed?
714 
715  if( @DEBUG )
716  {
717    "SSComputeLeadingSyzygyTerms::Output: ";
718    S;
719  } 
720
721  attrib(S, "isSB", 1);
722
723  return (S);
724}
725
726/// Compute Syz(L), where L is a monomial (leading) module
727proc SSCompute2LeadingSyzygyTerms(def L, int @TAILREDSYZ)
728{
729  if( typeof( attrib(basering, "DEBUG") ) == "int" )
730  {
731    int @DEBUG = attrib(basering, "DEBUG");
732  } else
733  {
734    int @DEBUG = !system("with", "ndebug");
735  }
736
737  if( typeof( attrib(basering, "SYZCHECK") ) == "int" )
738  {
739    int @SYZCHECK = attrib(basering, "SYZCHECK");
740  } else
741  {
742    int @SYZCHECK = @DEBUG;
743  }
744 
745
746  if( @DEBUG )
747  {
748    "SSCompute2LeadingSyzygyTerms::Input: ";
749    L;
750    "@TAILREDSYZ: ", @TAILREDSYZ;
751  }
752
753  int i, j, r;
754  int N = ncols(L);
755  def a, b;
756
757  poly aa, bb, @lcm;
758
759  bigint c;
760
761  module M;
762
763  module S = 0;
764
765  for(i = 1; i <= N; i++)
766  {
767    a = L[i];
768//    "a: ", a;
769    c = leadcomp(a);
770    r = int(c);
771
772    aa = leadmonomial(a);
773
774    M = 0;
775
776    for(j = i-1; j > 0; j--)
777    {
778      b = L[j];
779//      "b: ", b;
780
781      if( leadcomp(b) == c )
782      {
783        bb = leadmonomial(b);
784        @lcm = lcm(aa, bb);
785
786        M[j] = (@lcm / aa)* gen(i) - (@lcm / bb)* gen(j);
787      }
788    }
789   
790    M = simplify(M, 2);
791
792    // TODO: add quotient relations here...
793    S = S, M;
794  }
795
796  if( @TAILREDSYZ )
797  {
798    // Make sure that 2nd syzygy terms are not reducible by 1st
799    def opts = option(get);
800    option(redSB); option(redTail);
801    S = std(S); // binomial module
802    option(set, opts);
803    //  kill opts;
804  } else
805  {
806    S = simplify(S, 2 + 32);
807  }
808
809  S = sort(S, "ds", 1)[1]; // ,1 => reversed!
810
811  if( @DEBUG )
812  {
813    "SSCompute2LeadingSyzygyTerms::Syz(LEAD): "; S;
814  }
815
816  if( @SYZCHECK )
817  {
818    if( size(S) > 0 and size(L) > 0 )
819    {
820      if( size(module(transpose( transpose(S) * transpose(L) ))) > 0 )
821      {
822        transpose( transpose(S) * transpose(L) );
823        "ERROR: transpose( transpose(S) * transpose(L) ) != 0!!!";
824        $
825      }
826    }
827  }
828
829  module S2 = Tail(S);
830  S = lead(S); // (C,lp) on base ring!
831             
832  if( @DEBUG )
833  {
834    "SSCompute2LeadingSyzygyTerms::Output: "; S; S2;
835  } 
836 
837  attrib(S, "isSB", 1);
838
839  return (S, S2);
840}
841
842// -------------------------------------------------------- //
843
844/// TODO: save shortcut LM(m) * "t" -> ?
845proc SSReduceTerm(poly m, def t, def L, def T, list #)
846{
847  if( typeof( attrib(basering, "DEBUG") ) == "int" )
848  {
849    int @DEBUG = attrib(basering, "DEBUG");
850  } else
851  {
852    int @DEBUG = !system("with", "ndebug");
853  }
854
855  if( @DEBUG )
856  {
857    "SSReduce::Input: ";
858
859    "mult: ", m;
860    "term: ", t;
861    "L: ", L;
862    "T: ", T;
863    if( size(#) > 0 )
864    {
865      "LSyz: ", #;
866    }
867//    "attrib(LS, 'isSB')", attrib(LS, "isSB");
868  }
869 
870  vector s = 0;
871
872  if( t == 0 )
873  {
874    return (s);
875  }
876
877  def product = m * t;
878
879  bigint c = leadcomp(t);
880  int r = int(c);
881 
882  def a, b, nf, bb;
883
884  // looking for an appropriate reducer
885  for( int k = ncols(L); k > 0; k-- )
886  {
887    a = L[k];
888    // with the same mod. component
889    if( leadcomp(a) == c )
890    {
891      b = - (leadmonomial(product) / leadmonomial(L[k]));
892     
893      // which divides the product
894      if( b != 0 )
895      {
896//        "b: ", b;
897        bb = b * gen(k);
898        nf = bb;
899
900        if( size(#) > 0 )
901        {
902          if( typeof(#[1]) == "module" )
903          {
904            nf = NF(bb, #[1]);
905//        "NF: ", nf;
906          }
907        }
908
909        // while the complement (the fraction) is not reducible by leading syzygies
910        if( nf != 0 )
911        {
912          /// TODO: save shortcut LM(m) * T[i] -> ?
913
914          // choose ANY such reduction... (with the biggest index?)
915          s = bb + SSTraverseTail(b, T[k], L, T, #);
916          break;
917        }
918      }
919    }
920  } 
921  if( @DEBUG )
922  {
923    "SSReduceTerm::Output: ", s;
924  }
925  return (s);
926}
927
928// TODO: store m * @tail -.-^-.-^-.--> ?
929proc SSTraverseTail(poly m, def @tail, def L, def T, list #)
930{
931  if( typeof( attrib(basering, "DEBUG") ) == "int" )
932  {
933    int @DEBUG = attrib(basering, "DEBUG");
934  } else
935  {
936    int @DEBUG = !system("with", "ndebug");
937  }
938
939  if( @DEBUG )
940  {
941    "SSTraverse::Input: ";
942
943    "mult: ", m;
944    "tail: ", @tail; // T[i];
945
946    if( size(#) > 0 )
947    {
948      "LSyz: "; #[1];
949    }
950  }
951
952  vector s = 0;
953
954  def @l;
955
956  // iterate tail-terms in ANY order!
957  while( size(@tail) > 0 )
958  {
959    @l = lead(@tail);
960    s = s + SSReduceTerm(m, @l, L, T, #);
961    @tail = @tail - @l;
962  }
963
964  if( @DEBUG )
965  {
966    "SSTraverseTail::Output: ", s;
967  }
968  return (s);
969}
970
971// -------------------------------------------------------- //
972
973// module (N, LL, TT) = SSComputeSyzygy(L, T);
974// Compute Syz(L ++ T) = N = LL ++ TT
975proc SSComputeSyzygy(def L, def T)
976{
977  if( typeof( attrib(basering, "DEBUG") ) == "int" )
978  {
979    int @DEBUG = attrib(basering, "DEBUG");
980  } else
981  {
982    int @DEBUG = !system("with", "ndebug");
983  }
984 
985  if( @DEBUG )
986  {
987    "SSComputeSyzygy::Input";
988    "basering: ", basering; attrib(basering);
989//    DetailedPrint(basering);
990
991//    "iCompShift: ", iCompShift;
992
993    "L: "; L;
994    "T: "; T;
995  }
996
997  def a; bigint c; int r, k; poly aa;
998
999  int @LEAD2SYZ = 0;
1000  if( typeof( attrib(basering, "LEAD2SYZ") ) == "int" )
1001  {
1002    @LEAD2SYZ = attrib(basering, "LEAD2SYZ");
1003  }
1004
1005  int @TAILREDSYZ = 1;
1006  if( typeof( attrib(basering, "TAILREDSYZ") ) == "int" )
1007  {
1008    @TAILREDSYZ = attrib(basering, "TAILREDSYZ");
1009//    @TAILREDSYZ;
1010  }
1011 
1012  /// Get the critical leading syzygy terms
1013  if( @LEAD2SYZ ) // & 2nd syz. term
1014  {
1015    def a2; int r2; poly aa2; 
1016    module LL, LL2;
1017    (LL, LL2) = SSCompute2LeadingSyzygyTerms(L, @TAILREDSYZ); // ++
1018  } else
1019  {
1020    module LL = SSComputeLeadingSyzygyTerms(L);
1021  }
1022
1023  module TT, SYZ;
1024
1025  if( size(LL) > 0 )
1026  {
1027    list LS;
1028
1029    if( @TAILREDSYZ )
1030    {
1031      LS = list(LL);
1032    }
1033   
1034    vector @tail;
1035
1036    for(k = ncols(LL); k > 0; k-- )
1037    {
1038      // leading syz. term:
1039      a = LL[k]; c = leadcomp(a); r = int(c); aa = leadmonomial(a);
1040      //    "A: ", a, " --->>>> ", aa, " **** [", r, "]: ";
1041
1042      /// TODO: save shortcut (aa) * T[r] -> ?
1043      @tail = SSTraverseTail(aa, T[r], L, T, LS);
1044             
1045      // get the 2nd syzygy term...
1046     
1047      if( @LEAD2SYZ ) // with the 2nd syz. term:
1048      {     
1049        a2 = LL2[k]; c = leadcomp(a2); r2 = int(c); aa2 = leadmonomial(a2);
1050        @tail = @tail +
1051             /// TODO: save shortcut (aa2) * T[r2] -> ?
1052             a2 + SSTraverseTail(aa2, T[r2], L, T, LS);
1053      } else
1054      {
1055        @tail = @tail + SSReduceTerm(aa, L[r], L, T, LS);
1056      }
1057     
1058     
1059      TT[k] = @tail;
1060      SYZ[k] = a + @tail;
1061    }
1062  }
1063
1064/*
1065  def opts = option(get); option(redSB); option(redTail);
1066  module SYZ = std(syz(M));
1067  option(set, opts); kill opts;
1068 
1069  module LL = lead(SYZ); // TODO: WRONG ORDERING!!!!!!!!
1070  module TT = Tail(SYZ);
1071*/
1072 
1073  if( @DEBUG )
1074  {
1075    "SSComputeSyzygy::Output";
1076
1077    "SYZ: "; SYZ;
1078    "LL: "; LL;
1079    "TT: "; TT;
1080  }
1081
1082  return (SYZ, LL, TT);
1083}
1084
1085// resolution/syzygy step:
1086static proc SSstep()
1087{
1088  if( typeof( attrib(basering, "DEBUG") ) == "int" )
1089  {
1090    int @DEBUG = attrib(basering, "DEBUG");
1091  } else
1092  {
1093    int @DEBUG = !system("with", "ndebug");
1094  }
1095
1096
1097  if( typeof( attrib(basering, "SYZCHECK") ) == "int" )
1098  {
1099    int @SYZCHECK = attrib(basering, "SYZCHECK");
1100  } else
1101  {
1102    int @SYZCHECK = @DEBUG;
1103  }
1104
1105  if( @DEBUG )
1106  {
1107    "SSstep::NextInducedRing";
1108    "basering: ", basering; attrib(basering);
1109  }
1110
1111/*
1112  // is initial weights are all zeroes!
1113  def L =  lead(M);
1114  intvec @V = deg(M[1..ncols(M)]);  @W;  @V;  @W = @V;  attrib(L, "isHomog", @W); 
1115  SetInducedReferrence(L, @RANK, 0);
1116*/
1117
1118//  def L =  lead(MRES);
1119//  @W = @W, @V;
1120//  attrib(L, "isHomog", @W); 
1121
1122
1123  // General setting:
1124//  SetInducedReferrence(MRES, 0, 0); // limit: 0!
1125  int @l = size(RES);
1126
1127  def M =  RES[@l];
1128
1129  def L = LRES[@l];
1130  def T = TRES[@l];
1131
1132
1133  //// TODO: wrong !!!!!
1134  int @RANK = ncols(MRES) - ncols(M); // nrows(M); // what if M is zero?!
1135
1136 
1137
1138/*
1139  if( @RANK !=  nrows(M) )
1140  {
1141    type(MRES);
1142    @RANK;
1143    type(M);
1144    pause();
1145  }
1146*/
1147 
1148  intvec @W = attrib(M, "isHomog"); intvec @V = attrib(M, "degrees"); @V = @W, @V;
1149   
1150  if( @DEBUG )
1151  {
1152    "Sstep::NextInput: ";
1153    M;
1154    L;
1155    @V;
1156    @RANK;
1157//    DetailedPrint(MRES);
1158    attrib(MRES, "isHomog");
1159  }
1160
1161     
1162  // TODO: N  = SYZ( M )!!!
1163  module N, LL, TT;
1164  (N, LL, TT) = SSComputeSyzygy(/*M, */L, T/*, @RANK*/);
1165
1166  // shift syz.comp by @RANK:
1167  module Z;
1168  Z = 0; Z[@RANK] = 0; Z = Z, transpose(LL);   LL = transpose(Z);
1169  Z = 0; Z[@RANK] = 0; Z = Z, transpose(TT);   TT = transpose(Z);
1170  Z = 0; Z[@RANK] = 0; Z = Z, transpose(N);     N = transpose(Z);
1171
1172
1173  if( @SYZCHECK )
1174  {
1175    if( size(N) > 0 )
1176    {
1177      // next syz. property
1178      if( size(module(transpose( transpose(N) * transpose(MRES) ))) > 0 )
1179      {
1180        "MRES", MRES;
1181
1182        "N: "; N; // DetailedPrint(N, 2);
1183
1184        "LL:"; LL; // DetailedPrint(LL, 1);
1185        "TT:"; TT; // DetailedPrint(TT, 10);
1186
1187        "RANKS: ", @RANK;
1188
1189        "transpose( transpose(N) * transpose(MRES) ) != 0!!!";
1190        transpose( transpose(N) * transpose(MRES) );
1191
1192        "transpose(N) * transpose(MRES): ";
1193        transpose(N) * transpose(MRES);
1194        // DetailedPrint(module(_), 2);
1195        $
1196      }
1197    }
1198  }
1199
1200  attrib(N, "isHomog", @V);
1201
1202  // TODO: correct the following:
1203  intvec @DEGS = deg(N[1..ncols(N)]); // no mod. comp. weights :(
1204
1205 
1206  attrib(N, "degrees", @DEGS);
1207 
1208   RES[@l + 1] = N; // list of all syzygy modules
1209  LRES[@l + 1] = LL; // list of all syzygy modules
1210  TRES[@l + 1] = TT; // list of all syzygy modules
1211
1212  MRES = MRES, N;
1213 
1214  attrib(MRES, "isHomog", @V);
1215
1216//  L = L, lead(N);  attrib(basering, "InducionLeads", L);
1217
1218  if( @DEBUG )
1219  {
1220    "SSstep::NextSyzOutput: ";
1221    N;
1222//    DetailedPrint(N);
1223    attrib(N);
1224  }
1225
1226}
1227
1228proc SScontinue(int l)
1229"USAGE:  SScontinue(l)
1230RETURN:  nothing, instead it changes RES and MRES variables in the current ring
1231PURPOSE: computes further (at most l) syzygies
1232NOTE:    must be used within a ring returned by Sres or Ssyz. RES and MRES are
1233         explained in Sres
1234EXAMPLE: example Scontinue; shows an example
1235"
1236{
1237
1238  /// TODO!
1239//  def data = GetInducedData();
1240
1241  if( (!defined(RES)) || (!defined(MRES)) ) /* || (typeof(data) != "list") || (size(data) != 2) */
1242  {
1243    ERROR("Sorry, but basering does not seem to be returned by Sres or Ssyz");
1244  }
1245  for (;  (l != 0) && (size(RES[size(RES)]) > 0); l-- )
1246  {
1247    SSstep();
1248  }
1249}
1250example
1251{ "EXAMPLE:"; echo = 2;
1252  ring r;
1253  module M = maxideal(1); M;
1254  def S = SSsyz(M); setring S; S;
1255  "Only the first syzygy: ";
1256  RES; MRES;
1257  "More syzygies: ";
1258  SScontinue(10);
1259  RES; MRES;
1260}
1261
1262proc SSsyz(def M)
1263"USAGE:  SSsyz(M)
1264RETURN:  ring, containing a list of modules RES and a module MRES
1265PURPOSE: computes the first syzygy module of M (wrt some Schreyer ordering)?
1266NOTE:    The output is explained in Sres
1267EXAMPLE: example Ssyz; shows an example
1268"
1269{
1270  if( (typeof(M) != "module") && (typeof(M) != "ideal") )
1271  {
1272    ERROR("Sorry: need an ideal or a module for input");
1273  }
1274
1275  def SS = SSinit(M); setring SS;
1276 
1277  SSstep(); // NOTE: what if M is zero?
1278
1279  return (SS);
1280}
1281example
1282{ "EXAMPLE:"; echo = 2;
1283  ring r;
1284
1285/*  ideal M = 0;
1286  def S = SSsyz(M); setring S; S;
1287  "Only the first syzygy: ";
1288  RES; LRES; TRES;
1289  MRES;
1290 
1291  kill S; setring r; kill M;
1292*/ 
1293
1294  module M = maxideal(1); M;
1295  def S = SSres(M, 0); setring S; S;
1296  MRES;
1297  RES;
1298  "";
1299  LRES;
1300  "";
1301  TRES;
1302
1303  kill S; setring r; kill M;
1304
1305  kill r;
1306 
1307  ring R = 0, (w, x, y, z), dp;
1308  ideal M = w^2 - x*z,  w*x - y*z,  x^2 - w*y, x*y - z^2, y^2 - w*z;
1309 
1310  def S = SSres(M, 0); setring S; S;
1311  MRES;
1312  RES;
1313  "";
1314  LRES;
1315  "";
1316  TRES;
1317}
1318
1319proc SSres(def M, int l)
1320"USAGE:  SSres(I, l)
1321RETURN:  ring, containing a list of modules RES and a module MRES
1322PURPOSE: computes (at most l) syzygy modules of M wrt the classical Schreyer
1323         induced ordering with gen(i) > gen(j) if i > j, provided both gens
1324         are from the same syzygy level.???
1325NOTE:    RES contains the images of maps subsituting the beginning of the
1326         Schreyer free resolution of baseRing^r/M, while MRES is a sum of
1327         these images in a big free sum, containing all the syzygy modules.
1328         The syzygy modules are shifted so that gen(i) correspons to MRES[i].
1329         The leading zero module RES[0] indicates the fact that coker of the
1330         first map is zero. The number of zeroes inducates the rank of input.
1331NOTE:    If l == 0 then l is set to be nvars(basering) + 1
1332EXAMPLE: example SSres; shows an example
1333"
1334{
1335  if( (typeof(M) != "module") && (typeof(M) != "ideal") )
1336  {
1337    ERROR("Sorry: need an ideal or a module for input");
1338  }
1339
1340  def SS = SSinit(M); setring SS;
1341
1342  if (l == 0)
1343  {
1344    l = nvars(basering) + 1; // not really an estimate...?!
1345  }
1346
1347  SSstep(); l = l - 1;
1348
1349  SScontinue(l);
1350
1351  return (SS);
1352}
1353example
1354{ "EXAMPLE:"; echo = 2;
1355  ring r;
1356  module M = maxideal(1); M;
1357  def S = SSres(M, 0); setring S; S;
1358  RES;
1359  MRES;
1360  kill S;
1361  setring r; kill M;
1362
1363  def A = nc_algebra(-1,0); setring A;
1364  ideal Q = var(1)^2, var(2)^2, var(3)^2;
1365  qring SCA = twostd(Q);
1366  basering;
1367
1368  module M = maxideal(1);
1369  def S = SSres(M, 2); setring S; S;
1370  RES;
1371  MRES;
1372}
1373
1374
1375
1376static proc loadme()
1377{
1378  int @DEBUG = !system("with", "ndebug");
1379
1380  if( @DEBUG )
1381  {
1382   
1383    "ndebug?: ", system("with", "ndebug");
1384    "om_ndebug?: ", system("with", "om_ndebug");
1385
1386    listvar(Top);
1387    listvar(Schreyer);
1388  }
1389//  listvar(Syzextra); listvar(Syzextra_g);
1390
1391  if( !defined(DetailedPrint) )
1392  {
1393    if( !@DEBUG )
1394    {
1395
1396      if( @DEBUG )
1397      {
1398        "Loading the Release version!";
1399      }
1400      load("syzextra.so");
1401
1402      if( @DEBUG )
1403      {
1404        listvar(Syzextra);
1405      }
1406
1407      exportto(Top, Syzextra::ClearContent);
1408      exportto(Top, Syzextra::ClearDenominators);
1409     
1410//      export Syzextra;
1411
1412//      exportto(Schreyer, Syzextra::noop);
1413      exportto(Schreyer, Syzextra::DetailedPrint);
1414      exportto(Schreyer, Syzextra::leadmonomial);
1415      exportto(Schreyer, Syzextra::leadcomp);
1416//      exportto(Schreyer, Syzextra::leadrawexp);
1417//      exportto(Schreyer, Syzextra::ISUpdateComponents);
1418      exportto(Schreyer, Syzextra::SetInducedReferrence);
1419      exportto(Schreyer, Syzextra::GetInducedData);
1420//      exportto(Schreyer, Syzextra::GetAMData);
1421//      exportto(Schreyer, Syzextra::SetSyzComp);
1422      exportto(Schreyer, Syzextra::MakeInducedSchreyerOrdering);
1423//      exportto(Schreyer, Syzextra::MakeSyzCompOrdering);
1424      exportto(Schreyer, Syzextra::idPrepare);
1425//      exportto(Schreyer, Syzextra::reduce_syz);
1426//      exportto(Schreyer, Syzextra::p_Content);
1427
1428      exportto(Schreyer, Syzextra::Tail);
1429    }
1430    else
1431    {
1432      if( @DEBUG )
1433      {
1434        "Loading the Debug version!";
1435      }
1436
1437      load("syzextra_g.so");
1438
1439      if( @DEBUG )
1440      {     
1441        listvar(Syzextra_g);
1442      }
1443     
1444      exportto(Top, Syzextra_g::ClearContent);
1445      exportto(Top, Syzextra_g::ClearDenominators);
1446
1447//      export Syzextra_g;
1448//      exportto(Schreyer, Syzextra_g::noop);
1449      exportto(Schreyer, Syzextra_g::DetailedPrint);
1450      exportto(Schreyer, Syzextra_g::leadmonomial);
1451      exportto(Schreyer, Syzextra_g::leadcomp);
1452//      exportto(Schreyer, Syzextra_g::leadrawexp);
1453//      exportto(Schreyer, Syzextra_g::ISUpdateComponents);
1454      exportto(Schreyer, Syzextra_g::SetInducedReferrence);
1455      exportto(Schreyer, Syzextra_g::GetInducedData);
1456//      exportto(Schreyer, Syzextra_g::GetAMData);
1457//      exportto(Schreyer, Syzextra_g::SetSyzComp);
1458      exportto(Schreyer, Syzextra_g::MakeInducedSchreyerOrdering);
1459//      exportto(Schreyer, Syzextra_g::MakeSyzCompOrdering);
1460      exportto(Schreyer, Syzextra_g::idPrepare);
1461//      exportto(Schreyer, Syzextra_g::reduce_syz);
1462//      exportto(Schreyer, Syzextra_g::p_Content);
1463
1464      exportto(Schreyer, Syzextra_g::Tail);
1465     
1466    }
1467
1468    exportto(Top, DetailedPrint);
1469    exportto(Top, GetInducedData);
1470
1471    if( @DEBUG )
1472    {
1473      listvar(Top);
1474      listvar(Schreyer);
1475    }
1476  }
1477 
1478  if( !defined(GetInducedData) )
1479  {
1480    ERROR("Sorry but we are missing the dynamic module (syzextra(_g)?.so)...");
1481  }
1482
1483}
1484
1485static proc mod_init()
1486{
1487  loadme();
1488}
1489
1490
1491proc testallSexamples()
1492{
1493  example Ssyz;
1494  example Scontinue;
1495  example Sres; 
1496}
1497
1498proc testallSSexamples()
1499{
1500  example SSsyz;
1501  example SScontinue;
1502  example SSres; 
1503}
1504
1505example
1506{ "EXAMPLE:"; echo = 2;
1507  testallSexamples();
1508  testallSSexamples();
1509}
Note: See TracBrowser for help on using the repository browser.