# source:git/Singular/LIB/schreyer.lib@33694c

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