source: git/kernel/walkSupport.cc @ 762407

spielwiese
Last change on this file since 762407 was 762407, checked in by Oleksandr Motsak <motsak@…>, 12 years ago
config.h is for sources files only FIX: config.h should only be used by source (not from inside kernel/mod2.h!) NOTE: each source file should better include mod2.h right after config.h, while headers should better not include mod2.h.
  • Property mode set to 100644
File size: 31.9 KB
Line 
1#include <string.h>
2#include "config.h"
3#include <kernel/mod2.h>
4#include <misc/intvec.h>
5#include <misc/int64vec.h>
6#include <kernel/polys.h>
7#include <kernel/ideals.h>
8#include <polys/monomials/ring.h>
9#include <kernel/walkSupport.h>
10#include <polys/prCopy.h>
11#include <kernel/kstd1.h>
12#include <polys/matpol.h>
13#include <math.h>
14
15extern 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 containg the exponet 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  pGetExpV(p,e);
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  pGetExpV(p,e);
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 reslut 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 temp;
1027  int i=0;
1028  while(r->order[i]!=0 && pos2<n)
1029  {
1030    pos2=pos2+r->block1[i] - r->block0[i];
1031
1032    if(r->order[i]==ringorder_lp)
1033    {
1034      temp=pos1;
1035      for(int j=pos1; j<=pos2; j++)
1036        (*res)[j*n+j]=(int64)1;
1037    }
1038    else if(r->order[i]==ringorder_dp)
1039    {
1040      for(int j=pos1;j<=pos2;j++)
1041        (*res)[pos1*n+j]=(int64)1;
1042      for(int j=1;j<=(pos2-pos1);j++)
1043        (*res)[(pos1+j)*n+(pos2+1-j)]=(int64)-1;
1044    }
1045    else if(r->order[i]==ringorder_Dp)
1046    {
1047      for(int j=pos1;j<=pos2;j++)
1048        (*res)[pos1*n+j]=(int64)1;
1049      for(int j=1;j<=(pos2-pos1);j++)
1050        (*res)[(pos1+j)*n+(pos1+j-1)]=(int64)1;
1051    }
1052    else if(r->order[i]==ringorder_wp)
1053    {
1054      int* weights=r->wvhdl[i];
1055      for(int j=pos1;j<=pos2;j++)
1056        (*res)[pos1*n+j]=(int64)weights[j-pos1];
1057      for(int j=1;j<=(pos2-pos1);j++)
1058        (*res)[(pos1+j)*n+(pos2+1-j)]=(int64)-1;
1059    }
1060    else if(r->order[i]==ringorder_Wp)
1061    {
1062      int* weights=r->wvhdl[i];
1063      for(int j=pos1;j<=pos2;j++)
1064        (*res)[pos1*n+j]=(int64)weights[j-pos1];
1065      for(int j=1;j<=(pos2-pos1);j++)
1066        (*res)[(pos1+j)*n+(pos1+j-1)]=(int64)1;
1067    }
1068
1069    else if(r->order[0]==ringorder_M)
1070    {
1071      int* weights=r->wvhdl[0];
1072      for(int j=pos1;j<((pos2+1)*(pos2+1));j++)
1073        (*res)[j]=(int64)weights[j];
1074    }
1075
1076    pos1=pos2+1;
1077    pos2=pos2+1;
1078    i++;
1079  }
1080
1081  return(res);
1082}
1083
1084///////////////////////////////////////////////////////////////////
1085
1086
1087///////////////////////////////////////////////////////////////////
1088//rGetGlobalOrderWeightVec
1089///////////////////////////////////////////////////////////////////
1090//Description: returns a weight vector corresponding to the first
1091//row of a matrix corresponding to the order of r
1092///////////////////////////////////////////////////////////////////
1093//Uses: none
1094///////////////////////////////////////////////////////////////////
1095
1096int64vec* rGetGlobalOrderWeightVec(ring r)
1097{
1098  int n=rVar(r);
1099  int64vec* res=new int64vec(n);
1100
1101  if (rHasLocalOrMixedOrdering(r)) return res;
1102
1103  int length;
1104
1105  if(r->order[0]==ringorder_lp)
1106  {
1107    (*res)[0]=(int64)1;
1108  }
1109  else if( (r->order[0]==ringorder_dp) || (r->order[0]==ringorder_Dp) )
1110  {
1111    length=r->block1[0] - r->block0[0];
1112    for (int j=0;j<=length;j++)
1113      (*res)[j]=(int64)1;
1114  }
1115  else if( (r->order[0]==ringorder_wp) || (r->order[0]==ringorder_Wp) ||
1116      (r->order[0]==ringorder_a)  || (r->order[0]==ringorder_M)    )
1117  {
1118    int* weights=r->wvhdl[0];
1119    length=r->block1[0] - r->block0[0];
1120    for (int j=0;j<=length;j++)
1121      (*res)[j]=(int64)weights[j];
1122  }
1123  else if( (r->order[0]==ringorder_a64) )
1124  {
1125    int64* weights=(int64*)r->wvhdl[0];
1126    length=r->block1[0] - r->block0[0];
1127    for (int j=0;j<=length;j++)
1128      (*res)[j]=weights[j];
1129  }
1130
1131  return(res);
1132}
1133
1134///////////////////////////////////////////////////////////////////
1135
1136
1137///////////////////////////////////////////////////////////////////
1138//sortRedSB
1139///////////////////////////////////////////////////////////////////
1140//Description: sorts a reduced GB of an ideal after the leading
1141//terms of the polynomials with the smallest one first
1142///////////////////////////////////////////////////////////////////
1143//Assumes: that the given input is a minimal GB
1144///////////////////////////////////////////////////////////////////
1145//Uses:idealSize,idCopy,pLmCmp
1146///////////////////////////////////////////////////////////////////
1147
1148ideal sortRedSB(ideal G)
1149{
1150  int s=idealSize(G);
1151  poly* m=G->m;
1152  poly p,q;
1153  for (int i=0; i<(s-1); i++)
1154  {
1155    for (int j=0; j<((s-1)-i); j++)
1156    {
1157      p=m[j];
1158      q=m[j+1];
1159      if (pLmCmp(p,q)==1)
1160      {
1161        m[j+1]=p;
1162        m[j]=q;
1163      }
1164    }
1165  }
1166  return(G);
1167}
1168
1169///////////////////////////////////////////////////////////////////
1170
1171
1172///////////////////////////////////////////////////////////////////
1173//int64VecToIntVec
1174///////////////////////////////////////////////////////////////////
1175//Description: converts an int64vec to an intvec
1176// deletes the input
1177///////////////////////////////////////////////////////////////////
1178//Assumes: that the int64vec contains no entries larger than 2^32
1179///////////////////////////////////////////////////////////////////
1180//Uses: none
1181///////////////////////////////////////////////////////////////////
1182
1183intvec* int64VecToIntVec(int64vec* source)
1184{
1185  int r=source->rows();
1186  int c=source->cols();
1187  intvec* res=new intvec(r,c,0);
1188  for(int i=0;i<r;i++){
1189    for(int j=0;j<c;j++){
1190      (*res)[i*c+j]=(int)(*source)[i*c+j];
1191    }
1192  }
1193  delete source;
1194  return(res);
1195}
1196
1197///////////////////////////////////////////////////////////////////
Note: See TracBrowser for help on using the repository browser.