source: git/dyn_modules/syzextra/noro.sing @ 495328

spielwiese
Last change on this file since 495328 was 0917a96, checked in by Oleksandr Motsak <motsak@…>, 12 years ago
adaptation of dyn_modules/syzextra for spielwiese fix: fixed dyn_modules/syzextra to work with spielwiese chg: use n_Init for long -> bigint conversion NOTE: the resulting bigint may be negative for (implicitly) casted unsigned longs (e.g. for raw poly/vector exponents) chg: further minor adaptations
  • Property mode set to 100644
File size: 12.3 KB
Line 
1echo = 2;
2
3LIB("syzextra.so");
4
5noop();
6
7system("--min-time", "0.01");
8system("--ticks-per-sec", 100);
9
10LIB "poly.lib"; // for numerator & denominator
11
12option(redSB); option(redTail); // assumed for the results of kStd!
13option(prot);
14option(mem);
15option(notWarnSB);
16
17
18// option(noloadLib); option(noredefine);
19
20LIB "teachstd.lib";
21
22int reduce_time = 0;
23
24///////////////////////////////////////////////////////////////////////////////
25proc MyReduce(f, G)
26{
27  int t = timer;
28  def g = reduce_syz(f, G, 0);
29  int tt = timer;
30
31  reduce_time = reduce_time + (tt - t);
32
33  return(g);
34};
35
36proc separateSyzGB( module J, int c )
37{
38  module II, G; vector v; int i;
39
40  J = simplify(J, 2);
41
42  for( i = ncols(J); i > 0; i-- )
43  {
44    v = J[i];
45    if( leadcomp(v) > c )
46    {
47      II[i] = v;
48    } else
49    {
50      G[i] = v; // leave only gen(i): i <= c
51    }
52  }
53
54  II = simplify(II, 2);
55  G = simplify(G, 2);
56
57  return (list(G, II));
58}
59
60
61
62proc splitSyzGB( module J, int c )
63{
64  module JJ; vector v, vv; int i;
65
66  for( i = ncols(J); i > 0; i-- )
67  {
68    v = J[i];
69
70    vv = 0;
71   
72    while( leadcomp(v) <= c )
73    {
74      vv = vv + lead(v);
75      v  = v  - lead(v);
76    }
77
78    J[i] = vv;
79    JJ[i] = v;
80  }
81
82  J = simplify(J, 2);
83  JJ = simplify(JJ, 2);
84
85  return (list(J, JJ));
86}
87
88
89
90proc prepareSyz( module I, list # )
91{
92  int i;
93  int k = 0;
94  int r = nrows(I);
95  int c = ncols(I);
96
97
98  if( size(#) > 0 )
99  {
100    if( typeof(#[1]) == "int" || typeof(#[1]) == "bigint" )
101    {
102      k = #[1];
103    }
104  }
105
106  if( k < r )
107  {
108    "// *** Wrong k: ", k, " < nrows: ", r, " => setting k = r = ", r;
109    k = r;
110  }
111
112//   "k: ", k;  "c: ", c;   "I: ", I;
113
114  for( i = c; i > 0; i-- )
115  {
116    I[i] = I[i] + gen(k + i);
117  }
118
119//  DetailedPrint(I);
120
121  return(I);
122}
123
124
125/// is p - permissible?
126proc myPrimeTest(def I, int p)
127{
128  int i; def v; number c;
129
130  for(i = ncols(I); i > 0; i--)
131  {
132    v = I[i];
133    while(v != 0)
134    {
135      c = leadcoef(v);
136     
137      if( (numerator(c) mod p) == 0) { return(0); }
138      if( (denominator(c) mod p) == 0) { return(0); }
139
140      v = v - lead(v);
141    }
142  }
143  return(1);
144}
145
146proc myPrimeList(module I, int n, list #)
147{
148  intvec L;
149  int i,p;
150
151  if(size(#) == 0)
152  {
153    p = 2147483647;
154
155    while(myPrimeTest(I,p)==0)
156    {
157      p = prime(p-1);
158      if(p == 2) { ERROR("no more primes"); }
159    }
160    L[1] = p;
161  }
162  else
163  {
164    L = #[1];
165    p = prime(L[size(L)]-1);
166    while(!myPrimeTest(I,p))
167    {
168      p = prime(p-1);
169      if(p == 2) { ERROR("no more primes"); }
170    }
171    L[size(L)+1] = p;
172  }
173
174  if(p == 2) { ERROR("no more primes"); }
175  for(i = 2; i <= n; i++)
176  {
177    p = prime(p-1);
178    while(!myPrimeTest(I,p))
179    {
180      p = prime(p-1);
181      if(p == 2) { ERROR("no more primes"); }
182    }
183    L[size(L)+1] = p;
184  }
185
186  return(L);
187}
188
189
190
191////////////////////////////////////////////////////////////////////////////////
192
193// int try = 0;
194
195proc GBCandidate0( module M, def P)
196{
197
198  def Q = basering;
199
200  int iComp = attrib(basering, "SyzComp");
201 
202  setring P;
203  module M = imap(Q, M);
204
205  M = idPrepare(M, iComp);
206
207/* 
208  try  ++;
209 
210  if( try != 5 )
211  {
212    return (0);
213  }
214*/
215 
216  M = simplify(M, 2);
217 
218  setring Q;
219
220  M = imap(P, M);
221 
222  return (M);
223}
224
225
226
227proc Noro( module M )
228{
229  def save = basering;
230
231
232  if( char(save) != 0 )
233  {
234    ERROR("Noro: char != 0!");
235  }
236
237  int iComp = nrows(M);
238
239 
240 def Q = MakeSyzCompOrdering(); setring Q;
241//  def Q = MakeInducedSchreyerOrdering(); setring Q;
242
243
244  SetSyzComp(iComp);
245  attrib(Q, "SyzComp", iComp);
246  DetailedPrint(basering);
247
248  list L = ringlist(basering);
249
250  def M = imap(save, M);
251
252  M = prepareSyz(M, iComp);
253
254  DetailedPrint(basering);
255
256  intvec pp = myPrimeList(M, 10);
257 
258//  intvec pp = 99981793+1;
259 
260//  intvec pp = 32003 + 1; pp = myPrimeList(M, 10, pp);
261 
262  int p;
263  int i = 1;
264
265  int try = 0;
266 
267  while(1)
268  {
269    try ++;
270
271    // choose p:
272    p = pp[i];
273
274    "Current prime: ", p;
275     
276    if( typeof(L[1]) == "int" || typeof(L[1]) == "bigint" )
277    {
278      L[1] = p;
279    }
280    else
281    {
282      if( typeof(L[1][1]) == "int" || typeof(L[1][1]) == "bigint" )
283      {
284        L[1][1] = p;
285      }
286      else
287      {
288        ERROR("NORO: cannot create p-ring list: wrong input ring?");
289      }
290    }
291
292    // new ring over F_p
293    def P = ring(L); setring P;
294
295   SetSyzComp(iComp);
296//    DetailedPrint(basering);
297
298    setring Q;
299
300//    M;    print(M);
301
302    // Compute GB Candidate for M using the ring P / F_p
303    "Computing GB Candidate: ";
304    def result = GBCandidate(M, P);
305
306    if( defined(result) && (size(result) > 0) && (typeof(result) == "list") && (result[1] == "ok") )
307    {
308      "===================================================================";
309      "It is DONE!";
310      "try: ", try;
311
312      result[2];
313      print(result[2]);
314//      DetailedPrint(result[2]);
315
316      list GB_SYZ = separateSyzGB(result[2], iComp);
317
318      module GB  = GB_SYZ[1]; //splitSyzGB(, iComp)[1];
319      module SYZ = GB_SYZ[2];
320
321
322      GB;
323      SYZ;
324
325     
326
327      setring save;
328
329     
330      module GB  = imap(Q, GB);
331
332      if( size(GB) > 0 )
333      {
334        // need the top part!
335        GB = transpose(GB);
336        GB = GB[1..iComp];
337        GB = transpose(GB);
338      }
339
340     
341      module SYZ = imap(Q, SYZ);
342
343      if( size(SYZ) > 0 )
344      {
345        int r = nrows(SYZ);
346
347        // need the bottom part:
348        SYZ = transpose(SYZ);
349        SYZ = SYZ[iComp+1 .. r];
350        SYZ = transpose(SYZ);
351      }
352     
353      return (list((GB), (SYZ)));
354    }
355
356    kill P;
357    setring Q;
358
359    i ++;
360
361    if( i > size(pp) )
362    {
363      pp = myPrimeList(M, 10, pp);
364      pp;
365    }
366
367       
368    if( !defined(pp) )
369    {
370      ERROR("NORO: Sorry no more primes!");
371      return(0);
372    }
373  } // while(1)
374}
375   
376
377// basering is assumed to be Z[x_1, ..., x_n]
378
379// F is a list of vectors of the same length
380proc GBCandidate (module F, def Fp)
381{
382  int ss = 1; // symmetric S-polynomial
383 
384  // TODO: check input for errors
385  def br = basering;
386
387  int i, j;
388
389  int sizeF = size(F);
390
391  intvec C;
392  list D; // list of intvecs of length 2, as indices of F
393  for(i = 1; i <= sizeF; i++)
394  {
395    for(j = i+1; j <= sizeF; j++)
396    {
397      D = insert(D, intvec(i, j));
398    }
399  }
400
401  vector h, s;
402  number c;
403 
404  setring Fp;
405  vector h, s;
406  module F;
407 
408  int p = char(Fp);
409  // TODO: check that characteristic of Fp is permissible for F
410
411 
412
413  while(size(D) > 0)
414  {
415   
416    if(char(basering) != p )
417    {
418      "ERROR: wrong current ring!";
419      $$$
420    }
421   
422    C = D[size(D)];
423    D = delete(D, size(D));
424
425//    "Pair: (", C[1], ",", C[2], ")";
426   
427    F = fetch(br, F);
428    s = spoly(F[C[1]], F[C[2]], ss);
429    h = MyReduce(s, F);
430   
431    if(h != 0)
432    {
433      setring br;
434
435      s = spoly(F[C[1]], F[C[2]], ss);
436      h = MyReduce(s, F); // h \in Z_<p>[x] ???
437
438      c = leadcoef(h);
439
440      if(h != 0)
441      {
442        if( (denominator(c) mod p) == 0 )
443        {
444          "h: ", h;
445          "c: ", c;
446          "ERROR: GBCandidate: wrong coeff!";
447          $$
448        }
449      }
450     
451      if(h != 0 && (numerator(c) mod p) != 0)
452      {
453        for(i = 1; i <= sizeF; i++)
454        {
455          D = insert(D, intvec(i, sizeF+1));
456        }
457        sizeF++;
458        F[sizeF] = h;
459
460//        "new F[",sizeF,"] element: ", h;
461       
462        setring Fp;
463       
464        h = fetch(br, h);
465        F[sizeF] = h;
466      }
467      else
468      {
469        /*
470        "&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&";
471        matrix T; //  = lift(F, s);
472        basering; // /Q
473       
474        test(23+32);
475        option(teach);
476       
477        "i, j: ", C[1], C[2];
478       
479
480        F;
481        "f: ", F[C[1]];
482        "g: ", F[C[2]];
483        "s: ", s;
484        "ZERO::: h: ", h; // 0!!!
485
486        vector m = NFMora(s, F);
487        "MORA: ", m;
488
489//        matrix(s)-matrix(F)*T; print(T);
490
491        MyReduce(s, F);
492        "h-test: ", _ == h;
493       
494        "????";
495
496        setring Fp;
497
498        "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!";
499        basering; // /p
500       
501
502        F;
503        "F-proj-test: ", size(module(matrix(F) - matrix(fetch(br, F)))) == 0;
504       
505        "f: ", F[C[1]];
506        "g: ", F[C[2]];
507        "s: ", s;
508
509        "s-proj-test: ", s == fetch(br, s);
510       
511        "NON-ZERO::: h: ", h;
512
513//        matrix T = fetch(br, T);        def S = (matrix(F)*T);        print(S);        S;
514       
515        vector m = NFMora(s, F);
516        "MORA: ", m;
517
518        "m-proj-test: ", m == fetch(br, m);
519       
520       
521        MyReduce(s, F);
522        "h-test: ", _ == h;
523
524        "????";
525        $$$
526        */   
527        return(list("ng", module()));
528      }
529    }
530  }
531  setring br;
532  F = interred(F); // ?
533  return(list("ok", F));
534}
535
536
537
538
539echo = 1;
540
541
542ring R = 0, (w, x, y, z), dp;
543module I = w^2 - x*z,  w*x - y*z,  x^2 - w*y, x*y - z^2, y^2 - w*z;
544
545
546
547/*
548ring r = 0,(x,y,z),(c,dp);
549module I = [x+1, y, 1], [xy, z, z2]; // NO SYZ!
550*/
551
552
553/*
554// Test: kotsireas
555ring R = (0),(B, b, D, d, F, f),dp;
556module I = B*b-b*D-B*d+D*d-2*F+2,B*b+b*D-B*d-D*d-2*b*F+2*d*F-2*B+2*D,b^2-2*b*d+d^2-2*b-2*d+f+1,B^2*b^3-1,D^2*d^3-1,F^2*f^3-1;
557*/
558
559
560
561/*
562// Cohn3
563ring R = (0),(x1, y, z, t),dp;
564module I = -x1^3*y^2+2*x1^2*y^2*z-x1^2*y*z^2-144*x1^2*y^2-207*x1^2*y*z+288*x1*y^2*z+78*x1*y*z^2+x1*z^3-3456*x1^2*y-5184*x1*y^2-9504*x1*y*z-432*x1*z^2-248832*x1*y+62208*x1*z-2985984*x1,-x1^3*z*t^2+x1^2*z^2*t^2-6*x1^3*z*t+4*x1^2*z^2*t+32*x1^3*t^2-72*x1^2*z*t^2-87*x1*z^2*t^2-z^3*t^2-8*x1^3*z-432*x1^2*z*t-414*x1*z^2*t+2592*x1*z*t^2+864*z^2*t^2-1728*x1^2*z-20736*x1*z*t+3456*z^2*t-186624*z*t^2-124416*x1*z-1492992*z*t-2985984*z,x1^2*y*t^3-2*x1*y^2*t^3+y^3*t^3+8*x1^2*y*t^2-12*x1*y^2*t^2+4*y^3*t^2-24*x1*y*t^3+24*y^2*t^3+20*x1^2*y*t-20*x1*y^2*t-160*x1*y*t^2+96*y^2*t^2+128*x1*t^3+16*x1^2*y+96*x1*y*t+2304*x1*t^2+1152*x1*y+13824*x1*t+27648*x1,y^3*t^3-y^2*z*t^3+4*y^3*t^2-2*y^2*z*t^2+72*y^2*t^3+71*y*z*t^3+288*y^2*t^2+360*y*z*t^2+6*z^2*t^2+1728*y*t^3-464*z*t^3+432*y*z*t+8*z^2*t+6912*y*t^2-4320*z*t^2+13824*t^3+z^2-13824*z*t+55296*t^2-13824*z;
565*/
566
567
568
569/*
570// Cyclic 7: (done)
571ring R = (0),(a, b, c, d, e),dp;
572module I = a+b+c+d+e,a*b+a*e+b*c+c*d+d*e,a*b*c+a*b*e+a*d*e+b*c*d+c*d*e,a*b*c*d+a*b*c*e+a*b*d*e+a*c*d*e+b*c*d*e,a*b*c*d*e-1;
573*/
574
575/*
576// HCyclic 7: (running: 80% of 16GB mem while computing syz!... too slow... )
577ring R = (0),(x1, x2, x3, x4, x5, x6, x7, w),dp;
578module I = x1+x2+x3+x4+x5+x6+x7,x1*x2+x2*x3+x3*x4+x4*x5+x5*x6+x1*x7+x6*x7,x1*x2*x3+x2*x3*x4+x3*x4*x5+x4*x5*x6+x1*x2*x7+x1*x6*x7+x5*x6*x7,x1*x2*x3*x4+x2*x3*x4*x5+x3*x4*x5*x6+x1*x2*x3*x7+x1*x2*x6*x7+x1*x5*x6*x7+x4*x5*x6*x7,x1*x2*x3*x4*x5+x2*x3*x4*x5*x6+x1*x2*x3*x4*x7+x1*x2*x3*x6*x7+x1*x2*x5*x6*x7+x1*x4*x5*x6*x7+x3*x4*x5*x6*x7,x1*x2*x3*x4*x5*x6+x1*x2*x3*x4*x5*x7+x1*x2*x3*x4*x6*x7+x1*x2*x3*x5*x6*x7+x1*x2*x4*x5*x6*x7+x1*x3*x4*x5*x6*x7+x2*x3*x4*x5*x6*x7,x1*x2*x3*x4*x5*x6*x7-w^7;
579*/
580
581// syz(I); print(_); //
582
583/*
584module M;
585M[5]=abcde*gen(1)-gen(1)+gen(6);
586M[1]=a*gen(1)+b*gen(1)+c*gen(1)+d*gen(1)+e*gen(1)+gen(2);
587M[2]=ab*gen(1)+bc*gen(1)+cd*gen(1)+ae*gen(1)+de*gen(1)+gen(3);
588M[3]=abc*gen(1)+bcd*gen(1)+abe*gen(1)+ade*gen(1)+cde*gen(1)+gen(4);
589M[4]=abcd*gen(1)+abce*gen(1)+abde*gen(1)+acde*gen(1)+bcde*gen(1)+gen(5);
590
591ring P = (99981794),(a, b, c, d, e),dp;
592setring R;
593
594GBCandidate(M, P);$$
595*/
596
597int tG = timer;
598def G = groebner(I);
599int ttG = timer;
600
601"Time GB: ", ttG - tG;
602
603
604int tN = timer;
605def L = Noro( I );
606int ttN = timer;
607
608"Time Noro: ", ttN - tN;
609
610// All syzygies?
611int tS = timer;
612module S = syz(I);
613int ttS = timer;
614
615"Time Syz: ", ttS - tS;
616
617
618module GB = L[1];
619module SYZ = L[2];
620
621":::::::::::::::::::::::::::::::::::::::::::: GB :::::::::::::::::::::::::::::::::::: ";
622print(GB);
623
624// test GB:
625if( size(NF(I, GB)) > 0 )
626{
627  ERROR("NORO was wrong: I is bigger than GB!");
628}
629
630if( size(NF(GB, G)) > 0 )
631{
632  ERROR("NORO was wrong: GB is bigger than I!");
633}
634
635":::::::::::::::::::::::::::::::::::::::::::: SYZ :::::::::::::::::::::::::::::::::::: ";
636print(SYZ);
637
638
639if( size(SYZ) > 0 )
640{
641  // test syzygy
642  if( size( module(transpose(SYZ)*transpose(I)) ) > 0 )
643  {
644    ERROR("NORO was wrong: SYZ are NOT syzygies of I!");
645  }
646}
647
648
649
650// test SYZ:
651if( size(NF(SYZ, groebner(S))) > 0 )
652{
653  ERROR("NORO was wrong: too much syzygies found!!!");
654}
655
656if( size(NF(S, groebner(SYZ))) > 0 )
657{
658  ERROR("NORO was wrong: too few syzygies found!!!");
659}
660
661
662":::::::::::::::::::::::::::::::::::::::::::: GOOD :::::::::::::::::::::::::::::::::::: ";
663
664
665"Time GB: ", ttG - tG;
666"Time Syz: ", ttS - tS;
667"Time Noro: ", ttN - tN, " vs GB+SYZ: ", (ttS - tS) + (ttG - tG), "";
668
669"Time NF: ", reduce_time;
670
671
672
673number t = (ttN - tN);
674number tt = ((ttS - tS) + (ttG - tG));
675"Factor: ", t div tt;
676"Factor: ", (ttN - tN) div ((ttS - tS) + (ttG - tG));
677
678
679$$
Note: See TracBrowser for help on using the repository browser.