source: git/kernel/walkSupport.cc @ 07625cb

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