source: git/kernel/test.cc @ 5116a21

spielwiese
Last change on this file since 5116a21 was 5116a21, checked in by Oleksandr Motsak <motsak@…>, 12 years ago
test syzygy and all kinds of resolutions via make check ADD: test syz + sres/lres/nres/mres on the kernel level! ADD: print syStrategy ADD: get length + Betti numbers ADD: minimize resolutions TODO: conversion to list of ideals?
  • Property mode set to 100644
File size: 12.7 KB
Line 
1#include "mod2.h"
2
3#include <omalloc/omalloc.h>
4#include <misc/auxiliary.h>
5#ifdef HAVE_FACTORY
6#include <factory/factory.h> // :(
7#endif
8
9#include <misc/intvec.h>
10#include <misc/int64vec.h>
11#include <misc/mylimits.h>
12#include <misc/options.h>
13
14#include <reporter/reporter.h>
15
16#include <resources/feFopen.h>
17#include <resources/feResource.h>
18
19#include <coeffs/coeffs.h>
20
21#include <coeffs/si_gmp.h>
22
23#include <polys/kbuckets.h>
24#include <polys/matpol.h>
25#include <polys/mod_raw.h>
26#include <polys/prCopy.h>
27#include <polys/sbuckets.h>
28#include <polys/simpleideals.h>
29#include <polys/weight.h>
30
31#include <polys/monomials/maps.h>
32#include <polys/monomials/monomials.h>
33#include <polys/monomials/p_polys.h>
34#include <polys/monomials/ring.h>
35
36#include <polys/nc/nc.h>
37#include <polys/nc/ncSACache.h>
38#include <polys/nc/ncSAFormula.h>
39#include <polys/nc/ncSAMult.h>
40#include <polys/nc/sca.h>
41#include <polys/nc/summator.h>
42
43
44#include <polys/templates/p_MemAdd.h>
45#include <polys/templates/p_Procs.h>
46
47#include <polys/operations/pShallowCopyDelete.h>
48
49#include <polys/clapsing.h>
50
51
52#ifdef HAVE_FACTORY
53// The following are needed due to FACTORY (e.g. initCanonicalForm)
54int initializeGMP(){ return 1; } 
55int mmInit(void) {return 1; }
56#endif
57
58// // TODO: DUE to the use of HALT in npolygon.cc :(((
59extern "C" {void m2_end(int i){exit(i);}}
60
61// // TODO: DUE to its use in kutil.cc :(((
62char * showOption(){return NULL;}
63
64// // TODO: DUE to its use in feread.cc :(((
65char *iiArithGetCmd(int nPos){return NULL; }
66
67
68#include <coeffs/numbers.h>
69
70#include "mod2.h"
71#include "structs.h"
72
73
74// HEADERS:
75#include "hutil.h"
76//#include "idrec.h" // moved to Singular
77#include "stairc.h"
78#include "ideals.h"
79#include "syz.h"
80#include "fast_maps.h"
81#include "febase.h"
82#include "walkProc.h"
83#include "walkMain.h"
84#include "walkSupport.h"
85#include "khstd.h"
86/// #include "sparsmat.h" // TODO: install polys/this!
87//+
88
89#include "fglm.h"
90#include "kstd1.h"
91#include "fglmgauss.h"
92#include "fglmvec.h"
93#include "kstdfac.h"
94#include "kmatrix.h"
95#include "GMPrat.h"
96#include "multicnt.h"
97#include "npolygon.h"
98#include "semic.h"
99#include "spectrum.h"
100#include "splist.h"
101#include "multicnt.h"
102#include "eigenval.h"
103#include "units.h"
104#include "fegetopt.h"
105#include "ratgring.h"
106#include "shiftgb.h"
107// #include "mmalloc.h" // move to Singular!?
108#include "gfan.h"
109
110#include "kutil.h"
111
112// #include "dbm_sl.h" // TODO: needs si_link// already moved to Singular/!
113
114// #include "CCRing.h" // Too old!
115#include "digitech.h"
116#include "eigenval.h"
117#include "fast_maps.h"
118#include "fast_mult.h"
119#include "febase.h"
120#include "fegetopt.h"
121
122#include "fglmgauss.h"
123#include "fglm.h"
124#include "fglmvec.h"
125
126////////#include "F5cData.h"
127#include "f5c.h"
128#include "f5data.h"
129#include "f5gb.h"
130#include "f5lists.h"
131////////#include "F5cLists.h"
132
133#include "gfan.h"
134
135
136#include "GMPrat.h"
137
138// #include "htmlhelp.h" // For Windows //
139#include "hutil.h"
140// #include "Ideal.h" // Too old?
141
142
143#include "ideals.h"
144
145#include "kmatrix.h"
146#include "kstd1.h"
147#include "kstdfac.h"
148#include "khstd.h"
149
150#include "linearAlgebra.h"
151
152
153
154// #include "lplist.h" // Too old!
155#include "multicnt.h"
156#include "npolygon.h"
157// #include "Number.h" // Too old?
158#include "page.h"
159// #include "Poly.h" // Too old?
160// #include "PowerSeries.h" // Too old?
161
162#include "preimage.h"
163
164#include "nc.h"
165
166#include "ratgring.h"
167#include "ringgb.h"
168// #include "run.h" // moved to Singular
169#include "semic.h"
170#include "shiftgb.h"
171// #include "si_gmp.h" // remove it...
172#include "spectrum.h"
173#include "splist.h"
174#include "stairc.h"
175#include "structs.h"
176#include "syz.h"
177// #include "testpoly.h" // Too old?
178
179#include "tgbgauss.h"
180#include "tgb.h"
181
182#include "timer.h"
183
184#include "units.h"
185#include "walkMain.h"
186#include "walkProc.h"
187#include "walkSupport.h"
188
189
190// #include <polys/clapconv.h> // due to factory? :(
191// #include "tgb_internal.h" // :(
192// #include "F4.h" // uses tgb_internal // :(
193// #include "IIntvec.h" // :(
194
195
196
197
198// #include "fglmzero.cc" // looks like <factory/templates/ftmpl_list.h> must be installed!
199// TODO: looks like <coeffs/mpr_complex.h> must be installed!
200
201
202
203#include "polys.h"
204
205void TestGBEngine()
206{
207   
208  //  R = MPolynomialRing_polydict(QQ,5,'w,x,y,z,C', order='degrevlex')
209  //  J = (w*w - x*z, w*x - y*z, x*x - w*y, x*y - z*z, y*y - w*z)
210         
211  const short w = 1;
212  const short x = 2;
213  const short y = 3;
214  const short z = 4;
215
216  const short N = (z - w + 1);
217
218  char **n=(char**)omalloc(N*sizeof(char*));
219
220 
221  n[w-1]=omStrDup("w");
222  n[x-1]=omStrDup("x");
223  n[y-1]=omStrDup("y");
224  n[z-1]=omStrDup("z");
225
226
227  const int D = 3;
228  int *order = (int *) omAlloc0(D* sizeof(int));
229  int *block0 = (int *)omAlloc0(D * sizeof(int));
230  int *block1 = (int *)omAlloc0(D * sizeof(int));
231
232  order[0]  = ringorder_dp;
233  block0[0] = 1;
234  block1[0] = N;
235
236  order[1]  = ringorder_C;
237  block0[1] = 1;
238  block1[1] = N;
239
240  ring R = rDefault(0, N, n, D, order, block0, block1); 
241 
242//   ring R = rDefault(0, N, n);
243
244  rWrite(R); PrintLn();
245
246#ifdef RDEBUG
247  rDebugPrint(R);
248#endif
249
250  ideal I = idInit(5, 1);
251
252  int gen = 0;
253
254  {
255    // -xz
256    poly p = p_ISet(-1,R);
257   
258    p_SetExp(p, x, 1, R);
259    p_SetExp(p, z, 1, R);
260    p_Setm(p, R);
261
262    assume( p_GetExp(p, x, R) == 1 );
263    assume( p_GetExp(p, z, R) == 1 );
264    assume( p_GetExp(p, w, R) == 0 );
265    assume( p_GetExp(p, y, R) == 0 );
266
267    // +w2
268    poly lp = p_ISet(1,R);
269    p_SetExp(lp, w, 2, R);
270    p_Setm(lp, R);
271
272    assume( p_GetExp(lp, w, R) == 2 );
273    assume( p_GetExp(lp, x, R) == 0 );
274    assume( p_GetExp(lp, y, R) == 0 );
275    assume( p_GetExp(lp, z, R) == 0 );
276
277    MATELEM(I, 1, ++gen) = p_Add_q(lp, p, R); // w2 - xz
278  }
279
280  {
281    // -yz
282    poly p = p_ISet(-1,R);
283
284    p_SetExp(p, y, 1, R);
285    p_SetExp(p, z, 1, R);
286    p_Setm(p, R);
287
288    assume( p_GetExp(p, y, R) == 1 );
289    assume( p_GetExp(p, z, R) == 1 );
290    assume( p_GetExp(p, w, R) == 0 );
291    assume( p_GetExp(p, x, R) == 0 );
292
293    // +wx
294    poly lp = p_ISet(1,R);
295    p_SetExp(lp, w, 1, R);
296    p_SetExp(lp, x, 1, R);
297    p_Setm(lp, R);
298
299    assume( p_GetExp(lp, w, R) == 1 );
300    assume( p_GetExp(lp, x, R) == 1 );
301    assume( p_GetExp(lp, y, R) == 0 );
302    assume( p_GetExp(lp, z, R) == 0 );
303
304    MATELEM(I, 1, ++gen) = p_Add_q(lp, p, R); // wx - yz
305  }
306 
307
308  {
309    // -wy
310    poly p = p_ISet(-1,R);
311
312    p_SetExp(p, y, 1, R);
313    p_SetExp(p, w, 1, R);
314    p_Setm(p, R);
315
316    assume( p_GetExp(p, y, R) == 1 );
317    assume( p_GetExp(p, w, R) == 1 );
318    assume( p_GetExp(p, z, R) == 0 );
319    assume( p_GetExp(p, x, R) == 0 );
320
321    // +x2
322    poly lp = p_ISet(1,R);
323    p_SetExp(lp, x, 2, R);
324    p_Setm(lp, R);
325
326    assume( p_GetExp(lp, w, R) == 0 );
327    assume( p_GetExp(lp, x, R) == 2 );
328    assume( p_GetExp(lp, y, R) == 0 );
329    assume( p_GetExp(lp, z, R) == 0 );
330
331    MATELEM(I, 1, ++gen) = p_Add_q(lp, p, R); // x2 - wy
332  }
333 
334
335  {
336    // -z2
337    poly p = p_ISet(-1,R);
338
339    p_SetExp(p, z, 2, R);
340    p_Setm(p, R);
341
342    assume( p_GetExp(p, y, R) == 0 );
343    assume( p_GetExp(p, w, R) == 0 );
344    assume( p_GetExp(p, z, R) == 2 );
345    assume( p_GetExp(p, x, R) == 0 );
346
347    // +xy
348    poly lp = p_ISet(1,R);
349    p_SetExp(lp, x, 1, R);
350    p_SetExp(lp, y, 1, R);
351    p_Setm(lp, R);
352
353    assume( p_GetExp(lp, w, R) == 0 );
354    assume( p_GetExp(lp, x, R) == 1 );
355    assume( p_GetExp(lp, y, R) == 1 );
356    assume( p_GetExp(lp, z, R) == 0 );
357
358    MATELEM(I, 1, ++gen) = p_Add_q(lp, p, R); // xy - z2
359  }
360 
361
362  {
363    // -wz
364    poly p = p_ISet(-1,R);
365
366    p_SetExp(p, w, 1, R);
367    p_SetExp(p, z, 1, R);
368    p_Setm(p, R);
369
370    assume( p_GetExp(p, y, R) == 0 );
371    assume( p_GetExp(p, w, R) == 1 );
372    assume( p_GetExp(p, z, R) == 1 );
373    assume( p_GetExp(p, x, R) == 0 );
374
375    // +y2
376    poly lp = p_ISet(1,R);
377    p_SetExp(lp, y, 2, R);
378    p_Setm(lp, R);
379
380    assume( p_GetExp(lp, w, R) == 0 );
381    assume( p_GetExp(lp, x, R) == 0 );
382    assume( p_GetExp(lp, y, R) == 2 );
383    assume( p_GetExp(lp, z, R) == 0 );
384
385    MATELEM(I, 1, ++gen) = p_Add_q(lp, p, R); // y2 - wz
386  }
387#ifdef PDEBUG
388  PrintS("I: ");
389  idShow(I, R, R, 0);
390#endif
391
392
393//  ideal kStd(ideal F, ideal Q, tHomog h, intvec ** mw,intvec *hilb=NULL,
394//             int syzComp=0,int newIdeal=0, intvec *vw=NULL);
395  // make R the default ring:
396  rChangeCurrRing(R);
397
398  {
399    ideal G = kStd(I, currQuotient, testHomog, NULL);
400
401#ifdef PDEBUG
402    PrintS("GB: ");
403    idShow(G, R, R, 0);
404#endif
405
406    idDelete( &G, R);
407  }
408
409  {
410    intvec *weights = NULL;
411    ideal SYZ = idSyzygies(I, testHomog, &weights);
412   
413#ifdef PDEBUG
414    PrintS("SYZ: ");
415    idShow(SYZ, R, R, 0);
416#endif
417
418    idDelete(&SYZ, R);
419    if (weights!=NULL) { PrintS("weights: "); weights->show(); delete weights; }
420  }
421
422
423  {
424    PrintS("\n**********************************\n");
425    PrintS("lres: \n");
426    int dummy;
427    syStrategy r = syLaScala3(I,&dummy);
428
429    intvec *b = syBettiOfComputation(r, FALSE);
430    PrintS("non-min. betti: \n");    b->show();    PrintLn();
431    delete b;
432
433    Print("length: %d\n", sySize(r));
434
435    syPrint(r, "R");
436
437    r =  syMinimize(r); // syzstr->references ++ ==> memory leak :(((
438
439    b = syBettiOfComputation(r, TRUE);
440    PrintS("min. betti: \n");    b->show();    PrintLn();
441    delete b;   
442
443    Print("length: %d\n", sySize(r));
444
445    syPrint(r, "R");
446   
447    syKillComputation(r, R);
448  }
449
450  {
451    PrintS("\n**********************************\n");
452    PrintS("sres: \n");
453    const int maxl = rVar(R)-1; // +2*(1);
454
455    syStrategy r = sySchreyer(I, rVar(R));
456
457    intvec *b = syBettiOfComputation(r, FALSE);
458    PrintS("non-min. betti: \n");    b->show();    PrintLn();
459    delete b;
460
461    Print("length: %d\n", sySize(r));
462
463    syPrint(r, "R");
464
465    r =  syMinimize(r); // syzstr->references ++ ==> memory leak :(((
466
467    b = syBettiOfComputation(r, TRUE);
468    PrintS("min. betti: \n");    b->show();    PrintLn();
469    delete b;   
470
471    Print("length: %d\n", sySize(r));
472
473    syPrint(r, "R");
474
475    syKillComputation(r, R);
476  }
477
478   
479
480  {
481    PrintS("\n**********************************\n");
482    PrintS("nres: \n");
483    intvec *weights=NULL;
484//    const int maxl = rVar(R)-1 + 2*(1);
485    syStrategy r = syResolution(I, rVar(R)-1, weights, FALSE/*iiOp==MRES_CMD*/);
486
487    intvec *b = syBettiOfComputation(r, FALSE);
488    PrintS("non-min. betti: \n");    b->show();    PrintLn();
489    delete b;
490
491    Print("length: %d\n", sySize(r));
492
493    syPrint(r, "R");
494
495    r =  syMinimize(r); // syzstr->references ++ ==> memory leak :(((
496
497    b = syBettiOfComputation(r, TRUE);
498    PrintS("min. betti: \n");    b->show();    PrintLn();
499    delete b;   
500
501    Print("length: %d\n", sySize(r));
502
503    syPrint(r, "R");
504
505    syKillComputation(r, R);
506  }
507 
508
509  {
510    PrintS("\n**********************************\n");
511    PrintS("mres: \n");
512    intvec *weights=NULL;
513//    const int maxl = rVar(R)-1 + 2*(1);
514    syStrategy r = syResolution(I, rVar(R)+1, weights, TRUE/*iiOp==MRES_CMD*/);
515
516    intvec *b = syBettiOfComputation(r, FALSE);
517    PrintS("non-min. betti: \n");    b->show();    PrintLn();
518    delete b;
519
520    Print("length: %d\n", sySize(r));
521
522    syPrint(r, "R");
523
524    r =  syMinimize(r); // syzstr->references ++ ==> memory leak :(((
525
526    b = syBettiOfComputation(r, TRUE);
527    PrintS("min. betti: \n");    b->show();    PrintLn();
528    delete b;   
529
530    Print("length: %d\n", sySize(r));
531
532    syPrint(r, "R");
533
534    syKillComputation(r, R);
535  }
536 
537
538
539 
540  idDelete( &I, R); 
541  rDelete(R); // should cleanup every belonging polynomial, right!?
542   
543}
544
545
546
547void TestSimpleRingArithmetcs()
548{
549  // Libpolys tests:
550
551  // construct the ring Z/32003[x,y,z]
552  // the variable names
553  char **n=(char**)omalloc(3*sizeof(char*));
554  n[0]=omStrDup("x");
555  n[1]=omStrDup("y");
556  n[2]=omStrDup("z2");
557
558  ring R = rDefault(32003,3,n); //  ring R = rDefault(0,3,n);
559
560  rWrite(R); PrintLn();
561
562#ifdef RDEBUG
563  rDebugPrint(R);
564#endif
565
566
567  poly p = p_ISet(1,R); p_SetExp(p,1,1, R); p_Setm(p, R);
568
569  assume( p_GetExp(p,1, R) == 1 );
570
571  poly pp = pp_Mult_qq( p, p, R);
572
573  Print("p: "); p_Write0(p, R); Print(", deg(p): %d", p_Totaldegree(p, R)); assume( 1 == p_Totaldegree(p, R) );
574
575  Print("; p*p : "); p_Write0(pp, R); Print("deg(pp): %d\n", p_Totaldegree(pp, R)); assume( 2 == p_Totaldegree(pp, R) );
576
577
578  p_Delete(&p, R);
579
580  assume( p_GetExp(pp,1, R) == 2 );
581
582  p_Delete(&pp, R);
583
584
585//  rDelete(R);
586
587  // make R the default ring:
588  rChangeCurrRing(R);
589
590  // create the polynomial 1
591  poly p1=pISet(1);
592
593  // create tthe polynomial 2*x^3*z^2
594  poly p2=p_ISet(2,R);
595  pSetExp(p2,1,3);
596  pSetExp(p2,3,2);
597  pSetm(p2);
598
599  // print p1 + p2
600  Print("p1: "); pWrite0(p1); 
601  Print(" + p2: "); pWrite0(p2); 
602  Print("  ---- >>>> ");
603
604  // compute p1+p2
605  p1=p_Add_q(p1,p2,R); p2=NULL;
606  pWrite(p1); 
607
608  // clean up:
609//  pDelete(&p1);
610
611  rDelete(R); // should cleanup every belonging polynomial, right!?
612}
613
614
615int main( int, char *argv[] ) 
616{
617  feInitResources(argv[0]);
618
619  StringSetS("ressources in use (as reported by feStringAppendResources(0):\n");
620  feStringAppendResources(0);
621  PrintS(StringAppendS("\n"));
622
623  TestGBEngine();
624  TestSimpleRingArithmetcs();
625   
626  return 0;
627}
Note: See TracBrowser for help on using the repository browser.