source: git/kernel/walkSupport.cc @ fbc7cb

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