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