My Project
mpr_base.cc
Go to the documentation of this file.
1 /****************************************
2 * Computer Algebra System SINGULAR *
3 ****************************************/
4 
5 /*
6  * ABSTRACT - multipolynomial resultants - resultant matrices
7  * ( sparse, dense, u-resultant solver )
8  */
9 
10 //-> includes
11 
12 
13 
14 #include "kernel/mod2.h"
15 
16 #include "misc/mylimits.h"
17 #include "misc/options.h"
18 #include "misc/intvec.h"
19 #include "misc/sirandom.h"
20 
21 #include "coeffs/numbers.h"
22 #include "coeffs/mpr_global.h"
23 
24 #include "polys/matpol.h"
25 #include "polys/sparsmat.h"
26 
27 #include "polys/clapsing.h"
28 
29 #include "kernel/polys.h"
30 #include "kernel/ideals.h"
31 
32 #include "mpr_base.h"
33 #include "mpr_numeric.h"
34 
35 #include <cmath>
36 //<-
37 
38 //%s
39 //-----------------------------------------------------------------------------
40 //-------------- sparse resultant matrix --------------------------------------
41 //-----------------------------------------------------------------------------
42 
43 //-> definitions
44 
45 //#define mprTEST
46 //#define mprMINKSUM
47 
48 #define MAXPOINTS 10000
49 #define MAXINITELEMS 256
50 #define LIFT_COOR 50000 // siRand() % LIFT_COOR gives random lift value
51 #define SCALEDOWN 100.0 // lift value scale down for linear program
52 #define MINVDIST 0.0
53 #define RVMULT 0.0001 // multiplicator for random shift vector
54 #define MAXRVVAL 50000
55 #define MAXVARS 100
56 //<-
57 
58 //-> sparse resultant matrix
59 
60 /* set of points */
61 class pointSet;
62 
63 
64 
65 /* sparse resultant matrix class */
66 class resMatrixSparse : virtual public resMatrixBase
67 {
68 public:
69  resMatrixSparse( const ideal _gls, const int special = SNONE );
71 
72  // public interface according to base class resMatrixBase
73  ideal getMatrix();
74 
75  /** Fills in resMat[][] with evpoint[] and gets determinant
76  * uRPos[i][1]: row of matrix
77  * uRPos[i][idelem+1]: col of u(0)
78  * uRPos[i][2..idelem]: col of u(1) .. u(n)
79  * i= 1 .. numSet0
80  */
81  number getDetAt( const number* evpoint );
82 
83  poly getUDet( const number* evpoint );
84 
85 private:
87 
88  void randomVector( const int dim, mprfloat shift[] );
89 
90  /** Row Content Function
91  * Finds the largest i such that F[i] is a point, F[i]= a[ij] in A[i] for some j.
92  * Returns -1 iff the point vert does not lie in a cell
93  */
94  int RC( pointSet **pQ, pointSet *E, int vert, mprfloat shift[] );
95 
96  /* Remaps a result of LP to the according point set Qi.
97  * Returns false iff remaping was not possible, otherwise true.
98  */
99  bool remapXiToPoint( const int indx, pointSet **pQ, int *set, int *vtx );
100 
101  /** create coeff matrix
102  * uRPos[i][1]: row of matrix
103  * uRPos[i][idelem+1]: col of u(0)
104  * uRPos[i][2..idelem]: col of u(1) .. u(n)
105  * i= 1 .. numSet0
106  * Returns the dimension of the matrix or -1 in case of an error
107  */
108  int createMatrix( pointSet *E );
109 
110  pointSet * minkSumAll( pointSet **pQ, int numq, int dim );
111  pointSet * minkSumTwo( pointSet *Q1, pointSet *Q2, int dim );
112 
113 private:
114  ideal gls;
115 
116  int n, idelem; // number of variables, polynoms
117  int numSet0; // number of elements in S0
118  int msize; // size of matrix
119 
121 
122  ideal rmat; // sparse matrix representation
123 
124  simplex * LP; // linear programming stuff
125 };
126 //<-
127 
128 //-> typedefs and structs
129 poly monomAt( poly p, int i );
130 
131 typedef unsigned int Coord_t;
132 
133 struct setID
134 {
135  int set;
136  int pnt;
137 };
138 
139 struct onePoint
140 {
141  Coord_t * point; // point[0] is unused, maxial dimension is MAXVARS+1
142  setID rc; // filled in by Row Content Function
143  struct onePoint * rcPnt; // filled in by Row Content Function
144 };
145 
146 typedef struct onePoint * onePointP;
147 
148 /* sparse matrix entry */
149 struct _entry
150 {
151  number num;
152  int col;
153  struct _entry * next;
154 };
155 
156 typedef struct _entry * entry;
157 //<-
158 
159 //-> class pointSet
160 class pointSet
161 {
162 private:
163  onePointP *points; // set of onePoint's, index [1..num], supports of monoms
164  bool lifted;
165 
166 public:
167  int num; // number of elements in points
168  int max; // maximal entries in points, i.e. allocated mem
169  int dim; // dimension, i.e. valid coord entries in point
170  int index; // should hold unique identifier of point set
171 
172  pointSet( const int _dim, const int _index= 0, const int count= MAXINITELEMS );
173  ~pointSet();
174 
175  // pointSet.points[i] equals pointSet[i]
176  inline onePointP operator[] ( const int index );
177 
178  /** Adds a point to pointSet, copy vert[0,...,dim] ot point[num+1][0,...,dim].
179  * Returns false, iff additional memory was allocated ( i.e. num >= max )
180  * else returns true
181  */
182  bool addPoint( const onePointP vert );
183 
184  /** Adds a point to pointSet, copy vert[0,...,dim] ot point[num+1][0,...,dim].
185  * Returns false, iff additional memory was allocated ( i.e. num >= max )
186  * else returns true
187  */
188  bool addPoint( const int * vert );
189 
190  /** Adds a point to pointSet, copy vert[0,...,dim] ot point[num+1][0,...,dim].
191  * Returns false, iff additional memory was allocated ( i.e. num >= max )
192  * else returns true
193  */
194  bool addPoint( const Coord_t * vert );
195 
196  /* Removes the point at intex indx */
197  bool removePoint( const int indx );
198 
199  /** Adds point to pointSet, iff pointSet \cap point = \emptyset.
200  * Returns true, iff added, else false.
201  */
202  bool mergeWithExp( const onePointP vert );
203 
204  /** Adds point to pointSet, iff pointSet \cap point = \emptyset.
205  * Returns true, iff added, else false.
206  */
207  bool mergeWithExp( const int * vert );
208 
209  /* Adds support of poly p to pointSet, iff pointSet \cap point = \emptyset. */
210  void mergeWithPoly( const poly p );
211 
212  /* Returns the row polynom multiplicator in vert[] */
213  void getRowMP( const int indx, int * vert );
214 
215  /* Returns index of supp(LT(p)) in pointSet. */
216  int getExpPos( const poly p );
217 
218  /** sort lex
219  */
220  void sort();
221 
222  /** Lifts the point set using sufficiently generic linear lifting
223  * homogeneous forms l[1]..l[dim] in Z. Every l[i] is of the form
224  * L1x1+...+Lnxn, for generic L1..Ln in Z.
225  *
226  * Lifting raises dimension by one!
227  */
228  void lift( int *l= NULL ); // !! increments dim by 1
229  void unlift() { dim--; lifted= false; }
230 
231 private:
232  pointSet( const pointSet & );
233 
234  /** points[a] < points[b] ? */
235  inline bool smaller( int, int );
236 
237  /** points[a] > points[b] ? */
238  inline bool larger( int, int );
239 
240  /** Checks, if more mem is needed ( i.e. num >= max ),
241  * returns false, if more mem was allocated, else true
242  */
243  inline bool checkMem();
244 };
245 //<-
246 
247 //-> class convexHull
248 /* Compute convex hull of given exponent set */
250 {
251 public:
252  convexHull( simplex * _pLP ) : pLP(_pLP) {}
254 
255  /** Computes the point sets of the convex hulls of the supports given
256  * by the polynoms in gls.
257  * Returns Q[].
258  */
259  pointSet ** newtonPolytopesP( const ideal gls );
260  ideal newtonPolytopesI( const ideal gls );
261 
262 private:
263  /** Returns true iff the support of poly pointPoly is inside the
264  * convex hull of all points given by the support of poly p.
265  */
266  bool inHull(poly p, poly pointPoly, int m, int site);
267 
268 private:
270  int n;
272 };
273 //<-
274 
275 //-> class mayanPyramidAlg
276 /* Compute all lattice points in a given convex hull */
278 {
279 public:
280  mayanPyramidAlg( simplex * _pLP ) : n((currRing->N)), pLP(_pLP) {}
282 
283  /** Drive Mayan Pyramid Algorithm.
284  * The Alg computes conv(Qi[]+shift[]).
285  */
286  pointSet * getInnerPoints( pointSet **_q_i, mprfloat _shift[] );
287 
288 private:
289 
290  /** Recursive Mayan Pyramid algorithm for directly computing MinkowskiSum
291  * lattice points for (n+1)-fold MinkowskiSum of given point sets Qi[].
292  * Recursively for range of dim: dim in [0..n); acoords[0..var) fixed.
293  * Stores only MinkowskiSum points of udist > 0: done by storeMinkowskiSumPoints.
294  */
295  void runMayanPyramid( int dim );
296 
297  /** Compute v-distance via Linear Programing
298  * Linear Program finds the v-distance of the point in accords[].
299  * The v-distance is the distance along the direction v to boundary of
300  * Minkowski Sum of Qi (here vector v is represented by shift[]).
301  * Returns the v-distance or -1.0 if an error occurred.
302  */
304 
305  /** LP for finding min/max coord in MinkowskiSum, given previous coors.
306  * Assume MinkowskiSum in non-negative quadrants
307  * coor in [0,n); fixed coords in acoords[0..coor)
308  */
309  void mn_mx_MinkowskiSum( int dim, Coord_t *minR, Coord_t *maxR );
310 
311  /** Stores point in E->points[pt], iff v-distance != 0
312  * Returns true iff point was stored, else flase
313  */
314  bool storeMinkowskiSumPoint();
315 
316 private:
320 
321  int n,idelem;
322 
324 
326 };
327 //<-
328 
329 //-> debug output stuff
330 #if defined(mprDEBUG_PROT) || defined(mprDEBUG_ALL)
331 void print_mat(mprfloat **a, int maxrow, int maxcol)
332 {
333  int i, j;
334 
335  for (i = 1; i <= maxrow; i++)
336  {
337  PrintS("[");
338  for (j = 1; j <= maxcol; j++) Print("% 7.2f, ", a[i][j]);
339  PrintS("],\n");
340  }
341 }
342 void print_bmat(mprfloat **a, int nrows, int ncols, int N, int *iposv)
343 {
344  int i, j;
345 
346  printf("Output matrix from LinProg");
347  for (i = 1; i <= nrows; i++)
348  {
349  printf("\n[ ");
350  if (i == 1) printf(" ");
351  else if (iposv[i-1] <= N) printf("X%d", iposv[i-1]);
352  else printf("Y%d", iposv[i-1]-N+1);
353  for (j = 1; j <= ncols; j++) printf(" %7.2f ",(double)a[i][j]);
354  printf(" ]");
355  } printf("\n");
356  fflush(stdout);
357 }
358 
359 void print_exp( const onePointP vert, int n )
360 {
361  int i;
362  for ( i= 1; i <= n; i++ )
363  {
364  Print(" %d",vert->point[i] );
365 #ifdef LONG_OUTPUT
366  if ( i < n ) PrintS(", ");
367 #endif
368  }
369 }
370 void print_matrix( matrix omat )
371 {
372  int i,j;
373  int val;
374  Print(" matrix m[%d][%d]=(\n",MATROWS( omat ),MATCOLS( omat ));
375  for ( i= 1; i <= MATROWS( omat ); i++ )
376  {
377  for ( j= 1; j <= MATCOLS( omat ); j++ )
378  {
379  if ( (MATELEM( omat, i, j)!=NULL)
380  && (!nIsZero(pGetCoeff( MATELEM( omat, i, j)))))
381  {
382  val= n_Int(pGetCoeff( MATELEM( omat, i, j) ), currRing->cf);
383  if ( i==MATROWS(omat) && j==MATCOLS(omat) )
384  {
385  Print("%d ",val);
386  }
387  else
388  {
389  Print("%d, ",val);
390  }
391  }
392  else
393  {
394  if ( i==MATROWS(omat) && j==MATCOLS(omat) )
395  {
396  PrintS(" 0");
397  }
398  else
399  {
400  PrintS(" 0, ");
401  }
402  }
403  }
404  PrintLn();
405  }
406  PrintS(");\n");
407 }
408 #endif
409 //<-
410 
411 //-> pointSet::*
412 pointSet::pointSet( const int _dim, const int _index, const int count )
413  : num(0), max(count), dim(_dim), index(_index)
414 {
415  int i;
416  points = (onePointP *)omAlloc( (count+1) * sizeof(onePointP) );
417  for ( i= 0; i <= max; i++ )
418  {
419  points[i]= (onePointP)omAlloc( sizeof(onePoint) );
420  points[i]->point= (Coord_t *)omAlloc0( (dim+2) * sizeof(Coord_t) );
421  }
422  lifted= false;
423 }
424 
426 {
427  int i;
428  int fdim= lifted ? dim+1 : dim+2;
429  for ( i= 0; i <= max; i++ )
430  {
431  omFreeSize( (void *) points[i]->point, fdim * sizeof(Coord_t) );
432  omFreeSize( (void *) points[i], sizeof(onePoint) );
433  }
434  omFreeSize( (void *) points, (max+1) * sizeof(onePointP) );
435 }
436 
437 inline onePointP pointSet::operator[] ( const int index_i )
438 {
439  assume( index_i > 0 && index_i <= num );
440  return points[index_i];
441 }
442 
443 inline bool pointSet::checkMem()
444 {
445  if ( num >= max )
446  {
447  int i;
448  int fdim= lifted ? dim+1 : dim+2;
449  points= (onePointP*)omReallocSize( points,
450  (max+1) * sizeof(onePointP),
451  (2*max + 1) * sizeof(onePointP) );
452  for ( i= max+1; i <= max*2; i++ )
453  {
454  points[i]= (onePointP)omAlloc( sizeof(struct onePoint) );
455  points[i]->point= (Coord_t *)omAlloc0( fdim * sizeof(Coord_t) );
456  }
457  max*= 2;
459  return false;
460  }
461  return true;
462 }
463 
464 bool pointSet::addPoint( const onePointP vert )
465 {
466  int i;
467  bool ret;
468  num++;
469  ret= checkMem();
470  points[num]->rcPnt= NULL;
471  for ( i= 1; i <= dim; i++ ) points[num]->point[i]= vert->point[i];
472  return ret;
473 }
474 
475 bool pointSet::addPoint( const int * vert )
476 {
477  int i;
478  bool ret;
479  num++;
480  ret= checkMem();
481  points[num]->rcPnt= NULL;
482  for ( i= 1; i <= dim; i++ ) points[num]->point[i]= (Coord_t) vert[i];
483  return ret;
484 }
485 
486 bool pointSet::addPoint( const Coord_t * vert )
487 {
488  int i;
489  bool ret;
490  num++;
491  ret= checkMem();
492  points[num]->rcPnt= NULL;
493  for ( i= 0; i < dim; i++ ) points[num]->point[i+1]= vert[i];
494  return ret;
495 }
496 
497 bool pointSet::removePoint( const int indx )
498 {
499  assume( indx > 0 && indx <= num );
500  if ( indx != num )
501  {
502  onePointP tmp;
503  tmp= points[indx];
504  points[indx]= points[num];
505  points[num]= tmp;
506  }
507  num--;
508 
509  return true;
510 }
511 
512 bool pointSet::mergeWithExp( const onePointP vert )
513 {
514  int i,j;
515 
516  for ( i= 1; i <= num; i++ )
517  {
518  for ( j= 1; j <= dim; j++ )
519  if ( points[i]->point[j] != vert->point[j] ) break;
520  if ( j > dim ) break;
521  }
522 
523  if ( i > num )
524  {
525  addPoint( vert );
526  return true;
527  }
528  return false;
529 }
530 
531 bool pointSet::mergeWithExp( const int * vert )
532 {
533  int i,j;
534 
535  for ( i= 1; i <= num; i++ )
536  {
537  for ( j= 1; j <= dim; j++ )
538  if ( points[i]->point[j] != (Coord_t) vert[j] ) break;
539  if ( j > dim ) break;
540  }
541 
542  if ( i > num )
543  {
544  addPoint( vert );
545  return true;
546  }
547  return false;
548 }
549 
550 void pointSet::mergeWithPoly( const poly p )
551 {
552  int i,j;
553  poly piter= p;
554  int * vert;
555  vert= (int *)omAlloc( (dim+1) * sizeof(int) );
556 
557  while ( piter )
558  {
559  p_GetExpV( piter, vert, currRing );
560 
561  for ( i= 1; i <= num; i++ )
562  {
563  for ( j= 1; j <= dim; j++ )
564  if ( points[i]->point[j] != (Coord_t) vert[j] ) break;
565  if ( j > dim ) break;
566  }
567 
568  if ( i > num )
569  {
570  addPoint( vert );
571  }
572 
573  pIter( piter );
574  }
575  omFreeSize( (void *) vert, (dim+1) * sizeof(int) );
576 }
577 
578 int pointSet::getExpPos( const poly p )
579 {
580  int * vert;
581  int i,j;
582 
583  // hier unschoen...
584  vert= (int *)omAlloc( (dim+1) * sizeof(int) );
585 
586  p_GetExpV( p, vert, currRing );
587  for ( i= 1; i <= num; i++ )
588  {
589  for ( j= 1; j <= dim; j++ )
590  if ( points[i]->point[j] != (Coord_t) vert[j] ) break;
591  if ( j > dim ) break;
592  }
593  omFreeSize( (void *) vert, (dim+1) * sizeof(int) );
594 
595  if ( i > num ) return 0;
596  else return i;
597 }
598 
599 void pointSet::getRowMP( const int indx, int * vert )
600 {
601  assume( indx > 0 && indx <= num && points[indx]->rcPnt );
602  int i;
603 
604  vert[0]= 0;
605  for ( i= 1; i <= dim; i++ )
606  vert[i]= (int)(points[indx]->point[i] - points[indx]->rcPnt->point[i]);
607 }
608 
609 inline bool pointSet::smaller( int a, int b )
610 {
611  int i;
612 
613  for ( i= 1; i <= dim; i++ )
614  {
615  if ( points[a]->point[i] > points[b]->point[i] )
616  {
617  return false;
618  }
619  if ( points[a]->point[i] < points[b]->point[i] )
620  {
621  return true;
622  }
623  }
624 
625  return false; // they are equal
626 }
627 
628 inline bool pointSet::larger( int a, int b )
629 {
630  int i;
631 
632  for ( i= 1; i <= dim; i++ )
633  {
634  if ( points[a]->point[i] < points[b]->point[i] )
635  {
636  return false;
637  }
638  if ( points[a]->point[i] > points[b]->point[i] )
639  {
640  return true;
641  }
642  }
643 
644  return false; // they are equal
645 }
646 
648 {
649  int i;
650  bool found= true;
651  onePointP tmp;
652 
653  while ( found )
654  {
655  found= false;
656  for ( i= 1; i < num; i++ )
657  {
658  if ( larger( i, i+1 ) )
659  {
660  tmp= points[i];
661  points[i]= points[i+1];
662  points[i+1]= tmp;
663 
664  found= true;
665  }
666  }
667  }
668 }
669 
670 void pointSet::lift( int l[] )
671 {
672  bool outerL= true;
673  int i, j;
674  int sum;
675 
676  dim++;
677 
678  if ( l==NULL )
679  {
680  outerL= false;
681  l= (int *)omAlloc( (dim+1) * sizeof(int) ); // [1..dim-1]
682 
683  for(i = 1; i < dim; i++)
684  {
685  l[i]= 1 + siRand() % LIFT_COOR;
686  }
687  }
688  for ( j=1; j <= num; j++ )
689  {
690  sum= 0;
691  for ( i=1; i < dim; i++ )
692  {
693  sum += (int)points[j]->point[i] * l[i];
694  }
695  points[j]->point[dim]= sum;
696  }
697 
698 #ifdef mprDEBUG_ALL
699  PrintS(" lift vector: ");
700  for ( j=1; j < dim; j++ ) Print(" %d ",l[j] );
701  PrintLn();
702 #ifdef mprDEBUG_ALL
703  PrintS(" lifted points: \n");
704  for ( j=1; j <= num; j++ )
705  {
706  Print("%d: <",j);print_exp(points[j],dim);PrintS(">\n");
707  }
708  PrintLn();
709 #endif
710 #endif
711 
712  lifted= true;
713 
714  if ( !outerL ) omFreeSize( (void *) l, (dim+1) * sizeof(int) );
715 }
716 //<-
717 
718 //-> global functions
719 // Returns the monom at pos i in poly p
720 poly monomAt( poly p, int i )
721 {
722  assume( i > 0 );
723  poly iter= p;
724  for ( int j= 1; (j < i) && (iter!=NULL); j++ ) pIter(iter);
725  return iter;
726 }
727 //<-
728 
729 //-> convexHull::*
730 bool convexHull::inHull(poly p, poly pointPoly, int m, int site)
731 {
732  int i, j, col;
733 
734  pLP->m = n+1;
735  pLP->n = m; // this includes col of cts
736 
737  pLP->LiPM[1][1] = +0.0;
738  pLP->LiPM[1][2] = +1.0; // optimize (arbitrary) var
739  pLP->LiPM[2][1] = +1.0;
740  pLP->LiPM[2][2] = -1.0; // lambda vars sum up to 1
741 
742  for ( j=3; j <= pLP->n; j++)
743  {
744  pLP->LiPM[1][j] = +0.0;
745  pLP->LiPM[2][j] = -1.0;
746  }
747 
748  for( i= 1; i <= n; i++) { // each row constraints one coor
749  pLP->LiPM[i+2][1] = (mprfloat)pGetExp(pointPoly,i);
750  col = 2;
751  for( j= 1; j <= m; j++ )
752  {
753  if( j != site )
754  {
755  pLP->LiPM[i+2][col] = -(mprfloat)pGetExp( monomAt(p,j), i );
756  col++;
757  }
758  }
759  }
760 
761 #ifdef mprDEBUG_ALL
762  PrintS("Matrix of Linear Programming\n");
763  print_mat( pLP->LiPM, pLP->m+1,pLP->n);
764 #endif
765 
766  pLP->m3= pLP->m;
767 
768  pLP->compute();
769 
770  return (pLP->icase == 0);
771 }
772 
773 // mprSTICKYPROT:
774 // ST_SPARSE_VADD: new vertex of convex hull added
775 // ST_SPARSE_VREJ: point rejected (-> inside hull)
777 {
778  int i, j, k;
779  int m; // Anzahl der Exponentvektoren im i-ten Polynom (gls->m)[i] des Ideals gls
780  int idelem= IDELEMS(gls);
781  int * vert;
782 
783  n= (currRing->N);
784  vert= (int *)omAlloc( (idelem+1) * sizeof(int) );
785 
786  Q = (pointSet **)omAlloc( idelem * sizeof(pointSet*) ); // support hulls
787  for ( i= 0; i < idelem; i++ )
788  Q[i] = new pointSet( (currRing->N), i+1, pLength((gls->m)[i])+1 );
789 
790  for( i= 0; i < idelem; i++ )
791  {
792  k=1;
793  m = pLength( (gls->m)[i] );
794 
795  poly p= (gls->m)[i];
796  for( j= 1; j <= m; j++) { // für jeden Exponentvektor
797  if( !inHull( (gls->m)[i], p, m, j ) )
798  {
799  p_GetExpV( p, vert, currRing );
800  Q[i]->addPoint( vert );
801  k++;
803  }
804  else
805  {
807  }
808  pIter( p );
809  } // j
810  mprSTICKYPROT("\n");
811  } // i
812 
813  omFreeSize( (void *) vert, (idelem+1) * sizeof(int) );
814 
815 #ifdef mprDEBUG_PROT
816  PrintLn();
817  for( i= 0; i < idelem; i++ )
818  {
819  Print(" \\Conv(Qi[%d]): #%d\n", i,Q[i]->num );
820  for ( j=1; j <= Q[i]->num; j++ )
821  {
822  Print("%d: <",j);print_exp( (*Q[i])[j] , (currRing->N) );PrintS(">\n");
823  }
824  PrintLn();
825  }
826 #endif
827 
828  return Q;
829 }
830 
831 // mprSTICKYPROT:
832 // ST_SPARSE_VADD: new vertex of convex hull added
833 // ST_SPARSE_VREJ: point rejected (-> inside hull)
834 ideal convexHull::newtonPolytopesI( const ideal gls )
835 {
836  int i, j;
837  int m; // Anzahl der Exponentvektoren im i-ten Polynom (gls->m)[i] des Ideals gls
838  int idelem= IDELEMS(gls);
839  ideal id;
840  poly p,pid;
841  int * vert;
842 
843  n= (currRing->N);
844  vert= (int *)omAlloc( (idelem+1) * sizeof(int) );
845  id= idInit( idelem, 1 );
846 
847  for( i= 0; i < idelem; i++ )
848  {
849  m = pLength( (gls->m)[i] );
850 
851  p= (gls->m)[i];
852  for( j= 1; j <= m; j++) { // für jeden Exponentvektor
853  if( !inHull( (gls->m)[i], p, m, j ) )
854  {
855  if ( (id->m)[i] == NULL )
856  {
857  (id->m)[i]= pHead(p);
858  pid=(id->m)[i];
859  }
860  else
861  {
862  pNext(pid)= pHead(p);
863  pIter(pid);
864  pNext(pid)= NULL;
865  }
867  }
868  else
869  {
871  }
872  pIter( p );
873  } // j
874  mprSTICKYPROT("\n");
875  } // i
876 
877  omFreeSize( (void *) vert, (idelem+1) * sizeof(int) );
878 
879 #ifdef mprDEBUG_PROT
880  PrintLn();
881  for( i= 0; i < idelem; i++ )
882  {
883  }
884 #endif
885 
886  return id;
887 }
888 //<-
889 
890 //-> mayanPyramidAlg::*
892 {
893  int i;
894 
895  Qi= _q_i;
896  shift= _shift;
897 
898  E= new pointSet( Qi[0]->dim ); // E has same dim as Qi[...]
899 
900  for ( i= 0; i < MAXVARS+2; i++ ) acoords[i]= 0;
901 
902  runMayanPyramid(0);
903 
904  mprSTICKYPROT("\n");
905 
906  return E;
907 }
908 
910 {
911  int i, ii, j, k, col, r;
912  int numverts, cols;
913 
914  numverts = 0;
915  for( i=0; i<=n; i++)
916  {
917  numverts += Qi[i]->num;
918  }
919  cols = numverts + 2;
920 
921  //if( dim < 1 || dim > n )
922  // WerrorS("mayanPyramidAlg::vDistance: Known coords dim off range");
923 
924  pLP->LiPM[1][1] = 0.0;
925  pLP->LiPM[1][2] = 1.0; // maximize
926  for( j=3; j<=cols; j++) pLP->LiPM[1][j] = 0.0;
927 
928  for( i=0; i <= n; i++ )
929  {
930  pLP->LiPM[i+2][1] = 1.0;
931  pLP->LiPM[i+2][2] = 0.0;
932  }
933  for( i=1; i<=dim; i++)
934  {
935  pLP->LiPM[n+2+i][1] = (mprfloat)(acoords_a[i-1]);
936  pLP->LiPM[n+2+i][2] = -shift[i];
937  }
938 
939  ii = -1;
940  col = 2;
941  for ( i= 0; i <= n; i++ )
942  {
943  ii++;
944  for( k= 1; k <= Qi[ii]->num; k++ )
945  {
946  col++;
947  for ( r= 0; r <= n; r++ )
948  {
949  if ( r == i ) pLP->LiPM[r+2][col] = -1.0;
950  else pLP->LiPM[r+2][col] = 0.0;
951  }
952  for( r= 1; r <= dim; r++ )
953  pLP->LiPM[r+n+2][col] = -(mprfloat)((*Qi[ii])[k]->point[r]);
954  }
955  }
956 
957  if( col != cols)
958  Werror("mayanPyramidAlg::vDistance:"
959  "setting up matrix for udist: col %d != cols %d",col,cols);
960 
961  pLP->m = n+dim+1;
962  pLP->m3= pLP->m;
963  pLP->n=cols-1;
964 
965 #ifdef mprDEBUG_ALL
966  Print("vDistance LP, known koords dim=%d, constr %d, cols %d, acoords= ",
967  dim,pLP->m,cols);
968  for( i= 0; i < dim; i++ )
969  Print(" %d",acoords_a[i]);
970  PrintLn();
971  print_mat( pLP->LiPM, pLP->m+1, cols);
972 #endif
973 
974  pLP->compute();
975 
976 #ifdef mprDEBUG_ALL
977  PrintS("LP returns matrix\n");
978  print_bmat( pLP->LiPM, pLP->m+1, cols+1-pLP->m, cols, pLP->iposv);
979 #endif
980 
981  if( pLP->icase != 0 )
982  { // check for errors
983  WerrorS("mayanPyramidAlg::vDistance:");
984  if( pLP->icase == 1 )
985  WerrorS(" Unbounded v-distance: probably 1st v-coor=0");
986  else if( pLP->icase == -1 )
987  WerrorS(" Infeasible v-distance");
988  else
989  WerrorS(" Unknown error");
990  return -1.0;
991  }
992 
993  return pLP->LiPM[1][1];
994 }
995 
997 {
998  int i, j, k, cols, cons;
999  int la_cons_row;
1000 
1001  cons = n+dim+2;
1002 
1003  // first, compute minimum
1004  //
1005 
1006  // common part of the matrix
1007  pLP->LiPM[1][1] = 0.0;
1008  for( i=2; i<=n+2; i++)
1009  {
1010  pLP->LiPM[i][1] = 1.0; // 1st col
1011  pLP->LiPM[i][2] = 0.0; // 2nd col
1012  }
1013 
1014  la_cons_row = 1;
1015  cols = 2;
1016  for( i=0; i<=n; i++)
1017  {
1018  la_cons_row++;
1019  for( j=1; j<= Qi[i]->num; j++)
1020  {
1021  cols++;
1022  pLP->LiPM[1][cols] = 0.0; // set 1st row 0
1023  for( k=2; k<=n+2; k++)
1024  { // lambdas sum up to 1
1025  if( k != la_cons_row) pLP->LiPM[k][cols] = 0.0;
1026  else pLP->LiPM[k][cols] = -1.0;
1027  }
1028  for( k=1; k<=n; k++)
1029  pLP->LiPM[k+n+2][cols] = -(mprfloat)((*Qi[i])[j]->point[k]);
1030  } // j
1031  } // i
1032 
1033  for( i= 0; i < dim; i++ )
1034  { // fixed coords
1035  pLP->LiPM[i+n+3][1] = acoords[i];
1036  pLP->LiPM[i+n+3][2] = 0.0;
1037  }
1038  pLP->LiPM[dim+n+3][1] = 0.0;
1039 
1040 
1041  pLP->LiPM[1][2] = -1.0; // minimize
1042  pLP->LiPM[dim+n+3][2] = 1.0;
1043 
1044 #ifdef mprDEBUG_ALL
1045  Print("\nThats the matrix for minR, dim= %d, acoords= ",dim);
1046  for( i= 0; i < dim; i++ )
1047  Print(" %d",acoords[i]);
1048  PrintLn();
1049  print_mat( pLP->LiPM, cons+1, cols);
1050 #endif
1051 
1052  // simplx finds MIN for obj.fnc, puts it in [1,1]
1053  pLP->m= cons;
1054  pLP->n= cols-1;
1055  pLP->m3= cons;
1056 
1057  pLP->compute();
1058 
1059  if ( pLP->icase != 0 )
1060  { // check for errors
1061  if( pLP->icase < 0)
1062  WerrorS(" mn_mx_MinkowskiSum: LinearProgram: minR: infeasible");
1063  else if( pLP->icase > 0)
1064  WerrorS(" mn_mx_MinkowskiSum: LinearProgram: minR: unbounded");
1065  }
1066 
1067  *minR = (Coord_t)( -pLP->LiPM[1][1] + 1.0 - SIMPLEX_EPS );
1068 
1069  // now compute maximum
1070  //
1071 
1072  // common part of the matrix again
1073  pLP->LiPM[1][1] = 0.0;
1074  for( i=2; i<=n+2; i++)
1075  {
1076  pLP->LiPM[i][1] = 1.0;
1077  pLP->LiPM[i][2] = 0.0;
1078  }
1079  la_cons_row = 1;
1080  cols = 2;
1081  for( i=0; i<=n; i++)
1082  {
1083  la_cons_row++;
1084  for( j=1; j<=Qi[i]->num; j++)
1085  {
1086  cols++;
1087  pLP->LiPM[1][cols] = 0.0;
1088  for( k=2; k<=n+2; k++)
1089  {
1090  if( k != la_cons_row) pLP->LiPM[k][cols] = 0.0;
1091  else pLP->LiPM[k][cols] = -1.0;
1092  }
1093  for( k=1; k<=n; k++)
1094  pLP->LiPM[k+n+2][cols] = -(mprfloat)((*Qi[i])[j]->point[k]);
1095  } // j
1096  } // i
1097 
1098  for( i= 0; i < dim; i++ )
1099  { // fixed coords
1100  pLP->LiPM[i+n+3][1] = acoords[i];
1101  pLP->LiPM[i+n+3][2] = 0.0;
1102  }
1103  pLP->LiPM[dim+n+3][1] = 0.0;
1104 
1105  pLP->LiPM[1][2] = 1.0; // maximize
1106  pLP->LiPM[dim+n+3][2] = 1.0; // var = sum of pnt coords
1107 
1108 #ifdef mprDEBUG_ALL
1109  Print("\nThats the matrix for maxR, dim= %d\n",dim);
1110  print_mat( pLP->LiPM, cons+1, cols);
1111 #endif
1112 
1113  pLP->m= cons;
1114  pLP->n= cols-1;
1115  pLP->m3= cons;
1116 
1117  // simplx finds MAX for obj.fnc, puts it in [1,1]
1118  pLP->compute();
1119 
1120  if ( pLP->icase != 0 )
1121  {
1122  if( pLP->icase < 0)
1123  WerrorS(" mn_mx_MinkowskiSum: LinearProgram: maxR: infeasible");
1124  else if( pLP->icase > 0)
1125  WerrorS(" mn_mx_MinkowskiSum: LinearProgram: maxR: unbounded");
1126  }
1127 
1128  *maxR = (Coord_t)( pLP->LiPM[1][1] + SIMPLEX_EPS );
1129 
1130 #ifdef mprDEBUG_ALL
1131  Print(" Range for dim=%d: [%d,%d]\n", dim, *minR, *maxR);
1132 #endif
1133 }
1134 
1135 // mprSTICKYPROT:
1136 // ST_SPARSE_VREJ: rejected point
1137 // ST_SPARSE_VADD: added point to set
1139 {
1140  mprfloat dist;
1141 
1142  // determine v-distance of point pt
1143  dist= vDistance( &(acoords[0]), n );
1144 
1145  // store only points with v-distance > minVdist
1146  if( dist <= MINVDIST + SIMPLEX_EPS )
1147  {
1149  return false;
1150  }
1151 
1152  E->addPoint( &(acoords[0]) );
1154 
1155  return true;
1156 }
1157 
1158 // mprSTICKYPROT:
1159 // ST_SPARSE_MREC1: recurse
1160 // ST_SPARSE_MREC2: recurse with extra points
1161 // ST_SPARSE_MPEND: end
1163 {
1164  Coord_t minR, maxR;
1165  mprfloat dist;
1166 
1167  // step 3
1168  mn_mx_MinkowskiSum( dim, &minR, &maxR );
1169 
1170 #ifdef mprDEBUG_ALL
1171  int i;
1172  for( i=0; i <= dim; i++) Print("acoords[%d]=%d ",i,(int)acoords[i]);
1173  Print(":: [%d,%d]\n", minR, maxR);
1174 #endif
1175 
1176  // step 5 -> terminate
1177  if( dim == n-1 )
1178  {
1179  int lastKilled = 0;
1180  // insert points
1181  acoords[dim] = minR;
1182  while( acoords[dim] <= maxR )
1183  {
1184  if( !storeMinkowskiSumPoint() )
1185  lastKilled++;
1186  acoords[dim]++;
1187  }
1189  return;
1190  }
1191 
1192  // step 4 -> recurse at step 3
1193  acoords[dim] = minR;
1194  while ( acoords[dim] <= maxR )
1195  {
1196  if ( (acoords[dim] > minR) && (acoords[dim] <= maxR) )
1197  { // acoords[dim] >= minR ??
1199  runMayanPyramid( dim + 1 ); // recurse with higer dimension
1200  }
1201  else
1202  {
1203  // get v-distance of pt
1204  dist= vDistance( &(acoords[0]), dim + 1 );// dim+1 == known coordinates
1205 
1206  if( dist >= SIMPLEX_EPS )
1207  {
1209  runMayanPyramid( dim + 1 ); // recurse with higer dimension
1210  }
1211  }
1212  acoords[dim]++;
1213  } // while
1214 }
1215 //<-
1216 
1217 //-> resMatrixSparse::*
1218 bool resMatrixSparse::remapXiToPoint( const int indx, pointSet **pQ, int *set, int *pnt )
1219 {
1220  int i,nn= (currRing->N);
1221  int loffset= 0;
1222  for ( i= 0; i <= nn; i++ )
1223  {
1224  if ( (loffset < indx) && (indx <= pQ[i]->num + loffset) )
1225  {
1226  *set= i;
1227  *pnt= indx-loffset;
1228  return true;
1229  }
1230  else loffset+= pQ[i]->num;
1231  }
1232  return false;
1233 }
1234 
1235 // mprSTICKYPROT
1236 // ST_SPARSE_RC: point added
1237 int resMatrixSparse::RC( pointSet **pQ, pointSet *E, int vert, mprfloat shift[] )
1238 {
1239  int i, j, k,c ;
1240  int size;
1241  bool found= true;
1242  mprfloat cd;
1243  int onum;
1244  int bucket[MAXVARS+2];
1245  setID *optSum;
1246 
1247  LP->n = 1;
1248  LP->m = n + n + 1; // number of constrains
1249 
1250  // fill in LP matrix
1251  for ( i= 0; i <= n; i++ )
1252  {
1253  size= pQ[i]->num;
1254  for ( k= 1; k <= size; k++ )
1255  {
1256  LP->n++;
1257 
1258  // objective funtion, minimize
1259  LP->LiPM[1][LP->n] = - ( (mprfloat) (*pQ[i])[k]->point[pQ[i]->dim] / SCALEDOWN );
1260 
1261  // lambdas sum up to 1
1262  for ( j = 0; j <= n; j++ )
1263  {
1264  if ( i==j )
1265  LP->LiPM[j+2][LP->n] = -1.0;
1266  else
1267  LP->LiPM[j+2][LP->n] = 0.0;
1268  }
1269 
1270  // the points
1271  for ( j = 1; j <= n; j++ )
1272  {
1273  LP->LiPM[j+n+2][LP->n] = - ( (mprfloat) (*pQ[i])[k]->point[j] );
1274  }
1275  }
1276  }
1277 
1278  for ( j = 0; j <= n; j++ ) LP->LiPM[j+2][1] = 1.0;
1279  for ( j= 1; j <= n; j++ )
1280  {
1281  LP->LiPM[j+n+2][1]= (mprfloat)(*E)[vert]->point[j] - shift[j];
1282  }
1283  LP->n--;
1284 
1285  LP->LiPM[1][1] = 0.0;
1286 
1287 #ifdef mprDEBUG_ALL
1288  PrintLn();
1289  Print(" n= %d, LP->m=M= %d, LP->n=N= %d\n",n,LP->m,LP->n);
1290  print_mat(LP->LiPM, LP->m+1, LP->n+1);
1291 #endif
1292 
1293  LP->m3= LP->m;
1294 
1295  LP->compute();
1296 
1297  if ( LP->icase < 0 )
1298  {
1299  // infeasibility: the point does not lie in a cell -> remove it
1300  return -1;
1301  }
1302 
1303  // store result
1304  (*E)[vert]->point[E->dim]= (int)(-LP->LiPM[1][1] * SCALEDOWN);
1305 
1306 #ifdef mprDEBUG_ALL
1307  Print(" simplx returned %d, Objective value = %f\n", LP->icase, LP->LiPM[1][1]);
1308  //print_bmat(LP->LiPM, NumCons + 1, LP->n+1-NumCons, LP->n+1, LP->iposv); // ( rows= M+1, cols= N+1-m3 )
1309  //print_mat(LP->LiPM, NumCons+1, LP->n);
1310 #endif
1311 
1312 #if 1
1313  // sort LP results
1314  while (found)
1315  {
1316  found=false;
1317  for ( i= 1; i < LP->m; i++ )
1318  {
1319  if ( LP->iposv[i] > LP->iposv[i+1] )
1320  {
1321 
1322  c= LP->iposv[i];
1323  LP->iposv[i]=LP->iposv[i+1];
1324  LP->iposv[i+1]=c;
1325 
1326  cd=LP->LiPM[i+1][1];
1327  LP->LiPM[i+1][1]=LP->LiPM[i+2][1];
1328  LP->LiPM[i+2][1]=cd;
1329 
1330  found= true;
1331  }
1332  }
1333  }
1334 #endif
1335 
1336 #ifdef mprDEBUG_ALL
1337  print_bmat(LP->LiPM, LP->m + 1, LP->n+1-LP->m, LP->n+1, LP->iposv);
1338  PrintS(" now split into sets\n");
1339 #endif
1340 
1341 
1342  // init bucket
1343  for ( i= 0; i <= E->dim; i++ ) bucket[i]= 0;
1344  // remap results of LP to sets Qi
1345  c=0;
1346  optSum= (setID*)omAlloc( (LP->m) * sizeof(struct setID) );
1347  for ( i= 0; i < LP->m; i++ )
1348  {
1349  //Print("% .15f\n",LP->LiPM[i+2][1]);
1350  if ( LP->LiPM[i+2][1] > 1e-12 )
1351  {
1352  if ( !remapXiToPoint( LP->iposv[i+1], pQ, &(optSum[c].set), &(optSum[c].pnt) ) )
1353  {
1354  Werror(" resMatrixSparse::RC: Found bad solution in LP: %d!",LP->iposv[i+1]);
1355  WerrorS(" resMatrixSparse::RC: remapXiToPoint failed!");
1356  return -1;
1357  }
1358  bucket[optSum[c].set]++;
1359  c++;
1360  }
1361  }
1362 
1363  onum= c;
1364  // find last min in bucket[]: maximum i such that Fi is a point
1365  c= 0;
1366  for ( i= 1; i < E->dim; i++ )
1367  {
1368  if ( bucket[c] >= bucket[i] )
1369  {
1370  c= i;
1371  }
1372  }
1373  // find matching point set
1374  for ( i= onum - 1; i >= 0; i-- )
1375  {
1376  if ( optSum[i].set == c )
1377  break;
1378  }
1379  // store
1380  (*E)[vert]->rc.set= c;
1381  (*E)[vert]->rc.pnt= optSum[i].pnt;
1382  (*E)[vert]->rcPnt= (*pQ[c])[optSum[i].pnt];
1383  // count
1384  if ( (*E)[vert]->rc.set == linPolyS ) numSet0++;
1385 
1386 #ifdef mprDEBUG_PROT
1387  Print("\n Point E[%d] was <",vert);print_exp((*E)[vert],E->dim-1);Print(">, bucket={");
1388  for ( j= 0; j < E->dim; j++ )
1389  {
1390  Print(" %d",bucket[j]);
1391  }
1392  PrintS(" }\n optimal Sum: Qi ");
1393  for ( j= 0; j < LP->m; j++ )
1394  {
1395  Print(" [ %d, %d ]",optSum[j].set,optSum[j].pnt);
1396  }
1397  Print(" -> i= %d, j = %d\n",(*E)[vert]->rc.set,optSum[i].pnt);
1398 #endif
1399 
1400  // clean up
1401  omFreeSize( (void *) optSum, (LP->m) * sizeof(struct setID) );
1402 
1404 
1405  return (int)(-LP->LiPM[1][1] * SCALEDOWN);
1406 }
1407 
1408 // create coeff matrix
1410 {
1411  // sparse matrix
1412  // uRPos[i][1]: row of matrix
1413  // uRPos[i][idelem+1]: col of u(0)
1414  // uRPos[i][2..idelem]: col of u(1) .. u(n)
1415  // i= 1 .. numSet0
1416  int i,epos;
1417  int rp,cp;
1418  poly rowp,epp;
1419  poly iterp;
1420  int *epp_mon, *eexp;
1421 
1422  epp_mon= (int *)omAlloc( (n+2) * sizeof(int) );
1423  eexp= (int *)omAlloc0(((currRing->N)+1)*sizeof(int));
1424 
1425  totDeg= numSet0;
1426 
1427  mprSTICKYPROT2(" size of matrix: %d\n", E->num);
1428  mprSTICKYPROT2(" resultant deg: %d\n", numSet0);
1429 
1430  uRPos= new intvec( numSet0, pLength((gls->m)[0])+1, 0 );
1431 
1432  // sparse Matrix represented as a module where
1433  // each poly is column vector ( pSetComp(p,k) gives the row )
1434  rmat= idInit( E->num, E->num ); // cols, rank= number of rows
1435  msize= E->num;
1436 
1437  rp= 1;
1438  rowp= NULL;
1439  epp= pOne();
1440  for ( i= 1; i <= E->num; i++ )
1441  { // for every row
1442  E->getRowMP( i, epp_mon ); // compute (p-a[ij]), (i,j) = RC(p)
1443  pSetExpV( epp, epp_mon );
1444 
1445  //
1446  rowp= ppMult_qq( epp, (gls->m)[(*E)[i]->rc.set] ); // x^(p-a[ij]) * f(i)
1447 
1448  cp= 2;
1449  // get column for every monomial in rowp and store it
1450  iterp= rowp;
1451  while ( iterp!=NULL )
1452  {
1453  epos= E->getExpPos( iterp );
1454  if ( epos == 0 )
1455  {
1456  // this can happen, if the shift vektor or the lift funktions
1457  // are not generically chosen.
1458  Werror("resMatrixSparse::createMatrix: Found exponent not in E, id %d, set [%d, %d]!",
1459  i,(*E)[i]->rc.set,(*E)[i]->rc.pnt);
1460  return i;
1461  }
1462  pSetExpV(iterp,eexp);
1463  pSetComp(iterp, epos );
1464  pSetm(iterp);
1465  if ( (*E)[i]->rc.set == linPolyS )
1466  { // store coeff positions
1467  IMATELEM(*uRPos,rp,cp)= epos;
1468  cp++;
1469  }
1470  pIter( iterp );
1471  } // while
1472  if ( (*E)[i]->rc.set == linPolyS )
1473  { // store row
1474  IMATELEM(*uRPos,rp,1)= i-1;
1475  rp++;
1476  }
1477  (rmat->m)[i-1]= rowp;
1478  } // for
1479 
1480  pDelete( &epp );
1481  omFreeSize( (void *) epp_mon, (n+2) * sizeof(int) );
1482  omFreeSize( (void *) eexp, ((currRing->N)+1)*sizeof(int));
1483 
1484 #ifdef mprDEBUG_ALL
1485  if ( E->num <= 40 )
1486  {
1487  matrix mout= idModule2Matrix( idCopy(rmat) );
1488  print_matrix(mout);
1489  }
1490  for ( i= 1; i <= numSet0; i++ )
1491  {
1492  Print(" row %d contains coeffs of f_%d\n",IMATELEM(*uRPos,i,1),linPolyS);
1493  }
1494  PrintS(" Sparse Matrix done\n");
1495 #endif
1496 
1497  return E->num;
1498 }
1499 
1500 // find a sufficiently generic and small vector
1501 void resMatrixSparse::randomVector( const int dim, mprfloat shift[] )
1502 {
1503  int i,j;
1504  i= 1;
1505 
1506  while ( i <= dim )
1507  {
1508  shift[i]= (mprfloat) (RVMULT*(siRand()%MAXRVVAL)/(mprfloat)MAXRVVAL);
1509  i++;
1510  for ( j= 1; j < i-1; j++ )
1511  {
1512  if ( (shift[j] < shift[i-1] + SIMPLEX_EPS) && (shift[j] > shift[i-1] - SIMPLEX_EPS) )
1513  {
1514  i--;
1515  break;
1516  }
1517  }
1518  }
1519 }
1520 
1522 {
1523  pointSet *vs;
1524  onePoint vert;
1525  int j,k,l;
1526 
1527  vert.point=(Coord_t*)omAlloc( ((currRing->N)+2) * sizeof(Coord_t) );
1528 
1529  vs= new pointSet( dim );
1530 
1531  for ( j= 1; j <= Q1->num; j++ )
1532  {
1533  for ( k= 1; k <= Q2->num; k++ )
1534  {
1535  for ( l= 1; l <= dim; l++ )
1536  {
1537  vert.point[l]= (*Q1)[j]->point[l] + (*Q2)[k]->point[l];
1538  }
1539  vs->mergeWithExp( &vert );
1540  //vs->addPoint( &vert );
1541  }
1542  }
1543 
1544  omFreeSize( (void *) vert.point, ((currRing->N)+2) * sizeof(Coord_t) );
1545 
1546  return vs;
1547 }
1548 
1550 {
1551  pointSet *vs,*vs_old;
1552  int j;
1553 
1554  vs= new pointSet( dim );
1555 
1556  for ( j= 1; j <= pQ[0]->num; j++ ) vs->addPoint( (*pQ[0])[j] );
1557 
1558  for ( j= 1; j < numq; j++ )
1559  {
1560  vs_old= vs;
1561  vs= minkSumTwo( vs_old, pQ[j], dim );
1562 
1563  delete vs_old;
1564  }
1565 
1566  return vs;
1567 }
1568 
1569 //----------------------------------------------------------------------------------------
1570 
1571 resMatrixSparse::resMatrixSparse( const ideal _gls, const int special )
1572  : resMatrixBase(), gls( _gls )
1573 {
1574  pointSet **Qi; // vertices sets of Conv(Supp(f_i)), i=0..idelem
1575  pointSet *E; // all integer lattice points of the minkowski sum of Q0...Qn
1576  int i,k;
1577  int pnt;
1578  int totverts; // total number of exponent vectors in ideal gls
1579  mprfloat shift[MAXVARS+2]; // shiftvector delta, index [1..dim]
1580 
1581  if ( (currRing->N) > MAXVARS )
1582  {
1583  WerrorS("resMatrixSparse::resMatrixSparse: Too many variables!");
1584  return;
1585  }
1586 
1587  rmat= NULL;
1588  numSet0= 0;
1589 
1590  if ( special == SNONE ) linPolyS= 0;
1591  else linPolyS= special;
1592 
1594 
1595  n= (currRing->N);
1596  idelem= IDELEMS(gls); // should be n+1
1597 
1598  // prepare matrix LP->LiPM for Linear Programming
1599  totverts = 0;
1600  for( i=0; i < idelem; i++) totverts += pLength( (gls->m)[i] );
1601 
1602  LP = new simplex( idelem+totverts*2+5, totverts+5 ); // rows, cols
1603 
1604  // get shift vector
1605 #ifdef mprTEST
1606  shift[0]=0.005; shift[1]=0.003; shift[2]=0.008; shift[3]=0.005; shift[4]=0.002;
1607  shift[5]=0.1; shift[6]=0.3; shift[7]=0.2; shift[8]=0.4; shift[9]=0.2;
1608 #else
1609  randomVector( idelem, shift );
1610 #endif
1611 #ifdef mprDEBUG_PROT
1612  PrintS(" shift vector: ");
1613  for ( i= 1; i <= idelem; i++ ) Print(" %.12f ",(double)shift[i]);
1614  PrintLn();
1615 #endif
1616 
1617  // evaluate convex hull for supports of gls
1618  convexHull chnp( LP );
1619  Qi= chnp.newtonPolytopesP( gls );
1620 
1621 #ifdef mprMINKSUM
1622  E= minkSumAll( Qi, n+1, n);
1623 #else
1624  // get inner points
1625  mayanPyramidAlg mpa( LP );
1626  E= mpa.getInnerPoints( Qi, shift );
1627 #endif
1628 
1629 #ifdef mprDEBUG_PROT
1630 #ifdef mprMINKSUM
1631  PrintS("(MinkSum)");
1632 #endif
1633  PrintS("\n E = (Q_0 + ... + Q_n) \\cap \\N :\n");
1634  for ( pnt= 1; pnt <= E->num; pnt++ )
1635  {
1636  Print("%d: <",pnt);print_exp( (*E)[pnt], E->dim );PrintS(">\n");
1637  }
1638  PrintLn();
1639 #endif
1640 
1641 #ifdef mprTEST
1642  int lift[5][5];
1643  lift[0][1]=3; lift[0][2]=4; lift[0][3]=8; lift[0][4]=2;
1644  lift[1][1]=6; lift[1][2]=1; lift[1][3]=7; lift[1][4]=4;
1645  lift[2][1]=2; lift[2][2]=5; lift[2][3]=9; lift[2][4]=6;
1646  lift[3][1]=2; lift[3][2]=1; lift[3][3]=9; lift[3][4]=5;
1647  lift[4][1]=3; lift[4][2]=7; lift[4][3]=1; lift[4][4]=5;
1648  // now lift everything
1649  for ( i= 0; i <= n; i++ ) Qi[i]->lift( lift[i] );
1650 #else
1651  // now lift everything
1652  for ( i= 0; i <= n; i++ ) Qi[i]->lift();
1653 #endif
1654  E->dim++;
1655 
1656  // run Row Content Function for every point in E
1657  for ( pnt= 1; pnt <= E->num; pnt++ )
1658  {
1659  RC( Qi, E, pnt, shift );
1660  }
1661 
1662  // remove points not in cells
1663  k= E->num;
1664  for ( pnt= k; pnt > 0; pnt-- )
1665  {
1666  if ( (*E)[pnt]->rcPnt == NULL )
1667  {
1668  E->removePoint(pnt);
1670  }
1671  }
1672  mprSTICKYPROT("\n");
1673 
1674 #ifdef mprDEBUG_PROT
1675  PrintS(" points which lie in a cell:\n");
1676  for ( pnt= 1; pnt <= E->num; pnt++ )
1677  {
1678  Print("%d: <",pnt);print_exp( (*E)[pnt], E->dim );PrintS(">\n");
1679  }
1680  PrintLn();
1681 #endif
1682 
1683  // unlift to old dimension, sort
1684  for ( i= 0; i <= n; i++ ) Qi[i]->unlift();
1685  E->unlift();
1686  E->sort();
1687 
1688 #ifdef mprDEBUG_PROT
1689  Print(" points with a[ij] (%d):\n",E->num);
1690  for ( pnt= 1; pnt <= E->num; pnt++ )
1691  {
1692  Print("Punkt p \\in E[%d]: <",pnt);print_exp( (*E)[pnt], E->dim );
1693  Print(">, RC(p) = (i:%d, j:%d), a[i,j] = <",(*E)[pnt]->rc.set,(*E)[pnt]->rc.pnt);
1694  //print_exp( (Qi[(*E)[pnt]->rc.set])[(*E)[pnt]->rc.pnt], E->dim );PrintS("> = <");
1695  print_exp( (*E)[pnt]->rcPnt, E->dim );PrintS(">\n");
1696  }
1697 #endif
1698 
1699  // now create matrix
1700  if (E->num <1)
1701  {
1702  WerrorS("could not handle a degenerate situation: no inner points found");
1703  goto theEnd;
1704  }
1705  if ( createMatrix( E ) != E->num )
1706  {
1707  // this can happen if the shiftvector shift is to large or not generic
1709  WerrorS("resMatrixSparse::resMatrixSparse: Error in resMatrixSparse::createMatrix!");
1710  goto theEnd;
1711  }
1712 
1713  theEnd:
1714  // clean up
1715  for ( i= 0; i < idelem; i++ )
1716  {
1717  delete Qi[i];
1718  }
1719  omFreeSize( (void *) Qi, idelem * sizeof(pointSet*) );
1720 
1721  delete E;
1722 
1723  delete LP;
1724 }
1725 
1726 //----------------------------------------------------------------------------------------
1727 
1729 {
1730  delete uRPos;
1731  idDelete( &rmat );
1732 }
1733 
1735 {
1736  int i,/*j,*/cp;
1737  poly pp,phelp,piter,pgls;
1738 
1739  // copy original sparse res matrix
1740  ideal rmat_out= idCopy(rmat);
1741 
1742  // now fill in coeffs of f0
1743  for ( i= 1; i <= numSet0; i++ )
1744  {
1745 
1746  pgls= (gls->m)[0]; // f0
1747 
1748  // get matrix row and delete it
1749  pp= (rmat_out->m)[IMATELEM(*uRPos,i,1)];
1750  pDelete( &pp );
1751  pp= NULL;
1752  phelp= pp;
1753  piter= NULL;
1754 
1755  // u_1,..,u_k
1756  cp=2;
1757  while ( pNext(pgls)!=NULL )
1758  {
1759  phelp= pOne();
1760  pSetCoeff( phelp, nCopy(pGetCoeff(pgls)) );
1761  pSetComp( phelp, IMATELEM(*uRPos,i,cp) );
1762  pSetmComp( phelp );
1763  if ( piter!=NULL )
1764  {
1765  pNext(piter)= phelp;
1766  piter= phelp;
1767  }
1768  else
1769  {
1770  pp= phelp;
1771  piter= phelp;
1772  }
1773  cp++;
1774  pIter( pgls );
1775  }
1776  // u0, now pgls points to last monom
1777  phelp= pOne();
1778  pSetCoeff( phelp, nCopy(pGetCoeff(pgls)) );
1779  //pSetComp( phelp, IMATELEM(*uRPos,i,idelem+1) );
1780  pSetComp( phelp, IMATELEM(*uRPos,i,pLength((gls->m)[0])+1) );
1781  pSetmComp( phelp );
1782  if (piter!=NULL) pNext(piter)= phelp;
1783  else pp= phelp;
1784  (rmat_out->m)[IMATELEM(*uRPos,i,1)]= pp;
1785  }
1786 
1787  return rmat_out;
1788 }
1789 
1790 // Fills in resMat[][] with evpoint[] and gets determinant
1791 // uRPos[i][1]: row of matrix
1792 // uRPos[i][idelem+1]: col of u(0)
1793 // uRPos[i][2..idelem]: col of u(1) .. u(n)
1794 // i= 1 .. numSet0
1795 number resMatrixSparse::getDetAt( const number* evpoint )
1796 {
1797  int i,cp;
1798  poly pp,phelp,piter;
1799 
1800  mprPROTnl("smCallDet");
1801 
1802  for ( i= 1; i <= numSet0; i++ )
1803  {
1804  pp= (rmat->m)[IMATELEM(*uRPos,i,1)];
1805  pDelete( &pp );
1806  pp= NULL;
1807  phelp= pp;
1808  piter= NULL;
1809  // u_1,..,u_n
1810  for ( cp= 2; cp <= idelem; cp++ )
1811  {
1812  if ( !nIsZero(evpoint[cp-1]) )
1813  {
1814  phelp= pOne();
1815  pSetCoeff( phelp, nCopy(evpoint[cp-1]) );
1816  pSetComp( phelp, IMATELEM(*uRPos,i,cp) );
1817  pSetmComp( phelp );
1818  if ( piter )
1819  {
1820  pNext(piter)= phelp;
1821  piter= phelp;
1822  }
1823  else
1824  {
1825  pp= phelp;
1826  piter= phelp;
1827  }
1828  }
1829  }
1830  // u0
1831  phelp= pOne();
1832  pSetCoeff( phelp, nCopy(evpoint[0]) );
1833  pSetComp( phelp, IMATELEM(*uRPos,i,idelem+1) );
1834  pSetmComp( phelp );
1835  pNext(piter)= phelp;
1836  (rmat->m)[IMATELEM(*uRPos,i,1)]= pp;
1837  }
1838 
1839  mprSTICKYPROT(ST__DET); // 1
1840 
1841  poly pres= sm_CallDet( rmat, currRing );
1842  number numres= nCopy( pGetCoeff( pres ) );
1843  pDelete( &pres );
1844 
1845  mprSTICKYPROT(ST__DET); // 2
1846 
1847  return ( numres );
1848 }
1849 
1850 // Fills in resMat[][] with evpoint[] and gets determinant
1851 // uRPos[i][1]: row of matrix
1852 // uRPos[i][idelem+1]: col of u(0)
1853 // uRPos[i][2..idelem]: col of u(1) .. u(n)
1854 // i= 1 .. numSet0
1855 poly resMatrixSparse::getUDet( const number* evpoint )
1856 {
1857  int i,cp;
1858  poly pp,phelp/*,piter*/;
1859 
1860  mprPROTnl("smCallDet");
1861 
1862  for ( i= 1; i <= numSet0; i++ )
1863  {
1864  pp= (rmat->m)[IMATELEM(*uRPos,i,1)];
1865  pDelete( &pp );
1866  phelp= NULL;
1867  // piter= NULL;
1868  for ( cp= 2; cp <= idelem; cp++ )
1869  { // u1 .. un
1870  if ( !nIsZero(evpoint[cp-1]) )
1871  {
1872  phelp= pOne();
1873  pSetCoeff( phelp, nCopy(evpoint[cp-1]) );
1874  pSetComp( phelp, IMATELEM(*uRPos,i,cp) );
1875  //pSetmComp( phelp );
1876  pSetm( phelp );
1877  //Print("comp %d\n",IMATELEM(*uRPos,i,cp));
1878  #if 0
1879  if ( piter!=NULL )
1880  {
1881  pNext(piter)= phelp;
1882  piter= phelp;
1883  }
1884  else
1885  {
1886  pp= phelp;
1887  piter= phelp;
1888  }
1889  #else
1890  pp=pAdd(pp,phelp);
1891  #endif
1892  }
1893  }
1894  // u0
1895  phelp= pOne();
1896  pSetExp(phelp,1,1);
1897  pSetComp( phelp, IMATELEM(*uRPos,i,idelem+1) );
1898  // Print("comp %d\n",IMATELEM(*uRPos,i,idelem+1));
1899  pSetm( phelp );
1900  #if 0
1901  pNext(piter)= phelp;
1902  #else
1903  pp=pAdd(pp,phelp);
1904  #endif
1905  pTest(pp);
1906  (rmat->m)[IMATELEM(*uRPos,i,1)]= pp;
1907  }
1908 
1909  mprSTICKYPROT(ST__DET); // 1
1910 
1911  poly pres= sm_CallDet( rmat, currRing );
1912 
1913  mprSTICKYPROT(ST__DET); // 2
1914 
1915  return ( pres );
1916 }
1917 //<-
1918 
1919 //-----------------------------------------------------------------------------
1920 //-------------- dense resultant matrix ---------------------------------------
1921 //-----------------------------------------------------------------------------
1922 
1923 //-> dense resultant matrix
1924 //
1925 struct resVector;
1926 
1927 /* dense resultant matrix */
1928 class resMatrixDense : virtual public resMatrixBase
1929 {
1930 public:
1931  /**
1932  * _gls: system of multivariate polynoms
1933  * special: -1 -> resMatrixDense is a symbolic matrix
1934  * 0,1, ... -> resMatrixDense ist eine u-Resultante, wobei special das
1935  * lineare u-Polynom angibt
1936  */
1937  resMatrixDense( const ideal _gls, const int special = SNONE );
1938  ~resMatrixDense();
1939 
1940  /** column vector of matrix, index von 0 ... numVectors-1 */
1941  resVector *getMVector( const int i );
1942 
1943  /** Returns the matrix M in an usable presentation */
1944  ideal getMatrix();
1945 
1946  /** Returns the submatrix M' of M in an usable presentation */
1947  ideal getSubMatrix();
1948 
1949  /** Evaluate the determinant of the matrix M at the point evpoint
1950  * where the ui's are replaced by the components of evpoint.
1951  * Uses singclap_det from factory.
1952  */
1953  number getDetAt( const number* evpoint );
1954 
1955  /** Evaluates the determinant of the submatrix M'.
1956  * Since the matrix is numerically, no evaluation point is needed.
1957  * Uses singclap_det from factory.
1958  */
1959  number getSubDet();
1960 
1961 private:
1962  /** deactivated copy constructor */
1964 
1965  /** Generate the "matrix" M. Each column is presented by a resVector
1966  * holding all entries for this column.
1967  */
1968  void generateBaseData();
1969 
1970  /** Generates needed set of monoms, split them into sets S0, ... Sn and
1971  * check if reduced/nonreduced and calculate size of submatrix.
1972  */
1973  void generateMonomData( int deg, intvec* polyDegs , intvec* iVO );
1974 
1975  /** Recursively generate all homogeneous monoms of
1976  * (currRing->N) variables of degree deg.
1977  */
1978  void generateMonoms( poly m, int var, int deg );
1979 
1980  /** Creates quadratic matrix M of size numVectors for later use.
1981  * u0, u1, ...,un are replaced by 0.
1982  * Entries equal to 0 are not initialized ( == NULL)
1983  */
1984  void createMatrix();
1985 
1986 private:
1988 
1992  int subSize;
1993 
1995 };
1996 //<-
1997 
1998 //-> struct resVector
1999 /* Holds a row vector of the dense resultant matrix */
2001 {
2002 public:
2003  void init()
2004  {
2005  isReduced = FALSE;
2006  elementOfS = SFREE;
2007  mon = NULL;
2008  }
2009  void init( const poly m )
2010  {
2011  isReduced = FALSE;
2012  elementOfS = SFREE;
2013  mon = m;
2014  }
2015 
2016  /** index von 0 ... numVectors-1 */
2017  poly getElem( const int i );
2018 
2019  /** index von 0 ... numVectors-1 */
2020  number getElemNum( const int i );
2021 
2022  // variables
2023  poly mon;
2026 
2027  /** number of the set S mon is element of */
2029 
2030  /** holds the index of u0, u1, ..., un, if (elementOfS == linPolyS)
2031  * the size is given by (currRing->N)
2032  */
2034 
2035  /** holds the column vector if (elementOfS != linPolyS) */
2036  number *numColVector;
2037 
2038  /** size of numColVector */
2040 
2041  number *numColVecCopy;
2042 };
2043 //<-
2044 
2045 //-> resVector::*
2046 poly resVector::getElem( const int i ) // inline ???
2047 {
2048  assume( 0 < i || i > numColVectorSize );
2049  poly out= pOne();
2050  pSetCoeff( out, numColVector[i] );
2051  pTest( out );
2052  return( out );
2053 }
2054 
2055 number resVector::getElemNum( const int i ) // inline ??
2056 {
2057  assume( i >= 0 && i < numColVectorSize );
2058  return( numColVector[i] );
2059 }
2060 //<-
2061 
2062 //-> resMatrixDense::*
2063 resMatrixDense::resMatrixDense( const ideal _gls, const int special )
2064  : resMatrixBase()
2065 {
2066  int i;
2067 
2069  gls= idCopy( _gls );
2070  linPolyS= special;
2071  m=NULL;
2072 
2073  // init all
2074  generateBaseData();
2075 
2076  totDeg= 1;
2077  for ( i= 0; i < IDELEMS(gls); i++ )
2078  {
2079  totDeg*=pTotaldegree( (gls->m)[i] );
2080  }
2081 
2082  mprSTICKYPROT2(" resultant deg: %d\n",totDeg);
2083 
2085 }
2086 
2088 {
2089  int i,j;
2090  for (i=0; i < numVectors; i++)
2091  {
2092  pDelete( &resVectorList[i].mon );
2093  pDelete( &resVectorList[i].dividedBy );
2094  for ( j=0; j < resVectorList[i].numColVectorSize; j++ )
2095  {
2096  nDelete( resVectorList[i].numColVector+j );
2097  }
2098  // OB: ????? (solve_s.tst)
2099  if (resVectorList[i].numColVector!=NULL)
2100  omfreeSize( (void *)resVectorList[i].numColVector,
2101  numVectors * sizeof( number ) );
2102  if (resVectorList[i].numColParNr!=NULL)
2103  omfreeSize( (void *)resVectorList[i].numColParNr,
2104  ((currRing->N)+1) * sizeof(int) );
2105  }
2106 
2107  omFreeSize( (void *)resVectorList, veclistmax*sizeof( resVector ) );
2108 
2109  // free matrix m
2110  if ( m != NULL )
2111  {
2112  idDelete((ideal *)&m);
2113  }
2114 }
2115 
2116 // mprSTICKYPROT:
2117 // ST_DENSE_FR: found row S0
2118 // ST_DENSE_NR: normal row
2120 {
2121  int k,i,j;
2122  resVector *vecp;
2123 
2125 
2126  for ( i= 1; i <= MATROWS( m ); i++ )
2127  for ( j= 1; j <= MATCOLS( m ); j++ )
2128  {
2129  MATELEM(m,i,j)= pInit();
2130  pSetCoeff0( MATELEM(m,i,j), nInit(0) );
2131  }
2132 
2133 
2134  for ( k= 0; k <= numVectors - 1; k++ )
2135  {
2136  if ( linPolyS == getMVector(k)->elementOfS )
2137  {
2139  for ( i= 0; i < (currRing->N); i++ )
2140  {
2141  MATELEM(m,numVectors-k,numVectors-(getMVector(k)->numColParNr)[i])= pInit();
2142  }
2143  }
2144  else
2145  {
2147  vecp= getMVector(k);
2148  for ( i= 0; i < numVectors; i++)
2149  {
2150  if ( !nIsZero( vecp->getElemNum(i) ) )
2151  {
2152  MATELEM(m,numVectors - k,i + 1)= pInit();
2153  pSetCoeff0( MATELEM(m,numVectors - k,i + 1), nCopy(vecp->getElemNum(i)) );
2154  }
2155  }
2156  }
2157  } // for
2158  mprSTICKYPROT("\n");
2159 
2160 #ifdef mprDEBUG_ALL
2161  for ( k= numVectors - 1; k >= 0; k-- )
2162  {
2163  if ( linPolyS == getMVector(k)->elementOfS )
2164  {
2165  for ( i=0; i < (currRing->N); i++ )
2166  {
2167  Print(" %d ",(getMVector(k)->numColParNr)[i]);
2168  }
2169  PrintLn();
2170  }
2171  }
2172  for (i=1; i <= numVectors; i++)
2173  {
2174  for (j=1; j <= numVectors; j++ )
2175  {
2176  pWrite0(MATELEM(m,i,j));PrintS(" ");
2177  }
2178  PrintLn();
2179  }
2180 #endif
2181 }
2182 
2183 // mprSTICKYPROT:
2184 // ST_DENSE_MEM: more mem allocated
2185 // ST_DENSE_NMON: new monom added
2186 void resMatrixDense::generateMonoms( poly mm, int var, int deg )
2187 {
2188  if ( deg == 0 )
2189  {
2190  poly mon = pCopy( mm );
2191 
2192  if ( numVectors == veclistmax )
2193  {
2195  (veclistmax) * sizeof( resVector ),
2196  (veclistmax + veclistblock) * sizeof( resVector ) );
2197  int k;
2198  for ( k= veclistmax; k < (veclistmax + veclistblock); k++ )
2199  resVectorList[k].init();
2202 
2203  }
2204  resVectorList[numVectors].init( mon );
2205  numVectors++;
2207  return;
2208  }
2209  else
2210  {
2211  if ( var == (currRing->N)+1 ) return;
2212  poly newm = pCopy( mm );
2213  while ( deg >= 0 )
2214  {
2215  generateMonoms( newm, var+1, deg );
2216  pIncrExp( newm, var );
2217  pSetm( newm );
2218  deg--;
2219  }
2220  pDelete( & newm );
2221  }
2222 
2223  return;
2224 }
2225 
2226 void resMatrixDense::generateMonomData( int deg, intvec* polyDegs , intvec* iVO )
2227 {
2228  int i,j,k;
2229 
2230  // init monomData
2231  veclistblock= 512;
2234 
2235  // Init resVector()s
2236  for ( j= veclistmax - 1; j >= 0; j-- ) resVectorList[j].init();
2237  numVectors= 0;
2238 
2239  // Generate all monoms of degree deg
2240  poly start= pOne();
2241  generateMonoms( start, 1, deg );
2242  pDelete( & start );
2243 
2244  mprSTICKYPROT("\n");
2245 
2246  // Check for reduced monoms
2247  // First generate polyDegs.rows() monoms
2248  // x(k)^(polyDegs[k]), 0 <= k < polyDegs.rows()
2249  ideal pDegDiv= idInit( polyDegs->rows(), 1 );
2250  for ( k= 0; k < polyDegs->rows(); k++ )
2251  {
2252  poly p= pOne();
2253  pSetExp( p, k + 1, (*polyDegs)[k] );
2254  pSetm( p );
2255  (pDegDiv->m)[k]= p;
2256  }
2257 
2258  // Now check each monom if it is reduced.
2259  // A monom monom is called reduced if there exists
2260  // exactly one x(k)^(polyDegs[k]) that divides the monom.
2261  int divCount;
2262  for ( j= numVectors - 1; j >= 0; j-- )
2263  {
2264  divCount= 0;
2265  for ( k= 0; k < IDELEMS(pDegDiv); k++ )
2266  if ( pLmDivisibleByNoComp( (pDegDiv->m)[k], resVectorList[j].mon ) )
2267  divCount++;
2268  resVectorList[j].isReduced= (divCount == 1);
2269  }
2270 
2271  // create the sets S(k)s
2272  // a monom x(i)^deg, deg given, is element of the set S(i)
2273  // if all x(0)^(polyDegs[0]) ... x(i-1)^(polyDegs[i-1]) DONT divide
2274  // x(i)^deg and only x(i)^(polyDegs[i]) divides x(i)^deg
2275  bool doInsert;
2276  for ( k= 0; k < iVO->rows(); k++)
2277  {
2278  //mprPROTInl(" ------------ var:",(*iVO)[k]);
2279  for ( j= numVectors - 1; j >= 0; j-- )
2280  {
2281  //mprPROTPnl("testing monom",resVectorList[j].mon);
2282  if ( resVectorList[j].elementOfS == SFREE )
2283  {
2284  //mprPROTnl("\tfree");
2285  if ( pLmDivisibleByNoComp( (pDegDiv->m)[ (*iVO)[k] ], resVectorList[j].mon ) )
2286  {
2287  //mprPROTPnl("\tdivisible by ",(pDegDiv->m)[ (*iVO)[k] ]);
2288  doInsert=TRUE;
2289  for ( i= 0; i < k; i++ )
2290  {
2291  //mprPROTPnl("\tchecking db ",(pDegDiv->m)[ (*iVO)[i] ]);
2292  if ( pLmDivisibleByNoComp( (pDegDiv->m)[ (*iVO)[i] ], resVectorList[j].mon ) )
2293  {
2294  //mprPROTPnl("\t and divisible by",(pDegDiv->m)[ (*iVO)[i] ]);
2295  doInsert=FALSE;
2296  break;
2297  }
2298  }
2299  if ( doInsert )
2300  {
2301  //mprPROTInl("\t------------------> S ",(*iVO)[k]);
2302  resVectorList[j].elementOfS= (*iVO)[k];
2303  resVectorList[j].dividedBy= pCopy( (pDegDiv->m)[ (*iVO)[i] ] );
2304  }
2305  }
2306  }
2307  }
2308  }
2309 
2310  // size of submatrix M', equal to number of nonreduced monoms
2311  // (size of matrix M is equal to number of monoms=numVectors)
2312  subSize= 0;
2313  int sub;
2314  for ( i= 0; i < polyDegs->rows(); i++ )
2315  {
2316  sub= 1;
2317  for ( k= 0; k < polyDegs->rows(); k++ )
2318  if ( i != k ) sub*= (*polyDegs)[k];
2319  subSize+= sub;
2320  }
2322 
2323  // pDegDiv wieder freigeben!
2324  idDelete( &pDegDiv );
2325 
2326 #ifdef mprDEBUG_ALL
2327  // Print a list of monoms and their properties
2328  PrintS("// \n");
2329  for ( j= numVectors - 1; j >= 0; j-- )
2330  {
2331  Print("// %s, S(%d), db ",
2332  resVectorList[j].isReduced?"reduced":"nonreduced",
2333  resVectorList[j].elementOfS);
2334  pWrite0(resVectorList[j].dividedBy);
2335  PrintS(" monom ");
2336  pWrite(resVectorList[j].mon);
2337  }
2338  Print("// size: %d, subSize: %d\n",numVectors,subSize);
2339 #endif
2340 }
2341 
2343 {
2344  int k,j,i;
2345  number matEntry;
2346  poly pmatchPos;
2347  poly pi,factor,pmp;
2348 
2349  // holds the degrees of F0, F1, ..., Fn
2350  intvec polyDegs( IDELEMS(gls) );
2351  for ( k= 0; k < IDELEMS(gls); k++ )
2352  polyDegs[k]= pTotaldegree( (gls->m)[k] );
2353 
2354  // the internal Variable Ordering
2355  // make sure that the homogenization variable goes last!
2356  intvec iVO( (currRing->N) );
2357  if ( linPolyS != SNONE )
2358  {
2359  iVO[(currRing->N) - 1]= linPolyS;
2360  int p=0;
2361  for ( k= (currRing->N) - 1; k >= 0; k-- )
2362  {
2363  if ( k != linPolyS )
2364  {
2365  iVO[p]= k;
2366  p++;
2367  }
2368  }
2369  }
2370  else
2371  {
2372  linPolyS= 0;
2373  for ( k= 0; k < (currRing->N); k++ )
2374  iVO[k]= (currRing->N) - k - 1;
2375  }
2376 
2377  // the critical degree d= sum( deg(Fi) ) - n
2378  int sumDeg= 0;
2379  for ( k= 0; k < polyDegs.rows(); k++ )
2380  sumDeg+= polyDegs[k];
2381  sumDeg-= polyDegs.rows() - 1;
2382 
2383  // generate the base data
2384  generateMonomData( sumDeg, &polyDegs, &iVO );
2385 
2386  // generate "matrix"
2387  for ( k= numVectors - 1; k >= 0; k-- )
2388  {
2389  if ( resVectorList[k].elementOfS != linPolyS )
2390  {
2391  // column k is a normal column with numerical or symbolic entries
2392  // init stuff
2395  resVectorList[k].numColVector= (number *)omAlloc( numVectors*sizeof( number ) );
2396  for ( i= 0; i < numVectors; i++ ) resVectorList[k].numColVector[i]= nInit(0);
2397 
2398  // compute row poly
2399  poly pi= ppMult_qq( (gls->m)[ resVectorList[k].elementOfS ] , resVectorList[k].mon );
2400  pi= pp_DivideM( pi, resVectorList[k].dividedBy, currRing );
2401 
2402  // fill in "matrix"
2403  while ( pi != NULL )
2404  {
2405  matEntry= nCopy(pGetCoeff(pi));
2406  pmatchPos= pLmInit( pi );
2407  pSetCoeff0( pmatchPos, nInit(1) );
2408 
2409  for ( i= 0; i < numVectors; i++)
2410  if ( pLmEqual( pmatchPos, resVectorList[i].mon ) )
2411  break;
2412 
2413  resVectorList[k].numColVector[numVectors - i - 1] = nCopy(matEntry);
2414 
2415  pDelete( &pmatchPos );
2416  nDelete( &matEntry );
2417 
2418  pIter( pi );
2419  }
2420  pDelete( &pi );
2421  }
2422  else
2423  {
2424  // column is a special column, i.e. is generated by S0 and F0
2425  // safe only the positions of the ui's in the column
2426  //mprPROTInl(" setup of numColParNr ",k);
2429  resVectorList[k].numColParNr= (int *)omAlloc0( ((currRing->N)+1) * sizeof(int) );
2430 
2431  pi= (gls->m)[ resVectorList[k].elementOfS ];
2432  factor= pp_DivideM( resVectorList[k].mon, resVectorList[k].dividedBy, currRing );
2433 
2434  j=0;
2435  while ( pi != NULL )
2436  { // fill in "matrix"
2437  pmp= pMult( pCopy( factor ), pHead( pi ) );
2438  pTest( pmp );
2439 
2440  for ( i= 0; i < numVectors; i++)
2441  if ( pLmEqual( pmp, resVectorList[i].mon ) )
2442  break;
2443 
2445  pDelete( &pmp );
2446  pIter( pi );
2447  j++;
2448  }
2449  pDelete( &pi );
2450  pDelete( &factor );
2451  }
2452  } // for ( k= numVectors - 1; k >= 0; k-- )
2453 
2454  mprSTICKYPROT2(" size of matrix: %d\n",numVectors);
2455  mprSTICKYPROT2(" size of submatrix: %d\n",subSize);
2456 
2457  // create the matrix M
2458  createMatrix();
2459 
2460 }
2461 
2463 {
2464  assume( i >= 0 && i < numVectors );
2465  return &resVectorList[i];
2466 }
2467 
2469 {
2470  int i,j;
2471 
2472  // copy matrix
2473  matrix resmat= mpNew(numVectors,numVectors);
2474  poly p;
2475  for (i=1; i <= numVectors; i++)
2476  {
2477  for (j=1; j <= numVectors; j++ )
2478  {
2479  p=MATELEM(m,i,j);
2480  if (( p!=NULL)
2481  && (!nIsZero(pGetCoeff(p)))
2482  && (pGetCoeff(p)!=NULL)
2483  )
2484  {
2485  MATELEM(resmat,i,j)= pCopy( p );
2486  }
2487  }
2488  }
2489  for (i=0; i < numVectors; i++)
2490  {
2491  if ( resVectorList[i].elementOfS == linPolyS )
2492  {
2493  for (j=1; j <= (currRing->N); j++ )
2494  {
2495  if ( MATELEM(resmat,numVectors-i,
2496  numVectors-resVectorList[i].numColParNr[j-1])!=NULL )
2497  pDelete( &MATELEM(resmat,numVectors-i,numVectors-resVectorList[i].numColParNr[j-1]) );
2498  MATELEM(resmat,numVectors-i,numVectors-resVectorList[i].numColParNr[j-1])= pOne();
2499  // FIX ME
2500  if ( FALSE )
2501  {
2502  pSetCoeff( MATELEM(resmat,numVectors-i,numVectors-resVectorList[i].numColParNr[j-1]), n_Param(j,currRing) );
2503  }
2504  else
2505  {
2506  pSetExp( MATELEM(resmat,numVectors-i,numVectors-resVectorList[i].numColParNr[j-1]), j, 1 );
2507  pSetm(MATELEM(resmat,numVectors-i,numVectors-resVectorList[i].numColParNr[j-1]));
2508  }
2509  }
2510  }
2511  }
2512 
2513  // obachman: idMatrix2Module frees resmat !!
2514  ideal resmod= id_Matrix2Module(resmat,currRing);
2515  return resmod;
2516 }
2517 
2519 {
2520  int k,i,j,l;
2521  resVector *vecp;
2522 
2523  // generate quadratic matrix resmat of size subSize
2524  matrix resmat= mpNew( subSize, subSize );
2525 
2526  j=1;
2527  for ( k= numVectors - 1; k >= 0; k-- )
2528  {
2529  vecp= getMVector(k);
2530  if ( vecp->isReduced ) continue;
2531  l=1;
2532  for ( i= numVectors - 1; i >= 0; i-- )
2533  {
2534  if ( getMVector(i)->isReduced ) continue;
2535  if ( !nIsZero(vecp->getElemNum(numVectors - i - 1)) )
2536  {
2537  MATELEM(resmat,j,l)= pCopy( vecp->getElem(numVectors-i-1) );
2538  }
2539  l++;
2540  }
2541  j++;
2542  }
2543 
2544  // obachman: idMatrix2Module frees resmat !!
2545  ideal resmod= id_Matrix2Module(resmat,currRing);
2546  return resmod;
2547 }
2548 
2549 number resMatrixDense::getDetAt( const number* evpoint )
2550 {
2551  int k,i;
2552 
2553  // copy evaluation point into matrix
2554  // p0, p1, ..., pn replace u0, u1, ..., un
2555  for ( k= numVectors - 1; k >= 0; k-- )
2556  {
2557  if ( linPolyS == getMVector(k)->elementOfS )
2558  {
2559  for ( i= 0; i < (currRing->N); i++ )
2560  {
2561  number np=pGetCoeff(MATELEM(m,numVectors-k,numVectors-(getMVector(k)->numColParNr)[i]));
2562  if (np!=NULL) nDelete(&np);
2563  pSetCoeff0( MATELEM(m,numVectors-k,numVectors-(getMVector(k)->numColParNr)[i]),
2564  nCopy(evpoint[i]) );
2565  }
2566  }
2567  }
2568 
2570 
2571  // evaluate determinant of matrix m using factory singclap_det
2572  poly res= singclap_det( m, currRing );
2573 
2574  // avoid errors for det==0
2575  number numres;
2576  if ( (res!=NULL) && (!nIsZero(pGetCoeff( res ))) )
2577  {
2578  numres= nCopy( pGetCoeff( res ) );
2579  }
2580  else
2581  {
2582  numres= nInit(0);
2583  mprPROT("0");
2584  }
2585  pDelete( &res );
2586 
2588 
2589  return( numres );
2590 }
2591 
2593 {
2594  int k,i,j,l;
2595  resVector *vecp;
2596 
2597  // generate quadratic matrix mat of size subSize
2598  matrix mat= mpNew( subSize, subSize );
2599 
2600  for ( i= 1; i <= MATROWS( mat ); i++ )
2601  {
2602  for ( j= 1; j <= MATCOLS( mat ); j++ )
2603  {
2604  MATELEM(mat,i,j)= pInit();
2605  pSetCoeff0( MATELEM(mat,i,j), nInit(0) );
2606  }
2607  }
2608  j=1;
2609  for ( k= numVectors - 1; k >= 0; k-- )
2610  {
2611  vecp= getMVector(k);
2612  if ( vecp->isReduced ) continue;
2613  l=1;
2614  for ( i= numVectors - 1; i >= 0; i-- )
2615  {
2616  if ( getMVector(i)->isReduced ) continue;
2617  if ( vecp->getElemNum(numVectors - i - 1) && !nIsZero(vecp->getElemNum(numVectors - i - 1)) )
2618  {
2619  pSetCoeff(MATELEM(mat, j , l ), nCopy(vecp->getElemNum(numVectors - i - 1)));
2620  }
2621  /* else
2622  {
2623  MATELEM(mat, j , l )= pOne();
2624  pSetCoeff(MATELEM(mat, j , l ), nInit(0) );
2625  }
2626  */
2627  l++;
2628  }
2629  j++;
2630  }
2631 
2632  poly res= singclap_det( mat, currRing );
2633 
2634  number numres;
2635  if ((res != NULL) && (!nIsZero(pGetCoeff( res ))) )
2636  {
2637  numres= nCopy(pGetCoeff( res ));
2638  }
2639  else
2640  {
2641  numres= nInit(0);
2642  }
2643  pDelete( &res );
2644  return numres;
2645 }
2646 //<--
2647 
2648 //-----------------------------------------------------------------------------
2649 //-------------- uResultant ---------------------------------------------------
2650 //-----------------------------------------------------------------------------
2651 
2652 #define MAXEVPOINT 1000000 // 0x7fffffff
2653 //#define MPR_MASI
2654 
2655 //-> unsigned long over(unsigned long n,unsigned long d)
2656 // Calculates (n+d \over d) using gmp functionality
2657 //
2658 unsigned long over( const unsigned long n , const unsigned long d )
2659 { // (d+n)! / ( d! n! )
2660  mpz_t res;
2661  mpz_init(res);
2662  mpz_t m,md,mn;
2663  mpz_init(m);mpz_set_ui(m,1);
2664  mpz_init(md);mpz_set_ui(md,1);
2665  mpz_init(mn);mpz_set_ui(mn,1);
2666 
2667  mpz_fac_ui(m,n+d);
2668  mpz_fac_ui(md,d);
2669  mpz_fac_ui(mn,n);
2670 
2671  mpz_mul(res,md,mn);
2672  mpz_tdiv_q(res,m,res);
2673 
2674  mpz_clear(m);mpz_clear(md);mpz_clear(mn);
2675 
2676  unsigned long result = mpz_get_ui(res);
2677  mpz_clear(res);
2678 
2679  return result;
2680 }
2681 //<-
2682 
2683 //-> uResultant::*
2684 uResultant::uResultant( const ideal _gls, const resMatType _rmt, BOOLEAN extIdeal )
2685  : rmt( _rmt )
2686 {
2687  if ( extIdeal )
2688  {
2689  // extend given ideal by linear poly F0=u0x0 + u1x1 +...+ unxn
2690  gls= extendIdeal( _gls, linearPoly( rmt ), rmt );
2691  n= IDELEMS( gls );
2692  }
2693  else
2694  gls= idCopy( _gls );
2695 
2696  switch ( rmt )
2697  {
2698  case sparseResMat:
2699  resMat= new resMatrixSparse( gls );
2700  break;
2701  case denseResMat:
2702  resMat= new resMatrixDense( gls );
2703  break;
2704  default:
2705  WerrorS("uResultant::uResultant: Unknown chosen resultant matrix type!");
2706  }
2707 }
2708 
2710 {
2711  delete resMat;
2712 }
2713 
2714 ideal uResultant::extendIdeal( const ideal igls, poly linPoly, const resMatType rrmt )
2715 {
2716  ideal newGls= idCopy( igls );
2717  newGls->m= (poly *)omReallocSize( newGls->m,
2718  IDELEMS(igls) * sizeof(poly),
2719  (IDELEMS(igls) + 1) * sizeof(poly) );
2720  IDELEMS(newGls)++;
2721 
2722  switch ( rrmt )
2723  {
2724  case sparseResMat:
2725  case denseResMat:
2726  {
2727  int i;
2728  for ( i= IDELEMS(newGls)-1; i > 0; i-- )
2729  {
2730  newGls->m[i]= newGls->m[i-1];
2731  }
2732  newGls->m[0]= linPoly;
2733  }
2734  break;
2735  default:
2736  WerrorS("uResultant::extendIdeal: Unknown chosen resultant matrix type!");
2737  }
2738 
2739  return( newGls );
2740 }
2741 
2743 {
2744  int i;
2745 
2746  poly newlp= pOne();
2747  poly actlp, rootlp= newlp;
2748 
2749  for ( i= 1; i <= (currRing->N); i++ )
2750  {
2751  actlp= newlp;
2752  pSetExp( actlp, i, 1 );
2753  pSetm( actlp );
2754  newlp= pOne();
2755  actlp->next= newlp;
2756  }
2757  actlp->next= NULL;
2758  pDelete( &newlp );
2759 
2760  if ( rrmt == sparseResMat )
2761  {
2762  newlp= pOne();
2763  actlp->next= newlp;
2764  newlp->next= NULL;
2765  }
2766  return ( rootlp );
2767 }
2768 
2769 poly uResultant::interpolateDense( const number subDetVal )
2770 {
2771  int i,j,p;
2772  long tdg;
2773 
2774  // D is a Polynom homogeneous in the coeffs of F0 and of degree tdg = d0*d1*...*dn
2775  tdg= resMat->getDetDeg();
2776 
2777  // maximum number of terms in polynom D (homogeneous, of degree tdg)
2778  // long mdg= (facul(tdg+n-1) / facul( tdg )) / facul( n - 1 );
2779  long mdg= over( n-1, tdg );
2780 
2781  // maximal number of terms in a polynom of degree tdg
2782  long l=(long)pow( (double)(tdg+1), n );
2783 
2784 #ifdef mprDEBUG_PROT
2785  Print("// total deg of D: tdg %ld\n",tdg);
2786  Print("// maximum number of terms in D: mdg: %ld\n",mdg);
2787  Print("// maximum number of terms in polynom of deg tdg: l %ld\n",l);
2788 #endif
2789 
2790  // we need mdg results of D(p0,p1,...,pn)
2791  number *presults;
2792  presults= (number *)omAlloc( mdg * sizeof( number ) );
2793  for (i=0; i < mdg; i++) presults[i]= nInit(0);
2794 
2795  number *pevpoint= (number *)omAlloc( n * sizeof( number ) );
2796  number *pev= (number *)omAlloc( n * sizeof( number ) );
2797  for (i=0; i < n; i++) pev[i]= nInit(0);
2798 
2799  mprPROTnl("// initial evaluation point: ");
2800  // initial evaluatoin point
2801  p=1;
2802  for (i=0; i < n; i++)
2803  {
2804  // init pevpoint with primes 3,5,7,11, ...
2805  p= nextPrime( p );
2806  pevpoint[i]=nInit( p );
2807  nTest(pevpoint[i]);
2808  mprPROTNnl(" ",pevpoint[i]);
2809  }
2810 
2811  // evaluate the determinant in the points pev^0, pev^1, ..., pev^mdg
2812  mprPROTnl("// evaluating:");
2813  for ( i=0; i < mdg; i++ )
2814  {
2815  for (j=0; j < n; j++)
2816  {
2817  nDelete( &pev[j] );
2818  nPower(pevpoint[j],i,&pev[j]);
2819  mprPROTN(" ",pev[j]);
2820  }
2821  mprPROTnl("");
2822 
2823  nDelete( &presults[i] );
2824  presults[i]=resMat->getDetAt( pev );
2825 
2827  }
2828  mprSTICKYPROT("\n");
2829 
2830  // now interpolate using vandermode interpolation
2831  mprPROTnl("// interpolating:");
2832  number *ncpoly;
2833  {
2834  vandermonde vm( mdg, n, tdg, pevpoint );
2835  ncpoly= vm.interpolateDense( presults );
2836  }
2837 
2838  if ( subDetVal != NULL )
2839  { // divide by common factor
2840  number detdiv;
2841  for ( i= 0; i <= mdg; i++ )
2842  {
2843  detdiv= nDiv( ncpoly[i], subDetVal );
2844  nNormalize( detdiv );
2845  nDelete( &ncpoly[i] );
2846  ncpoly[i]= detdiv;
2847  }
2848  }
2849 
2850 #ifdef mprDEBUG_ALL
2851  PrintLn();
2852  for ( i=0; i < mdg; i++ )
2853  {
2854  nPrint(ncpoly[i]); PrintS(" --- ");
2855  }
2856  PrintLn();
2857 #endif
2858 
2859  // prepare ncpoly for later use
2860  number nn=nInit(0);
2861  for ( i=0; i < mdg; i++ )
2862  {
2863  if ( nEqual(ncpoly[i],nn) )
2864  {
2865  nDelete( &ncpoly[i] );
2866  ncpoly[i]=NULL;
2867  }
2868  }
2869  nDelete( &nn );
2870 
2871  // create poly presenting the determinat of the uResultant
2872  intvec exp( n );
2873  for ( i= 0; i < n; i++ ) exp[i]=0;
2874 
2875  poly result= NULL;
2876 
2877  long sum=0;
2878  long c=0;
2879 
2880  for ( i=0; i < l; i++ )
2881  {
2882  if ( sum == tdg )
2883  {
2884  if ( !nIsZero(ncpoly[c]) )
2885  {
2886  poly p= pOne();
2887  if ( rmt == denseResMat )
2888  {
2889  for ( j= 0; j < n; j++ ) pSetExp( p, j+1, exp[j] );
2890  }
2891  else if ( rmt == sparseResMat )
2892  {
2893  for ( j= 1; j < n; j++ ) pSetExp( p, j, exp[j] );
2894  }
2895  pSetCoeff( p, ncpoly[c] );
2896  pSetm( p );
2897  if (result!=NULL) result= pAdd( result, p );
2898  else result= p;
2899  }
2900  c++;
2901  }
2902  sum=0;
2903  exp[0]++;
2904  for ( j= 0; j < n - 1; j++ )
2905  {
2906  if ( exp[j] > tdg )
2907  {
2908  exp[j]= 0;
2909  exp[j + 1]++;
2910  }
2911  sum+=exp[j];
2912  }
2913  sum+=exp[n-1];
2914  }
2915 
2916  pTest( result );
2917 
2918  return result;
2919 }
2920 
2921 rootContainer ** uResultant::interpolateDenseSP( BOOLEAN matchUp, const number subDetVal )
2922 {
2923  int i,p,uvar;
2924  long tdg;
2925  int loops= (matchUp?n-2:n-1);
2926 
2927  mprPROTnl("uResultant::interpolateDenseSP");
2928 
2929  tdg= resMat->getDetDeg();
2930 
2931  // evaluate D in tdg+1 distinct points, so
2932  // we need tdg+1 results of D(p0,1,0,...,0) =
2933  // c(0)*u0^tdg + c(1)*u0^tdg-1 + ... + c(tdg-1)*u0 + c(tdg)
2934  number *presults;
2935  presults= (number *)omAlloc( (tdg + 1) * sizeof( number ) );
2936  for ( i=0; i <= tdg; i++ ) presults[i]= nInit(0);
2937 
2938  rootContainer ** roots;
2939  roots= (rootContainer **) omAlloc( loops * sizeof(rootContainer*) );
2940  for ( i=0; i < loops; i++ ) roots[i]= new rootContainer(); // 0..n-2
2941 
2942  number *pevpoint= (number *)omAlloc( n * sizeof( number ) );
2943  for (i=0; i < n; i++) pevpoint[i]= nInit(0);
2944 
2945  number *pev= (number *)omAlloc( n * sizeof( number ) );
2946  for (i=0; i < n; i++) pev[i]= nInit(0);
2947 
2948  // now we evaluate D(u0,-1,0,...0), D(u0,0,-1,0,...,0), ..., D(u0,0,..,0,-1)
2949  // or D(u0,k1,k2,0,...,0), D(u0,k1,k2,k3,0,...,0), ..., D(u0,k1,k2,k3,...,kn)
2950  // this gives us n-1 evaluations
2951  p=3;
2952  for ( uvar= 0; uvar < loops; uvar++ )
2953  {
2954  // generate initial evaluation point
2955  if ( matchUp )
2956  {
2957  for (i=0; i < n; i++)
2958  {
2959  // prime(random number) between 1 and MAXEVPOINT
2960  nDelete( &pevpoint[i] );
2961  if ( i == 0 )
2962  {
2963  //p= nextPrime( p );
2964  pevpoint[i]= nInit( p );
2965  }
2966  else if ( i <= uvar + 2 )
2967  {
2968  pevpoint[i]=nInit(1+siRand()%MAXEVPOINT);
2969  //pevpoint[i]=nInit(383);
2970  }
2971  else
2972  pevpoint[i]=nInit(0);
2973  mprPROTNnl(" ",pevpoint[i]);
2974  }
2975  }
2976  else
2977  {
2978  for (i=0; i < n; i++)
2979  {
2980  // init pevpoint with prime,0,...0,1,0,...,0
2981  nDelete( &pevpoint[i] );
2982  if ( i == 0 )
2983  {
2984  //p=nextPrime( p );
2985  pevpoint[i]=nInit( p );
2986  }
2987  else
2988  {
2989  if ( i == (uvar + 1) ) pevpoint[i]= nInit(-1);
2990  else pevpoint[i]= nInit(0);
2991  }
2992  mprPROTNnl(" ",pevpoint[i]);
2993  }
2994  }
2995 
2996  // prepare aktual evaluation point
2997  for (i=0; i < n; i++)
2998  {
2999  nDelete( &pev[i] );
3000  pev[i]= nCopy( pevpoint[i] );
3001  }
3002  // evaluate the determinant in the points pev^0, pev^1, ..., pev^tdg
3003  for ( i=0; i <= tdg; i++ )
3004  {
3005  nDelete( &pev[0] );
3006  nPower(pevpoint[0],i,&pev[0]); // new evpoint
3007 
3008  nDelete( &presults[i] );
3009  presults[i]=resMat->getDetAt( pev ); // evaluate det at point evpoint
3010 
3011  mprPROTNnl("",presults[i]);
3012 
3014  mprPROTL("",tdg-i);
3015  }
3016  mprSTICKYPROT("\n");
3017 
3018  // now interpolate
3019  vandermonde vm( tdg + 1, 1, tdg, pevpoint, FALSE );
3020  number *ncpoly= vm.interpolateDense( presults );
3021 
3022  if ( subDetVal != NULL )
3023  { // divide by common factor
3024  number detdiv;
3025  for ( i= 0; i <= tdg; i++ )
3026  {
3027  detdiv= nDiv( ncpoly[i], subDetVal );
3028  nNormalize( detdiv );
3029  nDelete( &ncpoly[i] );
3030  ncpoly[i]= detdiv;
3031  }
3032  }
3033 
3034 #ifdef mprDEBUG_ALL
3035  PrintLn();
3036  for ( i=0; i <= tdg; i++ )
3037  {
3038  nPrint(ncpoly[i]); PrintS(" --- ");
3039  }
3040  PrintLn();
3041 #endif
3042 
3043  // save results
3044  roots[uvar]->fillContainer( ncpoly, pevpoint, uvar+1, tdg,
3046  loops );
3047  }
3048 
3049  // free some stuff: pev, presult
3050  for ( i=0; i < n; i++ ) nDelete( pev + i );
3051  omFreeSize( (void *)pev, n * sizeof( number ) );
3052 
3053  for ( i=0; i <= tdg; i++ ) nDelete( presults+i );
3054  omFreeSize( (void *)presults, (tdg + 1) * sizeof( number ) );
3055 
3056  return roots;
3057 }
3058 
3059 rootContainer ** uResultant::specializeInU( BOOLEAN matchUp, const number subDetVal )
3060 {
3061  int i,/*p,*/uvar;
3062  long tdg;
3063  poly pures,piter;
3064  int loops=(matchUp?n-2:n-1);
3065  int nn=n;
3066  if (loops==0) { loops=1;nn++;}
3067 
3068  mprPROTnl("uResultant::specializeInU");
3069 
3070  tdg= resMat->getDetDeg();
3071 
3072  rootContainer ** roots;
3073  roots= (rootContainer **) omAlloc( loops * sizeof(rootContainer*) );
3074  for ( i=0; i < loops; i++ ) roots[i]= new rootContainer(); // 0..n-2
3075 
3076  number *pevpoint= (number *)omAlloc( nn * sizeof( number ) );
3077  for (i=0; i < nn; i++) pevpoint[i]= nInit(0);
3078 
3079  // now we evaluate D(u0,-1,0,...0), D(u0,0,-1,0,...,0), ..., D(u0,0,..,0,-1)
3080  // or D(u0,k1,k2,0,...,0), D(u0,k1,k2,k3,0,...,0), ..., D(u0,k1,k2,k3,...,kn)
3081  // p=3;
3082  for ( uvar= 0; uvar < loops; uvar++ )
3083  {
3084  // generate initial evaluation point
3085  if ( matchUp )
3086  {
3087  for (i=0; i < n; i++)
3088  {
3089  // prime(random number) between 1 and MAXEVPOINT
3090  nDelete( &pevpoint[i] );
3091  if ( i <= uvar + 2 )
3092  {
3093  pevpoint[i]=nInit(1+siRand()%MAXEVPOINT);
3094  //pevpoint[i]=nInit(383);
3095  }
3096  else pevpoint[i]=nInit(0);
3097  mprPROTNnl(" ",pevpoint[i]);
3098  }
3099  }
3100  else
3101  {
3102  for (i=0; i < n; i++)
3103  {
3104  // init pevpoint with prime,0,...0,-1,0,...,0
3105  nDelete( &(pevpoint[i]) );
3106  if ( i == (uvar + 1) ) pevpoint[i]= nInit(-1);
3107  else pevpoint[i]= nInit(0);
3108  mprPROTNnl(" ",pevpoint[i]);
3109  }
3110  }
3111 
3112  pures= resMat->getUDet( pevpoint );
3113 
3114  number *ncpoly= (number *)omAlloc( (tdg+1) * sizeof(number) );
3115 
3116 #ifdef MPR_MASI
3117  BOOLEAN masi=true;
3118 #endif
3119 
3120  piter= pures;
3121  for ( i= tdg; i >= 0; i-- )
3122  {
3123  if ( piter && pTotaldegree(piter) == i )
3124  {
3125  ncpoly[i]= nCopy( pGetCoeff( piter ) );
3126  pIter( piter );
3127 #ifdef MPR_MASI
3128  masi=false;
3129 #endif
3130  }
3131  else
3132  {
3133  ncpoly[i]= nInit(0);
3134  }
3135  mprPROTNnl("", ncpoly[i] );
3136  }
3137 #ifdef MPR_MASI
3138  if ( masi ) mprSTICKYPROT("MASI");
3139 #endif
3140 
3141  mprSTICKYPROT(ST_BASE_EV); // .
3142 
3143  if ( subDetVal != NULL ) // divide by common factor
3144  {
3145  number detdiv;
3146  for ( i= 0; i <= tdg; i++ )
3147  {
3148  detdiv= nDiv( ncpoly[i], subDetVal );
3149  nNormalize( detdiv );
3150  nDelete( &ncpoly[i] );
3151  ncpoly[i]= detdiv;
3152  }
3153  }
3154 
3155  pDelete( &pures );
3156 
3157  // save results
3158  roots[uvar]->fillContainer( ncpoly, pevpoint, uvar+1, tdg,
3160  loops );
3161  }
3162 
3163  mprSTICKYPROT("\n");
3164 
3165  // free some stuff: pev, presult
3166  for ( i=0; i < n; i++ ) nDelete( pevpoint + i );
3167  omFreeSize( (void *)pevpoint, n * sizeof( number ) );
3168 
3169  return roots;
3170 }
3171 
3172 int uResultant::nextPrime( const int i )
3173 {
3174  int init=i;
3175  int ii=i+2;
3176  extern int IsPrime(int p); // from Singular/ipshell.{h,cc}
3177  int j= IsPrime( ii );
3178  while ( j <= init )
3179  {
3180  ii+=2;
3181  j= IsPrime( ii );
3182  }
3183  return j;
3184 }
3185 //<-
3186 
3187 //-----------------------------------------------------------------------------
3188 
3189 //-> loNewtonPolytope(...)
3190 ideal loNewtonPolytope( const ideal id )
3191 {
3192  simplex * LP;
3193  int i;
3194  int /*n,*/totverts,idelem;
3195  ideal idr;
3196 
3197  // n= (currRing->N);
3198  idelem= IDELEMS(id); // should be n+1
3199 
3200  totverts = 0;
3201  for( i=0; i < idelem; i++) totverts += pLength( (id->m)[i] );
3202 
3203  LP = new simplex( idelem+totverts*2+5, totverts+5 ); // rows, cols
3204 
3205  // evaluate convex hull for supports of id
3206  convexHull chnp( LP );
3207  idr = chnp.newtonPolytopesI( id );
3208 
3209  delete LP;
3210 
3211  return idr;
3212 }
3213 //<-
3214 
3215 //%e
3216 
3217 //-----------------------------------------------------------------------------
3218 
3219 // local Variables: ***
3220 // folded-file: t ***
3221 // compile-command-1: "make installg" ***
3222 // compile-command-2: "make install" ***
3223 // End: ***
3224 
3225 // in folding: C-c x
3226 // leave fold: C-c y
3227 // foldmode: F10
Rational pow(const Rational &a, int e)
Definition: GMPrat.cc:411
int BOOLEAN
Definition: auxiliary.h:87
#define TRUE
Definition: auxiliary.h:100
#define FALSE
Definition: auxiliary.h:96
int size(const CanonicalForm &f, const Variable &v)
int size ( const CanonicalForm & f, const Variable & v )
Definition: cf_ops.cc:600
CanonicalForm FACTORY_PUBLIC pp(const CanonicalForm &)
CanonicalForm pp ( const CanonicalForm & f )
Definition: cf_gcd.cc:676
CanonicalForm num(const CanonicalForm &f)
const CanonicalForm CFMap CFMap & N
Definition: cfEzgcd.cc:56
int l
Definition: cfEzgcd.cc:100
int m
Definition: cfEzgcd.cc:128
for(int i=0;i<=n;i++) degsf[i]
Definition: cfEzgcd.cc:72
int i
Definition: cfEzgcd.cc:132
int k
Definition: cfEzgcd.cc:99
int p
Definition: cfModGcd.cc:4078
CanonicalForm cd(bCommonDen(FF))
Definition: cfModGcd.cc:4089
CanonicalForm b
Definition: cfModGcd.cc:4103
int int ncols
Definition: cf_linsys.cc:32
int nrows
Definition: cf_linsys.cc:32
poly singclap_det(const matrix m, const ring s)
Definition: clapsing.cc:1757
pointSet ** newtonPolytopesP(const ideal gls)
Computes the point sets of the convex hulls of the supports given by the polynoms in gls.
Definition: mpr_base.cc:776
ideal newtonPolytopesI(const ideal gls)
Definition: mpr_base.cc:834
pointSet ** Q
Definition: mpr_base.cc:269
bool inHull(poly p, poly pointPoly, int m, int site)
Returns true iff the support of poly pointPoly is inside the convex hull of all points given by the s...
Definition: mpr_base.cc:730
convexHull(simplex *_pLP)
Definition: mpr_base.cc:252
simplex * pLP
Definition: mpr_base.cc:271
Definition: intvec.h:23
int rows() const
Definition: intvec.h:96
pointSet * getInnerPoints(pointSet **_q_i, mprfloat _shift[])
Drive Mayan Pyramid Algorithm.
Definition: mpr_base.cc:891
pointSet * E
Definition: mpr_base.cc:318
bool storeMinkowskiSumPoint()
Stores point in E->points[pt], iff v-distance != 0 Returns true iff point was stored,...
Definition: mpr_base.cc:1138
void runMayanPyramid(int dim)
Recursive Mayan Pyramid algorithm for directly computing MinkowskiSum lattice points for (n+1)-fold M...
Definition: mpr_base.cc:1162
simplex * pLP
Definition: mpr_base.cc:325
mayanPyramidAlg(simplex *_pLP)
Definition: mpr_base.cc:280
mprfloat vDistance(Coord_t *acoords, int dim)
Compute v-distance via Linear Programing Linear Program finds the v-distance of the point in accords[...
Definition: mpr_base.cc:909
void mn_mx_MinkowskiSum(int dim, Coord_t *minR, Coord_t *maxR)
LP for finding min/max coord in MinkowskiSum, given previous coors.
Definition: mpr_base.cc:996
pointSet ** Qi
Definition: mpr_base.cc:317
mprfloat * shift
Definition: mpr_base.cc:319
Coord_t acoords[MAXVARS+2]
Definition: mpr_base.cc:323
bool addPoint(const onePointP vert)
Adds a point to pointSet, copy vert[0,...,dim] ot point[num+1][0,...,dim].
Definition: mpr_base.cc:464
bool lifted
Definition: mpr_base.cc:164
onePointP operator[](const int index)
Definition: mpr_base.cc:437
void mergeWithPoly(const poly p)
Definition: mpr_base.cc:550
void unlift()
Definition: mpr_base.cc:229
bool larger(int, int)
points[a] > points[b] ?
Definition: mpr_base.cc:628
void lift(int *l=NULL)
Lifts the point set using sufficiently generic linear lifting homogeneous forms l[1]....
Definition: mpr_base.cc:670
void sort()
sort lex
Definition: mpr_base.cc:647
int num
Definition: mpr_base.cc:167
pointSet(const pointSet &)
int dim
Definition: mpr_base.cc:169
int getExpPos(const poly p)
Definition: mpr_base.cc:578
int index
Definition: mpr_base.cc:170
onePointP * points
Definition: mpr_base.cc:163
bool mergeWithExp(const onePointP vert)
Adds point to pointSet, iff pointSet \cap point = \emptyset.
Definition: mpr_base.cc:512
void getRowMP(const int indx, int *vert)
Definition: mpr_base.cc:599
bool removePoint(const int indx)
Definition: mpr_base.cc:497
int max
Definition: mpr_base.cc:168
~pointSet()
Definition: mpr_base.cc:425
bool smaller(int, int)
points[a] < points[b] ?
Definition: mpr_base.cc:609
pointSet(const int _dim, const int _index=0, const int count=MAXINITELEMS)
Definition: mpr_base.cc:412
bool checkMem()
Checks, if more mem is needed ( i.e.
Definition: mpr_base.cc:443
Base class for sparse and dense u-Resultant computation.
Definition: mpr_base.h:23
ideal gls
Definition: mpr_base.h:46
virtual poly getUDet(const number *)
Definition: mpr_base.h:34
ring sourceRing
Definition: mpr_base.h:48
virtual number getDetAt(const number *)
Definition: mpr_base.h:36
int linPolyS
Definition: mpr_base.h:47
virtual long getDetDeg()
Definition: mpr_base.h:39
IStateType istate
Definition: mpr_base.h:44
number getSubDet()
Evaluates the determinant of the submatrix M'.
Definition: mpr_base.cc:2592
void generateBaseData()
Generate the "matrix" M.
Definition: mpr_base.cc:2342
ideal getSubMatrix()
Returns the submatrix M' of M in an usable presentation.
Definition: mpr_base.cc:2518
void generateMonoms(poly m, int var, int deg)
Recursively generate all homogeneous monoms of (currRing->N) variables of degree deg.
Definition: mpr_base.cc:2186
resMatrixDense(const ideal _gls, const int special=SNONE)
_gls: system of multivariate polynoms special: -1 -> resMatrixDense is a symbolic matrix 0,...
Definition: mpr_base.cc:2063
number getDetAt(const number *evpoint)
Evaluate the determinant of the matrix M at the point evpoint where the ui's are replaced by the comp...
Definition: mpr_base.cc:2549
resVector * getMVector(const int i)
column vector of matrix, index von 0 ...
Definition: mpr_base.cc:2462
void createMatrix()
Creates quadratic matrix M of size numVectors for later use.
Definition: mpr_base.cc:2119
ideal getMatrix()
Returns the matrix M in an usable presentation.
Definition: mpr_base.cc:2468
resMatrixDense(const resMatrixDense &)
deactivated copy constructor
void generateMonomData(int deg, intvec *polyDegs, intvec *iVO)
Generates needed set of monoms, split them into sets S0, ...
Definition: mpr_base.cc:2226
resVector * resVectorList
Definition: mpr_base.cc:1987
bool remapXiToPoint(const int indx, pointSet **pQ, int *set, int *vtx)
Definition: mpr_base.cc:1218
int RC(pointSet **pQ, pointSet *E, int vert, mprfloat shift[])
Row Content Function Finds the largest i such that F[i] is a point, F[i]= a[ij] in A[i] for some j.
Definition: mpr_base.cc:1237
resMatrixSparse(const ideal _gls, const int special=SNONE)
Definition: mpr_base.cc:1571
void randomVector(const int dim, mprfloat shift[])
Definition: mpr_base.cc:1501
pointSet * minkSumTwo(pointSet *Q1, pointSet *Q2, int dim)
Definition: mpr_base.cc:1521
simplex * LP
Definition: mpr_base.cc:124
ideal getMatrix()
Definition: mpr_base.cc:1734
int createMatrix(pointSet *E)
create coeff matrix uRPos[i][1]: row of matrix uRPos[i][idelem+1]: col of u(0) uRPos[i][2....
Definition: mpr_base.cc:1409
resMatrixSparse(const resMatrixSparse &)
intvec * uRPos
Definition: mpr_base.cc:120
poly getUDet(const number *evpoint)
Definition: mpr_base.cc:1855
pointSet * minkSumAll(pointSet **pQ, int numq, int dim)
Definition: mpr_base.cc:1549
number getDetAt(const number *evpoint)
Fills in resMat[][] with evpoint[] and gets determinant uRPos[i][1]: row of matrix uRPos[i][idelem+1]...
Definition: mpr_base.cc:1795
complex root finder for univariate polynomials based on laguers algorithm
Definition: mpr_numeric.h:66
void fillContainer(number *_coeffs, number *_ievpoint, const int _var, const int _tdg, const rootType _rt, const int _anz)
Definition: mpr_numeric.cc:300
Linear Programming / Linear Optimization using Simplex - Algorithm.
Definition: mpr_numeric.h:195
mprfloat ** LiPM
Definition: mpr_numeric.h:205
int * iposv
Definition: mpr_numeric.h:203
int icase
Definition: mpr_numeric.h:201
void compute()
poly linearPoly(const resMatType rmt)
Definition: mpr_base.cc:2742
resMatType rmt
Definition: mpr_base.h:91
rootContainer ** specializeInU(BOOLEAN matchUp=false, const number subDetVal=NULL)
Definition: mpr_base.cc:3059
uResultant(const ideal _gls, const resMatType _rmt=sparseResMat, BOOLEAN extIdeal=true)
Definition: mpr_base.cc:2684
resMatrixBase * resMat
Definition: mpr_base.h:92
ideal extendIdeal(const ideal gls, poly linPoly, const resMatType rmt)
Definition: mpr_base.cc:2714
rootContainer ** interpolateDenseSP(BOOLEAN matchUp=false, const number subDetVal=NULL)
Definition: mpr_base.cc:2921
@ denseResMat
Definition: mpr_base.h:65
@ sparseResMat
Definition: mpr_base.h:65
ideal gls
Definition: mpr_base.h:88
poly interpolateDense(const number subDetVal=NULL)
Definition: mpr_base.cc:2769
int nextPrime(const int p)
Definition: mpr_base.cc:3172
vandermonde system solver for interpolating polynomials from their values
Definition: mpr_numeric.h:29
number * interpolateDense(const number *q)
Solves the Vandermode linear system \sum_{i=1}^{n} x_i^k-1 w_i = q_k, k=1,..,n.
Definition: mpr_numeric.cc:146
static FORCE_INLINE number n_Param(const int iParameter, const coeffs r)
return the (iParameter^th) parameter as a NEW number NOTE: parameter numbering: 1....
Definition: coeffs.h:783
static FORCE_INLINE long n_Int(number &n, const coeffs r)
conversion of n to an int; 0 if not possible in Z/pZ: the representing int lying in (-p/2 ....
Definition: coeffs.h:547
#define Print
Definition: emacs.cc:80
CFFListIterator iter
Definition: facAbsBiFact.cc:53
return result
Definition: facAbsBiFact.cc:75
CanonicalForm res
Definition: facAbsFact.cc:60
REvaluation E(1, terms.length(), IntRandom(25))
CanonicalForm factor
Definition: facAbsFact.cc:97
bool found
Definition: facFactorize.cc:55
long isReduced(const mat_zz_p &M)
Definition: facFqBivar.cc:1468
int j
Definition: facHensel.cc:110
static int max(int a, int b)
Definition: fast_mult.cc:264
void WerrorS(const char *s)
Definition: feFopen.cc:24
#define idDelete(H)
delete an ideal
Definition: ideals.h:29
ideal idCopy(ideal A)
Definition: ideals.h:60
#define IMATELEM(M, I, J)
Definition: intvec.h:85
#define pi
Definition: libparse.cc:1145
#define phelp
Definition: libparse.cc:1242
ideal lift(const ideal J, const ring r, const ideal inI, const ring s)
Definition: lift.cc:26
matrix mpNew(int r, int c)
create a r x c zero-matrix
Definition: matpol.cc:37
#define MATELEM(mat, i, j)
1-based access to matrix
Definition: matpol.h:29
#define MATROWS(i)
Definition: matpol.h:26
#define MATCOLS(i)
Definition: matpol.h:27
#define assume(x)
Definition: mod2.h:389
#define pIter(p)
Definition: monomials.h:37
#define pNext(p)
Definition: monomials.h:36
static number & pGetCoeff(poly p)
return an alias to the leading coefficient of p assumes that p != NULL NOTE: not copy
Definition: monomials.h:44
#define pSetCoeff0(p, n)
Definition: monomials.h:59
#define MAXVARS
Definition: mpr_base.cc:55
unsigned int Coord_t
Definition: mpr_base.cc:131
poly monomAt(poly p, int i)
Definition: mpr_base.cc:720
unsigned long over(const unsigned long n, const unsigned long d)
Definition: mpr_base.cc:2658
int col
Definition: mpr_base.cc:152
setID rc
Definition: mpr_base.cc:142
#define MAXEVPOINT
Definition: mpr_base.cc:2652
#define RVMULT
Definition: mpr_base.cc:53
#define SCALEDOWN
Definition: mpr_base.cc:51
#define MAXRVVAL
Definition: mpr_base.cc:54
#define MINVDIST
Definition: mpr_base.cc:52
ideal loNewtonPolytope(const ideal id)
Definition: mpr_base.cc:3190
number num
Definition: mpr_base.cc:151
#define MAXINITELEMS
Definition: mpr_base.cc:49
Coord_t * point
Definition: mpr_base.cc:141
struct _entry * next
Definition: mpr_base.cc:153
int set
Definition: mpr_base.cc:135
struct onePoint * rcPnt
Definition: mpr_base.cc:143
int pnt
Definition: mpr_base.cc:136
#define LIFT_COOR
Definition: mpr_base.cc:50
Definition: mpr_base.cc:150
#define SNONE
Definition: mpr_base.h:14
#define SFREE
Definition: mpr_base.h:15
gmp_float exp(const gmp_float &a)
Definition: mpr_complex.cc:357
#define mprSTICKYPROT(msg)
Definition: mpr_global.h:54
#define ST_SPARSE_MREC2
Definition: mpr_global.h:74
#define mprSTICKYPROT2(msg, arg)
Definition: mpr_global.h:55
#define ST_SPARSE_VADD
Definition: mpr_global.h:70
#define ST_SPARSE_RCRJ
Definition: mpr_global.h:76
#define ST_SPARSE_RC
Definition: mpr_global.h:75
#define ST_DENSE_MEM
Definition: mpr_global.h:66
#define ST__DET
Definition: mpr_global.h:78
#define ST_SPARSE_VREJ
Definition: mpr_global.h:71
#define mprPROTL(msg, intval)
Definition: mpr_global.h:46
#define ST_DENSE_NR
Definition: mpr_global.h:65
#define mprPROTN(msg, nval)
Definition: mpr_global.h:48
#define mprPROT(msg)
Definition: mpr_global.h:41
#define ST_SPARSE_MPEND
Definition: mpr_global.h:72
#define ST_BASE_EV
Definition: mpr_global.h:62
#define mprPROTnl(msg)
Definition: mpr_global.h:42
#define ST_DENSE_FR
Definition: mpr_global.h:64
#define ST_SPARSE_MREC1
Definition: mpr_global.h:73
#define mprPROTNnl(msg, nval)
Definition: mpr_global.h:49
#define ST_DENSE_NMON
Definition: mpr_global.h:67
double mprfloat
Definition: mpr_global.h:17
#define ST_SPARSE_MEM
Definition: mpr_global.h:69
#define SIMPLEX_EPS
Definition: mpr_numeric.h:181
void init()
Definition: lintree.cc:864
#define nDiv(a, b)
Definition: numbers.h:32
#define nDelete(n)
Definition: numbers.h:16
#define nIsZero(n)
Definition: numbers.h:19
#define nEqual(n1, n2)
Definition: numbers.h:20
#define nCopy(n)
Definition: numbers.h:15
#define nPrint(a)
only for debug, over any initalized currRing
Definition: numbers.h:46
#define nNormalize(n)
Definition: numbers.h:30
#define nInit(i)
Definition: numbers.h:24
#define nTest(a)
Definition: numbers.h:35
#define nPower(a, b, res)
Definition: numbers.h:38
#define omFreeSize(addr, size)
Definition: omAllocDecl.h:260
#define omAlloc(size)
Definition: omAllocDecl.h:210
#define omReallocSize(addr, o_size, size)
Definition: omAllocDecl.h:220
#define omAlloc0(size)
Definition: omAllocDecl.h:211
#define omfreeSize(addr, size)
Definition: omAllocDecl.h:236
#define NULL
Definition: omList.c:12
static int index(p_Length length, p_Ord ord)
Definition: p_Procs_Impl.h:592
poly pp_DivideM(poly a, poly b, const ring r)
Definition: p_polys.cc:1629
static unsigned pLength(poly a)
Definition: p_polys.h:191
static void p_GetExpV(poly p, int *ev, const ring r)
Definition: p_polys.h:1522
VAR ring currRing
Widely used global variable which specifies the current polynomial ring for Singular interpreter and ...
Definition: polys.cc:13
Compatiblity layer for legacy polynomial operations (over currRing)
#define pAdd(p, q)
Definition: polys.h:203
static long pTotaldegree(poly p)
Definition: polys.h:282
#define pTest(p)
Definition: polys.h:415
#define pDelete(p_ptr)
Definition: polys.h:186
#define pHead(p)
returns newly allocated copy of Lm(p), coef is copied, next=NULL, p might be NULL
Definition: polys.h:67
#define pSetm(p)
Definition: polys.h:271
#define pLmEqual(p1, p2)
Definition: polys.h:111
#define pSetCoeff(p, n)
deletes old coeff before setting the new one
Definition: polys.h:31
#define ppMult_qq(p, q)
Definition: polys.h:208
#define pSetExpV(p, e)
Definition: polys.h:97
#define pLmInit(p)
like pInit, except that expvector is initialized to that of p, p must be != NULL
Definition: polys.h:64
#define pSetComp(p, v)
Definition: polys.h:38
void pWrite0(poly p)
Definition: polys.h:309
#define pIncrExp(p, i)
Definition: polys.h:43
#define pMult(p, q)
Definition: polys.h:207
void pWrite(poly p)
Definition: polys.h:308
#define pGetExp(p, i)
Exponent.
Definition: polys.h:41
#define pSetmComp(p)
TODO:
Definition: polys.h:273
#define pInit()
allocates a new monomial and initializes everything to 0
Definition: polys.h:61
#define pSetExp(p, i, v)
Definition: polys.h:42
#define pCopy(p)
return a copy of the poly
Definition: polys.h:185
#define pOne()
Definition: polys.h:315
#define pLmDivisibleByNoComp(a, b)
like pLmDivisibleBy, does not check components
Definition: polys.h:142
int IsPrime(int p)
Definition: prime.cc:61
void PrintS(const char *s)
Definition: reporter.cc:284
void PrintLn()
Definition: reporter.cc:310
void Werror(const char *fmt,...)
Definition: reporter.cc:189
int status int void size_t count
Definition: si_signals.h:59
ideal idInit(int idsize, int rank)
initialise an ideal / module
Definition: simpleideals.cc:35
ideal id_Matrix2Module(matrix mat, const ring R)
converts mat to module, destroys mat
#define IDELEMS(i)
Definition: simpleideals.h:23
int siRand()
Definition: sirandom.c:42
poly sm_CallDet(ideal I, const ring R)
Definition: sparsmat.cc:302
poly getElem(const int i)
index von 0 ...
Definition: mpr_base.cc:2046
poly dividedBy
Definition: mpr_base.cc:2024
poly mon
Definition: mpr_base.cc:2023
void init()
Definition: mpr_base.cc:2003
int numColVectorSize
size of numColVector
Definition: mpr_base.cc:2039
int elementOfS
number of the set S mon is element of
Definition: mpr_base.cc:2028
bool isReduced
Definition: mpr_base.cc:2025
void init(const poly m)
Definition: mpr_base.cc:2009
number * numColVecCopy
Definition: mpr_base.cc:2041
number getElemNum(const int i)
index von 0 ...
Definition: mpr_base.cc:2055
int * numColParNr
holds the index of u0, u1, ..., un, if (elementOfS == linPolyS) the size is given by (currRing->N)
Definition: mpr_base.cc:2033
number * numColVector
holds the column vector if (elementOfS != linPolyS)
Definition: mpr_base.cc:2036
int dim(ideal I, ring r)