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