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"
12#include "kernel/groebner_walk/walkSupport.h"
13#include "kernel/GBEngine/kstd1.h"
14
15EXTERN_VAR BOOLEAN overflow_error;
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
35int 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
54int 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
80int getMaxPosOfNthRow(intvec *v,int n)
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
109int64 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
141int 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
165intvec* getNthRow(intvec *v, int n)
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
181int64vec* getNthRow64(intvec *v, int n)
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
209void 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
266static inline int64 scalarProduct64(int64vec* a, int64vec* b)
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
299ideal 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
350BOOLEAN currwOnBorder64(ideal G, int64vec* currw64)
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
380BOOLEAN noPolysWithMoreThanTwoTerms(ideal Gw)
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
407int 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
435intvec* 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
481void 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
560void 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
604int64vec* nextw64(int64vec* currw, int64vec* targw,
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///////////////////////////////////////////////////////////////////
686poly 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 containing the exponent vector of p
742///////////////////////////////////////////////////////////////////
743//Uses: sizeof,omAlloc,omFree
744///////////////////////////////////////////////////////////////////
745
746intvec* leadExp(poly p)
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
769int64vec* leadExp64(poly p)
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
811static 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
836int 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
864int64 gcd64(int64 a, int64 b)
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
914lists makeTaunList64(int64vec *iv64,int64 i64)
915{
916 lists l=(lists)omAllocBin(slists_bin);
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
938ideal 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
958ideal 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 result as lift in Singular
974///////////////////////////////////////////////////////////////////
975//Uses: idLift,idModule2formatedMatrix
976///////////////////////////////////////////////////////////////////
977
978matrix 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);
983 matrix res=id_Module2formatedMatrix(Mtmp,rows,cols,currRing);
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
1000void rCopyAndChangeA(int64vec* w)
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
1019int64vec* rGetGlobalOrderMatrix(ring r)
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
1094int64vec* rGetGlobalOrderWeightVec(ring r)
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
1146ideal 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
1181intvec* int64VecToIntVec(int64vec* source)
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///////////////////////////////////////////////////////////////////
abs
Rational abs(const Rational &a)
Definition: GMPrat.cc:436
int64
long int64
Definition: auxiliary.h:68
BOOLEAN
int BOOLEAN
Definition: auxiliary.h:87
TRUE
#define TRUE
Definition: auxiliary.h:100
FALSE
#define FALSE
Definition: auxiliary.h:96
N
const CanonicalForm CFMap CFMap & N
Definition: cfEzgcd.cc:56
l
int l
Definition: cfEzgcd.cc:100
m
int m
Definition: cfEzgcd.cc:128
for
for(int i=0;i<=n;i++) degsf[i]
Definition: cfEzgcd.cc:72
i
int i
Definition: cfEzgcd.cc:132
p
int p
Definition: cfModGcd.cc:4078
g
g
Definition: cfModGcd.cc:4090
GG
const CanonicalForm & GG
Definition: cfModGcd.cc:4076
b
CanonicalForm b
Definition: cfModGcd.cc:4103
int64vec::length
int length() const
Definition: int64vec.h:64
int64vec::rows
int rows() const
Definition: int64vec.h:66
int64vec::cols
int cols() const
Definition: int64vec.h:65
intvec
Definition: intvec.h:23
intvec::length
int length() const
Definition: intvec.h:94
intvec::rows
int rows() const
Definition: intvec.h:96
ip_smatrix
Definition: matpol.h:15
slists
Definition: lists.h:24
result
return result
Definition: facAbsBiFact.cc:75
s
const CanonicalForm int s
Definition: facAbsFact.cc:51
res
CanonicalForm res
Definition: facAbsFact.cc:60
w
const CanonicalForm & w
Definition: facAbsFact.cc:51
v
const Variable & v
< [in] a sqrfree bivariate poly
Definition: facBivar.h:39
j
int j
Definition: facHensel.cc:110
EXTERN_VAR
#define EXTERN_VAR
Definition: globaldefs.h:6
idLift
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
idDelete
#define idDelete(H)
delete an ideal
Definition: ideals.h:29
iv64Add
int64vec * iv64Add(int64vec *a, int64vec *b)
Definition: int64vec.cc:172
iv64Sub
int64vec * iv64Sub(int64vec *a, int64vec *b)
Definition: int64vec.cc:202
iv64Copy
int64vec * iv64Copy(int64vec *o)
Definition: int64vec.h:84
length
static BOOLEAN length(leftv result, leftv arg)
Definition: interval.cc:257
ivSub
intvec * ivSub(intvec *a, intvec *b)
Definition: intvec.cc:297
IMATELEM
#define IMATELEM(M, I, J)
Definition: intvec.h:85
G
STATIC_VAR TreeM * G
Definition: janet.cc:31
kInterRedOld
ideal kInterRedOld(ideal F, ideal Q)
Definition: kstd1.cc:3409
kStd
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:2447
slists_bin
VAR omBin slists_bin
Definition: lists.cc:23
assume
#define assume(x)
Definition: mod2.h:389
pIter
#define pIter(p)
Definition: monomials.h:37
pNext
#define pNext(p)
Definition: monomials.h:36
lists
slists * lists
Definition: mpr_numeric.h:146
omAlloc
#define omAlloc(size)
Definition: omAllocDecl.h:210
omAllocBin
#define omAllocBin(bin)
Definition: omAllocDecl.h:205
omFree
#define omFree(addr)
Definition: omAllocDecl.h:261
NULL
#define NULL
Definition: omList.c:12
pLength
static int pLength(poly a)
Definition: p_polys.h:188
p_GetExpV
static void p_GetExpV(poly p, int *ev, const ring r)
Definition: p_polys.h:1518
rChangeCurrRing
void rChangeCurrRing(ring r)
Definition: polys.cc:15
currRing
VAR ring currRing
Widely used global variable which specifies the current polynomial ring for Singular interpreter and ...
Definition: polys.cc:13
polys.h
Compatibility layer for legacy polynomial operations (over currRing)
pAdd
#define pAdd(p, q)
Definition: polys.h:203
pTotaldegree
static long pTotaldegree(poly p)
Definition: polys.h:282
pDelete
#define pDelete(p_ptr)
Definition: polys.h:186
pHead
#define pHead(p)
returns newly allocated copy of Lm(p), coef is copied, next=NULL, p might be NULL
Definition: polys.h:67
pLmCmp
#define pLmCmp(p, q)
returns 0|1|-1 if p=q|p>q|p<q w.r.t monomial ordering
Definition: polys.h:105
rSetWeightVec
void rSetWeightVec(ring r, int64 *wv)
Definition: ring.cc:5232
rComplete
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:3450
rCopy0
ring rCopy0(const ring r, BOOLEAN copy_qideal, BOOLEAN copy_ordering)
Definition: ring.cc:1421
ringorder_lp
@ ringorder_lp
Definition: ring.h:77
ringorder_a
@ ringorder_a
Definition: ring.h:70
ringorder_a64
@ ringorder_a64
for int64 weights
Definition: ring.h:71
ringorder_Dp
@ ringorder_Dp
Definition: ring.h:80
ringorder_dp
@ ringorder_dp
Definition: ring.h:78
ringorder_Wp
@ ringorder_Wp
Definition: ring.h:82
ringorder_wp
@ ringorder_wp
Definition: ring.h:81
ringorder_M
@ ringorder_M
Definition: ring.h:74
rVar
static short rVar(const ring r)
#define rVar(r) (r->N)
Definition: ring.h:592
rHasLocalOrMixedOrdering
BOOLEAN rHasLocalOrMixedOrdering(const ring r)
Definition: ring.h:760
idInit
ideal idInit(int idsize, int rank)
initialise an ideal / module
Definition: simpleideals.cc:35
id_Module2formatedMatrix
matrix id_Module2formatedMatrix(ideal mod, int rows, int cols, const ring R)
Definition: simpleideals.cc:1522
idSkipZeroes
void idSkipZeroes(ideal ide)
gives an ideal/module the minimal possible size
Definition: simpleideals.cc:181
IDELEMS
#define IDELEMS(i)
Definition: simpleideals.h:23
M
#define M
Definition: sirandom.c:25
testHomog
@ testHomog
Definition: structs.h:38
INTVEC_CMD
@ INTVEC_CMD
Definition: tok.h:101
INT_CMD
@ INT_CMD
Definition: tok.h:96
rGetGlobalOrderMatrix
int64vec * rGetGlobalOrderMatrix(ring r)
Definition: walkSupport.cc:1019
nextt64
void nextt64(ideal G, int64vec *currw64, int64vec *targw64, int64 &tvec0, int64 &tvec1)
Definition: walkSupport.cc:560
getNthPolyOfId
poly getNthPolyOfId(ideal I, int n)
Definition: walkSupport.cc:686
matIdLift
matrix matIdLift(ideal Gomega, ideal M)
Definition: walkSupport.cc:978
invEpsOk64
int invEpsOk64(ideal I, intvec *targm, int pertdeg, int64 inveps64)
Definition: walkSupport.cc:141
overflow_error
EXTERN_VAR BOOLEAN overflow_error
Definition: walkSupport.cc:15
getTaun64
void getTaun64(ideal G, intvec *targm, int pertdeg, int64vec **v64, int64 &i64)
Definition: walkSupport.cc:209
getMaxTdeg
int getMaxTdeg(ideal I)
Definition: walkSupport.cc:54
DIFF
intvec * DIFF(ideal G)
Definition: walkSupport.cc:435
getMaxPosOfNthRow
int getMaxPosOfNthRow(intvec *v, int n)
Definition: walkSupport.cc:80
getNthRow
intvec * getNthRow(intvec *v, int n)
Definition: walkSupport.cc:165
currwOnBorder64
BOOLEAN currwOnBorder64(ideal G, int64vec *currw64)
Definition: walkSupport.cc:350
nextw64
int64vec * nextw64(int64vec *currw, int64vec *targw, int64 nexttvec0, int64 nexttvec1)
Definition: walkSupport.cc:604
sortRedSB
ideal sortRedSB(ideal G)
Definition: walkSupport.cc:1146
scalarProduct
static long scalarProduct(intvec *a, intvec *b)
Definition: walkSupport.cc:811
noPolysWithMoreThanTwoTerms
BOOLEAN noPolysWithMoreThanTwoTerms(ideal Gw)
Definition: walkSupport.cc:380
rCopyAndChangeA
void rCopyAndChangeA(int64vec *w)
Definition: walkSupport.cc:1000
tdeg
int tdeg(poly p)
Definition: walkSupport.cc:35
int64VecToIntVec
intvec * int64VecToIntVec(int64vec *source)
Definition: walkSupport.cc:1181
gcd64
int64 gcd64(int64 a, int64 b)
Definition: walkSupport.cc:864
init64
ideal init64(ideal G, int64vec *currw64)
Definition: walkSupport.cc:299
gett64
void gett64(intvec *listw, int64vec *currw64, int64vec *targw64, int64 &tvec0, int64 &tvec1)
Definition: walkSupport.cc:481
getNthRow64
int64vec * getNthRow64(intvec *v, int n)
Definition: walkSupport.cc:181
leadExp64
int64vec * leadExp64(poly p)
Definition: walkSupport.cc:769
leadExp
intvec * leadExp(poly p)
Definition: walkSupport.cc:746
DIFFspy
int DIFFspy(ideal G)
Definition: walkSupport.cc:407
idInterRed
ideal idInterRed(ideal G)
Definition: walkSupport.cc:958
scalarProduct64
static int64 scalarProduct64(int64vec *a, int64vec *b)
Definition: walkSupport.cc:266
getInvEps64
int64 getInvEps64(ideal G, intvec *targm, int pertdeg)
Definition: walkSupport.cc:109
rGetGlobalOrderWeightVec
int64vec * rGetGlobalOrderWeightVec(ring r)
Definition: walkSupport.cc:1094
gcd
int gcd(int a, int b)
Definition: walkSupport.cc:836
idStd
ideal idStd(ideal G)
Definition: walkSupport.cc:938
iv64Size
int iv64Size(int64vec *v)
Definition: walkSupport.h:37
ivSize
int ivSize(intvec *v)
Definition: walkSupport.h:36
idealSize
#define idealSize(I)
Definition: walkSupport.h:35
abs64
int64 abs64(int64 i)
Definition: walkSupport.h:44
Generated on Thu Apr 27 2023 08:16:02 for My Project by   doxygen 1.9.4 for Singular