walkSupport.cc
Go to the documentation of this file.
1 #include <kernel/mod2.h>
2 
3 #include <misc/intvec.h>
4 #include <misc/int64vec.h>
5 
6 #include <polys/monomials/ring.h>
7 #include <polys/prCopy.h>
8 #include <polys/matpol.h>
9 
10 #include <kernel/polys.h>
11 #include <kernel/ideals.h>
13 #include <kernel/GBEngine/kstd1.h>
14 
15 #include <string.h>
16 #include <math.h>
17 
18 extern BOOLEAN overflow_error;
19 
20 ///////////////////////////////////////////////////////////////////
21 //Support functions for Groebner Walk and Fractal Walk
22 ///////////////////////////////////////////////////////////////////
23 //v1.2 2004-06-17
24 ///////////////////////////////////////////////////////////////////
25 //implemented by Henrik Strohmayer
26 ///////////////////////////////////////////////////////////////////
27 
28 
29 
30 ///////////////////////////////////////////////////////////////////
31 //tdeg
32 ///////////////////////////////////////////////////////////////////
33 //Description: returns the normal degree of the input polynomial
34 ///////////////////////////////////////////////////////////////////
35 //Uses: pTotaldegree
36 ///////////////////////////////////////////////////////////////////
37 
38 int tdeg(poly p)
39 {
40  int res=0;
41  if(p!=NULL) res=pTotaldegree(p);
42  return(res);
43 }
44 
45 ///////////////////////////////////////////////////////////////////
46 
47 
48 ///////////////////////////////////////////////////////////////////
49 //getMaxTdeg
50 ///////////////////////////////////////////////////////////////////
51 //Description: returns the maximum normal degree of the
52 //polynomials of the input ideal
53 ///////////////////////////////////////////////////////////////////
54 //Uses: pTotaldegree
55 ///////////////////////////////////////////////////////////////////
56 
57 int getMaxTdeg(ideal I)
58 {
59  int res=-1;
60  int length=(int)I->ncols;
61  for(int j=length-1;j>=0;j--)
62  {
63  if ((I->m)[j]!=NULL)
64  {
65  int temp=pTotaldegree((I->m)[j]);
66  if(temp>res) {res=temp;}
67  }
68  }
69  return(res);
70 }
71 ///////////////////////////////////////////////////////////////////
72 
73 
74 ///////////////////////////////////////////////////////////////////
75 //getMaxPosOfNthRow
76 ///////////////////////////////////////////////////////////////////
77 //Description: returns the greatest integer of row n of the
78 //input intvec
79 ///////////////////////////////////////////////////////////////////
80 //Uses: none
81 ///////////////////////////////////////////////////////////////////
82 
84 {
85  int res=0;
86  assume( (0<n) && (n<=(v->rows())) );
87  {
88  int c=v->cols();
89  int cc=(n-1)*c;
90  res=abs((*v)[0+cc /*(n-1)*c*/]);
91  for (int i=c-1;i>=0;i--)
92  {
93  int temp=abs((*v)[i+cc /*(n-1)*c*/]);
94  if(temp>res){res=temp;}
95  }
96  }
97  return(res);
98 }
99 
100 ///////////////////////////////////////////////////////////////////
101 
102 
103 ///////////////////////////////////////////////////////////////////
104 //getInvEps64
105 ///////////////////////////////////////////////////////////////////
106 //Description: returns the inverse of epsilon (see relevant
107 //part of thesis for definition of epsilon)
108 ///////////////////////////////////////////////////////////////////
109 //Uses: getMaxPosOfNthRow
110 ///////////////////////////////////////////////////////////////////
111 
112 int64 getInvEps64(ideal G,intvec *targm,int pertdeg)
113 {
114  int n;
115  int64 temp64;
116  int64 sum64=0;
117 //think n=2 is enough (instead of n=1)
118  for (n=pertdeg; n>1; n--)
119  {
120  temp64=getMaxPosOfNthRow(targm,n);
121  sum64 += temp64;
122  }
123  int64 inveps64=getMaxTdeg(G)*sum64+1;
124 
125  //overflow test
126  if( sum64!=0 && (((inveps64-1)/sum64)!=getMaxTdeg(G)) )
127  overflow_error=11;
128 
129  return(inveps64);
130 }
131 
132 ///////////////////////////////////////////////////////////////////
133 
134 
135 ///////////////////////////////////////////////////////////////////
136 //invEpsOk64
137 ///////////////////////////////////////////////////////////////////
138 //Description: checks whether the input inveps64 is the same
139 //as the one you get from the current ideal I
140 ///////////////////////////////////////////////////////////////////
141 //Uses: getInvEps64
142 ///////////////////////////////////////////////////////////////////
143 
144 int invEpsOk64(ideal I, intvec *targm, int pertdeg, int64 inveps64)
145 {
146  int64 temp64=getInvEps64(I,targm,pertdeg);
147  if (inveps64>=temp64)
148  {
149  return(1);
150  }
151  else
152  {
153  return(0);
154  }
155 }
156 
157 ///////////////////////////////////////////////////////////////////
158 
159 
160 ///////////////////////////////////////////////////////////////////
161 //getNthRow
162 ///////////////////////////////////////////////////////////////////
163 //Description: returns an intvec containing the nth row of v
164 ///////////////////////////////////////////////////////////////////
165 //Uses: none
166 ///////////////////////////////////////////////////////////////////
167 
169 {
170  int r=v->rows();
171  int c=v->cols();
172  intvec *res=new intvec(c);
173  if((0<n) && (n<=r))
174  {
175  int cc=(n-1)*c;
176  for (int i=0; i<c; i++)
177  {
178  (*res)[i]=(*v)[i+cc /*(n-1)*c*/ ];
179  }
180  }
181  return(res);
182 }
183 
185 {
186  int r=v->rows();
187  int c=v->cols();
188  int64vec *res=new int64vec(c);
189  if((0<n) && (n<=r))
190  {
191  int cc=(n-1)*c;
192  for (int i=0; i<c; i++)
193  {
194  (*res)[i]=(int64)(*v)[i+cc /*(n-1)*c*/ ];
195  }
196  }
197  return(res);
198 }
199 
200 ///////////////////////////////////////////////////////////////////
201 
202 
203 ///////////////////////////////////////////////////////////////////
204 //getTaun64
205 ///////////////////////////////////////////////////////////////////
206 //Description: returns a list containing the nth perturbed vector
207 //and the inveps64 used to calculate the vector
208 ///////////////////////////////////////////////////////////////////
209 //Uses: getNthRow,getInvEps64
210 ///////////////////////////////////////////////////////////////////
211 
212 void getTaun64(ideal G,intvec* targm,int pertdeg, int64vec** v64, int64 & i64)
213 {
214  int64vec* taun64=getNthRow64(targm,1);
215  int64vec *temp64,*add64;
216  int64 inveps64=1;
217  if (pertdeg>1) inveps64=getInvEps64(G,targm,pertdeg);
218 
219  int n;
220  //temp64 is used to check for overflow:
221  for (n=2; n<=pertdeg; n++)
222  {
223  if (inveps64!=1)
224  {
225  temp64=iv64Copy(taun64);
226  (*taun64)*=inveps64;
227  for(int i=0; i<rVar(currRing);i++)
228  {
229  if((*temp64)[i]!=0 && (((*taun64)[i])/((*temp64)[i]))!=inveps64)
230  overflow_error=12;
231  }
232  delete temp64;
233  }
234  temp64=iv64Copy(taun64);
235  add64=getNthRow64(targm,n);
236  taun64=iv64Add(add64,taun64);
237  for(int i=0; i<rVar(currRing);i++)
238  {
239  if( ( ((*temp64)[i]) > 0 ) && ( ((*add64)[i]) > 0 ) )
240  {
241  if( ((*taun64)[i]) < ((*temp64)[i]) )
242  overflow_error=13;
243  }
244  if( ( ((*temp64)[i]) < 0 ) && ( ((*add64)[i]) < 0 ) )
245  {
246  if( ((*taun64)[i]) > ((*temp64)[i]) )
247  overflow_error=13;
248  }
249  }
250  delete temp64;
251  }
252 
253  //lists taunlist64=makeTaunList64(taun64,inveps64);
254  //return(taunlist64);
255  *v64=taun64;
256  i64=inveps64;
257 }
258 
259 ///////////////////////////////////////////////////////////////////
260 //scalarProduct64
261 ///////////////////////////////////////////////////////////////////
262 //Description: returns the scalar product of int64vecs a and b
263 ///////////////////////////////////////////////////////////////////
264 //Assume: Overflow tests assumes input has nonnegative entries
265 ///////////////////////////////////////////////////////////////////
266 //Uses: none
267 ///////////////////////////////////////////////////////////////////
268 
270 {
271  assume( a->length() == b->length());
272  int i, n = a->length();
273  int64 result = 0;
274  int64 temp1,temp2;
275  for(i=n-1; i>=0; i--)
276  {
277  temp1=(*a)[i] * (*b)[i];
278  if((*a)[i]!=0 && (temp1/(*a)[i])!=(*b)[i]) overflow_error=1;
279 
280  temp2=result;
281  result += temp1;
282 
283  //since both vectors always have nonnegative entries in init64
284  //result should be >=temp2
285  if(temp2>result) overflow_error=2;
286  }
287 
288  return result;
289 }
290 ///////////////////////////////////////////////////////////////////
291 
292 
293 ///////////////////////////////////////////////////////////////////
294 //init64
295 ///////////////////////////////////////////////////////////////////
296 //Description: returns the initial form ideal of the input ideal
297 //I w.r.t. the weight vector currw64
298 ///////////////////////////////////////////////////////////////////
299 //Uses: idealSize,getNthPolyOfId,leadExp64,pLength
300 ///////////////////////////////////////////////////////////////////
301 
302 ideal init64(ideal G,int64vec *currw64)
303 {
304  int length=IDELEMS(G);
305  ideal I=idInit(length,G->rank);
306  int j;
307  int64 leadingweight,templeadingweight;
308  poly p=NULL;
309  poly leadp=NULL;
310  for (j=1; j<=length; j++)
311  {
312  p=getNthPolyOfId(G,j);
313  int64vec *tt=leadExp64(p);
314  leadingweight=scalarProduct64(currw64,tt /*leadExp64(p)*/);
315  delete tt;
316  while (p!=NULL)
317  {
318  tt=leadExp64(p);
319  templeadingweight=scalarProduct64(currw64,tt /*leadExp64(p)*/);
320  delete tt;
321  if(templeadingweight==leadingweight)
322  {
323  leadp=pAdd(leadp,pHead(p));
324  }
325  if(templeadingweight>leadingweight)
326  {
327  pDelete(&leadp);
328  leadp=pHead(p);
329  leadingweight=templeadingweight;
330  }
331  pIter(p);
332  }
333  (I->m)[j-1]=leadp;
334  p=NULL;
335  leadp=NULL;
336  }
337  return(I);
338 }
339 
340 ///////////////////////////////////////////////////////////////////
341 
342 
343 ///////////////////////////////////////////////////////////////////
344 //currwOnBorder64
345 ///////////////////////////////////////////////////////////////////
346 //Description: returns TRUE if currw64 lies on the border of the
347 //groebner cone of the order w.r.t. which the reduced GB G
348 //is calculated, otherwise FALSE
349 ///////////////////////////////////////////////////////////////////
350 //Uses: init64
351 ///////////////////////////////////////////////////////////////////
352 
354 {
355  ideal J=init64(G,currw64);
356  int length=idealSize(J);
357  BOOLEAN res=FALSE;
358  for(int i=length; i>0; i--)
359  {
360  //if(pLength(getNthPolyOfId(J,i))>1)
361  poly p=getNthPolyOfId(J,i);
362  if ((p!=NULL) && (pNext(p)!=NULL))
363  {
364  res=TRUE;break;
365  }
366  }
367  idDelete(&J);
368  return res;
369 }
370 
371 ///////////////////////////////////////////////////////////////////
372 
373 
374 ///////////////////////////////////////////////////////////////////
375 //noPolysWithMoreThanTwoTerms
376 ///////////////////////////////////////////////////////////////////
377 //Description: returns TRUE if all polynomials of Gw are of length
378 //<=2, otherwise FALSE
379 ///////////////////////////////////////////////////////////////////
380 //Uses: idealSize, getNthPolyOfId
381 ///////////////////////////////////////////////////////////////////
382 
384 {
385  int length=idealSize(Gw);
386  for(int i=length; i>0; i--)
387  {
388  //if(pLength(getNthPolyOfId(Gw,i))>2)
389  poly p=getNthPolyOfId(Gw,i);
390  if ((p!=NULL) && (pNext(p)!=NULL) && (pNext(pNext(p))!=NULL))
391  {
392  return FALSE;
393  }
394  }
395  return TRUE;
396 }
397 
398 ///////////////////////////////////////////////////////////////////
399 
400 
401 ///////////////////////////////////////////////////////////////////
402 //DIFFspy
403 ///////////////////////////////////////////////////////////////////
404 //Description: returns the length of the list to be created in
405 //DIFF
406 ///////////////////////////////////////////////////////////////////
407 //Uses: idealSize,pLength,getNthPolyOfId
408 ///////////////////////////////////////////////////////////////////
409 
410 int DIFFspy(ideal G)
411 {
412  int s=idealSize(G);
413  int j;
414  int temp;
415  int sum=0;
416  poly p;
417  for (j=1; j<=s; j++)
418  {
419  p=getNthPolyOfId(G,j);
420  if((temp=pLength(p))>0) {sum += (temp-1);}
421  }
422  return(sum);
423 }
424 
425 ///////////////////////////////////////////////////////////////////
426 
427 
428 ///////////////////////////////////////////////////////////////////
429 //DIFF
430 ///////////////////////////////////////////////////////////////////
431 //Description: returns a list of all differences of leading
432 //exponents and nonleading exponents of elements of the current
433 //GB (see relevant part of thesis for further details)
434 ///////////////////////////////////////////////////////////////////
435 //Uses: DIFFspy,idealSize,getNthPolyOfId,leadExp,pLength
436 ///////////////////////////////////////////////////////////////////
437 
438 intvec* DIFF(ideal G)
439 {
440  intvec *v,*w;
441  poly p;
442  int s=idealSize(G);
443  int n=rVar(currRing);
444  int m=DIFFspy(G);
445  intvec* diffm=new intvec(m,n,0);
446  int j,l;
447  int inc=0;
448  for (j=1; j<=s; j++)
449  {
450  p=getNthPolyOfId(G,j);
451  v=leadExp(p);
452  pIter(p);
453  while(p!=NULL)
454  {
455  inc++;
456  intvec *lep=leadExp(p);
457  w=ivSub(v,lep /*leadExp(p)*/);
458  delete lep;
459  pIter(p);
460  for (l=1; l<=n; l++)
461  {
462  // setPosOfIM(diffm,inc,l,(*w)[l-1]);
463  IMATELEM(*diffm,inc,l) =(*w)[l-1] ;
464  }
465  delete w;
466  }
467  delete v;
468  }
469  return(diffm);
470 }
471 
472 ///////////////////////////////////////////////////////////////////
473 
474 
475 ///////////////////////////////////////////////////////////////////
476 //gett64
477 ///////////////////////////////////////////////////////////////////
478 //Description: returns the t corresponding to the vector listw
479 //which contains a vector from the list returned by DIFF
480 ///////////////////////////////////////////////////////////////////
481 //Uses: ivSize
482 ///////////////////////////////////////////////////////////////////
483 
484 void gett64(intvec* listw, int64vec* currw64, int64vec* targw64, int64 &tvec0, int64 &tvec1)
485 {
486  int s=ivSize(listw);
487  int j;
488  int64 zaehler64=0;
489  int64 nenner64=0;
490  int64 temp1,temp2,temp3,temp4; //overflowstuff
491  for(j=1; j<=s; j++)
492  {
493 
494  temp3=zaehler64;
495  temp1=((int64)((*listw)[j-1]));
496  temp2=((*currw64)[j-1]);
497  temp4=temp1*temp2;
498  zaehler64=temp3-temp4;
499 
500  //overflow test
501  if(temp1!=0 && (temp4/temp1)!=temp2) overflow_error=3;
502 
503  if( ( temp3<0 && temp4>0 ) || ( temp3>0 && temp4<0 ) )
504  {
505  int64 abs_t3=abs64(temp3);
506  if( (abs_t3+abs64(temp4))<abs_t3 ) overflow_error=4;
507  }
508 
509  //overflow test
510  temp1=((*targw64)[j-1])-((*currw64)[j-1]);
511  //this subtraction can never yield an overflow since both number
512  //will always be positive
513  temp2=((int64)((*listw)[j-1]));
514  temp3=nenner64;
515  temp4=temp1*temp2;
516  nenner64=temp3+temp4;
517 
518  //overflow test
519  if(temp1!=0 && ((temp1*temp2)/temp1)!=temp2) overflow_error=5;
520 
521  if( (temp3>0 && temp4>0) ||
522  (temp3<0 && temp4<0) )
523  {
524  int64 abs_t3=abs64(temp3);
525  if( (abs_t3+abs64(temp4))<abs_t3 )
526  {
527  overflow_error=6;
528  }
529  }
530  }
531 
532  if (nenner64==0)
533  {
534  zaehler64=2;
535  }
536  else
537  {
538  if ( (zaehler64<=0) && (nenner64<0) )
539  {
540  zaehler64=-zaehler64;
541  nenner64=-nenner64;
542  }
543  }
544 
545  int64 g=gcd64(zaehler64,nenner64);
546 
547  tvec0=zaehler64/g;
548  tvec1=nenner64/g;
549 
550 }
551 
552 ///////////////////////////////////////////////////////////////////
553 
554 
555 ///////////////////////////////////////////////////////////////////
556 //nextt64
557 ///////////////////////////////////////////////////////////////////
558 //Description: returns the t determining the next weight vector
559 ///////////////////////////////////////////////////////////////////
560 //Uses:
561 ///////////////////////////////////////////////////////////////////
562 
563 void nextt64(ideal G,int64vec* currw64,int64vec* targw64, int64 & tvec0, int64 & tvec1)
564 {
565  intvec* diffm=DIFF(G);
566  int s=diffm->rows();
567  tvec0=(int64)2;
568  tvec1=(int64)0;
569  intvec *tt;
570  for(int j=1; j<=s; j++)
571  {
572  tt=getNthRow(diffm,j);
573  int64 temptvec0, temptvec1;
574  gett64(tt,currw64,targw64,temptvec0, temptvec1);
575  delete tt;
576 
577  //if tempt>0 both parts will be>0
578  if ( (temptvec1!=0) //that tempt is defined
579  &&
580  (temptvec0>0) && (temptvec1>0) //that tempt>0
581  )
582  {
583  if( ( (temptvec0) <= (temptvec1) ) //that tempt<=1
584  &&
585  ( ( (temptvec0) * (tvec1) ) <
586  ( (temptvec1) * (tvec0) ) )
587  )
588  { //that tempt<t
589  tvec0=temptvec0;
590  tvec1=temptvec1;
591  }
592  }
593  }
594  delete diffm;
595  return;
596 }
597 
598 ///////////////////////////////////////////////////////////////////
599 
600 
601 ///////////////////////////////////////////////////////////////////
602 //nextw64
603 ///////////////////////////////////////////////////////////////////
604 //Uses:iv64Size,gcd,
605 ///////////////////////////////////////////////////////////////////
606 
608  int64 nexttvec0, int64 nexttvec1)
609 {
610  //to do (targw-currw)*tvec[0]+currw*tvec[1]
611  int64vec* tempv;
612  int64vec* nextweight;
613  int64vec* a=iv64Sub(targw,currw);
614  //no overflow can occur since both are>=0
615 
616  //to test overflow
617  tempv=iv64Copy(a);
618  *a *= (nexttvec0);
619  for(int i=0; i<rVar(currRing); i++)
620  {
621  if( (nexttvec0) !=0 &&
622  (((*a)[i])/(nexttvec0))!=((*tempv)[i]) )
623  {
624  overflow_error=7;
625  break;
626  }
627  }
628  delete tempv;
629  int64vec* b=currw;
630  tempv=iv64Copy(b);
631  *b *= (nexttvec1);
632  for(int i=0; i<rVar(currRing); i++)
633  {
634  if( (nexttvec1) !=0 &&
635  (((*b)[i])/(nexttvec1))!=((*tempv)[i]) )
636  {
637  overflow_error=8;
638  break;
639  }
640  }
641  delete tempv;
642  nextweight=iv64Add(a,b);
643 
644  for(int i=0; i<rVar(currRing); i++)
645  {
646  if( (((*a)[i])>=0 && ((*b)[i])>=0) ||
647  (((*a)[i])<0 && ((*b)[i])<0) )
648  {
649  if( (abs64((*a)[i]))>abs64((*nextweight)[i]) ||
650  (abs64((*b)[i]))>abs64((*nextweight)[i])
651  )
652  {
653  overflow_error=9;
654  break;
655  }
656  }
657  }
658 
659  //to reduce common factors of nextweight
660  int s=iv64Size(nextweight);
661  int64 g,temp;
662  g=(*nextweight)[0];
663  for (int i=1; i<s; i++)
664  {
665  temp=(*nextweight)[i];
666  g=gcd64(g,temp);
667  if (g==1) break;
668  }
669 
670  if (g!=1) *nextweight /= g;
671 
672  return(nextweight);
673 }
674 
675 ///////////////////////////////////////////////////////////////////
676 
677 
678 
679 ///////////////////////////////////////////////////////////////////
680 //FUNCTIONS NOT ORIGINATING FROM THE SINGULAR IMPLEMENTATION CODE
681 ///////////////////////////////////////////////////////////////////
682 
683 
684 ///////////////////////////////////////////////////////////////////
685 //getNthPolyOfId
686 ///////////////////////////////////////////////////////////////////
687 //Description: returns the nth poly of ideal I
688 ///////////////////////////////////////////////////////////////////
689 poly getNthPolyOfId(ideal I,int n)
690 {
691  if(0<n && n<=((int)I->ncols))
692  {
693  return (I->m)[n-1];
694  }
695  else
696  {
697  return(NULL);
698  }
699 }
700 
701 ///////////////////////////////////////////////////////////////////
702 
703 
704 ///////////////////////////////////////////////////////////////////
705 //idealSize
706 ///////////////////////////////////////////////////////////////////
707 //Description: returns the number of generator of input ideal I
708 ///////////////////////////////////////////////////////////////////
709 //Uses: none
710 ///////////////////////////////////////////////////////////////////
711 // #define idealSize(I) IDELEMS(I)
712 ///////////////////////////////////////////////////////////////////
713 
714 
715 ///////////////////////////////////////////////////////////////////
716 //ivSize
717 ///////////////////////////////////////////////////////////////////
718 //Description: returns the number of entries of v
719 ///////////////////////////////////////////////////////////////////
720 //Uses: none
721 ///////////////////////////////////////////////////////////////////
722 
723 // inline int ivSize(intvec* v){ return((v->rows())*(v->cols())); }
724 
725 ///////////////////////////////////////////////////////////////////
726 
727 
728 ///////////////////////////////////////////////////////////////////
729 //iv64Size
730 ///////////////////////////////////////////////////////////////////
731 //Description: returns the number of entries of v
732 ///////////////////////////////////////////////////////////////////
733 //Uses: none
734 ///////////////////////////////////////////////////////////////////
735 
736 // int iv64Size(int64vec* v){ return((v->rows())*(v->cols())); }
737 
738 ///////////////////////////////////////////////////////////////////
739 
740 
741 ///////////////////////////////////////////////////////////////////
742 //leadExp
743 ///////////////////////////////////////////////////////////////////
744 //Description: returns an intvec containg the exponet vector of p
745 ///////////////////////////////////////////////////////////////////
746 //Uses: sizeof,omAlloc,omFree
747 ///////////////////////////////////////////////////////////////////
748 
750 {
751  int N=rVar(currRing);
752  int *e=(int*)omAlloc((N+1)*sizeof(int));
753  pGetExpV(p,e);
754  intvec* iv=new intvec(N);
755  for(int i=N;i>0;i--) { (*iv)[i-1]=e[i];}
756  omFree(e);
757  return(iv);
758 }
759 
760 ///////////////////////////////////////////////////////////////////
761 
762 
763 ///////////////////////////////////////////////////////////////////
764 //leadExp64
765 ///////////////////////////////////////////////////////////////////
766 //Description: returns an int64vec containing the exponent
767 //vector of p
768 ///////////////////////////////////////////////////////////////////
769 //Uses: sizeof,omAlloc,omFree
770 ///////////////////////////////////////////////////////////////////
771 
773 {
774  int N=rVar(currRing);
775  int *e=(int*)omAlloc((N+1)*sizeof(int));
776  pGetExpV(p,e);
777  int64vec* iv64=new int64vec(N);
778  for(int i=N;i>0;i--) { (*iv64)[i-1]=(int64)e[i];}
779  omFree(e);
780  return(iv64);
781 }
782 
783 ///////////////////////////////////////////////////////////////////
784 
785 
786 ///////////////////////////////////////////////////////////////////
787 //setPosOfIM
788 ///////////////////////////////////////////////////////////////////
789 //Description: sets entry i,j of im to val
790 ///////////////////////////////////////////////////////////////////
791 //Uses: none
792 ///////////////////////////////////////////////////////////////////
793 
794 //void setPosOfIM(intvec* im,int i,int j,int val){
795 // int r=im->rows();
796 // int c=im->cols();
797 // if( (0<i && i<=r) && (0<j && j<=c) ){
798 // (*im)[(i-1)*c+j-1]=val;
799 // }
800 // return;
801 //}
802 
803 ///////////////////////////////////////////////////////////////////
804 
805 
806 ///////////////////////////////////////////////////////////////////
807 //scalarProduct
808 ///////////////////////////////////////////////////////////////////
809 //Description: returns the scalar product of intvecs a and b
810 ///////////////////////////////////////////////////////////////////
811 //Uses: none
812 ///////////////////////////////////////////////////////////////////
813 
814 static inline long scalarProduct(intvec* a, intvec* b)
815 {
816  assume( a->length() == b->length());
817  int i, n = a->length();
818  long result = 0;
819  for(i=n-1; i>=0; i--)
820  result += (*a)[i] * (*b)[i];
821  return result;
822 }
823 
824 ///////////////////////////////////////////////////////////////////
825 
826 
827 
828 ///////////////////////////////////////////////////////////////////
829 
830 
831 ///////////////////////////////////////////////////////////////////
832 //gcd
833 ///////////////////////////////////////////////////////////////////
834 //Description: returns the gcd of a and b
835 ///////////////////////////////////////////////////////////////////
836 //Uses: none
837 ///////////////////////////////////////////////////////////////////
838 
839 int gcd(int a, int b)
840 {
841  int r, p0 = a, p1 = b;
842  if(p0 < 0)
843  p0 = -p0;
844 
845  if(p1 < 0)
846  p1 = -p1;
847  while(p1 != 0)
848  {
849  r = p0 % p1;
850  p0 = p1;
851  p1 = r;
852  }
853  return p0;
854 }
855 
856 ///////////////////////////////////////////////////////////////////
857 
858 
859 ///////////////////////////////////////////////////////////////////
860 //gcd64
861 ///////////////////////////////////////////////////////////////////
862 //Description: returns the gcd of a and b
863 ///////////////////////////////////////////////////////////////////
864 //Uses: none
865 ///////////////////////////////////////////////////////////////////
866 
868 {
869  int64 r, p0 = a, p1 = b;
870  if(p0 < 0)
871  p0 = -p0;
872 
873  if(p1 < 0)
874  p1 = -p1;
875 
876  while(p1 != ((int64)0) )
877  {
878  r = p0 % p1;
879  p0 = p1;
880  p1 = r;
881  }
882 
883  return p0;
884 }
885 
886 ///////////////////////////////////////////////////////////////////
887 
888 
889 ///////////////////////////////////////////////////////////////////
890 //abs64
891 ///////////////////////////////////////////////////////////////////
892 //Description: returns the absolute value of int64 i
893 ///////////////////////////////////////////////////////////////////
894 //Uses: none
895 ///////////////////////////////////////////////////////////////////
896 
897 //int64 abs64(int64 i)
898 //{
899 // if(i>=0)
900 // return(i);
901 //else
902 // return((-i));
903 //}
904 
905 ///////////////////////////////////////////////////////////////////
906 
907 
908 ///////////////////////////////////////////////////////////////////
909 //makeTaunList64
910 ///////////////////////////////////////////////////////////////////
911 //Description: makes a list of an int64vec and an int64
912 ///////////////////////////////////////////////////////////////////
913 //Uses: omAllocBin
914 ///////////////////////////////////////////////////////////////////
915 
916 #if 0
917 lists makeTaunList64(int64vec *iv64,int64 i64)
918 {
920  l->Init(2);
921  l->m[0].rtyp=INTVEC_CMD;
922  l->m[1].rtyp=INT_CMD;
923  l->m[0].data=(void *)iv64;
924  l->m[1].data=(void *)i64;
925  return l;
926 }
927 #endif
928 
929 ///////////////////////////////////////////////////////////////////
930 
931 
932 ///////////////////////////////////////////////////////////////////
933 //idStd
934 ///////////////////////////////////////////////////////////////////
935 //Description: returns the GB of G calculated w.r.t. the order of
936 //currRing
937 ///////////////////////////////////////////////////////////////////
938 //Uses: kStd,idSkipZeroes
939 ///////////////////////////////////////////////////////////////////
940 
941 ideal idStd(ideal G)
942 {
943  ideal GG = kStd(G, NULL, testHomog, NULL);
944  idSkipZeroes(GG);
945  return GG;
946 }
947 
948 ///////////////////////////////////////////////////////////////////
949 
950 
951 ///////////////////////////////////////////////////////////////////
952 //idInterRed
953 ///////////////////////////////////////////////////////////////////
954 //Description: returns the interreduction of G
955 ///////////////////////////////////////////////////////////////////
956 //Assumes: that the input is a GB
957 ///////////////////////////////////////////////////////////////////
958 //Uses: kInterRed,idSkipZeroes
959 ///////////////////////////////////////////////////////////////////
960 
961 ideal idInterRed(ideal G)
962 {
963  assume(G != NULL);
964 
965  ideal GG = kInterRedOld(G, NULL);
966  idDelete(&G);
967  return GG;
968 }
969 
970 ///////////////////////////////////////////////////////////////////
971 
972 
973 ///////////////////////////////////////////////////////////////////
974 //matIdLift
975 ///////////////////////////////////////////////////////////////////
976 //Description: yields the same reslut as lift in Singular
977 ///////////////////////////////////////////////////////////////////
978 //Uses: idLift,idModule2formatedMatrix
979 ///////////////////////////////////////////////////////////////////
980 
981 matrix matIdLift(ideal Gomega, ideal M)
982 {
983  ideal Mtmp = idLift(Gomega, M, NULL, FALSE, FALSE, TRUE, NULL);
984  int rows=IDELEMS(Gomega);
985  int cols=IDELEMS(Mtmp);
987  return res;
988 }
989 ///////////////////////////////////////////////////////////////////
990 
991 
992 ///////////////////////////////////////////////////////////////////
993 //rCopyAndChangeA
994 ///////////////////////////////////////////////////////////////////
995 //Description: changes currRing to a copy of currRing with the
996 //int64vec w instead of the old one
997 ///////////////////////////////////////////////////////////////////
998 //Assumes: that currRing is alterd as mentioned in rCopy0AndAddA
999 ///////////////////////////////////////////////////////////////////
1000 //Uses: rCopy0,rComplete,rChangeCurrRing,rSetWeightVec
1001 ///////////////////////////////////////////////////////////////////
1002 
1004 {
1005  ring rnew=rCopy0(currRing);
1006  rComplete(rnew);
1007  rSetWeightVec(rnew,w->iv64GetVec());
1008  rChangeCurrRing(rnew);
1009 }
1010 
1011 ///////////////////////////////////////////////////////////////////
1012 //rGetGlobalOrderMatrix
1013 ///////////////////////////////////////////////////////////////////
1014 //Description: returns a matrix corresponding to the order of r
1015 ///////////////////////////////////////////////////////////////////
1016 //Assumes: the order of r is a combination of the orders M,lp,dp,
1017 //Dp,wp,Wp,C
1018 ///////////////////////////////////////////////////////////////////
1019 //Uses: none
1020 ///////////////////////////////////////////////////////////////////
1021 
1023 {
1024  int n=rVar(r);
1025  int64vec* res=new int64vec(n,n,(int64)0);
1026  if (rHasLocalOrMixedOrdering(r)) return res;
1027  int pos1=0;
1028  int pos2=0;
1029  int temp;
1030  int i=0;
1031  while(r->order[i]!=0 && pos2<n)
1032  {
1033  pos2=pos2+r->block1[i] - r->block0[i];
1034 
1035  if(r->order[i]==ringorder_lp)
1036  {
1037  temp=pos1;
1038  for(int j=pos1; j<=pos2; j++)
1039  (*res)[j*n+j]=(int64)1;
1040  }
1041  else if(r->order[i]==ringorder_dp)
1042  {
1043  for(int j=pos1;j<=pos2;j++)
1044  (*res)[pos1*n+j]=(int64)1;
1045  for(int j=1;j<=(pos2-pos1);j++)
1046  (*res)[(pos1+j)*n+(pos2+1-j)]=(int64)-1;
1047  }
1048  else if(r->order[i]==ringorder_Dp)
1049  {
1050  for(int j=pos1;j<=pos2;j++)
1051  (*res)[pos1*n+j]=(int64)1;
1052  for(int j=1;j<=(pos2-pos1);j++)
1053  (*res)[(pos1+j)*n+(pos1+j-1)]=(int64)1;
1054  }
1055  else if(r->order[i]==ringorder_wp)
1056  {
1057  int* weights=r->wvhdl[i];
1058  for(int j=pos1;j<=pos2;j++)
1059  (*res)[pos1*n+j]=(int64)weights[j-pos1];
1060  for(int j=1;j<=(pos2-pos1);j++)
1061  (*res)[(pos1+j)*n+(pos2+1-j)]=(int64)-1;
1062  }
1063  else if(r->order[i]==ringorder_Wp)
1064  {
1065  int* weights=r->wvhdl[i];
1066  for(int j=pos1;j<=pos2;j++)
1067  (*res)[pos1*n+j]=(int64)weights[j-pos1];
1068  for(int j=1;j<=(pos2-pos1);j++)
1069  (*res)[(pos1+j)*n+(pos1+j-1)]=(int64)1;
1070  }
1071 
1072  else if(r->order[0]==ringorder_M)
1073  {
1074  int* weights=r->wvhdl[0];
1075  for(int j=pos1;j<((pos2+1)*(pos2+1));j++)
1076  (*res)[j]=(int64)weights[j];
1077  }
1078 
1079  pos1=pos2+1;
1080  pos2=pos2+1;
1081  i++;
1082  }
1083 
1084  return(res);
1085 }
1086 
1087 ///////////////////////////////////////////////////////////////////
1088 
1089 
1090 ///////////////////////////////////////////////////////////////////
1091 //rGetGlobalOrderWeightVec
1092 ///////////////////////////////////////////////////////////////////
1093 //Description: returns a weight vector corresponding to the first
1094 //row of a matrix corresponding to the order of r
1095 ///////////////////////////////////////////////////////////////////
1096 //Uses: none
1097 ///////////////////////////////////////////////////////////////////
1098 
1100 {
1101  int n=rVar(r);
1102  int64vec* res=new int64vec(n);
1103 
1104  if (rHasLocalOrMixedOrdering(r)) return res;
1105 
1106  int length;
1107 
1108  if(r->order[0]==ringorder_lp)
1109  {
1110  (*res)[0]=(int64)1;
1111  }
1112  else if( (r->order[0]==ringorder_dp) || (r->order[0]==ringorder_Dp) )
1113  {
1114  length=r->block1[0] - r->block0[0];
1115  for (int j=0;j<=length;j++)
1116  (*res)[j]=(int64)1;
1117  }
1118  else if( (r->order[0]==ringorder_wp) || (r->order[0]==ringorder_Wp) ||
1119  (r->order[0]==ringorder_a) || (r->order[0]==ringorder_M) )
1120  {
1121  int* weights=r->wvhdl[0];
1122  length=r->block1[0] - r->block0[0];
1123  for (int j=0;j<=length;j++)
1124  (*res)[j]=(int64)weights[j];
1125  }
1126  else if( /*(*/ r->order[0]==ringorder_a64 /*)*/ )
1127  {
1128  int64* weights=(int64*)r->wvhdl[0];
1129  length=r->block1[0] - r->block0[0];
1130  for (int j=0;j<=length;j++)
1131  (*res)[j]=weights[j];
1132  }
1133 
1134  return(res);
1135 }
1136 
1137 ///////////////////////////////////////////////////////////////////
1138 
1139 
1140 ///////////////////////////////////////////////////////////////////
1141 //sortRedSB
1142 ///////////////////////////////////////////////////////////////////
1143 //Description: sorts a reduced GB of an ideal after the leading
1144 //terms of the polynomials with the smallest one first
1145 ///////////////////////////////////////////////////////////////////
1146 //Assumes: that the given input is a minimal GB
1147 ///////////////////////////////////////////////////////////////////
1148 //Uses:idealSize,idCopy,pLmCmp
1149 ///////////////////////////////////////////////////////////////////
1150 
1151 ideal sortRedSB(ideal G)
1152 {
1153  int s=idealSize(G);
1154  poly* m=G->m;
1155  poly p,q;
1156  for (int i=0; i<(s-1); i++)
1157  {
1158  for (int j=0; j<((s-1)-i); j++)
1159  {
1160  p=m[j];
1161  q=m[j+1];
1162  if (pLmCmp(p,q)==1)
1163  {
1164  m[j+1]=p;
1165  m[j]=q;
1166  }
1167  }
1168  }
1169  return(G);
1170 }
1171 
1172 ///////////////////////////////////////////////////////////////////
1173 
1174 
1175 ///////////////////////////////////////////////////////////////////
1176 //int64VecToIntVec
1177 ///////////////////////////////////////////////////////////////////
1178 //Description: converts an int64vec to an intvec
1179 // deletes the input
1180 ///////////////////////////////////////////////////////////////////
1181 //Assumes: that the int64vec contains no entries larger than 2^32
1182 ///////////////////////////////////////////////////////////////////
1183 //Uses: none
1184 ///////////////////////////////////////////////////////////////////
1185 
1187 {
1188  int r=source->rows();
1189  int c=source->cols();
1190  intvec* res=new intvec(r,c,0);
1191  for(int i=0;i<r;i++){
1192  for(int j=0;j<c;j++){
1193  (*res)[i*c+j]=(int)(*source)[i*c+j];
1194  }
1195  }
1196  delete source;
1197  return(res);
1198 }
1199 
1200 ///////////////////////////////////////////////////////////////////
BOOLEAN rHasLocalOrMixedOrdering(const ring r)
Definition: ring.h:747
void nextt64(ideal G, int64vec *currw64, int64vec *targw64, int64 &tvec0, int64 &tvec1)
Definition: walkSupport.cc:563
int64vec * rGetGlobalOrderWeightVec(ring r)
ideal sortRedSB(ideal G)
#define omAllocBin(bin)
Definition: omAllocDecl.h:205
void rCopyAndChangeA(int64vec *w)
const CanonicalForm int s
Definition: facAbsFact.cc:55
sleftv * m
Definition: lists.h:45
for int64 weights
Definition: ring.h:79
const poly a
Definition: syzextra.cc:212
Definition: tok.h:95
#define pAdd(p, q)
Definition: polys.h:186
static long scalarProduct(intvec *a, intvec *b)
Definition: walkSupport.cc:814
#define idDelete(H)
delete an ideal
Definition: ideals.h:29
Definition: lists.h:22
#define FALSE
Definition: auxiliary.h:95
Compatiblity layer for legacy polynomial operations (over currRing)
return P p
Definition: myNF.cc:203
int64vec * iv64Copy(int64vec *o)
Definition: int64vec.h:75
#define pLmCmp(p, q)
returns 0|1|-1 if p=q|p>q|p<q w.r.t monomial ordering
Definition: polys.h:105
int rows() const
Definition: intvec.h:88
int cols() const
Definition: int64vec.h:56
int64vec * iv64Add(int64vec *a, int64vec *b)
Definition: int64vec.cc:173
static short rVar(const ring r)
#define rVar(r) (r->N)
Definition: ring.h:580
long int64
Definition: auxiliary.h:67
int64vec * nextw64(int64vec *currw, int64vec *targw, int64 nexttvec0, int64 nexttvec1)
Definition: walkSupport.cc:607
#define TRUE
Definition: auxiliary.h:99
ideal kStd(ideal F, ideal Q, tHomog h, intvec **w, intvec *hilb, int syzComp, int newIdeal, intvec *vw, s_poly_proc_t sp)
Definition: kstd1.cc:2231
intvec * leadExp(poly p)
Definition: walkSupport.cc:749
intvec * ivSub(intvec *a, intvec *b)
Definition: intvec.cc:280
int ivSize(intvec *v)
Definition: walkSupport.h:36
int64vec * getNthRow64(intvec *v, int n)
Definition: walkSupport.cc:184
g
Definition: cfModGcd.cc:4031
int64vec * leadExp64(poly p)
Definition: walkSupport.cc:772
static TreeM * G
Definition: janet.cc:38
#define omAlloc(size)
Definition: omAllocDecl.h:210
int tdeg(poly p)
Definition: walkSupport.cc:38
Rational abs(const Rational &a)
Definition: GMPrat.cc:443
const CanonicalForm & GG
Definition: cfModGcd.cc:4017
void * data
Definition: subexpr.h:89
#define pIter(p)
Definition: monomials.h:44
poly res
Definition: myNF.cc:322
int invEpsOk64(ideal I, intvec *targm, int pertdeg, int64 inveps64)
Definition: walkSupport.cc:144
#define M
Definition: sirandom.c:24
ring currRing
Widely used global variable which specifies the current polynomial ring for Singular interpreter and ...
Definition: polys.cc:10
int64 abs64(int64 i)
Definition: walkSupport.h:44
const ring r
Definition: syzextra.cc:208
BOOLEAN overflow_error
Definition: walkMain.cc:37
ideal idStd(ideal G)
Definition: walkSupport.cc:941
int length() const
Definition: int64vec.h:55
Definition: intvec.h:14
const CanonicalForm CFMap CFMap & N
Definition: cfEzgcd.cc:49
BOOLEAN rComplete(ring r, int force)
this needs to be called whenever a new ring is created: new fields in ring are created (like VarOffse...
Definition: ring.cc:3351
for(int i=0;i< R->ExpL_Size;i++) Print("%09lx "
Definition: cfEzgcd.cc:66
int j
Definition: myNF.cc:70
#define omFree(addr)
Definition: omAllocDecl.h:261
static long pTotaldegree(poly p)
Definition: polys.h:265
#define assume(x)
Definition: mod2.h:403
ideal init64(ideal G, int64vec *currw64)
Definition: walkSupport.cc:302
int rows() const
Definition: int64vec.h:57
ring rCopy0(const ring r, BOOLEAN copy_qideal, BOOLEAN copy_ordering)
Definition: ring.cc:1321
int64vec * iv64Sub(int64vec *a, int64vec *b)
Definition: int64vec.cc:203
ideal idInterRed(ideal G)
Definition: walkSupport.cc:961
void getTaun64(ideal G, intvec *targm, int pertdeg, int64vec **v64, int64 &i64)
Definition: walkSupport.cc:212
int m
Definition: cfEzgcd.cc:119
int getMaxTdeg(ideal I)
Definition: walkSupport.cc:57
int i
Definition: cfEzgcd.cc:123
ideal kInterRedOld(ideal F, ideal Q)
Definition: kstd1.cc:3177
int iv64Size(int64vec *v)
Definition: walkSupport.h:37
matrix matIdLift(ideal Gomega, ideal M)
Definition: walkSupport.cc:981
static unsigned pLength(poly a)
Definition: p_polys.h:189
intvec * DIFF(ideal G)
Definition: walkSupport.cc:438
#define pHead(p)
returns newly allocated copy of Lm(p), coef is copied, next=NULL, p might be NULL ...
Definition: polys.h:67
#define IDELEMS(i)
Definition: simpleideals.h:24
int64 gcd64(int64 a, int64 b)
Definition: walkSupport.cc:867
void idSkipZeroes(ideal ide)
gives an ideal/module the minimal possible size
void rChangeCurrRing(ring r)
Definition: polys.cc:12
INLINE_THIS void Init(int l=0)
ideal idInit(int idsize, int rank)
initialise an ideal / module
Definition: simpleideals.cc:38
void rSetWeightVec(ring r, int64 *wv)
Definition: ring.cc:5104
intvec * getNthRow(intvec *v, int n)
Definition: walkSupport.cc:168
BOOLEAN noPolysWithMoreThanTwoTerms(ideal Gw)
Definition: walkSupport.cc:383
const Variable & v
< [in] a sqrfree bivariate poly
Definition: facBivar.h:37
int getMaxPosOfNthRow(intvec *v, int n)
Definition: walkSupport.cc:83
#define NULL
Definition: omList.c:10
matrix id_Module2formatedMatrix(ideal mod, int rows, int cols, const ring R)
slists * lists
Definition: mpr_numeric.h:146
int length() const
Definition: intvec.h:86
int gcd(int a, int b)
Definition: walkSupport.cc:839
BOOLEAN currwOnBorder64(ideal G, int64vec *currw64)
Definition: walkSupport.cc:353
int DIFFspy(ideal G)
Definition: walkSupport.cc:410
const CanonicalForm & w
Definition: facAbsFact.cc:55
#define pDelete(p_ptr)
Definition: polys.h:169
int cols() const
Definition: intvec.h:87
#define pGetExpV(p, e)
Gets a copy of (resp. set) the exponent vector, where e is assumed to point to (r->N +1)*sizeof(long)...
Definition: polys.h:96
int rtyp
Definition: subexpr.h:92
int64vec * rGetGlobalOrderMatrix(ring r)
#define idealSize(I)
Definition: walkSupport.h:35
#define pNext(p)
Definition: monomials.h:43
ideal idLift(ideal mod, ideal submod, ideal *rest, BOOLEAN goodShape, BOOLEAN isSB, BOOLEAN divide, matrix *unit)
Definition: ideals.cc:891
omBin slists_bin
Definition: lists.cc:23
void gett64(intvec *listw, int64vec *currw64, int64vec *targw64, int64 &tvec0, int64 &tvec1)
Definition: walkSupport.cc:484
intvec * int64VecToIntVec(int64vec *source)
static int64 scalarProduct64(int64vec *a, int64vec *b)
Definition: walkSupport.cc:269
poly getNthPolyOfId(ideal I, int n)
Definition: walkSupport.cc:689
polyrec * poly
Definition: hilb.h:10
int BOOLEAN
Definition: auxiliary.h:86
#define IMATELEM(M, I, J)
Definition: intvec.h:77
const poly b
Definition: syzextra.cc:213
int64 * iv64GetVec()
Definition: int64vec.h:61
int l
Definition: cfEzgcd.cc:94
return result
Definition: facAbsBiFact.cc:76
int64 getInvEps64(ideal G, intvec *targm, int pertdeg)
Definition: walkSupport.cc:112