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