source: git/kernel/walkSupport.cc @ 06662e

spielwiese
Last change on this file since 06662e was 098f98f, checked in by Hans Schönemann <hannes@…>, 14 years ago
IDTYP etc -> ipid.h git-svn-id: file:///usr/local/Singular/svn/trunk@12409 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<rVar(currRing);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<rVar(currRing);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=rVar(currRing);
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<rVar(currRing); 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<rVar(currRing); 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<rVar(currRing); i++)
610  {
611    if( (((*a)[i])>=0 && ((*b)[i])>=0) ||
612        (((*a)[i])<0 && ((*b)[i])<0) )
613    {
614      if( (abs64((*a)[i]))>abs64((*nextweight)[i]) ||
615          (abs64((*b)[i]))>abs64((*nextweight)[i])
616        )
617      {
618        overflow_error=9;
619        break;
620      }
621    }
622  }
623
624  //to reduce common factors of nextweight
625  int s=iv64Size(nextweight);
626  int64 g,temp;
627  g=(*nextweight)[0];
628  for (int i=1; i<s; i++)
629  {
630    temp=(*nextweight)[i];
631    g=gcd64(g,temp);
632    if (g==1) break;
633  }
634
635  if (g!=1) *nextweight /= g;
636
637  return(nextweight);
638}
639
640///////////////////////////////////////////////////////////////////
641
642
643
644///////////////////////////////////////////////////////////////////
645//FUNCTIONS NOT ORIGINATING FROM THE SINGULAR IMPLEMENTATION CODE
646///////////////////////////////////////////////////////////////////
647
648
649///////////////////////////////////////////////////////////////////
650//getNthPolyOfId
651///////////////////////////////////////////////////////////////////
652//Description: returns the nth poly of ideal I
653///////////////////////////////////////////////////////////////////
654poly getNthPolyOfId(ideal I,int n)
655{
656  if(0<n && n<=((int)I->ncols))
657  {
658    return (I->m)[n-1];
659  }
660  else
661  {
662    return(NULL);
663  }
664}
665
666///////////////////////////////////////////////////////////////////
667
668
669///////////////////////////////////////////////////////////////////
670//idealSize
671///////////////////////////////////////////////////////////////////
672//Description: returns the number of generator of input ideal I
673///////////////////////////////////////////////////////////////////
674//Uses: none
675///////////////////////////////////////////////////////////////////
676// #define idealSize(I) IDELEMS(I)
677///////////////////////////////////////////////////////////////////
678
679
680///////////////////////////////////////////////////////////////////
681//ivSize
682///////////////////////////////////////////////////////////////////
683//Description: returns the number of entries of v
684///////////////////////////////////////////////////////////////////
685//Uses: none
686///////////////////////////////////////////////////////////////////
687
688// inline int ivSize(intvec* v){ return((v->rows())*(v->cols())); }
689
690///////////////////////////////////////////////////////////////////
691
692
693///////////////////////////////////////////////////////////////////
694//iv64Size
695///////////////////////////////////////////////////////////////////
696//Description: returns the number of entries of v
697///////////////////////////////////////////////////////////////////
698//Uses: none
699///////////////////////////////////////////////////////////////////
700
701// int iv64Size(int64vec* v){ return((v->rows())*(v->cols())); }
702
703///////////////////////////////////////////////////////////////////
704
705
706///////////////////////////////////////////////////////////////////
707//leadExp
708///////////////////////////////////////////////////////////////////
709//Description: returns an intvec containg the exponet vector of p
710///////////////////////////////////////////////////////////////////
711//Uses: sizeof,omAlloc,omFree
712///////////////////////////////////////////////////////////////////
713
714intvec* leadExp(poly p)
715{
716  int N=rVar(currRing);
717  int *e=(int*)omAlloc((N+1)*sizeof(int));
718  pGetExpV(p,e);
719  intvec* iv=new intvec(N);
720  for(int i=N;i>0;i--) { (*iv)[i-1]=e[i];}
721  omFree(e);
722  return(iv);
723}
724
725///////////////////////////////////////////////////////////////////
726
727
728///////////////////////////////////////////////////////////////////
729//leadExp64
730///////////////////////////////////////////////////////////////////
731//Description: returns an int64vec containing the exponent
732//vector of p
733///////////////////////////////////////////////////////////////////
734//Uses: sizeof,omAlloc,omFree
735///////////////////////////////////////////////////////////////////
736
737int64vec* leadExp64(poly p)
738{
739  int N=rVar(currRing);
740  int *e=(int*)omAlloc((N+1)*sizeof(int));
741  pGetExpV(p,e);
742  int64vec* iv64=new int64vec(N);
743  for(int i=N;i>0;i--) { (*iv64)[i-1]=(int64)e[i];}
744  omFree(e);
745  return(iv64);
746}
747
748///////////////////////////////////////////////////////////////////
749
750
751///////////////////////////////////////////////////////////////////
752//setPosOfIM
753///////////////////////////////////////////////////////////////////
754//Description: sets entry i,j of im to val
755///////////////////////////////////////////////////////////////////
756//Uses: none
757///////////////////////////////////////////////////////////////////
758
759//void setPosOfIM(intvec* im,int i,int j,int val){
760//  int r=im->rows();
761//  int c=im->cols();
762//  if( (0<i && i<=r) && (0<j && j<=c) ){
763//    (*im)[(i-1)*c+j-1]=val;
764//    }
765//  return;
766//}
767
768///////////////////////////////////////////////////////////////////
769
770
771///////////////////////////////////////////////////////////////////
772//scalarProduct
773///////////////////////////////////////////////////////////////////
774//Description: returns the scalar product of intvecs a and b
775///////////////////////////////////////////////////////////////////
776//Uses: none
777///////////////////////////////////////////////////////////////////
778
779static inline long  scalarProduct(intvec* a, intvec* b)
780{
781  assume( a->length() ==  b->length());
782  int i, n = a->length();
783  long result = 0;
784  for(i=n-1; i>=0; i--)
785    result += (*a)[i] * (*b)[i];
786  return result;
787}
788
789///////////////////////////////////////////////////////////////////
790
791
792///////////////////////////////////////////////////////////////////
793//scalarProduct64
794///////////////////////////////////////////////////////////////////
795//Description: returns the scalar product of int64vecs a and b
796///////////////////////////////////////////////////////////////////
797//Assume: Overflow tests assumes input has nonnegative entries
798///////////////////////////////////////////////////////////////////
799//Uses: none
800///////////////////////////////////////////////////////////////////
801
802static inline int64  scalarProduct64(int64vec* a, int64vec* b)
803{
804  assume( a->length() ==  b->length());
805  int i, n = a->length();
806  int64 result = 0;
807  int64 temp1,temp2;
808  for(i=n-1; i>=0; i--)
809  {
810    temp1=(*a)[i] * (*b)[i];
811    if((*a)[i]!=0 && (temp1/(*a)[i])!=(*b)[i]) overflow_error=1;
812
813    temp2=result;
814    result += temp1;
815
816    //since both vectors always have nonnegative entries in init64
817    //result should be >=temp2
818    if(temp2>result) overflow_error=2;
819  }
820
821  return result;
822}
823
824///////////////////////////////////////////////////////////////////
825
826
827///////////////////////////////////////////////////////////////////
828//gcd
829///////////////////////////////////////////////////////////////////
830//Description: returns the gcd of a and b
831///////////////////////////////////////////////////////////////////
832//Uses: none
833///////////////////////////////////////////////////////////////////
834
835int gcd(int a, int b)
836{
837  int r, p0 = a, p1 = b;
838  if(p0 < 0)
839    p0 = -p0;
840
841  if(p1 < 0)
842    p1 = -p1;
843  while(p1 != 0)
844  {
845    r = p0 % p1;
846    p0 = p1;
847    p1 = r;
848  }
849  return p0;
850}
851
852///////////////////////////////////////////////////////////////////
853
854
855///////////////////////////////////////////////////////////////////
856//gcd64
857///////////////////////////////////////////////////////////////////
858//Description: returns the gcd of a and b
859///////////////////////////////////////////////////////////////////
860//Uses: none
861///////////////////////////////////////////////////////////////////
862
863int64 gcd64(int64 a, int64 b)
864{
865  int64 r, p0 = a, p1 = b;
866  if(p0 < 0)
867    p0 = -p0;
868
869  if(p1 < 0)
870    p1 = -p1;
871
872  while(p1 != ((int64)0) )
873  {
874    r = p0 % p1;
875    p0 = p1;
876    p1 = r;
877  }
878
879  return p0;
880}
881
882///////////////////////////////////////////////////////////////////
883
884
885///////////////////////////////////////////////////////////////////
886//abs64
887///////////////////////////////////////////////////////////////////
888//Description: returns the absolute value of int64 i
889///////////////////////////////////////////////////////////////////
890//Uses: none
891///////////////////////////////////////////////////////////////////
892
893//int64 abs64(int64 i)
894//{
895//  if(i>=0)
896//  return(i);
897//else
898//  return((-i));
899//}
900
901///////////////////////////////////////////////////////////////////
902
903
904///////////////////////////////////////////////////////////////////
905//makeTaunList64
906///////////////////////////////////////////////////////////////////
907//Description: makes a list of an int64vec and an int64
908///////////////////////////////////////////////////////////////////
909//Uses: omAllocBin
910///////////////////////////////////////////////////////////////////
911
912#if 0
913lists makeTaunList64(int64vec *iv64,int64 i64)
914{
915  lists l=(lists)omAllocBin(slists_bin);
916  l->Init(2);
917  l->m[0].rtyp=INTVEC_CMD;
918  l->m[1].rtyp=INT_CMD;
919  l->m[0].data=(void *)iv64;
920  l->m[1].data=(void *)i64;
921  return l;
922}
923#endif
924
925///////////////////////////////////////////////////////////////////
926
927
928///////////////////////////////////////////////////////////////////
929//idStd
930///////////////////////////////////////////////////////////////////
931//Description: returns the GB of G calculated w.r.t. the order of
932//currRing
933///////////////////////////////////////////////////////////////////
934//Uses: kStd,idSkipZeroes
935///////////////////////////////////////////////////////////////////
936
937ideal idStd(ideal G)
938{
939  ideal GG = kStd(G, NULL, testHomog, NULL);
940  idSkipZeroes(GG);
941  return GG;
942}
943
944///////////////////////////////////////////////////////////////////
945
946
947///////////////////////////////////////////////////////////////////
948//idInterRed
949///////////////////////////////////////////////////////////////////
950//Description: returns the interreduction of G
951///////////////////////////////////////////////////////////////////
952//Assumes: that the input is a GB
953///////////////////////////////////////////////////////////////////
954//Uses: kInterRed,idSkipZeroes
955///////////////////////////////////////////////////////////////////
956
957ideal idInterRed(ideal G)
958{
959  assume(G != NULL);
960
961  ideal GG = kInterRedOld(G, NULL);
962  idDelete(&G);
963  return GG;
964}
965
966///////////////////////////////////////////////////////////////////
967
968
969///////////////////////////////////////////////////////////////////
970//matIdLift
971///////////////////////////////////////////////////////////////////
972//Description: yields the same reslut as lift in Singular
973///////////////////////////////////////////////////////////////////
974//Uses: idLift,idModule2formatedMatrix
975///////////////////////////////////////////////////////////////////
976
977matrix matIdLift(ideal Gomega, ideal M)
978{
979  ideal Mtmp = idLift(Gomega, M, NULL, FALSE, FALSE, TRUE, NULL);
980  int rows=IDELEMS(Gomega);
981  int cols=IDELEMS(Mtmp);
982  matrix res=idModule2formatedMatrix(Mtmp,rows,cols);
983  return res;
984}
985///////////////////////////////////////////////////////////////////
986
987
988///////////////////////////////////////////////////////////////////
989//rCopyAndChangeA
990///////////////////////////////////////////////////////////////////
991//Description: changes currRing to a copy of currRing with the
992//int64vec w instead of the old one
993///////////////////////////////////////////////////////////////////
994//Assumes: that currRing is alterd as mentioned in rCopy0AndAddA
995///////////////////////////////////////////////////////////////////
996//Uses: rCopy0,rComplete,rChangeCurrRing,rSetWeightVec
997///////////////////////////////////////////////////////////////////
998
999void rCopyAndChangeA(int64vec* w)
1000{
1001  ring rnew=rCopy0(currRing);
1002  rComplete(rnew);
1003  rSetWeightVec(rnew,w->iv64GetVec());
1004  rChangeCurrRing(rnew);
1005}
1006
1007///////////////////////////////////////////////////////////////////
1008
1009
1010///////////////////////////////////////////////////////////////////
1011//rCopy0AndAddA
1012///////////////////////////////////////////////////////////////////
1013//Description: returns a copy of currRing with the
1014//int64vec wv64 added in front of the current order
1015//i.e. if old order was O, then the new one will be (A(w),O)
1016///////////////////////////////////////////////////////////////////
1017//Uses: omAllocBin,memcpy4,sizeof,nCopy,omStrDup,omMemDup,
1018//idrCopyR_NoSort
1019///////////////////////////////////////////////////////////////////
1020
1021ring rCopy0AndAddA(ring r, int64vec *wv64, BOOLEAN copy_qideal,
1022                   BOOLEAN copy_ordering)
1023{
1024
1025  if (r == NULL) return NULL;
1026  int i,j;
1027  ring res=(ring)omAllocBin(sip_sring_bin);
1028
1029  memcpy4(res,r,sizeof(ip_sring));
1030  res->VarOffset = NULL;
1031  res->ref=0;
1032  if (r->algring!=NULL)
1033    r->algring->ref++;
1034  if (r->parameter!=NULL)
1035  {
1036    res->minpoly=nCopy(r->minpoly);
1037    int l=rPar(r);
1038    res->parameter=(char **)omAlloc(l*sizeof(char_ptr));
1039    int i;
1040    for(i=0;i<rPar(r);i++)
1041    {
1042      res->parameter[i]=omStrDup(r->parameter[i]);
1043    }
1044  }
1045
1046  i=rBlocks(r);
1047  res->wvhdl   = (int **)omAlloc((i+1) * sizeof(int_ptr));
1048  res->order   = (int *) omAlloc((i+1) * sizeof(int));
1049  res->block0  = (int *) omAlloc((i+1) * sizeof(int));
1050  res->block1  = (int *) omAlloc((i+1) * sizeof(int));
1051  for (j=0; j<i; j++)
1052  {
1053    if (r->wvhdl[j]!=NULL)
1054    {
1055      res->wvhdl[j+1] = (int*) omMemDup(r->wvhdl[j]);
1056    }
1057    else
1058      res->wvhdl[j+1]=NULL;
1059  }
1060  for (j=0;j<i;j++)
1061  {
1062    res->order[j+1]=r->order[j];
1063    res->block0[j+1]=r->block0[j];
1064    res->block1[j+1]=r->block1[j];
1065  }
1066
1067 //the added A
1068  res->order[0]=ringorder_a64;
1069  int length=wv64->rows();
1070  int64 *A=(int64 *)omAlloc(length*sizeof(int64));
1071  for(j=length-1;j>=0;j--)
1072  {
1073    A[j]=(*wv64)[j];
1074  }
1075  res->wvhdl[0]=(int *)A;
1076  res->block0[0]=1;
1077  res->block1[0]=length;
1078
1079  res->names   = (char **)omAlloc0(rVar(r) * sizeof(char_ptr));
1080  for (i=rVar(res)-1;i>=0; i--)
1081  {
1082    res->names[i] = omStrDup(r->names[i]);
1083  }
1084  res->idroot = NULL;
1085  if (r->qideal!=NULL)
1086  {
1087    if (copy_qideal) res->qideal= idrCopyR_NoSort(r->qideal, r);
1088    else res->qideal = NULL;
1089  }
1090
1091  return res;
1092}
1093
1094///////////////////////////////////////////////////////////////////
1095
1096
1097///////////////////////////////////////////////////////////////////
1098//rGetGlobalOrderMatrix
1099///////////////////////////////////////////////////////////////////
1100//Description: returns a matrix corresponding to the order of r
1101///////////////////////////////////////////////////////////////////
1102//Assumes: the order of r is a combination of the orders M,lp,dp,
1103//Dp,wp,Wp,C
1104///////////////////////////////////////////////////////////////////
1105//Uses: none
1106///////////////////////////////////////////////////////////////////
1107
1108int64vec* rGetGlobalOrderMatrix(ring r)
1109{
1110  int n=rVar(r);
1111  int64vec* res=new int64vec(n,n,(int64)0);
1112  if (rHasLocalOrMixedOrdering(r)) return res;
1113  int pos1=0;
1114  int pos2=0;
1115  int temp;
1116  int i=0;
1117  while(r->order[i]!=0 && pos2<n)
1118  {
1119    pos2=pos2+r->block1[i] - r->block0[i];
1120
1121    if(r->order[i]==ringorder_lp)
1122    {
1123      temp=pos1;
1124      for(int j=pos1; j<=pos2; j++)
1125        (*res)[j*n+j]=(int64)1;
1126    }
1127    else if(r->order[i]==ringorder_dp)
1128    {
1129      for(int j=pos1;j<=pos2;j++)
1130        (*res)[pos1*n+j]=(int64)1;
1131      for(int j=1;j<=(pos2-pos1);j++)
1132        (*res)[(pos1+j)*n+(pos2+1-j)]=(int64)-1;
1133    }
1134    else if(r->order[i]==ringorder_Dp)
1135    {
1136      for(int j=pos1;j<=pos2;j++)
1137        (*res)[pos1*n+j]=(int64)1;
1138      for(int j=1;j<=(pos2-pos1);j++)
1139        (*res)[(pos1+j)*n+(pos1+j-1)]=(int64)1;
1140    }
1141    else if(r->order[i]==ringorder_wp)
1142    {
1143      int* weights=r->wvhdl[i];
1144      for(int j=pos1;j<=pos2;j++)
1145        (*res)[pos1*n+j]=(int64)weights[j-pos1];
1146      for(int j=1;j<=(pos2-pos1);j++)
1147        (*res)[(pos1+j)*n+(pos2+1-j)]=(int64)-1;
1148    }
1149    else if(r->order[i]==ringorder_Wp)
1150    {
1151      int* weights=r->wvhdl[i];
1152      for(int j=pos1;j<=pos2;j++)
1153        (*res)[pos1*n+j]=(int64)weights[j-pos1];
1154      for(int j=1;j<=(pos2-pos1);j++)
1155        (*res)[(pos1+j)*n+(pos1+j-1)]=(int64)1;
1156    }
1157
1158    else if(r->order[0]==ringorder_M)
1159    {
1160      int* weights=r->wvhdl[0];
1161      for(int j=pos1;j<((pos2+1)*(pos2+1));j++)
1162        (*res)[j]=(int64)weights[j];
1163    }
1164
1165    pos1=pos2+1;
1166    pos2=pos2+1;
1167    i++;
1168  }
1169
1170  return(res);
1171}
1172
1173///////////////////////////////////////////////////////////////////
1174
1175
1176///////////////////////////////////////////////////////////////////
1177//rGetGlobalOrderWeightVec
1178///////////////////////////////////////////////////////////////////
1179//Description: returns a weight vector corresponding to the first
1180//row of a matrix corresponding to the order of r
1181///////////////////////////////////////////////////////////////////
1182//Uses: none
1183///////////////////////////////////////////////////////////////////
1184
1185int64vec* rGetGlobalOrderWeightVec(ring r)
1186{
1187  int n=rVar(r);
1188  int64vec* res=new int64vec(n);
1189
1190  if (rHasLocalOrMixedOrdering(r)) return res;
1191
1192  int length;
1193
1194  if(r->order[0]==ringorder_lp)
1195  {
1196    (*res)[0]=(int64)1;
1197  }
1198  else if( (r->order[0]==ringorder_dp) || (r->order[0]==ringorder_Dp) )
1199  {
1200    length=r->block1[0] - r->block0[0];
1201    for (int j=0;j<=length;j++)
1202      (*res)[j]=(int64)1;
1203  }
1204  else if( (r->order[0]==ringorder_wp) || (r->order[0]==ringorder_Wp) ||
1205      (r->order[0]==ringorder_a)  || (r->order[0]==ringorder_M)    )
1206  {
1207    int* weights=r->wvhdl[0];
1208    length=r->block1[0] - r->block0[0];
1209    for (int j=0;j<=length;j++)
1210      (*res)[j]=(int64)weights[j];
1211  }
1212  else if( (r->order[0]==ringorder_a64) )
1213  {
1214    int64* weights=(int64*)r->wvhdl[0];
1215    length=r->block1[0] - r->block0[0];
1216    for (int j=0;j<=length;j++)
1217      (*res)[j]=weights[j];
1218  }
1219
1220  return(res);
1221}
1222
1223///////////////////////////////////////////////////////////////////
1224
1225
1226///////////////////////////////////////////////////////////////////
1227//sortRedSB
1228///////////////////////////////////////////////////////////////////
1229//Description: sorts a reduced GB of an ideal after the leading
1230//terms of the polynomials with the smallest one first
1231///////////////////////////////////////////////////////////////////
1232//Assumes: that the given input is a minimal GB
1233///////////////////////////////////////////////////////////////////
1234//Uses:idealSize,idCopy,pLmCmp
1235///////////////////////////////////////////////////////////////////
1236
1237ideal sortRedSB(ideal G)
1238{
1239  int s=idealSize(G);
1240  poly* m=G->m;
1241  poly p,q;
1242  for (int i=0; i<(s-1); i++)
1243  {
1244    for (int j=0; j<((s-1)-i); j++)
1245    {
1246      p=m[j];
1247      q=m[j+1];
1248      if (pLmCmp(p,q)==1)
1249      {
1250        m[j+1]=p;
1251        m[j]=q;
1252      }
1253    }
1254  }
1255  return(G);
1256}
1257
1258///////////////////////////////////////////////////////////////////
1259
1260
1261///////////////////////////////////////////////////////////////////
1262//int64VecToIntVec
1263///////////////////////////////////////////////////////////////////
1264//Description: converts an int64vec to an intvec
1265// deletes the input
1266///////////////////////////////////////////////////////////////////
1267//Assumes: that the int64vec contains no entries larger than 2^32
1268///////////////////////////////////////////////////////////////////
1269//Uses: none
1270///////////////////////////////////////////////////////////////////
1271
1272intvec* int64VecToIntVec(int64vec* source)
1273{
1274  int r=source->rows();
1275  int c=source->cols();
1276  intvec* res=new intvec(r,c,0);
1277  for(int i=0;i<r;i++){
1278    for(int j=0;j<c;j++){
1279      (*res)[i*c+j]=(int)(*source)[i*c+j];
1280    }
1281  }
1282  delete source;
1283  return(res);
1284}
1285
1286///////////////////////////////////////////////////////////////////
Note: See TracBrowser for help on using the repository browser.