source: git/kernel/test.cc @ 8d1432e

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