My Project
hdegree.cc
Go to the documentation of this file.
1 /****************************************
2 * Computer Algebra System SINGULAR *
3 ****************************************/
4 /*
5 * ABSTRACT - dimension, multiplicity, HC, kbase
6 */
7 
8 #include "kernel/mod2.h"
9 
10 #include "misc/intvec.h"
11 #include "coeffs/numbers.h"
12 
13 #include "kernel/structs.h"
14 #include "kernel/ideals.h"
15 #include "kernel/polys.h"
16 
20 #include "reporter/reporter.h"
21 
22 #ifdef HAVE_SHIFTBBA
23 #include <vector>
24 #include "misc/options.h"
25 #endif
26 
27 VAR int hCo, hMu2;
28 VAR long hMu;
29 VAR omBin indlist_bin = omGetSpecBin(sizeof(indlist));
30 
31 /*0 implementation*/
32 
33 // dimension
34 
35 void hDimSolve(scmon pure, int Npure, scfmon rad, int Nrad,
36  varset var, int Nvar)
37 {
38  int dn, iv, rad0, b, c, x;
39  scmon pn;
40  scfmon rn;
41  if (Nrad < 2)
42  {
43  dn = Npure + Nrad;
44  if (dn < hCo)
45  hCo = dn;
46  return;
47  }
48  if (Npure+1 >= hCo)
49  return;
50  iv = Nvar;
51  while(pure[var[iv]]) iv--;
52  hStepR(rad, Nrad, var, iv, &rad0);
53  if (rad0!=0)
54  {
55  iv--;
56  if (rad0 < Nrad)
57  {
58  pn = hGetpure(pure);
59  rn = hGetmem(Nrad, rad, radmem[iv]);
60  hDimSolve(pn, Npure + 1, rn, rad0, var, iv);
61  b = rad0;
62  c = Nrad;
63  hElimR(rn, &rad0, b, c, var, iv);
64  hPure(rn, b, &c, var, iv, pn, &x);
65  hLex2R(rn, rad0, b, c, var, iv, hwork);
66  rad0 += (c - b);
67  hDimSolve(pn, Npure + x, rn, rad0, var, iv);
68  }
69  else
70  {
71  hDimSolve(pure, Npure, rad, Nrad, var, iv);
72  }
73  }
74  else
75  hCo = Npure + 1;
76 }
77 
78 int scDimInt(ideal S, ideal Q)
79 {
80  id_Test(S, currRing);
81  if( Q!=NULL ) id_Test(Q, currRing);
82 
83  int mc;
84  hexist = hInit(S, Q, &hNexist);
85  if (!hNexist)
86  return (currRing->N);
87  hwork = (scfmon)omAlloc(hNexist * sizeof(scmon));
88  hvar = (varset)omAlloc(((currRing->N) + 1) * sizeof(int));
89  hpure = (scmon)omAlloc((1 + ((currRing->N) * (currRing->N))) * sizeof(int));
90  mc = hisModule;
91  if (!mc)
92  {
93  hrad = hexist;
94  hNrad = hNexist;
95  }
96  else
97  hrad = (scfmon)omAlloc(hNexist * sizeof(scmon));
98  radmem = hCreate((currRing->N) - 1);
99  hCo = (currRing->N) + 1;
100  loop
101  {
102  if (mc)
103  hComp(hexist, hNexist, mc, hrad, &hNrad);
104  if (hNrad)
105  {
106  hNvar = (currRing->N);
107  hRadical(hrad, &hNrad, hNvar);
108  hSupp(hrad, hNrad, hvar, &hNvar);
109  if (hNvar)
110  {
111  memset(hpure, 0, ((currRing->N) + 1) * sizeof(int));
112  hPure(hrad, 0, &hNrad, hvar, hNvar, hpure, &hNpure);
113  hLexR(hrad, hNrad, hvar, hNvar);
115  }
116  }
117  else
118  {
119  hCo = 0;
120  break;
121  }
122  mc--;
123  if (mc <= 0)
124  break;
125  }
126  hKill(radmem, (currRing->N) - 1);
127  omFreeSize((ADDRESS)hpure, (1 + ((currRing->N) * (currRing->N))) * sizeof(int));
128  omFreeSize((ADDRESS)hvar, ((currRing->N) + 1) * sizeof(int));
129  omFreeSize((ADDRESS)hwork, hNexist * sizeof(scmon));
131  if (hisModule)
132  omFreeSize((ADDRESS)hrad, hNexist * sizeof(scmon));
133  return (currRing->N) - hCo;
134 }
135 
136 int scDimIntRing(ideal vid, ideal Q)
137 {
138 #ifdef HAVE_RINGS
140  {
141  int i = idPosConstant(vid);
142  if ((i != -1) && (n_IsUnit(pGetCoeff(vid->m[i]),currRing->cf)))
143  { /* ideal v contains unit; dim = -1 */
144  return(-1);
145  }
146  ideal vv = id_Head(vid,currRing);
147  idSkipZeroes(vv);
148  i = idPosConstant(vid);
149  int d;
150  if(i == -1)
151  {
152  d = scDimInt(vv, Q);
153  if(rField_is_Z(currRing))
154  d++;
155  }
156  else
157  {
158  if(n_IsUnit(pGetCoeff(vv->m[i]),currRing->cf))
159  d = -1;
160  else
161  d = scDimInt(vv, Q);
162  }
163  //Anne's Idea for std(4,2x) = 0 bug
164  int dcurr = d;
165  for(unsigned ii=0;ii<(unsigned)IDELEMS(vv);ii++)
166  {
167  if(vv->m[ii] != NULL && !n_IsUnit(pGetCoeff(vv->m[ii]),currRing->cf))
168  {
169  ideal vc = idCopy(vv);
170  poly c = pInit();
171  pSetCoeff0(c,nCopy(pGetCoeff(vv->m[ii])));
172  idInsertPoly(vc,c);
173  idSkipZeroes(vc);
174  for(unsigned jj = 0;jj<(unsigned)IDELEMS(vc)-1;jj++)
175  {
176  if((vc->m[jj]!=NULL)
177  && (n_DivBy(pGetCoeff(vc->m[jj]),pGetCoeff(c),currRing->cf)))
178  {
179  pDelete(&vc->m[jj]);
180  }
181  }
182  idSkipZeroes(vc);
183  i = idPosConstant(vc);
184  if (i != -1) pDelete(&vc->m[i]);
185  dcurr = scDimInt(vc, Q);
186  // the following assumes the ground rings to be either zero- or one-dimensional
187  if((i==-1) && rField_is_Z(currRing))
188  {
189  // should also be activated for other euclidean domains as groundfield
190  dcurr++;
191  }
192  idDelete(&vc);
193  }
194  if(dcurr > d)
195  d = dcurr;
196  }
197  idDelete(&vv);
198  return d;
199  }
200 #endif
201  return scDimInt(vid,Q);
202 }
203 
204 // independent set
206 
207 static void hIndSolve(scmon pure, int Npure, scfmon rad, int Nrad,
208  varset var, int Nvar)
209 {
210  int dn, iv, rad0, b, c, x;
211  scmon pn;
212  scfmon rn;
213  if (Nrad < 2)
214  {
215  dn = Npure + Nrad;
216  if (dn < hCo)
217  {
218  hCo = dn;
219  for (iv=(currRing->N); iv; iv--)
220  {
221  if (pure[iv])
222  hInd[iv] = 0;
223  else
224  hInd[iv] = 1;
225  }
226  if (Nrad)
227  {
228  pn = *rad;
229  iv = Nvar;
230  loop
231  {
232  x = var[iv];
233  if (pn[x])
234  {
235  hInd[x] = 0;
236  break;
237  }
238  iv--;
239  }
240  }
241  }
242  return;
243  }
244  if (Npure+1 >= hCo)
245  return;
246  iv = Nvar;
247  while(pure[var[iv]]) iv--;
248  hStepR(rad, Nrad, var, iv, &rad0);
249  if (rad0)
250  {
251  iv--;
252  if (rad0 < Nrad)
253  {
254  pn = hGetpure(pure);
255  rn = hGetmem(Nrad, rad, radmem[iv]);
256  pn[var[iv + 1]] = 1;
257  hIndSolve(pn, Npure + 1, rn, rad0, var, iv);
258  pn[var[iv + 1]] = 0;
259  b = rad0;
260  c = Nrad;
261  hElimR(rn, &rad0, b, c, var, iv);
262  hPure(rn, b, &c, var, iv, pn, &x);
263  hLex2R(rn, rad0, b, c, var, iv, hwork);
264  rad0 += (c - b);
265  hIndSolve(pn, Npure + x, rn, rad0, var, iv);
266  }
267  else
268  {
269  hIndSolve(pure, Npure, rad, Nrad, var, iv);
270  }
271  }
272  else
273  {
274  hCo = Npure + 1;
275  for (x=(currRing->N); x; x--)
276  {
277  if (pure[x])
278  hInd[x] = 0;
279  else
280  hInd[x] = 1;
281  }
282  hInd[var[iv]] = 0;
283  }
284 }
285 
286 intvec * scIndIntvec(ideal S, ideal Q)
287 {
288  id_Test(S, currRing);
289  if( Q!=NULL ) id_Test(Q, currRing);
290 
291  intvec *Set=new intvec((currRing->N));
292  int mc,i;
293  hexist = hInit(S, Q, &hNexist);
294  if (hNexist==0)
295  {
296  for(i=0; i<(currRing->N); i++)
297  (*Set)[i]=1;
298  return Set;
299  }
300  hwork = (scfmon)omAlloc(hNexist * sizeof(scmon));
301  hvar = (varset)omAlloc(((currRing->N) + 1) * sizeof(int));
302  hpure = (scmon)omAlloc((1 + ((currRing->N) * (currRing->N))) * sizeof(int));
303  hInd = (scmon)omAlloc0((1 + (currRing->N)) * sizeof(int));
304  mc = hisModule;
305  if (mc==0)
306  {
307  hrad = hexist;
308  hNrad = hNexist;
309  }
310  else
311  hrad = (scfmon)omAlloc(hNexist * sizeof(scmon));
312  radmem = hCreate((currRing->N) - 1);
313  hCo = (currRing->N) + 1;
314  loop
315  {
316  if (mc!=0)
317  hComp(hexist, hNexist, mc, hrad, &hNrad);
318  if (hNrad!=0)
319  {
320  hNvar = (currRing->N);
321  hRadical(hrad, &hNrad, hNvar);
322  hSupp(hrad, hNrad, hvar, &hNvar);
323  if (hNvar!=0)
324  {
325  memset(hpure, 0, ((currRing->N) + 1) * sizeof(int));
326  hPure(hrad, 0, &hNrad, hvar, hNvar, hpure, &hNpure);
327  hLexR(hrad, hNrad, hvar, hNvar);
329  }
330  }
331  else
332  {
333  hCo = 0;
334  break;
335  }
336  mc--;
337  if (mc <= 0)
338  break;
339  }
340  for(i=0; i<(currRing->N); i++)
341  (*Set)[i] = hInd[i+1];
342  hKill(radmem, (currRing->N) - 1);
343  omFreeSize((ADDRESS)hpure, (1 + ((currRing->N) * (currRing->N))) * sizeof(int));
344  omFreeSize((ADDRESS)hInd, (1 + (currRing->N)) * sizeof(int));
345  omFreeSize((ADDRESS)hvar, ((currRing->N) + 1) * sizeof(int));
346  omFreeSize((ADDRESS)hwork, hNexist * sizeof(scmon));
348  if (hisModule)
349  omFreeSize((ADDRESS)hrad, hNexist * sizeof(scmon));
350  return Set;
351 }
352 
354 
355 static BOOLEAN hNotZero(scfmon rad, int Nrad, varset var, int Nvar)
356 {
357  int k1, i;
358  k1 = var[Nvar];
359  i = 0;
360  loop
361  {
362  if (rad[i][k1]==0)
363  return FALSE;
364  i++;
365  if (i == Nrad)
366  return TRUE;
367  }
368 }
369 
370 static void hIndep(scmon pure)
371 {
372  int iv;
373  intvec *Set;
374 
375  Set = ISet->set = new intvec((currRing->N));
376  for (iv=(currRing->N); iv!=0 ; iv--)
377  {
378  (*Set)[iv-1] = (pure[iv]==0);
379  }
381  hMu++;
382 }
383 
384 void hIndMult(scmon pure, int Npure, scfmon rad, int Nrad,
385  varset var, int Nvar)
386 {
387  int dn, iv, rad0, b, c, x;
388  scmon pn;
389  scfmon rn;
390  if (Nrad < 2)
391  {
392  dn = Npure + Nrad;
393  if (dn == hCo)
394  {
395  if (Nrad==0)
396  hIndep(pure);
397  else
398  {
399  pn = *rad;
400  for (iv = Nvar; iv!=0; iv--)
401  {
402  x = var[iv];
403  if (pn[x])
404  {
405  pure[x] = 1;
406  hIndep(pure);
407  pure[x] = 0;
408  }
409  }
410  }
411  }
412  return;
413  }
414  iv = Nvar;
415  dn = Npure+1;
416  if (dn >= hCo)
417  {
418  if (dn > hCo)
419  return;
420  loop
421  {
422  if(!pure[var[iv]])
423  {
424  if(hNotZero(rad, Nrad, var, iv))
425  {
426  pure[var[iv]] = 1;
427  hIndep(pure);
428  pure[var[iv]] = 0;
429  }
430  }
431  iv--;
432  if (!iv)
433  return;
434  }
435  }
436  while(pure[var[iv]]) iv--;
437  hStepR(rad, Nrad, var, iv, &rad0);
438  iv--;
439  if (rad0 < Nrad)
440  {
441  pn = hGetpure(pure);
442  rn = hGetmem(Nrad, rad, radmem[iv]);
443  pn[var[iv + 1]] = 1;
444  hIndMult(pn, Npure + 1, rn, rad0, var, iv);
445  pn[var[iv + 1]] = 0;
446  b = rad0;
447  c = Nrad;
448  hElimR(rn, &rad0, b, c, var, iv);
449  hPure(rn, b, &c, var, iv, pn, &x);
450  hLex2R(rn, rad0, b, c, var, iv, hwork);
451  rad0 += (c - b);
452  hIndMult(pn, Npure + x, rn, rad0, var, iv);
453  }
454  else
455  {
456  hIndMult(pure, Npure, rad, Nrad, var, iv);
457  }
458 }
459 
460 /*3
461 * consider indset x := !pure
462 * (for all i) (if(sm(i) > x) return FALSE)
463 * else return TRUE
464 */
465 static BOOLEAN hCheck1(indset sm, scmon pure)
466 {
467  int iv;
468  intvec *Set;
469  while (sm->nx != NULL)
470  {
471  Set = sm->set;
472  iv=(currRing->N);
473  loop
474  {
475  if (((*Set)[iv-1] == 0) && (pure[iv] == 0))
476  break;
477  iv--;
478  if (iv == 0)
479  return FALSE;
480  }
481  sm = sm->nx;
482  }
483  return TRUE;
484 }
485 
486 /*3
487 * consider indset x := !pure
488 * (for all i) if(x > sm(i)) delete sm(i))
489 * return (place for x)
490 */
491 static indset hCheck2(indset sm, scmon pure)
492 {
493  int iv;
494  intvec *Set;
495  indset be, a1 = NULL;
496  while (sm->nx != NULL)
497  {
498  Set = sm->set;
499  iv=(currRing->N);
500  loop
501  {
502  if ((pure[iv] == 1) && ((*Set)[iv-1] == 1))
503  break;
504  iv--;
505  if (iv == 0)
506  {
507  if (a1 == NULL)
508  {
509  a1 = sm;
510  }
511  else
512  {
513  hMu2--;
514  be->nx = sm->nx;
515  delete Set;
517  sm = be;
518  }
519  break;
520  }
521  }
522  be = sm;
523  sm = sm->nx;
524  }
525  if (a1 != NULL)
526  {
527  return a1;
528  }
529  else
530  {
531  hMu2++;
532  sm->set = new intvec((currRing->N));
533  sm->nx = (indset)omAlloc0Bin(indlist_bin);
534  return sm;
535  }
536 }
537 
538 /*2
539 * definition x >= y
540 * x(i) == 0 => y(i) == 0
541 * > ex. j with x(j) == 1 and y(j) == 0
542 */
543 static void hCheckIndep(scmon pure)
544 {
545  intvec *Set;
546  indset res;
547  int iv;
548  if (hCheck1(ISet, pure))
549  {
550  if (hCheck1(JSet, pure))
551  {
552  res = hCheck2(JSet,pure);
553  if (res == NULL)
554  return;
555  Set = res->set;
556  for (iv=(currRing->N); iv; iv--)
557  {
558  (*Set)[iv-1] = (pure[iv]==0);
559  }
560  }
561  }
562 }
563 
564 void hIndAllMult(scmon pure, int Npure, scfmon rad, int Nrad,
565  varset var, int Nvar)
566 {
567  int dn, iv, rad0, b, c, x;
568  scmon pn;
569  scfmon rn;
570  if (Nrad < 2)
571  {
572  dn = Npure + Nrad;
573  if (dn > hCo)
574  {
575  if (!Nrad)
576  hCheckIndep(pure);
577  else
578  {
579  pn = *rad;
580  for (iv = Nvar; iv; iv--)
581  {
582  x = var[iv];
583  if (pn[x])
584  {
585  pure[x] = 1;
586  hCheckIndep(pure);
587  pure[x] = 0;
588  }
589  }
590  }
591  }
592  return;
593  }
594  iv = Nvar;
595  while(pure[var[iv]]) iv--;
596  hStepR(rad, Nrad, var, iv, &rad0);
597  iv--;
598  if (rad0 < Nrad)
599  {
600  pn = hGetpure(pure);
601  rn = hGetmem(Nrad, rad, radmem[iv]);
602  pn[var[iv + 1]] = 1;
603  hIndAllMult(pn, Npure + 1, rn, rad0, var, iv);
604  pn[var[iv + 1]] = 0;
605  b = rad0;
606  c = Nrad;
607  hElimR(rn, &rad0, b, c, var, iv);
608  hPure(rn, b, &c, var, iv, pn, &x);
609  hLex2R(rn, rad0, b, c, var, iv, hwork);
610  rad0 += (c - b);
611  hIndAllMult(pn, Npure + x, rn, rad0, var, iv);
612  }
613  else
614  {
615  hIndAllMult(pure, Npure, rad, Nrad, var, iv);
616  }
617 }
618 
619 // multiplicity
620 
621 static long hZeroMult(scmon pure, scfmon stc, int Nstc, varset var, int Nvar)
622 {
623  int iv = Nvar -1, a, a0, a1, b, i;
624  long sum;
625  int x, x0;
626  scmon pn;
627  scfmon sn;
628  if (!iv)
629  return pure[var[1]];
630  else if (!Nstc)
631  {
632  sum = 1;
633  for (i = Nvar; i; i--)
634  sum *= pure[var[i]];
635  return sum;
636  }
637  x = a = 0;
638  pn = hGetpure(pure);
639  sn = hGetmem(Nstc, stc, stcmem[iv]);
640  hStepS(sn, Nstc, var, Nvar, &a, &x);
641  if (a == Nstc)
642  {
643  #if SIZEOF_LONG==8
644  return (long)pure[var[Nvar]] * hZeroMult(pn, sn, a, var, iv);
645  #else
646  int64 t=hZeroMult(pn, sn, a, var, iv);
647  t *= pure[var[Nvar]];
648  if ((t>=INT_MIN)&&(t<=INT_MAX)) sum=t;
649  else if (!errorreported) WerrorS("int overflow in vdim 3");
650  return sum;
651  #endif
652  }
653  else
654  {
655  #if SIZEOF_LONG==8
656  sum = x * hZeroMult(pn, sn, a, var, iv);
657  #else
658  int64 t=hZeroMult(pn, sn, a, var, iv);
659  t *= x;
660  if ((t>=INT_MIN)&&(t<=INT_MAX)) sum=t;
661  else if (!errorreported) WerrorS("int overflow in vdim 4");
662  #endif
663  }
664  b = a;
665  loop
666  {
667  a0 = a;
668  x0 = x;
669  hStepS(sn, Nstc, var, Nvar, &a, &x);
670  hElimS(sn, &b, a0, a, var, iv);
671  a1 = a;
672  hPure(sn, a0, &a1, var, iv, pn, &i);
673  hLex2S(sn, b, a0, a1, var, iv, hwork);
674  b += (a1 - a0);
675  if (a < Nstc)
676  {
677  #if SIZEOF_LONG==8
678  sum += (long)(x - x0) * hZeroMult(pn, sn, b, var, iv);
679  #else
680  int64 t=hZeroMult(pn, sn, b, var, iv);
681  t *= (x-x0);
682  t += sum;
683  if ((t>=INT_MIN)&&(t<=INT_MAX)) sum=t;
684  else if (!errorreported) WerrorS("int overflow in vdim 1");
685  #endif
686  }
687  else
688  {
689  #if SIZEOF_LONG==8
690  sum += (long)(pure[var[Nvar]] - x0) * hZeroMult(pn, sn, b, var, iv);
691  #else
692  int64 t=hZeroMult(pn, sn, b, var, iv);
693  t *= (pure[var[Nvar]]-x0);
694  t += sum;
695  if ((t>=INT_MIN)&&(t<=INT_MAX)) sum=t;
696  else if (!errorreported) WerrorS("int overflow in vdim 2");
697  #endif
698  return sum;
699  }
700  }
701 }
702 
703 static void hProject(scmon pure, varset sel)
704 {
705  int i, i0, k;
706  i0 = 0;
707  for (i = 1; i <= (currRing->N); i++)
708  {
709  if (pure[i])
710  {
711  i0++;
712  sel[i0] = i;
713  }
714  }
715  i = hNstc;
716  memcpy(hwork, hstc, i * sizeof(scmon));
717  hStaircase(hwork, &i, sel, i0);
718  if ((i0 > 2) && (i > 10))
719  hOrdSupp(hwork, i, sel, i0);
720  memset(hpur0, 0, ((currRing->N) + 1) * sizeof(int));
721  hPure(hwork, 0, &i, sel, i0, hpur0, &k);
722  hLexS(hwork, i, sel, i0);
723  hMu += hZeroMult(hpur0, hwork, i, sel, i0);
724 }
725 
726 static void hDimMult(scmon pure, int Npure, scfmon rad, int Nrad,
727  varset var, int Nvar)
728 {
729  int dn, iv, rad0, b, c, x;
730  scmon pn;
731  scfmon rn;
732  if (Nrad < 2)
733  {
734  dn = Npure + Nrad;
735  if (dn == hCo)
736  {
737  if (!Nrad)
738  hProject(pure, hsel);
739  else
740  {
741  pn = *rad;
742  for (iv = Nvar; iv; iv--)
743  {
744  x = var[iv];
745  if (pn[x])
746  {
747  pure[x] = 1;
748  hProject(pure, hsel);
749  pure[x] = 0;
750  }
751  }
752  }
753  }
754  return;
755  }
756  iv = Nvar;
757  dn = Npure+1;
758  if (dn >= hCo)
759  {
760  if (dn > hCo)
761  return;
762  loop
763  {
764  if(!pure[var[iv]])
765  {
766  if(hNotZero(rad, Nrad, var, iv))
767  {
768  pure[var[iv]] = 1;
769  hProject(pure, hsel);
770  pure[var[iv]] = 0;
771  }
772  }
773  iv--;
774  if (!iv)
775  return;
776  }
777  }
778  while(pure[var[iv]]) iv--;
779  hStepR(rad, Nrad, var, iv, &rad0);
780  iv--;
781  if (rad0 < Nrad)
782  {
783  pn = hGetpure(pure);
784  rn = hGetmem(Nrad, rad, radmem[iv]);
785  pn[var[iv + 1]] = 1;
786  hDimMult(pn, Npure + 1, rn, rad0, var, iv);
787  pn[var[iv + 1]] = 0;
788  b = rad0;
789  c = Nrad;
790  hElimR(rn, &rad0, b, c, var, iv);
791  hPure(rn, b, &c, var, iv, pn, &x);
792  hLex2R(rn, rad0, b, c, var, iv, hwork);
793  rad0 += (c - b);
794  hDimMult(pn, Npure + x, rn, rad0, var, iv);
795  }
796  else
797  {
798  hDimMult(pure, Npure, rad, Nrad, var, iv);
799  }
800 }
801 
802 static void hDegree(ideal S, ideal Q)
803 {
804  id_Test(S, currRing);
805  if( Q!=NULL ) id_Test(Q, currRing);
806 
807  int di;
808  int mc;
809  hexist = hInit(S, Q, &hNexist);
810  if (!hNexist)
811  {
812  hCo = 0;
813  hMu = 1;
814  return;
815  }
816  //hWeight();
817  hwork = (scfmon)omAlloc(hNexist * sizeof(scmon));
818  hvar = (varset)omAlloc(((currRing->N) + 1) * sizeof(int));
819  hsel = (varset)omAlloc(((currRing->N) + 1) * sizeof(int));
820  hpure = (scmon)omAlloc((1 + ((currRing->N) * (currRing->N))) * sizeof(int));
821  hpur0 = (scmon)omAlloc((1 + ((currRing->N) * (currRing->N))) * sizeof(int));
822  mc = hisModule;
823  hrad = (scfmon)omAlloc(hNexist * sizeof(scmon));
824  if (!mc)
825  {
826  memcpy(hrad, hexist, hNexist * sizeof(scmon));
827  hstc = hexist;
828  hNrad = hNstc = hNexist;
829  }
830  else
831  hstc = (scfmon)omAlloc(hNexist * sizeof(scmon));
832  radmem = hCreate((currRing->N) - 1);
833  stcmem = hCreate((currRing->N) - 1);
834  hCo = (currRing->N) + 1;
835  di = hCo + 1;
836  loop
837  {
838  if (mc)
839  {
840  hComp(hexist, hNexist, mc, hrad, &hNrad);
841  hNstc = hNrad;
842  memcpy(hstc, hrad, hNrad * sizeof(scmon));
843  }
844  if (hNrad)
845  {
846  hNvar = (currRing->N);
847  hRadical(hrad, &hNrad, hNvar);
848  hSupp(hrad, hNrad, hvar, &hNvar);
849  if (hNvar)
850  {
851  hCo = hNvar;
852  memset(hpure, 0, ((currRing->N) + 1) * sizeof(int));
853  hPure(hrad, 0, &hNrad, hvar, hNvar, hpure, &hNpure);
854  hLexR(hrad, hNrad, hvar, hNvar);
856  }
857  }
858  else
859  {
860  hNvar = 1;
861  hCo = 0;
862  }
863  if (hCo < di)
864  {
865  di = hCo;
866  hMu = 0;
867  }
868  if (hNvar && (hCo == di))
869  {
870  if (di && (di < (currRing->N)))
872  else if (!di)
873  hMu++;
874  else
875  {
877  if ((hNvar > 2) && (hNstc > 10))
879  memset(hpur0, 0, ((currRing->N) + 1) * sizeof(int));
880  hPure(hstc, 0, &hNstc, hvar, hNvar, hpur0, &hNpure);
881  hLexS(hstc, hNstc, hvar, hNvar);
883  }
884  }
885  mc--;
886  if (mc <= 0)
887  break;
888  }
889  hCo = di;
890  hKill(stcmem, (currRing->N) - 1);
891  hKill(radmem, (currRing->N) - 1);
892  omFreeSize((ADDRESS)hpur0, (1 + ((currRing->N) * (currRing->N))) * sizeof(int));
893  omFreeSize((ADDRESS)hpure, (1 + ((currRing->N) * (currRing->N))) * sizeof(int));
894  omFreeSize((ADDRESS)hsel, ((currRing->N) + 1) * sizeof(int));
895  omFreeSize((ADDRESS)hvar, ((currRing->N) + 1) * sizeof(int));
896  omFreeSize((ADDRESS)hwork, hNexist * sizeof(scmon));
897  omFreeSize((ADDRESS)hrad, hNexist * sizeof(scmon));
899  if (hisModule)
900  omFreeSize((ADDRESS)hstc, hNexist * sizeof(scmon));
901 }
902 
903 int scMultInt(ideal S, ideal Q)
904 {
905  id_Test(S, currRing);
906  if( Q!=NULL ) id_Test(Q, currRing);
907 
908  hDegree(S, Q);
909  return hMu;
910 }
911 
912 void scPrintDegree(int co, int mu)
913 {
914  int di = (currRing->N)-co;
915  if (currRing->OrdSgn == 1)
916  {
917  if (di>0)
918  Print("// dimension (proj.) = %d\n// degree (proj.) = %d\n", di-1, mu);
919  else
920  Print("// dimension (affine) = 0\n// degree (affine) = %d\n", mu);
921  }
922  else
923  Print("// dimension (local) = %d\n// multiplicity = %d\n", di, mu);
924 }
925 
926 void scDegree(ideal S, intvec *modulweight, ideal Q)
927 {
928  id_Test(S, currRing);
929  if( Q!=NULL ) id_Test(Q, currRing);
930 
931  int co, mu, l;
932  intvec *hseries2;
933  intvec *hseries1 = hFirstSeries(S, modulweight, Q);
934  if (errorreported) return;
935  l = hseries1->length()-1;
936  if (l > 1)
937  hseries2 = hSecondSeries(hseries1);
938  else
939  hseries2 = hseries1;
940  hDegreeSeries(hseries1, hseries2, &co, &mu);
941  if ((l == 1) &&(mu == 0))
942  scPrintDegree((currRing->N)+1, 0);
943  else
944  scPrintDegree(co, mu);
945  if (l>1)
946  delete hseries1;
947  delete hseries2;
948 }
949 
950 long scMult0Int(ideal S, ideal Q)
951 {
952  id_LmTest(S, currRing);
953  if (Q!=NULL) id_LmTest(Q, currRing);
954 
955  int mc;
956  hexist = hInit(S, Q, &hNexist);
957  if (!hNexist)
958  {
959  hMu = -1;
960  return -1;
961  }
962  else
963  hMu = 0;
964 
965  const ring r = currRing;
966 
967  hwork = (scfmon)omAlloc(hNexist * sizeof(scmon));
968  hvar = (varset)omAlloc(((r->N) + 1) * sizeof(int));
969  hpur0 = (scmon)omAlloc((1 + ((r->N) * (r->N))) * sizeof(int));
970  mc = hisModule;
971  if (!mc)
972  {
973  hstc = hexist;
974  hNstc = hNexist;
975  }
976  else
977  hstc = (scfmon)omAlloc(hNexist * sizeof(scmon));
978  stcmem = hCreate((r->N) - 1);
979  loop
980  {
981  if (mc)
982  {
983  hComp(hexist, hNexist, mc, hstc, &hNstc);
984  if (!hNstc)
985  {
986  hMu = -1;
987  break;
988  }
989  }
990  hNvar = (r->N);
991  for (int i = hNvar; i; i--)
992  hvar[i] = i;
994  hSupp(hstc, hNstc, hvar, &hNvar);
995  if ((hNvar == (r->N)) && (hNstc >= (r->N)))
996  {
997  if ((hNvar > 2) && (hNstc > 10))
999  memset(hpur0, 0, ((r->N) + 1) * sizeof(int));
1000  hPure(hstc, 0, &hNstc, hvar, hNvar, hpur0, &hNpure);
1001  if (hNpure == hNvar)
1002  {
1003  hLexS(hstc, hNstc, hvar, hNvar);
1005  }
1006  else
1007  hMu = -1;
1008  }
1009  else if (hNvar)
1010  hMu = -1;
1011  mc--;
1012  if (mc <= 0 || hMu < 0)
1013  break;
1014  }
1015  hKill(stcmem, (r->N) - 1);
1016  omFreeSize((ADDRESS)hpur0, (1 + ((r->N) * (r->N))) * sizeof(int));
1017  omFreeSize((ADDRESS)hvar, ((r->N) + 1) * sizeof(int));
1018  omFreeSize((ADDRESS)hwork, hNexist * sizeof(scmon));
1020  if (hisModule)
1021  omFreeSize((ADDRESS)hstc, hNexist * sizeof(scmon));
1022  return hMu;
1023 }
1024 
1025 // HC
1026 
1028 
1029 static void hHedge(poly hEdge)
1030 {
1031  pSetm(pWork);
1032  if (pLmCmp(pWork, hEdge) == currRing->OrdSgn)
1033  {
1034  for (int i = hNvar; i>0; i--)
1035  pSetExp(hEdge,i, pGetExp(pWork,i));
1036  pSetm(hEdge);
1037  }
1038 }
1039 
1040 
1041 static void hHedgeStep(scmon pure, scfmon stc,
1042  int Nstc, varset var, int Nvar,poly hEdge)
1043 {
1044  int iv = Nvar -1, k = var[Nvar], a, a0, a1, b, i;
1045  int x/*, x0*/;
1046  scmon pn;
1047  scfmon sn;
1048  if (iv==0)
1049  {
1050  pSetExp(pWork, k, pure[k]);
1051  hHedge(hEdge);
1052  return;
1053  }
1054  else if (Nstc==0)
1055  {
1056  for (i = Nvar; i>0; i--)
1057  pSetExp(pWork, var[i], pure[var[i]]);
1058  hHedge(hEdge);
1059  return;
1060  }
1061  x = a = 0;
1062  pn = hGetpure(pure);
1063  sn = hGetmem(Nstc, stc, stcmem[iv]);
1064  hStepS(sn, Nstc, var, Nvar, &a, &x);
1065  if (a == Nstc)
1066  {
1067  pSetExp(pWork, k, pure[k]);
1068  hHedgeStep(pn, sn, a, var, iv,hEdge);
1069  return;
1070  }
1071  else
1072  {
1073  pSetExp(pWork, k, x);
1074  hHedgeStep(pn, sn, a, var, iv,hEdge);
1075  }
1076  b = a;
1077  loop
1078  {
1079  a0 = a;
1080  // x0 = x;
1081  hStepS(sn, Nstc, var, Nvar, &a, &x);
1082  hElimS(sn, &b, a0, a, var, iv);
1083  a1 = a;
1084  hPure(sn, a0, &a1, var, iv, pn, &i);
1085  hLex2S(sn, b, a0, a1, var, iv, hwork);
1086  b += (a1 - a0);
1087  if (a < Nstc)
1088  {
1089  pSetExp(pWork, k, x);
1090  hHedgeStep(pn, sn, b, var, iv,hEdge);
1091  }
1092  else
1093  {
1094  pSetExp(pWork, k, pure[k]);
1095  hHedgeStep(pn, sn, b, var, iv,hEdge);
1096  return;
1097  }
1098  }
1099 }
1100 
1101 void scComputeHC(ideal S, ideal Q, int ak, poly &hEdge)
1102 {
1103  id_LmTest(S, currRing);
1104  if (Q!=NULL) id_LmTest(Q, currRing);
1105 
1106  int i;
1107  int k = ak;
1108  #ifdef HAVE_RINGS
1109  if (rField_is_Ring(currRing) && (currRing->OrdSgn == -1))
1110  {
1111  //consider just monic generators (over rings with zero-divisors)
1112  ideal SS=id_Head(S,currRing);
1113  for(i=0;i<=idElem(S);i++)
1114  {
1115  if((SS->m[i]!=NULL)
1116  && ((p_IsPurePower(SS->m[i],currRing)==0)
1117  ||(!n_IsUnit(pGetCoeff(SS->m[i]), currRing->cf))))
1118  {
1119  p_Delete(&SS->m[i],currRing);
1120  }
1121  }
1122  S=id_Copy(SS,currRing);
1123  idSkipZeroes(S);
1124  }
1125  #if 0
1126  printf("\nThis is HC:\n");
1127  for(int ii=0;ii<=idElem(S);ii++)
1128  {
1129  pWrite(S->m[ii]);
1130  }
1131  //getchar();
1132  #endif
1133  #endif
1134  if(idElem(S) == 0)
1135  return;
1136  hNvar = (currRing->N);
1137  hexist = hInit(S, Q, &hNexist);
1138  if (k!=0)
1139  hComp(hexist, hNexist, k, hexist, &hNstc);
1140  else
1141  hNstc = hNexist;
1142  assume(hNexist > 0);
1143  hwork = (scfmon)omAlloc(hNexist * sizeof(scmon));
1144  hvar = (varset)omAlloc((hNvar + 1) * sizeof(int));
1145  hpure = (scmon)omAlloc((1 + (hNvar * hNvar)) * sizeof(int));
1146  stcmem = hCreate(hNvar - 1);
1147  for (i = hNvar; i>0; i--)
1148  hvar[i] = i;
1150  if ((hNvar > 2) && (hNstc > 10))
1152  memset(hpure, 0, (hNvar + 1) * sizeof(int));
1153  hPure(hexist, 0, &hNstc, hvar, hNvar, hpure, &hNpure);
1154  hLexS(hexist, hNstc, hvar, hNvar);
1155  if (hEdge!=NULL)
1156  pLmFree(hEdge);
1157  hEdge = pInit();
1158  pWork = pInit();
1159  hHedgeStep(hpure, hexist, hNstc, hvar, hNvar,hEdge);
1160  pSetComp(hEdge,ak);
1161  hKill(stcmem, hNvar - 1);
1162  omFreeSize((ADDRESS)hwork, hNexist * sizeof(scmon));
1163  omFreeSize((ADDRESS)hvar, (hNvar + 1) * sizeof(int));
1164  omFreeSize((ADDRESS)hpure, (1 + (hNvar * hNvar)) * sizeof(int));
1166  pLmFree(pWork);
1167 }
1168 
1169 
1170 
1171 // kbase
1172 
1175 
1176 static void scElKbase()
1177 {
1178  poly q = pInit();
1179  pSetCoeff0(q,nInit(1));
1180  pSetExpV(q,act);
1181  pNext(q) = NULL;
1182  last = pNext(last) = q;
1183 }
1184 
1185 static int scMax( int i, scfmon stc, int Nvar)
1186 {
1187  int x, y=stc[0][Nvar];
1188  for (; i;)
1189  {
1190  i--;
1191  x = stc[i][Nvar];
1192  if (x > y) y = x;
1193  }
1194  return y;
1195 }
1196 
1197 static int scMin( int i, scfmon stc, int Nvar)
1198 {
1199  int x, y=stc[0][Nvar];
1200  for (; i;)
1201  {
1202  i--;
1203  x = stc[i][Nvar];
1204  if (x < y) y = x;
1205  }
1206  return y;
1207 }
1208 
1209 static int scRestrict( int &Nstc, scfmon stc, int Nvar)
1210 {
1211  int x, y;
1212  int i, j, Istc = Nstc;
1213 
1214  y = MAX_INT_VAL;
1215  for (i=Nstc-1; i>=0; i--)
1216  {
1217  j = Nvar-1;
1218  loop
1219  {
1220  if(stc[i][j] != 0) break;
1221  j--;
1222  if (j == 0)
1223  {
1224  Istc--;
1225  x = stc[i][Nvar];
1226  if (x < y) y = x;
1227  stc[i] = NULL;
1228  break;
1229  }
1230  }
1231  }
1232  if (Istc < Nstc)
1233  {
1234  for (i=Nstc-1; i>=0; i--)
1235  {
1236  if (stc[i] && (stc[i][Nvar] >= y))
1237  {
1238  Istc--;
1239  stc[i] = NULL;
1240  }
1241  }
1242  j = 0;
1243  while (stc[j]) j++;
1244  i = j+1;
1245  for(; i<Nstc; i++)
1246  {
1247  if (stc[i])
1248  {
1249  stc[j] = stc[i];
1250  j++;
1251  }
1252  }
1253  Nstc = Istc;
1254  return y;
1255  }
1256  else
1257  return -1;
1258 }
1259 
1260 static void scAll( int Nvar, int deg)
1261 {
1262  int i;
1263  int d = deg;
1264  if (d == 0)
1265  {
1266  for (i=Nvar; i; i--) act[i] = 0;
1267  scElKbase();
1268  return;
1269  }
1270  if (Nvar == 1)
1271  {
1272  act[1] = d;
1273  scElKbase();
1274  return;
1275  }
1276  do
1277  {
1278  act[Nvar] = d;
1279  scAll(Nvar-1, deg-d);
1280  d--;
1281  } while (d >= 0);
1282 }
1283 
1284 static void scAllKbase( int Nvar, int ideg, int deg)
1285 {
1286  do
1287  {
1288  act[Nvar] = ideg;
1289  scAll(Nvar-1, deg-ideg);
1290  ideg--;
1291  } while (ideg >= 0);
1292 }
1293 
1294 static void scDegKbase( scfmon stc, int Nstc, int Nvar, int deg)
1295 {
1296  int Ivar, Istc, i, j;
1297  scfmon sn;
1298  int x, ideg;
1299 
1300  if (deg == 0)
1301  {
1302  for (i=Nstc-1; i>=0; i--)
1303  {
1304  for (j=Nvar;j;j--){ if(stc[i][j]) break; }
1305  if (j==0){return;}
1306  }
1307  for (i=Nvar; i; i--) act[i] = 0;
1308  scElKbase();
1309  return;
1310  }
1311  if (Nvar == 1)
1312  {
1313  for (i=Nstc-1; i>=0; i--) if(deg >= stc[i][1]) return;
1314  act[1] = deg;
1315  scElKbase();
1316  return;
1317  }
1318  Ivar = Nvar-1;
1319  sn = hGetmem(Nstc, stc, stcmem[Ivar]);
1320  x = scRestrict(Nstc, sn, Nvar);
1321  if (x <= 0)
1322  {
1323  if (x == 0) return;
1324  ideg = deg;
1325  }
1326  else
1327  {
1328  if (deg < x) ideg = deg;
1329  else ideg = x-1;
1330  if (Nstc == 0)
1331  {
1332  scAllKbase(Nvar, ideg, deg);
1333  return;
1334  }
1335  }
1336  loop
1337  {
1338  x = scMax(Nstc, sn, Nvar);
1339  while (ideg >= x)
1340  {
1341  act[Nvar] = ideg;
1342  scDegKbase(sn, Nstc, Ivar, deg-ideg);
1343  ideg--;
1344  }
1345  if (ideg < 0) return;
1346  Istc = Nstc;
1347  for (i=Nstc-1; i>=0; i--)
1348  {
1349  if (ideg < sn[i][Nvar])
1350  {
1351  Istc--;
1352  sn[i] = NULL;
1353  }
1354  }
1355  if (Istc == 0)
1356  {
1357  scAllKbase(Nvar, ideg, deg);
1358  return;
1359  }
1360  j = 0;
1361  while (sn[j]) j++;
1362  i = j+1;
1363  for (; i<Nstc; i++)
1364  {
1365  if (sn[i])
1366  {
1367  sn[j] = sn[i];
1368  j++;
1369  }
1370  }
1371  Nstc = Istc;
1372  }
1373 }
1374 
1375 static void scInKbase( scfmon stc, int Nstc, int Nvar)
1376 {
1377  int Ivar, Istc, i, j;
1378  scfmon sn;
1379  int x, ideg;
1380 
1381  if (Nvar == 1)
1382  {
1383  ideg = scMin(Nstc, stc, 1);
1384  while (ideg > 0)
1385  {
1386  ideg--;
1387  act[1] = ideg;
1388  scElKbase();
1389  }
1390  return;
1391  }
1392  Ivar = Nvar-1;
1393  sn = hGetmem(Nstc, stc, stcmem[Ivar]);
1394  x = scRestrict(Nstc, sn, Nvar);
1395  if (x == 0) return;
1396  ideg = x-1;
1397  loop
1398  {
1399  x = scMax(Nstc, sn, Nvar);
1400  while (ideg >= x)
1401  {
1402  act[Nvar] = ideg;
1403  scInKbase(sn, Nstc, Ivar);
1404  ideg--;
1405  }
1406  if (ideg < 0) return;
1407  Istc = Nstc;
1408  for (i=Nstc-1; i>=0; i--)
1409  {
1410  if (ideg < sn[i][Nvar])
1411  {
1412  Istc--;
1413  sn[i] = NULL;
1414  }
1415  }
1416  j = 0;
1417  while (sn[j]) j++;
1418  i = j+1;
1419  for (; i<Nstc; i++)
1420  {
1421  if (sn[i])
1422  {
1423  sn[j] = sn[i];
1424  j++;
1425  }
1426  }
1427  Nstc = Istc;
1428  }
1429 }
1430 
1431 static ideal scIdKbase(poly q, const int rank)
1432 {
1433  ideal res = idInit(pLength(q), rank);
1434  polyset mm = res->m;
1435  do
1436  {
1437  *mm = q; ++mm;
1438 
1439  const poly p = pNext(q);
1440  pNext(q) = NULL;
1441  q = p;
1442 
1443  } while (q!=NULL);
1444 
1445  id_Test(res, currRing); // WRONG RANK!!!???
1446  return res;
1447 }
1448 
1449 ideal scKBase(int deg, ideal s, ideal Q, intvec * mv)
1450 {
1451  if( Q!=NULL) id_Test(Q, currRing);
1452 
1453  int i, di;
1454  poly p;
1455 
1456  if (deg < 0)
1457  {
1458  di = scDimInt(s, Q);
1459  if (di != 0)
1460  {
1461  //Werror("KBase not finite");
1462  return idInit(1,s->rank);
1463  }
1464  }
1465  stcmem = hCreate((currRing->N) - 1);
1466  hexist = hInit(s, Q, &hNexist);
1467  p = last = pInit();
1468  /*pNext(p) = NULL;*/
1469  act = (scmon)omAlloc(((currRing->N) + 1) * sizeof(int));
1470  *act = 0;
1471  if (!hNexist)
1472  {
1473  scAll((currRing->N), deg);
1474  goto ende;
1475  }
1476  if (!hisModule)
1477  {
1478  if (deg < 0) scInKbase(hexist, hNexist, (currRing->N));
1479  else scDegKbase(hexist, hNexist, (currRing->N), deg);
1480  }
1481  else
1482  {
1483  hstc = (scfmon)omAlloc(hNexist * sizeof(scmon));
1484  for (i = 1; i <= hisModule; i++)
1485  {
1486  *act = i;
1487  hComp(hexist, hNexist, i, hstc, &hNstc);
1488  int deg_ei=deg;
1489  if (mv!=NULL) deg_ei -= (*mv)[i-1];
1490  if ((deg < 0) || (deg_ei>=0))
1491  {
1492  if (hNstc)
1493  {
1494  if (deg < 0) scInKbase(hstc, hNstc, (currRing->N));
1495  else scDegKbase(hstc, hNstc, (currRing->N), deg_ei);
1496  }
1497  else
1498  scAll((currRing->N), deg_ei);
1499  }
1500  }
1501  omFreeSize((ADDRESS)hstc, hNexist * sizeof(scmon));
1502  }
1503 ende:
1505  omFreeSize((ADDRESS)act, ((currRing->N) + 1) * sizeof(int));
1506  hKill(stcmem, (currRing->N) - 1);
1507  pLmFree(&p);
1508  if (p == NULL)
1509  return idInit(1,s->rank);
1510 
1511  last = p;
1512  return scIdKbase(p, s->rank);
1513 }
1514 
1515 #if 0 //-- alternative implementation of scComputeHC
1516 /*
1517 void scComputeHCw(ideal ss, ideal Q, int ak, poly &hEdge)
1518 {
1519  id_LmTest(ss, currRing);
1520  if (Q!=NULL) id_LmTest(Q, currRing);
1521 
1522  int i, di;
1523  poly p;
1524 
1525  if (hEdge!=NULL)
1526  pLmFree(hEdge);
1527 
1528  ideal s=idInit(IDELEMS(ss),ak);
1529  for(i=IDELEMS(ss)-1;i>=0;i--)
1530  {
1531  if (ss->m[i]!=NULL) s->m[i]=pHead(ss->m[i]);
1532  }
1533  di = scDimInt(s, Q);
1534  stcmem = hCreate((currRing->N) - 1);
1535  hexist = hInit(s, Q, &hNexist);
1536  p = last = pInit();
1537  // pNext(p) = NULL;
1538  act = (scmon)omAlloc(((currRing->N) + 1) * sizeof(int));
1539  *act = 0;
1540  if (!hNexist)
1541  {
1542  scAll((currRing->N), -1);
1543  goto ende;
1544  }
1545  if (!hisModule)
1546  {
1547  scInKbase(hexist, hNexist, (currRing->N));
1548  }
1549  else
1550  {
1551  hstc = (scfmon)omAlloc(hNexist * sizeof(scmon));
1552  for (i = 1; i <= hisModule; i++)
1553  {
1554  *act = i;
1555  hComp(hexist, hNexist, i, hstc, &hNstc);
1556  if (hNstc)
1557  {
1558  scInKbase(hstc, hNstc, (currRing->N));
1559  }
1560  else
1561  scAll((currRing->N), -1);
1562  }
1563  omFreeSize((ADDRESS)hstc, hNexist * sizeof(scmon));
1564  }
1565 ende:
1566  hDelete(hexist, hNexist);
1567  omFreeSize((ADDRESS)act, ((currRing->N) + 1) * sizeof(int));
1568  hKill(stcmem, (currRing->N) - 1);
1569  pDeleteLm(&p);
1570  idDelete(&s);
1571  if (p == NULL)
1572  {
1573  return; // no HEdge
1574  }
1575  else
1576  {
1577  last = p;
1578  ideal res=scIdKbase(p, ss->rank);
1579  poly p_ind=res->m[0]; int ind=0;
1580  for(i=IDELEMS(res)-1;i>0;i--)
1581  {
1582  if (pCmp(res->m[i],p_ind)==-1) { p_ind=res->m[i]; ind=i; }
1583  }
1584  assume(p_ind!=NULL);
1585  assume(res->m[ind]==p_ind);
1586  hEdge=p_ind;
1587  res->m[ind]=NULL;
1588  nDelete(&pGetCoeff(hEdge));
1589  pGetCoeff(hEdge)=NULL;
1590  for(i=(currRing->N);i>0;i--)
1591  pIncrExp(hEdge,i);
1592  pSetm(hEdge);
1593 
1594  idDelete(&res);
1595  return;
1596  }
1597 }
1598  */
1599 #endif
1600 
1601 #ifdef HAVE_SHIFTBBA
1602 
1603 /*
1604  * Computation of the Gel'fand-Kirillov Dimension
1605  */
1606 
1607 #include "polys/shiftop.h"
1608 #include <vector>
1609 
1610 static std::vector<int> countCycles(const intvec* _G, int v, std::vector<int> path, std::vector<BOOLEAN> visited, std::vector<BOOLEAN> cyclic, std::vector<int> cache)
1611 {
1612  intvec* G = ivCopy(_G); // modifications must be local
1613 
1614  if (cache[v] != -2) return cache; // value is already cached
1615 
1616  visited[v] = TRUE;
1617  path.push_back(v);
1618 
1619  int cycles = 0;
1620  for (int w = 0; w < G->cols(); w++)
1621  {
1622  if (IMATELEM(*G, v + 1, w + 1)) // edge v -> w exists in G
1623  {
1624  if (!visited[w])
1625  { // continue with w
1626  cache = countCycles(G, w, path, visited, cyclic, cache);
1627  if (cache[w] == -1)
1628  {
1629  cache[v] = -1;
1630  return cache;
1631  }
1632  cycles = si_max(cycles, cache[w]);
1633  }
1634  else
1635  { // found new cycle
1636  int pathIndexOfW = -1;
1637  for (int i = path.size() - 1; i >= 0; i--) {
1638  if (cyclic[path[i]] == 1) { // found an already cyclic vertex
1639  cache[v] = -1;
1640  return cache;
1641  }
1642  cyclic[path[i]] = TRUE;
1643 
1644  if (path[i] == w) { // end of the cycle
1645  assume(IMATELEM(*G, v + 1, w + 1) != 0);
1646  IMATELEM(*G, v + 1, w + 1) = 0; // remove edge v -> w
1647  pathIndexOfW = i;
1648  break;
1649  } else {
1650  assume(IMATELEM(*G, path[i - 1] + 1, path[i] + 1) != 0);
1651  IMATELEM(*G, path[i - 1] + 1, path[i] + 1) = 0; // remove edge vi-1 -> vi
1652  }
1653  }
1654  assume(pathIndexOfW != -1); // should never happen
1655  for (int i = path.size() - 1; i >= pathIndexOfW; i--) {
1656  cache = countCycles(G, path[i], path, visited, cyclic, cache);
1657  if (cache[path[i]] == -1)
1658  {
1659  cache[v] = -1;
1660  return cache;
1661  }
1662  cycles = si_max(cycles, cache[path[i]] + 1);
1663  }
1664  }
1665  }
1666  }
1667  cache[v] = cycles;
1668 
1669  delete G;
1670  return cache;
1671 }
1672 
1673 // -1 is infinity
1674 static int graphGrowth(const intvec* G)
1675 {
1676  // init
1677  int n = G->cols();
1678  std::vector<int> path;
1679  std::vector<BOOLEAN> visited;
1680  std::vector<BOOLEAN> cyclic;
1681  std::vector<int> cache;
1682  visited.resize(n, FALSE);
1683  cyclic.resize(n, FALSE);
1684  cache.resize(n, -2);
1685 
1686  // get max number of cycles
1687  int cycles = 0;
1688  for (int v = 0; v < n; v++)
1689  {
1690  cache = countCycles(G, v, path, visited, cyclic, cache);
1691  if (cache[v] == -1)
1692  return -1;
1693  cycles = si_max(cycles, cache[v]);
1694  }
1695  return cycles;
1696 }
1697 
1698 // ATTENTION:
1699 // - `words` contains the words normal modulo M of length n
1700 // - `numberOfNormalWords` contains the number of words normal modulo M of length 0 ... n
1701 static void _lp_computeNormalWords(ideal words, int& numberOfNormalWords, int length, ideal M, int minDeg, int& last)
1702 {
1703  if (length <= 0){
1704  poly one = pOne();
1705  if (p_LPDivisibleBy(M, one, currRing)) // 1 \in M => no normal words at all
1706  {
1707  pDelete(&one);
1708  last = -1;
1709  numberOfNormalWords = 0;
1710  }
1711  else
1712  {
1713  words->m[0] = one;
1714  last = 0;
1715  numberOfNormalWords = 1;
1716  }
1717  return;
1718  }
1719 
1720  _lp_computeNormalWords(words, numberOfNormalWords, length - 1, M, minDeg, last);
1721 
1722  int nVars = currRing->isLPring - currRing->LPncGenCount;
1723  int numberOfNewNormalWords = 0;
1724 
1725  for (int j = nVars - 1; j >= 0; j--)
1726  {
1727  for (int i = last; i >= 0; i--)
1728  {
1729  int index = (j * (last + 1)) + i;
1730 
1731  if (words->m[i] != NULL)
1732  {
1733  if (j > 0) {
1734  words->m[index] = pCopy(words->m[i]);
1735  }
1736 
1737  int varOffset = ((length - 1) * currRing->isLPring) + 1;
1738  pSetExp(words->m[index], varOffset + j, 1);
1739  pSetm(words->m[index]);
1740  pTest(words->m[index]);
1741 
1742  if (length >= minDeg && p_LPDivisibleBy(M, words->m[index], currRing))
1743  {
1744  pDelete(&words->m[index]);
1745  words->m[index] = NULL;
1746  }
1747  else
1748  {
1749  numberOfNewNormalWords++;
1750  }
1751  }
1752  }
1753  }
1754 
1755  last = nVars * last + nVars - 1;
1756 
1757  numberOfNormalWords += numberOfNewNormalWords;
1758 }
1759 
1760 static ideal lp_computeNormalWords(int length, ideal M)
1761 {
1762  long minDeg = IDELEMS(M) > 0 ? pTotaldegree(M->m[0]) : 0;
1763  for (int i = 1; i < IDELEMS(M); i++)
1764  {
1765  minDeg = si_min(minDeg, pTotaldegree(M->m[i]));
1766  }
1767 
1768  int nVars = currRing->isLPring - currRing->LPncGenCount;
1769 
1770  int maxElems = 1;
1771  for (int i = 0; i < length; i++) // maxElems = nVars^n
1772  maxElems *= nVars;
1773  ideal words = idInit(maxElems);
1774  int last, numberOfNormalWords;
1775  _lp_computeNormalWords(words, numberOfNormalWords, length, M, minDeg, last);
1776  idSkipZeroes(words);
1777  return words;
1778 }
1779 
1780 static int lp_countNormalWords(int upToLength, ideal M)
1781 {
1782  long minDeg = IDELEMS(M) > 0 ? pTotaldegree(M->m[0]) : 0;
1783  for (int i = 1; i < IDELEMS(M); i++)
1784  {
1785  minDeg = si_min(minDeg, pTotaldegree(M->m[i]));
1786  }
1787 
1788  int nVars = currRing->isLPring - currRing->LPncGenCount;
1789 
1790  int maxElems = 1;
1791  for (int i = 0; i < upToLength; i++) // maxElems = nVars^n
1792  maxElems *= nVars;
1793  ideal words = idInit(maxElems);
1794  int last, numberOfNormalWords;
1795  _lp_computeNormalWords(words, numberOfNormalWords, upToLength, M, minDeg, last);
1796  idDelete(&words);
1797  return numberOfNormalWords;
1798 }
1799 
1800 // NULL if graph is undefined
1801 intvec* lp_ufnarovskiGraph(ideal G, ideal &standardWords)
1802 {
1803  long l = 0;
1804  for (int i = 0; i < IDELEMS(G); i++)
1805  l = si_max(pTotaldegree(G->m[i]), l);
1806  l--;
1807  if (l <= 0)
1808  {
1809  WerrorS("Ufnarovski graph not implemented for l <= 0");
1810  return NULL;
1811  }
1812  int lV = currRing->isLPring;
1813 
1814  standardWords = lp_computeNormalWords(l, G);
1815 
1816  int n = IDELEMS(standardWords);
1817  intvec* UG = new intvec(n, n, 0);
1818  for (int i = 0; i < n; i++)
1819  {
1820  for (int j = 0; j < n; j++)
1821  {
1822  poly v = standardWords->m[i];
1823  poly w = standardWords->m[j];
1824 
1825  // check whether v*x1 = x2*w (overlap)
1826  bool overlap = true;
1827  for (int k = 1; k <= (l - 1) * lV; k++)
1828  {
1829  if (pGetExp(v, k + lV) != pGetExp(w, k)) {
1830  overlap = false;
1831  break;
1832  }
1833  }
1834 
1835  if (overlap)
1836  {
1837  // create the overlap
1838  poly p = pMult(pCopy(v), p_LPVarAt(w, l, currRing));
1839 
1840  // check whether the overlap is normal
1841  bool normal = true;
1842  for (int k = 0; k < IDELEMS(G); k++)
1843  {
1844  if (p_LPDivisibleBy(G->m[k], p, currRing))
1845  {
1846  normal = false;
1847  break;
1848  }
1849  }
1850 
1851  if (normal)
1852  {
1853  IMATELEM(*UG, i + 1, j + 1) = 1;
1854  }
1855  }
1856  }
1857  }
1858  return UG;
1859 }
1860 
1861 // -1 is infinity, -2 is error
1862 int lp_gkDim(const ideal _G)
1863 {
1864  id_Test(_G, currRing);
1865 
1866  if (rField_is_Ring(currRing)) {
1867  WerrorS("GK-Dim not implemented for rings");
1868  return -2;
1869  }
1870 
1871  for (int i=IDELEMS(_G)-1;i>=0; i--)
1872  {
1873  if (_G->m[i] != NULL)
1874  {
1875  if (pGetComp(_G->m[i]) != 0)
1876  {
1877  WerrorS("GK-Dim not implemented for modules");
1878  return -2;
1879  }
1880  if (pGetNCGen(_G->m[i]) != 0)
1881  {
1882  WerrorS("GK-Dim not implemented for bi-modules");
1883  return -2;
1884  }
1885  }
1886  }
1887 
1888  ideal G = id_Head(_G, currRing); // G = LM(G) (and copy)
1889  idSkipZeroes(G); // remove zeros
1890  id_DelLmEquals(G, currRing); // remove duplicates
1891 
1892  // check if G is the zero ideal
1893  if (IDELEMS(G) == 1 && G->m[0] == NULL)
1894  {
1895  // NOTE: this is needed because if the ideal is <0>, then idSkipZeroes keeps this element, and IDELEMS is still 1!
1896  int lV = currRing->isLPring;
1897  int ncGenCount = currRing->LPncGenCount;
1898  if (lV - ncGenCount == 0)
1899  {
1900  idDelete(&G);
1901  return 0;
1902  }
1903  if (lV - ncGenCount == 1)
1904  {
1905  idDelete(&G);
1906  return 1;
1907  }
1908  if (lV - ncGenCount >= 2)
1909  {
1910  idDelete(&G);
1911  return -1;
1912  }
1913  }
1914 
1915  // get the max deg
1916  long maxDeg = 0;
1917  for (int i = 0; i < IDELEMS(G); i++)
1918  {
1919  maxDeg = si_max(maxDeg, pTotaldegree(G->m[i]));
1920 
1921  // also check whether G = <1>
1922  if (pIsConstantComp(G->m[i]))
1923  {
1924  WerrorS("GK-Dim not defined for 0-ring");
1925  idDelete(&G);
1926  return -2;
1927  }
1928  }
1929 
1930  // early termination if G \subset X
1931  if (maxDeg <= 1)
1932  {
1933  int lV = currRing->isLPring;
1934  int ncGenCount = currRing->LPncGenCount;
1935  if (IDELEMS(G) == lV - ncGenCount) // V = {1} no edges
1936  {
1937  idDelete(&G);
1938  return 0;
1939  }
1940  if (IDELEMS(G) == lV - ncGenCount - 1) // V = {1} with loop
1941  {
1942  idDelete(&G);
1943  return 1;
1944  }
1945  if (IDELEMS(G) <= lV - ncGenCount - 2) // V = {1} with more than one loop
1946  {
1947  idDelete(&G);
1948  return -1;
1949  }
1950  }
1951 
1952  ideal standardWords;
1953  intvec* UG = lp_ufnarovskiGraph(G, standardWords);
1954  if (UG == NULL)
1955  {
1956  idDelete(&G);
1957  return -2;
1958  }
1959  if (errorreported)
1960  {
1961  delete UG;
1962  idDelete(&G);
1963  return -2;
1964  }
1965  int gkDim = graphGrowth(UG);
1966  delete UG;
1967  idDelete(&G);
1968  return gkDim;
1969 }
1970 
1971 // converts an intvec matrix to a vector<vector<int> >
1972 static std::vector<std::vector<int> > iv2vv(intvec* M)
1973 {
1974  int rows = M->rows();
1975  int cols = M->cols();
1976 
1977  std::vector<std::vector<int> > mat(rows, std::vector<int>(cols));
1978 
1979  for (int i = 0; i < rows; i++)
1980  {
1981  for (int j = 0; j < cols; j++)
1982  {
1983  mat[i][j] = IMATELEM(*M, i + 1, j + 1);
1984  }
1985  }
1986 
1987  return mat;
1988 }
1989 
1990 static void vvPrint(const std::vector<std::vector<int> >& mat)
1991 {
1992  for (int i = 0; i < mat.size(); i++)
1993  {
1994  for (int j = 0; j < mat[i].size(); j++)
1995  {
1996  Print("%d ", mat[i][j]);
1997  }
1998  PrintLn();
1999  }
2000 }
2001 
2002 static void vvTest(const std::vector<std::vector<int> >& mat)
2003 {
2004  if (mat.size() > 0)
2005  {
2006  int cols = mat[0].size();
2007  for (int i = 1; i < mat.size(); i++)
2008  {
2009  if (cols != mat[i].size())
2010  WerrorS("number of cols in matrix inconsistent");
2011  }
2012  }
2013 }
2014 
2015 static void vvDeleteRow(std::vector<std::vector<int> >& mat, int row)
2016 {
2017  mat.erase(mat.begin() + row);
2018 }
2019 
2020 static void vvDeleteColumn(std::vector<std::vector<int> >& mat, int col)
2021 {
2022  for (int i = 0; i < mat.size(); i++)
2023  {
2024  mat[i].erase(mat[i].begin() + col);
2025  }
2026 }
2027 
2028 static BOOLEAN vvIsRowZero(const std::vector<std::vector<int> >& mat, int row)
2029 {
2030  for (int i = 0; i < mat[row].size(); i++)
2031  {
2032  if (mat[row][i] != 0)
2033  return FALSE;
2034  }
2035  return TRUE;
2036 }
2037 
2038 static BOOLEAN vvIsColumnZero(const std::vector<std::vector<int> >& mat, int col)
2039 {
2040  for (int i = 0; i < mat.size(); i++)
2041  {
2042  if (mat[i][col] != 0)
2043  return FALSE;
2044  }
2045  return TRUE;
2046 }
2047 
2048 static BOOLEAN vvIsZero(const std::vector<std::vector<int> >& mat)
2049 {
2050  for (int i = 0; i < mat.size(); i++)
2051  {
2052  if (!vvIsRowZero(mat, i))
2053  return FALSE;
2054  }
2055  return TRUE;
2056 }
2057 
2058 static std::vector<std::vector<int> > vvMult(const std::vector<std::vector<int> >& a, const std::vector<std::vector<int> >& b)
2059 {
2060  int ra = a.size();
2061  int rb = b.size();
2062  int ca = a.size() > 0 ? a[0].size() : 0;
2063  int cb = b.size() > 0 ? b[0].size() : 0;
2064 
2065  if (ca != rb)
2066  {
2067  WerrorS("matrix dimensions do not match");
2068  return std::vector<std::vector<int> >();
2069  }
2070 
2071  std::vector<std::vector<int> > res(ra, std::vector<int>(cb));
2072  for (int i = 0; i < ra; i++)
2073  {
2074  for (int j = 0; j < cb; j++)
2075  {
2076  int sum = 0;
2077  for (int k = 0; k < ca; k++)
2078  sum += a[i][k] * b[k][j];
2079  res[i][j] = sum;
2080  }
2081  }
2082  return res;
2083 }
2084 
2085 static BOOLEAN isAcyclic(const intvec* G)
2086 {
2087  // init
2088  int n = G->cols();
2089  std::vector<int> path;
2090  std::vector<BOOLEAN> visited;
2091  std::vector<BOOLEAN> cyclic;
2092  std::vector<int> cache;
2093  visited.resize(n, FALSE);
2094  cyclic.resize(n, FALSE);
2095  cache.resize(n, -2);
2096 
2097  for (int v = 0; v < n; v++)
2098  {
2099  cache = countCycles(G, v, path, visited, cyclic, cache);
2100  // check that there are 0 cycles from v
2101  if (cache[v] != 0)
2102  return FALSE;
2103  }
2104  return TRUE;
2105 }
2106 
2107 /*
2108  * Computation of the K-Dimension
2109  */
2110 
2111 // -1 is infinity, -2 is error
2112 int lp_kDim(const ideal _G)
2113 {
2114  if (rField_is_Ring(currRing)) {
2115  WerrorS("K-Dim not implemented for rings");
2116  return -2;
2117  }
2118 
2119  for (int i=IDELEMS(_G)-1;i>=0; i--)
2120  {
2121  if (_G->m[i] != NULL)
2122  {
2123  if (pGetComp(_G->m[i]) != 0)
2124  {
2125  WerrorS("K-Dim not implemented for modules");
2126  return -2;
2127  }
2128  if (pGetNCGen(_G->m[i]) != 0)
2129  {
2130  WerrorS("K-Dim not implemented for bi-modules");
2131  return -2;
2132  }
2133  }
2134  }
2135 
2136  ideal G = id_Head(_G, currRing); // G = LM(G) (and copy)
2137  if (TEST_OPT_PROT)
2138  Print("%d original generators\n", IDELEMS(G));
2139  idSkipZeroes(G); // remove zeros
2140  id_DelLmEquals(G, currRing); // remove duplicates
2141  if (TEST_OPT_PROT)
2142  Print("%d non-zero unique generators\n", IDELEMS(G));
2143 
2144  // check if G is the zero ideal
2145  if (IDELEMS(G) == 1 && G->m[0] == NULL)
2146  {
2147  // NOTE: this is needed because if the ideal is <0>, then idSkipZeroes keeps this element, and IDELEMS is still 1!
2148  int lV = currRing->isLPring;
2149  int ncGenCount = currRing->LPncGenCount;
2150  if (lV - ncGenCount == 0)
2151  {
2152  idDelete(&G);
2153  return 1;
2154  }
2155  if (lV - ncGenCount == 1)
2156  {
2157  idDelete(&G);
2158  return -1;
2159  }
2160  if (lV - ncGenCount >= 2)
2161  {
2162  idDelete(&G);
2163  return -1;
2164  }
2165  }
2166 
2167  // get the max deg
2168  long maxDeg = 0;
2169  for (int i = 0; i < IDELEMS(G); i++)
2170  {
2171  maxDeg = si_max(maxDeg, pTotaldegree(G->m[i]));
2172 
2173  // also check whether G = <1>
2174  if (pIsConstantComp(G->m[i]))
2175  {
2176  WerrorS("K-Dim not defined for 0-ring"); // TODO is it minus infinity ?
2177  idDelete(&G);
2178  return -2;
2179  }
2180  }
2181  if (TEST_OPT_PROT)
2182  Print("max deg: %ld\n", maxDeg);
2183 
2184 
2185  // for normal words of length minDeg ... maxDeg-1
2186  // brute-force the normal words
2187  if (TEST_OPT_PROT)
2188  PrintS("Computing normal words normally...\n");
2189  long numberOfNormalWords = lp_countNormalWords(maxDeg - 1, G);
2190 
2191  if (TEST_OPT_PROT)
2192  Print("%ld normal words up to length %ld\n", numberOfNormalWords, maxDeg - 1);
2193 
2194  // early termination if G \subset X
2195  if (maxDeg <= 1)
2196  {
2197  int lV = currRing->isLPring;
2198  int ncGenCount = currRing->LPncGenCount;
2199  if (IDELEMS(G) == lV - ncGenCount) // V = {1} no edges
2200  {
2201  idDelete(&G);
2202  return numberOfNormalWords;
2203  }
2204  if (IDELEMS(G) == lV - ncGenCount - 1) // V = {1} with loop
2205  {
2206  idDelete(&G);
2207  return -1;
2208  }
2209  if (IDELEMS(G) <= lV - ncGenCount - 2) // V = {1} with more than one loop
2210  {
2211  idDelete(&G);
2212  return -1;
2213  }
2214  }
2215 
2216  if (TEST_OPT_PROT)
2217  PrintS("Computing Ufnarovski graph...\n");
2218 
2219  ideal standardWords;
2220  intvec* UG = lp_ufnarovskiGraph(G, standardWords);
2221  if (UG == NULL)
2222  {
2223  idDelete(&G);
2224  return -2;
2225  }
2226  if (errorreported)
2227  {
2228  delete UG;
2229  idDelete(&G);
2230  return -2;
2231  }
2232 
2233  if (TEST_OPT_PROT)
2234  Print("Ufnarovski graph is %dx%d.\n", UG->rows(), UG->cols());
2235 
2236  if (TEST_OPT_PROT)
2237  PrintS("Checking whether Ufnarovski graph is acyclic...\n");
2238 
2239  if (!isAcyclic(UG))
2240  {
2241  // in this case we have infinitely many normal words
2242  return -1;
2243  }
2244 
2245  std::vector<std::vector<int> > vvUG = iv2vv(UG);
2246  for (int i = 0; i < vvUG.size(); i++)
2247  {
2248  if (vvIsRowZero(vvUG, i) && vvIsColumnZero(vvUG, i)) // i is isolated vertex
2249  {
2250  vvDeleteRow(vvUG, i);
2251  vvDeleteColumn(vvUG, i);
2252  i--;
2253  }
2254  }
2255  if (TEST_OPT_PROT)
2256  Print("Simplified Ufnarovski graph to %dx%d.\n", (int)vvUG.size(), (int)vvUG.size());
2257 
2258  // for normal words of length >= maxDeg
2259  // use Ufnarovski graph
2260  if (TEST_OPT_PROT)
2261  PrintS("Computing normal words via Ufnarovski graph...\n");
2262  std::vector<std::vector<int> > UGpower = vvUG;
2263  long nUGpower = 1;
2264  while (!vvIsZero(UGpower))
2265  {
2266  if (TEST_OPT_PROT)
2267  PrintS("Start count graph entries.\n");
2268  for (int i = 0; i < UGpower.size(); i++)
2269  {
2270  for (int j = 0; j < UGpower[i].size(); j++)
2271  {
2272  numberOfNormalWords += UGpower[i][j];
2273  }
2274  }
2275 
2276  if (TEST_OPT_PROT)
2277  {
2278  PrintS("Done count graph entries.\n");
2279  Print("%ld normal words up to length %ld\n", numberOfNormalWords, maxDeg - 1 + nUGpower);
2280  }
2281 
2282  if (TEST_OPT_PROT)
2283  PrintS("Start mat mult.\n");
2284  UGpower = vvMult(UGpower, vvUG); // TODO: avoid creation of new intvec
2285  if (TEST_OPT_PROT)
2286  PrintS("Done mat mult.\n");
2287  nUGpower++;
2288  }
2289 
2290  delete UG;
2291  idDelete(&G);
2292  return numberOfNormalWords;
2293 }
2294 #endif
long int64
Definition: auxiliary.h:68
static int si_max(const int a, const int b)
Definition: auxiliary.h:124
int BOOLEAN
Definition: auxiliary.h:87
#define TRUE
Definition: auxiliary.h:100
#define FALSE
Definition: auxiliary.h:96
void * ADDRESS
Definition: auxiliary.h:119
static int si_min(const int a, const int b)
Definition: auxiliary.h:125
int size(const CanonicalForm &f, const Variable &v)
int size ( const CanonicalForm & f, const Variable & v )
Definition: cf_ops.cc:600
int l
Definition: cfEzgcd.cc:100
int i
Definition: cfEzgcd.cc:132
int k
Definition: cfEzgcd.cc:99
Variable x
Definition: cfModGcd.cc:4082
int p
Definition: cfModGcd.cc:4078
CanonicalForm b
Definition: cfModGcd.cc:4103
void mu(int **points, int sizePoints)
Definition: intvec.h:23
int length() const
Definition: intvec.h:94
int cols() const
Definition: intvec.h:95
int rows() const
Definition: intvec.h:96
static FORCE_INLINE BOOLEAN n_IsUnit(number n, const coeffs r)
TRUE iff n has a multiplicative inverse in the given coeff field/ring r.
Definition: coeffs.h:515
static FORCE_INLINE BOOLEAN n_DivBy(number a, number b, const coeffs r)
test whether 'a' is divisible 'b'; for r encoding a field: TRUE iff 'b' does not represent zero in Z:...
Definition: coeffs.h:753
#define Print
Definition: emacs.cc:80
const CanonicalForm int s
Definition: facAbsFact.cc:51
const CanonicalForm int const CFList const Variable & y
Definition: facAbsFact.cc:53
CanonicalForm res
Definition: facAbsFact.cc:60
const CanonicalForm & w
Definition: facAbsFact.cc:51
const Variable & v
< [in] a sqrfree bivariate poly
Definition: facBivar.h:39
int j
Definition: facHensel.cc:110
VAR short errorreported
Definition: feFopen.cc:23
void WerrorS(const char *s)
Definition: feFopen.cc:24
#define STATIC_VAR
Definition: globaldefs.h:7
#define VAR
Definition: globaldefs.h:5
static long hZeroMult(scmon pure, scfmon stc, int Nstc, varset var, int Nvar)
Definition: hdegree.cc:621
static ideal lp_computeNormalWords(int length, ideal M)
Definition: hdegree.cc:1760
void scComputeHC(ideal S, ideal Q, int ak, poly &hEdge)
Definition: hdegree.cc:1101
STATIC_VAR scmon hInd
Definition: hdegree.cc:205
static void hHedgeStep(scmon pure, scfmon stc, int Nstc, varset var, int Nvar, poly hEdge)
Definition: hdegree.cc:1041
static void hDimMult(scmon pure, int Npure, scfmon rad, int Nrad, varset var, int Nvar)
Definition: hdegree.cc:726
static std::vector< std::vector< int > > iv2vv(intvec *M)
Definition: hdegree.cc:1972
ideal scKBase(int deg, ideal s, ideal Q, intvec *mv)
Definition: hdegree.cc:1449
int scDimIntRing(ideal vid, ideal Q)
scDimInt for ring-coefficients
Definition: hdegree.cc:136
long scMult0Int(ideal S, ideal Q)
Definition: hdegree.cc:950
void hIndMult(scmon pure, int Npure, scfmon rad, int Nrad, varset var, int Nvar)
Definition: hdegree.cc:384
intvec * lp_ufnarovskiGraph(ideal G, ideal &standardWords)
Definition: hdegree.cc:1801
intvec * scIndIntvec(ideal S, ideal Q)
Definition: hdegree.cc:286
static int scMin(int i, scfmon stc, int Nvar)
Definition: hdegree.cc:1197
static void vvDeleteRow(std::vector< std::vector< int > > &mat, int row)
Definition: hdegree.cc:2015
static indset hCheck2(indset sm, scmon pure)
Definition: hdegree.cc:491
STATIC_VAR poly last
Definition: hdegree.cc:1173
static BOOLEAN hCheck1(indset sm, scmon pure)
Definition: hdegree.cc:465
static int graphGrowth(const intvec *G)
Definition: hdegree.cc:1674
static BOOLEAN vvIsColumnZero(const std::vector< std::vector< int > > &mat, int col)
Definition: hdegree.cc:2038
VAR omBin indlist_bin
Definition: hdegree.cc:29
STATIC_VAR poly pWork
Definition: hdegree.cc:1027
VAR int hMu2
Definition: hdegree.cc:27
static void hDegree(ideal S, ideal Q)
Definition: hdegree.cc:802
static void vvDeleteColumn(std::vector< std::vector< int > > &mat, int col)
Definition: hdegree.cc:2020
static BOOLEAN hNotZero(scfmon rad, int Nrad, varset var, int Nvar)
Definition: hdegree.cc:355
int lp_kDim(const ideal _G)
Definition: hdegree.cc:2112
static void scElKbase()
Definition: hdegree.cc:1176
static void hHedge(poly hEdge)
Definition: hdegree.cc:1029
static void hIndSolve(scmon pure, int Npure, scfmon rad, int Nrad, varset var, int Nvar)
Definition: hdegree.cc:207
VAR int hCo
Definition: hdegree.cc:27
static int scRestrict(int &Nstc, scfmon stc, int Nvar)
Definition: hdegree.cc:1209
int lp_gkDim(const ideal _G)
Definition: hdegree.cc:1862
VAR indset ISet
Definition: hdegree.cc:353
static void vvPrint(const std::vector< std::vector< int > > &mat)
Definition: hdegree.cc:1990
static void vvTest(const std::vector< std::vector< int > > &mat)
Definition: hdegree.cc:2002
static void scAllKbase(int Nvar, int ideg, int deg)
Definition: hdegree.cc:1284
VAR long hMu
Definition: hdegree.cc:28
static void scAll(int Nvar, int deg)
Definition: hdegree.cc:1260
int scMultInt(ideal S, ideal Q)
Definition: hdegree.cc:903
static void scDegKbase(scfmon stc, int Nstc, int Nvar, int deg)
Definition: hdegree.cc:1294
STATIC_VAR scmon act
Definition: hdegree.cc:1174
static void hCheckIndep(scmon pure)
Definition: hdegree.cc:543
void scPrintDegree(int co, int mu)
Definition: hdegree.cc:912
VAR indset JSet
Definition: hdegree.cc:353
static int lp_countNormalWords(int upToLength, ideal M)
Definition: hdegree.cc:1780
static BOOLEAN isAcyclic(const intvec *G)
Definition: hdegree.cc:2085
static int scMax(int i, scfmon stc, int Nvar)
Definition: hdegree.cc:1185
static ideal scIdKbase(poly q, const int rank)
Definition: hdegree.cc:1431
static void hIndep(scmon pure)
Definition: hdegree.cc:370
static void scInKbase(scfmon stc, int Nstc, int Nvar)
Definition: hdegree.cc:1375
static void hProject(scmon pure, varset sel)
Definition: hdegree.cc:703
void scDegree(ideal S, intvec *modulweight, ideal Q)
Definition: hdegree.cc:926
static BOOLEAN vvIsZero(const std::vector< std::vector< int > > &mat)
Definition: hdegree.cc:2048
int scDimInt(ideal S, ideal Q)
ideal dimension
Definition: hdegree.cc:78
static BOOLEAN vvIsRowZero(const std::vector< std::vector< int > > &mat, int row)
Definition: hdegree.cc:2028
static std::vector< std::vector< int > > vvMult(const std::vector< std::vector< int > > &a, const std::vector< std::vector< int > > &b)
Definition: hdegree.cc:2058
static std::vector< int > countCycles(const intvec *_G, int v, std::vector< int > path, std::vector< BOOLEAN > visited, std::vector< BOOLEAN > cyclic, std::vector< int > cache)
Definition: hdegree.cc:1610
static void _lp_computeNormalWords(ideal words, int &numberOfNormalWords, int length, ideal M, int minDeg, int &last)
Definition: hdegree.cc:1701
void hDimSolve(scmon pure, int Npure, scfmon rad, int Nrad, varset var, int Nvar)
Definition: hdegree.cc:35
void hIndAllMult(scmon pure, int Npure, scfmon rad, int Nrad, varset var, int Nvar)
Definition: hdegree.cc:564
intvec * hFirstSeries(ideal A, intvec *module_w, ideal Q, intvec *wdegree)
Definition: hilb.cc:1947
void hDegreeSeries(intvec *s1, intvec *s2, int *co, int *mu)
Definition: hilb.cc:732
intvec * hSecondSeries(intvec *hseries1)
Definition: hilb.cc:697
monf hCreate(int Nvar)
Definition: hutil.cc:996
void hComp(scfmon exist, int Nexist, int ak, scfmon stc, int *Nstc)
Definition: hutil.cc:154
void hLex2S(scfmon rad, int e1, int a2, int e2, varset var, int Nvar, scfmon w)
Definition: hutil.cc:812
VAR scfmon hstc
Definition: hutil.cc:16
VAR varset hvar
Definition: hutil.cc:18
void hKill(monf xmem, int Nvar)
Definition: hutil.cc:1010
VAR int hNexist
Definition: hutil.cc:19
void hElimS(scfmon stc, int *e1, int a2, int e2, varset var, int Nvar)
Definition: hutil.cc:672
void hLexS(scfmon stc, int Nstc, varset var, int Nvar)
Definition: hutil.cc:506
void hDelete(scfmon ev, int ev_length)
Definition: hutil.cc:140
VAR scmon hpur0
Definition: hutil.cc:17
VAR monf stcmem
Definition: hutil.cc:21
scfmon hGetmem(int lm, scfmon old, monp monmem)
Definition: hutil.cc:1023
void hPure(scfmon stc, int a, int *Nstc, varset var, int Nvar, scmon pure, int *Npure)
Definition: hutil.cc:621
VAR scfmon hwork
Definition: hutil.cc:16
void hSupp(scfmon stc, int Nstc, varset var, int *Nvar)
Definition: hutil.cc:174
void hLexR(scfmon rad, int Nrad, varset var, int Nvar)
Definition: hutil.cc:565
VAR scmon hpure
Definition: hutil.cc:17
void hStepR(scfmon rad, int Nrad, varset var, int Nvar, int *a)
Definition: hutil.cc:974
void hLex2R(scfmon rad, int e1, int a2, int e2, varset var, int Nvar, scfmon w)
Definition: hutil.cc:880
VAR scfmon hrad
Definition: hutil.cc:16
VAR int hisModule
Definition: hutil.cc:20
void hStepS(scfmon stc, int Nstc, varset var, int Nvar, int *a, int *x)
Definition: hutil.cc:949
void hStaircase(scfmon stc, int *Nstc, varset var, int Nvar)
Definition: hutil.cc:313
void hElimR(scfmon rad, int *e1, int a2, int e2, varset var, int Nvar)
Definition: hutil.cc:742
VAR monf radmem
Definition: hutil.cc:21
void hOrdSupp(scfmon stc, int Nstc, varset var, int Nvar)
Definition: hutil.cc:202
VAR varset hsel
Definition: hutil.cc:18
VAR int hNpure
Definition: hutil.cc:19
VAR int hNrad
Definition: hutil.cc:19
scfmon hInit(ideal S, ideal Q, int *Nexist)
Definition: hutil.cc:31
VAR scfmon hexist
Definition: hutil.cc:16
void hRadical(scfmon rad, int *Nrad, int Nvar)
Definition: hutil.cc:411
scmon hGetpure(scmon p)
Definition: hutil.cc:1052
VAR int hNstc
Definition: hutil.cc:19
VAR int hNvar
Definition: hutil.cc:19
scmon * scfmon
Definition: hutil.h:15
indlist * indset
Definition: hutil.h:28
int * varset
Definition: hutil.h:16
int * scmon
Definition: hutil.h:14
#define idDelete(H)
delete an ideal
Definition: ideals.h:29
BOOLEAN idInsertPoly(ideal h1, poly h2)
insert h2 into h1 (if h2 is not the zero polynomial) return TRUE iff h2 was indeed inserted
ideal id_Copy(ideal h1, const ring r)
copy an ideal
ideal idCopy(ideal A)
Definition: ideals.h:60
#define idPosConstant(I)
index of generator with leading term in ground ring (if any); otherwise -1
Definition: ideals.h:37
static BOOLEAN length(leftv result, leftv arg)
Definition: interval.cc:257
#define IMATELEM(M, I, J)
Definition: intvec.h:85
intvec * ivCopy(const intvec *o)
Definition: intvec.h:145
STATIC_VAR TreeM * G
Definition: janet.cc:31
STATIC_VAR jList * Q
Definition: janet.cc:30
#define assume(x)
Definition: mod2.h:389
#define pNext(p)
Definition: monomials.h:36
static number & pGetCoeff(poly p)
return an alias to the leading coefficient of p assumes that p != NULL NOTE: not copy
Definition: monomials.h:44
#define pSetCoeff0(p, n)
Definition: monomials.h:59
const int MAX_INT_VAL
Definition: mylimits.h:12
#define nCopy(n)
Definition: numbers.h:15
#define nInit(i)
Definition: numbers.h:24
#define omFreeSize(addr, size)
Definition: omAllocDecl.h:260
#define omAlloc(size)
Definition: omAllocDecl.h:210
#define omAlloc0Bin(bin)
Definition: omAllocDecl.h:206
#define omAlloc0(size)
Definition: omAllocDecl.h:211
#define omFreeBin(addr, bin)
Definition: omAllocDecl.h:259
#define omGetSpecBin(size)
Definition: omBin.h:11
#define NULL
Definition: omList.c:12
omBin_t * omBin
Definition: omStructs.h:12
#define TEST_OPT_PROT
Definition: options.h:104
static int index(p_Length length, p_Ord ord)
Definition: p_Procs_Impl.h:592
int p_IsPurePower(const poly p, const ring r)
return i, if head depends only on var(i)
Definition: p_polys.cc:1226
static void p_Delete(poly *p, const ring r)
Definition: p_polys.h:903
static unsigned pLength(poly a)
Definition: p_polys.h:191
VAR ring currRing
Widely used global variable which specifies the current polynomial ring for Singular interpreter and ...
Definition: polys.cc:13
Compatiblity layer for legacy polynomial operations (over currRing)
static long pTotaldegree(poly p)
Definition: polys.h:282
#define pTest(p)
Definition: polys.h:415
#define pDelete(p_ptr)
Definition: polys.h:186
#define pSetm(p)
Definition: polys.h:271
#define pGetComp(p)
Component.
Definition: polys.h:37
#define pIsConstantComp(p)
return true if p is either NULL, or if all exponents of p are 0, Comp of p might be !...
Definition: polys.h:236
#define pSetExpV(p, e)
Definition: polys.h:97
#define pSetComp(p, v)
Definition: polys.h:38
#define pMult(p, q)
Definition: polys.h:207
static void pLmFree(poly p)
frees the space of the monomial m, assumes m != NULL coef is not freed, m is not advanced
Definition: polys.h:70
void pWrite(poly p)
Definition: polys.h:308
#define pGetExp(p, i)
Exponent.
Definition: polys.h:41
#define pInit()
allocates a new monomial and initializes everything to 0
Definition: polys.h:61
#define pSetExp(p, i, v)
Definition: polys.h:42
#define pLmCmp(p, q)
returns 0|1|-1 if p=q|p>q|p<q w.r.t monomial ordering
Definition: polys.h:105
#define pCopy(p)
return a copy of the poly
Definition: polys.h:185
#define pOne()
Definition: polys.h:315
poly * polyset
Definition: polys.h:259
void PrintS(const char *s)
Definition: reporter.cc:284
void PrintLn()
Definition: reporter.cc:310
static BOOLEAN rField_is_Z(const ring r)
Definition: ring.h:510
#define rField_is_Ring(R)
Definition: ring.h:486
BOOLEAN p_LPDivisibleBy(poly a, poly b, const ring r)
Definition: shiftop.cc:776
poly p_LPVarAt(poly p, int pos, const ring r)
Definition: shiftop.cc:845
#define pGetNCGen(p)
Definition: shiftop.h:65
ideal idInit(int idsize, int rank)
initialise an ideal / module
Definition: simpleideals.cc:35
int idElem(const ideal F)
count non-zero elements
ideal id_Head(ideal h, const ring r)
returns the ideals of initial terms
void id_DelLmEquals(ideal id, const ring r)
Delete id[j], if Lm(j) == Lm(i) and both LC(j), LC(i) are units and j > i.
void idSkipZeroes(ideal ide)
gives an ideal/module the minimal possible size
#define IDELEMS(i)
Definition: simpleideals.h:23
#define id_Test(A, lR)
Definition: simpleideals.h:78
#define id_LmTest(A, lR)
Definition: simpleideals.h:79
#define M
Definition: sirandom.c:25
#define loop
Definition: structs.h:75