source: git/kernel/test.cc @ 576f5b

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