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 <omalloc/omalloc.h>
17 
18 #include <misc/mylimits.h>
19 #include <misc/options.h>
20 #include <misc/intvec.h>
21 #include <misc/sirandom.h>
22 
23 #include <coeffs/numbers.h>
24 #include <coeffs/mpr_global.h>
25 
26 #include <polys/matpol.h>
27 #include <polys/sparsmat.h>
28 
29 #include <polys/clapsing.h>
30 
31 #include <kernel/polys.h>
32 #include <kernel/ideals.h>
33 
34 #include "mpr_base.h"
35 #include "mpr_numeric.h"
36 
37 #include <math.h>
38 //<-
39 
40 //%s
41 //-----------------------------------------------------------------------------
42 //-------------- sparse resultant matrix --------------------------------------
43 //-----------------------------------------------------------------------------
44 
45 //-> definitions
46 
47 //#define mprTEST
48 //#define mprMINKSUM
49 
50 #define MAXPOINTS 10000
51 #define MAXINITELEMS 256
52 #define LIFT_COOR 50000 // siRand() % LIFT_COOR gives random lift value
53 #define SCALEDOWN 100.0 // lift value scale down for linear program
54 #define MINVDIST 0.0
55 #define RVMULT 0.0001 // multiplicator for random shift vector
56 #define MAXRVVAL 50000
57 #define MAXVARS 100
58 //<-
59 
60 //-> sparse resultant matrix
61 
62 /* set of points */
63 class pointSet;
64 
65 
66 
67 /* sparse resultant matrix class */
68 class resMatrixSparse : virtual public resMatrixBase
69 {
70 public:
71  resMatrixSparse( const ideal _gls, const int special = SNONE );
73 
74  // public interface according to base class resMatrixBase
75  ideal getMatrix();
76 
77  /** Fills in resMat[][] with evpoint[] and gets determinant
78  * uRPos[i][1]: row of matrix
79  * uRPos[i][idelem+1]: col of u(0)
80  * uRPos[i][2..idelem]: col of u(1) .. u(n)
81  * i= 1 .. numSet0
82  */
83  number getDetAt( const number* evpoint );
84 
85  poly getUDet( const number* evpoint );
86 
87 private:
89 
90  void randomVector( const int dim, mprfloat shift[] );
91 
92  /** Row Content Function
93  * Finds the largest i such that F[i] is a point, F[i]= a[ij] in A[i] for some j.
94  * Returns -1 iff the point vert does not lie in a cell
95  */
96  int RC( pointSet **pQ, pointSet *E, int vert, mprfloat shift[] );
97 
98  /* Remaps a result of LP to the according point set Qi.
99  * Returns false iff remaping was not possible, otherwise true.
100  */
101  bool remapXiToPoint( const int indx, pointSet **pQ, int *set, int *vtx );
102 
103  /** create coeff matrix
104  * uRPos[i][1]: row of matrix
105  * uRPos[i][idelem+1]: col of u(0)
106  * uRPos[i][2..idelem]: col of u(1) .. u(n)
107  * i= 1 .. numSet0
108  * Returns the dimension of the matrix or -1 in case of an error
109  */
110  int createMatrix( pointSet *E );
111 
112  pointSet * minkSumAll( pointSet **pQ, int numq, int dim );
113  pointSet * minkSumTwo( pointSet *Q1, pointSet *Q2, int dim );
114 
115 private:
116  ideal gls;
117 
118  int n, idelem; // number of variables, polynoms
119  int numSet0; // number of elements in S0
120  int msize; // size of matrix
121 
123 
124  ideal rmat; // sparse matrix representation
125 
126  simplex * LP; // linear programming stuff
127 };
128 //<-
129 
130 //-> typedefs and structs
131 poly monomAt( poly p, int i );
132 
133 typedef unsigned int Coord_t;
134 
135 struct setID
136 {
137  int set;
138  int pnt;
139 };
140 
141 struct onePoint
142 {
143  Coord_t * point; // point[0] is unused, maxial dimension is MAXVARS+1
144  setID rc; // filled in by Row Content Function
145  struct onePoint * rcPnt; // filled in by Row Content Function
146 };
147 
148 typedef struct onePoint * onePointP;
149 
150 /* sparse matrix entry */
151 struct _entry
152 {
153  number num;
154  int col;
155  struct _entry * next;
156 };
157 
158 typedef struct _entry * entry;
159 //<-
160 
161 //-> class pointSet
162 class pointSet
163 {
164 private:
165  onePointP *points; // set of onePoint's, index [1..num], supports of monoms
166  bool lifted;
167 
168 public:
169  int num; // number of elements in points
170  int max; // maximal entries in points, i.e. allocated mem
171  int dim; // dimension, i.e. valid coord entries in point
172  int index; // should hold unique identifier of point set
173 
174  pointSet( const int _dim, const int _index= 0, const int count= MAXINITELEMS );
175  ~pointSet();
176 
177  // pointSet.points[i] equals pointSet[i]
178  inline onePointP operator[] ( const int index );
179 
180  /** Adds a point to pointSet, copy vert[0,...,dim] ot point[num+1][0,...,dim].
181  * Returns false, iff additional memory was allocated ( i.e. num >= max )
182  * else returns true
183  */
184  bool addPoint( const onePointP vert );
185 
186  /** Adds a point to pointSet, copy vert[0,...,dim] ot point[num+1][0,...,dim].
187  * Returns false, iff additional memory was allocated ( i.e. num >= max )
188  * else returns true
189  */
190  bool addPoint( const int * vert );
191 
192  /** Adds a point to pointSet, copy vert[0,...,dim] ot point[num+1][0,...,dim].
193  * Returns false, iff additional memory was allocated ( i.e. num >= max )
194  * else returns true
195  */
196  bool addPoint( const Coord_t * vert );
197 
198  /* Removes the point at intex indx */
199  bool removePoint( const int indx );
200 
201  /** Adds point to pointSet, iff pointSet \cap point = \emptyset.
202  * Returns true, iff added, else false.
203  */
204  bool mergeWithExp( const onePointP vert );
205 
206  /** Adds point to pointSet, iff pointSet \cap point = \emptyset.
207  * Returns true, iff added, else false.
208  */
209  bool mergeWithExp( const int * vert );
210 
211  /* Adds support of poly p to pointSet, iff pointSet \cap point = \emptyset. */
212  void mergeWithPoly( const poly p );
213 
214  /* Returns the row polynom multiplicator in vert[] */
215  void getRowMP( const int indx, int * vert );
216 
217  /* Returns index of supp(LT(p)) in pointSet. */
218  int getExpPos( const poly p );
219 
220  /** sort lex
221  */
222  void sort();
223 
224  /** Lifts the point set using sufficiently generic linear lifting
225  * homogeneous forms l[1]..l[dim] in Z. Every l[i] is of the form
226  * L1x1+...+Lnxn, for generic L1..Ln in Z.
227  *
228  * Lifting raises dimension by one!
229  */
230  void lift( int *l= NULL ); // !! increments dim by 1
231  void unlift() { dim--; lifted= false; }
232 
233 private:
234  pointSet( const pointSet & );
235 
236  /** points[a] < points[b] ? */
237  inline bool smaller( int, int );
238 
239  /** points[a] > points[b] ? */
240  inline bool larger( int, int );
241 
242  /** Checks, if more mem is needed ( i.e. num >= max ),
243  * returns false, if more mem was allocated, else true
244  */
245  inline bool checkMem();
246 };
247 //<-
248 
249 //-> class convexHull
250 /* Compute convex hull of given exponent set */
252 {
253 public:
254  convexHull( simplex * _pLP ) : pLP(_pLP) {}
256 
257  /** Computes the point sets of the convex hulls of the supports given
258  * by the polynoms in gls.
259  * Returns Q[].
260  */
261  pointSet ** newtonPolytopesP( const ideal gls );
262  ideal newtonPolytopesI( const ideal gls );
263 
264 private:
265  /** Returns true iff the support of poly pointPoly is inside the
266  * convex hull of all points given by the support of poly p.
267  */
268  bool inHull(poly p, poly pointPoly, int m, int site);
269 
270 private:
272  int n;
274 };
275 //<-
276 
277 //-> class mayanPyramidAlg
278 /* Compute all lattice points in a given convex hull */
280 {
281 public:
282  mayanPyramidAlg( simplex * _pLP ) : n((currRing->N)), pLP(_pLP) {}
284 
285  /** Drive Mayan Pyramid Algorithm.
286  * The Alg computes conv(Qi[]+shift[]).
287  */
288  pointSet * getInnerPoints( pointSet **_q_i, mprfloat _shift[] );
289 
290 private:
291 
292  /** Recursive Mayan Pyramid algorithm for directly computing MinkowskiSum
293  * lattice points for (n+1)-fold MinkowskiSum of given point sets Qi[].
294  * Recursively for range of dim: dim in [0..n); acoords[0..var) fixed.
295  * Stores only MinkowskiSum points of udist > 0: done by storeMinkowskiSumPoints.
296  */
297  void runMayanPyramid( int dim );
298 
299  /** Compute v-distance via Linear Programing
300  * Linear Program finds the v-distance of the point in accords[].
301  * The v-distance is the distance along the direction v to boundary of
302  * Minkowski Sum of Qi (here vector v is represented by shift[]).
303  * Returns the v-distance or -1.0 if an error occurred.
304  */
305  mprfloat vDistance( Coord_t * acoords, int dim );
306 
307  /** LP for finding min/max coord in MinkowskiSum, given previous coors.
308  * Assume MinkowskiSum in non-negative quadrants
309  * coor in [0,n); fixed coords in acoords[0..coor)
310  */
311  void mn_mx_MinkowskiSum( int dim, Coord_t *minR, Coord_t *maxR );
312 
313  /** Stores point in E->points[pt], iff v-distance != 0
314  * Returns true iff point was stored, else flase
315  */
316  bool storeMinkowskiSumPoint();
317 
318 private:
322 
323  int n,idelem;
324 
325  Coord_t acoords[MAXVARS+2];
326 
328 };
329 //<-
330 
331 //-> debug output stuff
332 #if defined(mprDEBUG_PROT) || defined(mprDEBUG_ALL)
333 void print_mat(mprfloat **a, int maxrow, int maxcol)
334 {
335  int i, j;
336 
337  for (i = 1; i <= maxrow; i++)
338  {
339  PrintS("[");
340  for (j = 1; j <= maxcol; j++) Print("% 7.2f, ", a[i][j]);
341  PrintS("],\n");
342  }
343 }
344 void print_bmat(mprfloat **a, int nrows, int ncols, int N, int *iposv)
345 {
346  int i, j;
347 
348  printf("Output matrix from LinProg");
349  for (i = 1; i <= nrows; i++)
350  {
351  printf("\n[ ");
352  if (i == 1) printf(" ");
353  else if (iposv[i-1] <= N) printf("X%d", iposv[i-1]);
354  else printf("Y%d", iposv[i-1]-N+1);
355  for (j = 1; j <= ncols; j++) printf(" %7.2f ",(double)a[i][j]);
356  printf(" ]");
357  } printf("\n");
358  fflush(stdout);
359 }
360 
361 void print_exp( const onePointP vert, int n )
362 {
363  int i;
364  for ( i= 1; i <= n; i++ )
365  {
366  Print(" %d",vert->point[i] );
367 #ifdef LONG_OUTPUT
368  if ( i < n ) PrintS(", ");
369 #endif
370  }
371 }
372 void print_matrix( matrix omat )
373 {
374  int i,j;
375  int val;
376  Print(" matrix m[%d][%d]=(\n",MATROWS( omat ),MATCOLS( omat ));
377  for ( i= 1; i <= MATROWS( omat ); i++ )
378  {
379  for ( j= 1; j <= MATCOLS( omat ); j++ )
380  {
381  if ( (MATELEM( omat, i, j)!=NULL)
382  && (!nIsZero(pGetCoeff( MATELEM( omat, i, j)))))
383  {
384  val= n_Int(pGetCoeff( MATELEM( omat, i, j) ), currRing->cf);
385  if ( i==MATROWS(omat) && j==MATCOLS(omat) )
386  {
387  Print("%d ",val);
388  }
389  else
390  {
391  Print("%d, ",val);
392  }
393  }
394  else
395  {
396  if ( i==MATROWS(omat) && j==MATCOLS(omat) )
397  {
398  PrintS(" 0");
399  }
400  else
401  {
402  PrintS(" 0, ");
403  }
404  }
405  }
406  PrintLn();
407  }
408  PrintS(");\n");
409 }
410 #endif
411 //<-
412 
413 //-> pointSet::*
414 pointSet::pointSet( const int _dim, const int _index, const int count )
415  : num(0), max(count), dim(_dim), index(_index)
416 {
417  int i;
418  points = (onePointP *)omAlloc( (count+1) * sizeof(onePointP) );
419  for ( i= 0; i <= max; i++ )
420  {
421  points[i]= (onePointP)omAlloc( sizeof(onePoint) );
422  points[i]->point= (Coord_t *)omAlloc0( (dim+2) * sizeof(Coord_t) );
423  }
424  lifted= false;
425 }
426 
428 {
429  int i;
430  int fdim= lifted ? dim+1 : dim+2;
431  for ( i= 0; i <= max; i++ )
432  {
433  omFreeSize( (void *) points[i]->point, fdim * sizeof(Coord_t) );
434  omFreeSize( (void *) points[i], sizeof(onePoint) );
435  }
436  omFreeSize( (void *) points, (max+1) * sizeof(onePointP) );
437 }
438 
439 inline onePointP pointSet::operator[] ( const int index_i )
440 {
441  assume( index_i > 0 && index_i <= num );
442  return points[index_i];
443 }
444 
445 inline bool pointSet::checkMem()
446 {
447  if ( num >= max )
448  {
449  int i;
450  int fdim= lifted ? dim+1 : dim+2;
451  points= (onePointP*)omReallocSize( points,
452  (max+1) * sizeof(onePointP),
453  (2*max + 1) * sizeof(onePointP) );
454  for ( i= max+1; i <= max*2; i++ )
455  {
456  points[i]= (onePointP)omAlloc( sizeof(struct onePoint) );
457  points[i]->point= (Coord_t *)omAlloc0( fdim * sizeof(Coord_t) );
458  }
459  max*= 2;
461  return false;
462  }
463  return true;
464 }
465 
466 bool pointSet::addPoint( const onePointP vert )
467 {
468  int i;
469  bool ret;
470  num++;
471  ret= checkMem();
472  points[num]->rcPnt= NULL;
473  for ( i= 1; i <= dim; i++ ) points[num]->point[i]= vert->point[i];
474  return ret;
475 }
476 
477 bool pointSet::addPoint( const int * vert )
478 {
479  int i;
480  bool ret;
481  num++;
482  ret= checkMem();
483  points[num]->rcPnt= NULL;
484  for ( i= 1; i <= dim; i++ ) points[num]->point[i]= (Coord_t) vert[i];
485  return ret;
486 }
487 
488 bool pointSet::addPoint( const Coord_t * vert )
489 {
490  int i;
491  bool ret;
492  num++;
493  ret= checkMem();
494  points[num]->rcPnt= NULL;
495  for ( i= 0; i < dim; i++ ) points[num]->point[i+1]= vert[i];
496  return ret;
497 }
498 
499 bool pointSet::removePoint( const int indx )
500 {
501  assume( indx > 0 && indx <= num );
502  if ( indx != num )
503  {
504  onePointP tmp;
505  tmp= points[indx];
506  points[indx]= points[num];
507  points[num]= tmp;
508  }
509  num--;
510 
511  return true;
512 }
513 
514 bool pointSet::mergeWithExp( const onePointP vert )
515 {
516  int i,j;
517 
518  for ( i= 1; i <= num; i++ )
519  {
520  for ( j= 1; j <= dim; j++ )
521  if ( points[i]->point[j] != vert->point[j] ) break;
522  if ( j > dim ) break;
523  }
524 
525  if ( i > num )
526  {
527  addPoint( vert );
528  return true;
529  }
530  return false;
531 }
532 
533 bool pointSet::mergeWithExp( const int * vert )
534 {
535  int i,j;
536 
537  for ( i= 1; i <= num; i++ )
538  {
539  for ( j= 1; j <= dim; j++ )
540  if ( points[i]->point[j] != (Coord_t) vert[j] ) break;
541  if ( j > dim ) break;
542  }
543 
544  if ( i > num )
545  {
546  addPoint( vert );
547  return true;
548  }
549  return false;
550 }
551 
553 {
554  int i,j;
555  poly piter= p;
556  int * vert;
557  vert= (int *)omAlloc( (dim+1) * sizeof(int) );
558 
559  while ( piter )
560  {
561  pGetExpV( piter, vert );
562 
563  for ( i= 1; i <= num; i++ )
564  {
565  for ( j= 1; j <= dim; j++ )
566  if ( points[i]->point[j] != (Coord_t) vert[j] ) break;
567  if ( j > dim ) break;
568  }
569 
570  if ( i > num )
571  {
572  addPoint( vert );
573  }
574 
575  pIter( piter );
576  }
577  omFreeSize( (void *) vert, (dim+1) * sizeof(int) );
578 }
579 
581 {
582  int * vert;
583  int i,j;
584 
585  // hier unschoen...
586  vert= (int *)omAlloc( (dim+1) * sizeof(int) );
587 
588  pGetExpV( p, vert );
589  for ( i= 1; i <= num; i++ )
590  {
591  for ( j= 1; j <= dim; j++ )
592  if ( points[i]->point[j] != (Coord_t) vert[j] ) break;
593  if ( j > dim ) break;
594  }
595  omFreeSize( (void *) vert, (dim+1) * sizeof(int) );
596 
597  if ( i > num ) return 0;
598  else return i;
599 }
600 
601 void pointSet::getRowMP( const int indx, int * vert )
602 {
603  assume( indx > 0 && indx <= num && points[indx]->rcPnt );
604  int i;
605 
606  vert[0]= 0;
607  for ( i= 1; i <= dim; i++ )
608  vert[i]= (int)(points[indx]->point[i] - points[indx]->rcPnt->point[i]);
609 }
610 
611 inline bool pointSet::smaller( int a, int b )
612 {
613  int i;
614 
615  for ( i= 1; i <= dim; i++ )
616  {
617  if ( points[a]->point[i] > points[b]->point[i] )
618  {
619  return false;
620  }
621  if ( points[a]->point[i] < points[b]->point[i] )
622  {
623  return true;
624  }
625  }
626 
627  return false; // they are equal
628 }
629 
630 inline bool pointSet::larger( int a, int b )
631 {
632  int i;
633 
634  for ( i= 1; i <= dim; i++ )
635  {
636  if ( points[a]->point[i] < points[b]->point[i] )
637  {
638  return false;
639  }
640  if ( points[a]->point[i] > points[b]->point[i] )
641  {
642  return true;
643  }
644  }
645 
646  return false; // they are equal
647 }
648 
650 {
651  int i;
652  bool found= true;
653  onePointP tmp;
654 
655  while ( found )
656  {
657  found= false;
658  for ( i= 1; i < num; i++ )
659  {
660  if ( larger( i, i+1 ) )
661  {
662  tmp= points[i];
663  points[i]= points[i+1];
664  points[i+1]= tmp;
665 
666  found= true;
667  }
668  }
669  }
670 }
671 
672 void pointSet::lift( int l[] )
673 {
674  bool outerL= true;
675  int i, j;
676  int sum;
677 
678  dim++;
679 
680  if ( l==NULL )
681  {
682  outerL= false;
683  l= (int *)omAlloc( (dim+1) * sizeof(int) ); // [1..dim-1]
684 
685  for(i = 1; i < dim; i++)
686  {
687  l[i]= 1 + siRand() % LIFT_COOR;
688  }
689  }
690  for ( j=1; j <= num; j++ )
691  {
692  sum= 0;
693  for ( i=1; i < dim; i++ )
694  {
695  sum += (int)points[j]->point[i] * l[i];
696  }
697  points[j]->point[dim]= sum;
698  }
699 
700 #ifdef mprDEBUG_ALL
701  PrintS(" lift vector: ");
702  for ( j=1; j < dim; j++ ) Print(" %d ",l[j] );
703  PrintLn();
704 #ifdef mprDEBUG_ALL
705  PrintS(" lifted points: \n");
706  for ( j=1; j <= num; j++ )
707  {
708  Print("%d: <",j);print_exp(points[j],dim);PrintS(">\n");
709  }
710  PrintLn();
711 #endif
712 #endif
713 
714  lifted= true;
715 
716  if ( !outerL ) omFreeSize( (void *) l, (dim+1) * sizeof(int) );
717 }
718 //<-
719 
720 //-> global functions
721 // Returns the monom at pos i in poly p
722 poly monomAt( poly p, int i )
723 {
724  assume( i > 0 );
725  poly iter= p;
726  for ( int j= 1; (j < i) && (iter!=NULL); j++ ) pIter(iter);
727  return iter;
728 }
729 //<-
730 
731 //-> convexHull::*
732 bool convexHull::inHull(poly p, poly pointPoly, int m, int site)
733 {
734  int i, j, col;
735 
736  pLP->m = n+1;
737  pLP->n = m; // this includes col of cts
738 
739  pLP->LiPM[1][1] = +0.0;
740  pLP->LiPM[1][2] = +1.0; // optimize (arbitrary) var
741  pLP->LiPM[2][1] = +1.0;
742  pLP->LiPM[2][2] = -1.0; // lambda vars sum up to 1
743 
744  for ( j=3; j <= pLP->n; j++)
745  {
746  pLP->LiPM[1][j] = +0.0;
747  pLP->LiPM[2][j] = -1.0;
748  }
749 
750  for( i= 1; i <= n; i++) { // each row constraints one coor
751  pLP->LiPM[i+2][1] = (mprfloat)pGetExp(pointPoly,i);
752  col = 2;
753  for( j= 1; j <= m; j++ )
754  {
755  if( j != site )
756  {
757  pLP->LiPM[i+2][col] = -(mprfloat)pGetExp( monomAt(p,j), i );
758  col++;
759  }
760  }
761  }
762 
763 #ifdef mprDEBUG_ALL
764  PrintS("Matrix of Linear Programming\n");
765  print_mat( pLP->LiPM, pLP->m+1,pLP->n);
766 #endif
767 
768  pLP->m3= pLP->m;
769 
770  pLP->compute();
771 
772  return (pLP->icase == 0);
773 }
774 
775 // mprSTICKYPROT:
776 // ST_SPARSE_VADD: new vertex of convex hull added
777 // ST_SPARSE_VREJ: point rejected (-> inside hull)
779 {
780  int i, j, k;
781  int m; // Anzahl der Exponentvektoren im i-ten Polynom (gls->m)[i] des Ideals gls
782  int idelem= IDELEMS(gls);
783  int * vert;
784 
785  n= (currRing->N);
786  vert= (int *)omAlloc( (idelem+1) * sizeof(int) );
787 
788  Q = (pointSet **)omAlloc( idelem * sizeof(pointSet*) ); // support hulls
789  for ( i= 0; i < idelem; i++ )
790  Q[i] = new pointSet( (currRing->N), i+1, pLength((gls->m)[i])+1 );
791 
792  for( i= 0; i < idelem; i++ )
793  {
794  k=1;
795  m = pLength( (gls->m)[i] );
796 
797  poly p= (gls->m)[i];
798  for( j= 1; j <= m; j++) { // für jeden Exponentvektor
799  if( !inHull( (gls->m)[i], p, m, j ) )
800  {
801  pGetExpV( p, vert );
802  Q[i]->addPoint( vert );
803  k++;
805  }
806  else
807  {
809  }
810  pIter( p );
811  } // j
812  mprSTICKYPROT("\n");
813  } // i
814 
815  omFreeSize( (void *) vert, (idelem+1) * sizeof(int) );
816 
817 #ifdef mprDEBUG_PROT
818  PrintLn();
819  for( i= 0; i < idelem; i++ )
820  {
821  Print(" \\Conv(Qi[%d]): #%d\n", i,Q[i]->num );
822  for ( j=1; j <= Q[i]->num; j++ )
823  {
824  Print("%d: <",j);print_exp( (*Q[i])[j] , (currRing->N) );PrintS(">\n");
825  }
826  PrintLn();
827  }
828 #endif
829 
830  return Q;
831 }
832 
833 // mprSTICKYPROT:
834 // ST_SPARSE_VADD: new vertex of convex hull added
835 // ST_SPARSE_VREJ: point rejected (-> inside hull)
836 ideal convexHull::newtonPolytopesI( const ideal gls )
837 {
838  int i, j;
839  int m; // Anzahl der Exponentvektoren im i-ten Polynom (gls->m)[i] des Ideals gls
840  int idelem= IDELEMS(gls);
841  ideal id;
842  poly p,pid;
843  int * vert;
844 
845  n= (currRing->N);
846  vert= (int *)omAlloc( (idelem+1) * sizeof(int) );
847  id= idInit( idelem, 1 );
848 
849  for( i= 0; i < idelem; i++ )
850  {
851  m = pLength( (gls->m)[i] );
852 
853  p= (gls->m)[i];
854  for( j= 1; j <= m; j++) { // für jeden Exponentvektor
855  if( !inHull( (gls->m)[i], p, m, j ) )
856  {
857  if ( (id->m)[i] == NULL )
858  {
859  (id->m)[i]= pHead(p);
860  pid=(id->m)[i];
861  }
862  else
863  {
864  pNext(pid)= pHead(p);
865  pIter(pid);
866  pNext(pid)= NULL;
867  }
869  }
870  else
871  {
873  }
874  pIter( p );
875  } // j
876  mprSTICKYPROT("\n");
877  } // i
878 
879  omFreeSize( (void *) vert, (idelem+1) * sizeof(int) );
880 
881 #ifdef mprDEBUG_PROT
882  PrintLn();
883  for( i= 0; i < idelem; i++ )
884  {
885  }
886 #endif
887 
888  return id;
889 }
890 //<-
891 
892 //-> mayanPyramidAlg::*
894 {
895  int i;
896 
897  Qi= _q_i;
898  shift= _shift;
899 
900  E= new pointSet( Qi[0]->dim ); // E has same dim as Qi[...]
901 
902  for ( i= 0; i < MAXVARS+2; i++ ) acoords[i]= 0;
903 
904  runMayanPyramid(0);
905 
906  mprSTICKYPROT("\n");
907 
908  return E;
909 }
910 
912 {
913  int i, ii, j, k, col, r;
914  int numverts, cols;
915 
916  numverts = 0;
917  for( i=0; i<=n; i++)
918  {
919  numverts += Qi[i]->num;
920  }
921  cols = numverts + 2;
922 
923  //if( dim < 1 || dim > n )
924  // WerrorS("mayanPyramidAlg::vDistance: Known coords dim off range");
925 
926  pLP->LiPM[1][1] = 0.0;
927  pLP->LiPM[1][2] = 1.0; // maximize
928  for( j=3; j<=cols; j++) pLP->LiPM[1][j] = 0.0;
929 
930  for( i=0; i <= n; i++ )
931  {
932  pLP->LiPM[i+2][1] = 1.0;
933  pLP->LiPM[i+2][2] = 0.0;
934  }
935  for( i=1; i<=dim; i++)
936  {
937  pLP->LiPM[n+2+i][1] = (mprfloat)(acoords_a[i-1]);
938  pLP->LiPM[n+2+i][2] = -shift[i];
939  }
940 
941  ii = -1;
942  col = 2;
943  for ( i= 0; i <= n; i++ )
944  {
945  ii++;
946  for( k= 1; k <= Qi[ii]->num; k++ )
947  {
948  col++;
949  for ( r= 0; r <= n; r++ )
950  {
951  if ( r == i ) pLP->LiPM[r+2][col] = -1.0;
952  else pLP->LiPM[r+2][col] = 0.0;
953  }
954  for( r= 1; r <= dim; r++ )
955  pLP->LiPM[r+n+2][col] = -(mprfloat)((*Qi[ii])[k]->point[r]);
956  }
957  }
958 
959  if( col != cols)
960  Werror("mayanPyramidAlg::vDistance:"
961  "setting up matrix for udist: col %d != cols %d",col,cols);
962 
963  pLP->m = n+dim+1;
964  pLP->m3= pLP->m;
965  pLP->n=cols-1;
966 
967 #ifdef mprDEBUG_ALL
968  Print("vDistance LP, known koords dim=%d, constr %d, cols %d, acoords= ",
969  dim,pLP->m,cols);
970  for( i= 0; i < dim; i++ )
971  Print(" %d",acoords_a[i]);
972  PrintLn();
973  print_mat( pLP->LiPM, pLP->m+1, cols);
974 #endif
975 
976  pLP->compute();
977 
978 #ifdef mprDEBUG_ALL
979  PrintS("LP returns matrix\n");
980  print_bmat( pLP->LiPM, pLP->m+1, cols+1-pLP->m, cols, pLP->iposv);
981 #endif
982 
983  if( pLP->icase != 0 )
984  { // check for errors
985  WerrorS("mayanPyramidAlg::vDistance:");
986  if( pLP->icase == 1 )
987  WerrorS(" Unbounded v-distance: probably 1st v-coor=0");
988  else if( pLP->icase == -1 )
989  WerrorS(" Infeasible v-distance");
990  else
991  WerrorS(" Unknown error");
992  return -1.0;
993  }
994 
995  return pLP->LiPM[1][1];
996 }
997 
999 {
1000  int i, j, k, cols, cons;
1001  int la_cons_row;
1002 
1003  cons = n+dim+2;
1004 
1005  // first, compute minimum
1006  //
1007 
1008  // common part of the matrix
1009  pLP->LiPM[1][1] = 0.0;
1010  for( i=2; i<=n+2; i++)
1011  {
1012  pLP->LiPM[i][1] = 1.0; // 1st col
1013  pLP->LiPM[i][2] = 0.0; // 2nd col
1014  }
1015 
1016  la_cons_row = 1;
1017  cols = 2;
1018  for( i=0; i<=n; i++)
1019  {
1020  la_cons_row++;
1021  for( j=1; j<= Qi[i]->num; j++)
1022  {
1023  cols++;
1024  pLP->LiPM[1][cols] = 0.0; // set 1st row 0
1025  for( k=2; k<=n+2; k++)
1026  { // lambdas sum up to 1
1027  if( k != la_cons_row) pLP->LiPM[k][cols] = 0.0;
1028  else pLP->LiPM[k][cols] = -1.0;
1029  }
1030  for( k=1; k<=n; k++)
1031  pLP->LiPM[k+n+2][cols] = -(mprfloat)((*Qi[i])[j]->point[k]);
1032  } // j
1033  } // i
1034 
1035  for( i= 0; i < dim; i++ )
1036  { // fixed coords
1037  pLP->LiPM[i+n+3][1] = acoords[i];
1038  pLP->LiPM[i+n+3][2] = 0.0;
1039  }
1040  pLP->LiPM[dim+n+3][1] = 0.0;
1041 
1042 
1043  pLP->LiPM[1][2] = -1.0; // minimize
1044  pLP->LiPM[dim+n+3][2] = 1.0;
1045 
1046 #ifdef mprDEBUG_ALL
1047  Print("\nThats the matrix for minR, dim= %d, acoords= ",dim);
1048  for( i= 0; i < dim; i++ )
1049  Print(" %d",acoords[i]);
1050  PrintLn();
1051  print_mat( pLP->LiPM, cons+1, cols);
1052 #endif
1053 
1054  // simplx finds MIN for obj.fnc, puts it in [1,1]
1055  pLP->m= cons;
1056  pLP->n= cols-1;
1057  pLP->m3= cons;
1058 
1059  pLP->compute();
1060 
1061  if ( pLP->icase != 0 )
1062  { // check for errors
1063  if( pLP->icase < 0)
1064  WerrorS(" mn_mx_MinkowskiSum: LinearProgram: minR: infeasible");
1065  else if( pLP->icase > 0)
1066  WerrorS(" mn_mx_MinkowskiSum: LinearProgram: minR: unbounded");
1067  }
1068 
1069  *minR = (Coord_t)( -pLP->LiPM[1][1] + 1.0 - SIMPLEX_EPS );
1070 
1071  // now compute maximum
1072  //
1073 
1074  // common part of the matrix again
1075  pLP->LiPM[1][1] = 0.0;
1076  for( i=2; i<=n+2; i++)
1077  {
1078  pLP->LiPM[i][1] = 1.0;
1079  pLP->LiPM[i][2] = 0.0;
1080  }
1081  la_cons_row = 1;
1082  cols = 2;
1083  for( i=0; i<=n; i++)
1084  {
1085  la_cons_row++;
1086  for( j=1; j<=Qi[i]->num; j++)
1087  {
1088  cols++;
1089  pLP->LiPM[1][cols] = 0.0;
1090  for( k=2; k<=n+2; k++)
1091  {
1092  if( k != la_cons_row) pLP->LiPM[k][cols] = 0.0;
1093  else pLP->LiPM[k][cols] = -1.0;
1094  }
1095  for( k=1; k<=n; k++)
1096  pLP->LiPM[k+n+2][cols] = -(mprfloat)((*Qi[i])[j]->point[k]);
1097  } // j
1098  } // i
1099 
1100  for( i= 0; i < dim; i++ )
1101  { // fixed coords
1102  pLP->LiPM[i+n+3][1] = acoords[i];
1103  pLP->LiPM[i+n+3][2] = 0.0;
1104  }
1105  pLP->LiPM[dim+n+3][1] = 0.0;
1106 
1107  pLP->LiPM[1][2] = 1.0; // maximize
1108  pLP->LiPM[dim+n+3][2] = 1.0; // var = sum of pnt coords
1109 
1110 #ifdef mprDEBUG_ALL
1111  Print("\nThats the matrix for maxR, dim= %d\n",dim);
1112  print_mat( pLP->LiPM, cons+1, cols);
1113 #endif
1114 
1115  pLP->m= cons;
1116  pLP->n= cols-1;
1117  pLP->m3= cons;
1118 
1119  // simplx finds MAX for obj.fnc, puts it in [1,1]
1120  pLP->compute();
1121 
1122  if ( pLP->icase != 0 )
1123  {
1124  if( pLP->icase < 0)
1125  WerrorS(" mn_mx_MinkowskiSum: LinearProgram: maxR: infeasible");
1126  else if( pLP->icase > 0)
1127  WerrorS(" mn_mx_MinkowskiSum: LinearProgram: maxR: unbounded");
1128  }
1129 
1130  *maxR = (Coord_t)( pLP->LiPM[1][1] + SIMPLEX_EPS );
1131 
1132 #ifdef mprDEBUG_ALL
1133  Print(" Range for dim=%d: [%d,%d]\n", dim, *minR, *maxR);
1134 #endif
1135 }
1136 
1137 // mprSTICKYPROT:
1138 // ST_SPARSE_VREJ: rejected point
1139 // ST_SPARSE_VADD: added point to set
1141 {
1142  mprfloat dist;
1143 
1144  // determine v-distance of point pt
1145  dist= vDistance( &(acoords[0]), n );
1146 
1147  // store only points with v-distance > minVdist
1148  if( dist <= MINVDIST + SIMPLEX_EPS )
1149  {
1151  return false;
1152  }
1153 
1154  E->addPoint( &(acoords[0]) );
1156 
1157  return true;
1158 }
1159 
1160 // mprSTICKYPROT:
1161 // ST_SPARSE_MREC1: recurse
1162 // ST_SPARSE_MREC2: recurse with extra points
1163 // ST_SPARSE_MPEND: end
1165 {
1166  Coord_t minR, maxR;
1167  mprfloat dist;
1168 
1169  // step 3
1170  mn_mx_MinkowskiSum( dim, &minR, &maxR );
1171 
1172 #ifdef mprDEBUG_ALL
1173  int i;
1174  for( i=0; i <= dim; i++) Print("acoords[%d]=%d ",i,(int)acoords[i]);
1175  Print(":: [%d,%d]\n", minR, maxR);
1176 #endif
1177 
1178  // step 5 -> terminate
1179  if( dim == n-1 )
1180  {
1181  int lastKilled = 0;
1182  // insert points
1183  acoords[dim] = minR;
1184  while( acoords[dim] <= maxR )
1185  {
1186  if( !storeMinkowskiSumPoint() )
1187  lastKilled++;
1188  acoords[dim]++;
1189  }
1191  return;
1192  }
1193 
1194  // step 4 -> recurse at step 3
1195  acoords[dim] = minR;
1196  while ( acoords[dim] <= maxR )
1197  {
1198  if ( (acoords[dim] > minR) && (acoords[dim] <= maxR) )
1199  { // acoords[dim] >= minR ??
1201  runMayanPyramid( dim + 1 ); // recurse with higer dimension
1202  }
1203  else
1204  {
1205  // get v-distance of pt
1206  dist= vDistance( &(acoords[0]), dim + 1 );// dim+1 == known coordinates
1207 
1208  if( dist >= SIMPLEX_EPS )
1209  {
1211  runMayanPyramid( dim + 1 ); // recurse with higer dimension
1212  }
1213  }
1214  acoords[dim]++;
1215  } // while
1216 }
1217 //<-
1218 
1219 //-> resMatrixSparse::*
1220 bool resMatrixSparse::remapXiToPoint( const int indx, pointSet **pQ, int *set, int *pnt )
1221 {
1222  int i,nn= (currRing->N);
1223  int loffset= 0;
1224  for ( i= 0; i <= nn; i++ )
1225  {
1226  if ( (loffset < indx) && (indx <= pQ[i]->num + loffset) )
1227  {
1228  *set= i;
1229  *pnt= indx-loffset;
1230  return true;
1231  }
1232  else loffset+= pQ[i]->num;
1233  }
1234  return false;
1235 }
1236 
1237 // mprSTICKYPROT
1238 // ST_SPARSE_RC: point added
1239 int resMatrixSparse::RC( pointSet **pQ, pointSet *E, int vert, mprfloat shift[] )
1240 {
1241  int i, j, k,c ;
1242  int size;
1243  bool found= true;
1244  mprfloat cd;
1245  int onum;
1246  int bucket[MAXVARS+2];
1247  setID *optSum;
1248 
1249  LP->n = 1;
1250  LP->m = n + n + 1; // number of constrains
1251 
1252  // fill in LP matrix
1253  for ( i= 0; i <= n; i++ )
1254  {
1255  size= pQ[i]->num;
1256  for ( k= 1; k <= size; k++ )
1257  {
1258  LP->n++;
1259 
1260  // objective funtion, minimize
1261  LP->LiPM[1][LP->n] = - ( (mprfloat) (*pQ[i])[k]->point[pQ[i]->dim] / SCALEDOWN );
1262 
1263  // lambdas sum up to 1
1264  for ( j = 0; j <= n; j++ )
1265  {
1266  if ( i==j )
1267  LP->LiPM[j+2][LP->n] = -1.0;
1268  else
1269  LP->LiPM[j+2][LP->n] = 0.0;
1270  }
1271 
1272  // the points
1273  for ( j = 1; j <= n; j++ )
1274  {
1275  LP->LiPM[j+n+2][LP->n] = - ( (mprfloat) (*pQ[i])[k]->point[j] );
1276  }
1277  }
1278  }
1279 
1280  for ( j = 0; j <= n; j++ ) LP->LiPM[j+2][1] = 1.0;
1281  for ( j= 1; j <= n; j++ )
1282  {
1283  LP->LiPM[j+n+2][1]= (mprfloat)(*E)[vert]->point[j] - shift[j];
1284  }
1285  LP->n--;
1286 
1287  LP->LiPM[1][1] = 0.0;
1288 
1289 #ifdef mprDEBUG_ALL
1290  PrintLn();
1291  Print(" n= %d, LP->m=M= %d, LP->n=N= %d\n",n,LP->m,LP->n);
1292  print_mat(LP->LiPM, LP->m+1, LP->n+1);
1293 #endif
1294 
1295  LP->m3= LP->m;
1296 
1297  LP->compute();
1298 
1299  if ( LP->icase < 0 )
1300  {
1301  // infeasibility: the point does not lie in a cell -> remove it
1302  return -1;
1303  }
1304 
1305  // store result
1306  (*E)[vert]->point[E->dim]= (int)(-LP->LiPM[1][1] * SCALEDOWN);
1307 
1308 #ifdef mprDEBUG_ALL
1309  Print(" simplx returned %d, Objective value = %f\n", LP->icase, LP->LiPM[1][1]);
1310  //print_bmat(LP->LiPM, NumCons + 1, LP->n+1-NumCons, LP->n+1, LP->iposv); // ( rows= M+1, cols= N+1-m3 )
1311  //print_mat(LP->LiPM, NumCons+1, LP->n);
1312 #endif
1313 
1314 #if 1
1315  // sort LP results
1316  while (found)
1317  {
1318  found=false;
1319  for ( i= 1; i < LP->m; i++ )
1320  {
1321  if ( LP->iposv[i] > LP->iposv[i+1] )
1322  {
1323 
1324  c= LP->iposv[i];
1325  LP->iposv[i]=LP->iposv[i+1];
1326  LP->iposv[i+1]=c;
1327 
1328  cd=LP->LiPM[i+1][1];
1329  LP->LiPM[i+1][1]=LP->LiPM[i+2][1];
1330  LP->LiPM[i+2][1]=cd;
1331 
1332  found= true;
1333  }
1334  }
1335  }
1336 #endif
1337 
1338 #ifdef mprDEBUG_ALL
1339  print_bmat(LP->LiPM, LP->m + 1, LP->n+1-LP->m, LP->n+1, LP->iposv);
1340  PrintS(" now split into sets\n");
1341 #endif
1342 
1343 
1344  // init bucket
1345  for ( i= 0; i <= E->dim; i++ ) bucket[i]= 0;
1346  // remap results of LP to sets Qi
1347  c=0;
1348  optSum= (setID*)omAlloc( (LP->m) * sizeof(struct setID) );
1349  for ( i= 0; i < LP->m; i++ )
1350  {
1351  //Print("% .15f\n",LP->LiPM[i+2][1]);
1352  if ( LP->LiPM[i+2][1] > 1e-12 )
1353  {
1354  if ( !remapXiToPoint( LP->iposv[i+1], pQ, &(optSum[c].set), &(optSum[c].pnt) ) )
1355  {
1356  Werror(" resMatrixSparse::RC: Found bad solution in LP: %d!",LP->iposv[i+1]);
1357  WerrorS(" resMatrixSparse::RC: remapXiToPoint failed!");
1358  return -1;
1359  }
1360  bucket[optSum[c].set]++;
1361  c++;
1362  }
1363  }
1364 
1365  onum= c;
1366  // find last min in bucket[]: maximum i such that Fi is a point
1367  c= 0;
1368  for ( i= 1; i < E->dim; i++ )
1369  {
1370  if ( bucket[c] >= bucket[i] )
1371  {
1372  c= i;
1373  }
1374  }
1375  // find matching point set
1376  for ( i= onum - 1; i >= 0; i-- )
1377  {
1378  if ( optSum[i].set == c )
1379  break;
1380  }
1381  // store
1382  (*E)[vert]->rc.set= c;
1383  (*E)[vert]->rc.pnt= optSum[i].pnt;
1384  (*E)[vert]->rcPnt= (*pQ[c])[optSum[i].pnt];
1385  // count
1386  if ( (*E)[vert]->rc.set == linPolyS ) numSet0++;
1387 
1388 #ifdef mprDEBUG_PROT
1389  Print("\n Point E[%d] was <",vert);print_exp((*E)[vert],E->dim-1);Print(">, bucket={");
1390  for ( j= 0; j < E->dim; j++ )
1391  {
1392  Print(" %d",bucket[j]);
1393  }
1394  PrintS(" }\n optimal Sum: Qi ");
1395  for ( j= 0; j < LP->m; j++ )
1396  {
1397  Print(" [ %d, %d ]",optSum[j].set,optSum[j].pnt);
1398  }
1399  Print(" -> i= %d, j = %d\n",(*E)[vert]->rc.set,optSum[i].pnt);
1400 #endif
1401 
1402  // clean up
1403  omFreeSize( (void *) optSum, (LP->m) * sizeof(struct setID) );
1404 
1406 
1407  return (int)(-LP->LiPM[1][1] * SCALEDOWN);
1408 }
1409 
1410 // create coeff matrix
1412 {
1413  // sparse matrix
1414  // uRPos[i][1]: row of matrix
1415  // uRPos[i][idelem+1]: col of u(0)
1416  // uRPos[i][2..idelem]: col of u(1) .. u(n)
1417  // i= 1 .. numSet0
1418  int i,epos;
1419  int rp,cp;
1420  poly rowp,epp;
1421  poly iterp;
1422  int *epp_mon, *eexp;
1423 
1424  epp_mon= (int *)omAlloc( (n+2) * sizeof(int) );
1425  eexp= (int *)omAlloc0(((currRing->N)+1)*sizeof(int));
1426 
1427  totDeg= numSet0;
1428 
1429  mprSTICKYPROT2(" size of matrix: %d\n", E->num);
1430  mprSTICKYPROT2(" resultant deg: %d\n", numSet0);
1431 
1432  uRPos= new intvec( numSet0, pLength((gls->m)[0])+1, 0 );
1433 
1434  // sparse Matrix represented as a module where
1435  // each poly is column vector ( pSetComp(p,k) gives the row )
1436  rmat= idInit( E->num, E->num ); // cols, rank= number of rows
1437  msize= E->num;
1438 
1439  rp= 1;
1440  rowp= NULL;
1441  epp= pOne();
1442  for ( i= 1; i <= E->num; i++ )
1443  { // for every row
1444  E->getRowMP( i, epp_mon ); // compute (p-a[ij]), (i,j) = RC(p)
1445  pSetExpV( epp, epp_mon );
1446 
1447  //
1448  rowp= ppMult_qq( epp, (gls->m)[(*E)[i]->rc.set] ); // x^(p-a[ij]) * f(i)
1449 
1450  cp= 2;
1451  // get column for every monomial in rowp and store it
1452  iterp= rowp;
1453  while ( iterp!=NULL )
1454  {
1455  epos= E->getExpPos( iterp );
1456  if ( epos == 0 )
1457  {
1458  // this can happen, if the shift vektor or the lift funktions
1459  // are not generically chosen.
1460  Werror("resMatrixSparse::createMatrix: Found exponent not in E, id %d, set [%d, %d]!",
1461  i,(*E)[i]->rc.set,(*E)[i]->rc.pnt);
1462  return i;
1463  }
1464  pSetExpV(iterp,eexp);
1465  pSetComp(iterp, epos );
1466  pSetm(iterp);
1467  if ( (*E)[i]->rc.set == linPolyS )
1468  { // store coeff positions
1469  IMATELEM(*uRPos,rp,cp)= epos;
1470  cp++;
1471  }
1472  pIter( iterp );
1473  } // while
1474  if ( (*E)[i]->rc.set == linPolyS )
1475  { // store row
1476  IMATELEM(*uRPos,rp,1)= i-1;
1477  rp++;
1478  }
1479  (rmat->m)[i-1]= rowp;
1480  } // for
1481 
1482  pDelete( &epp );
1483  omFreeSize( (void *) epp_mon, (n+2) * sizeof(int) );
1484  omFreeSize( (void *) eexp, ((currRing->N)+1)*sizeof(int));
1485 
1486 #ifdef mprDEBUG_ALL
1487  if ( E->num <= 40 )
1488  {
1489  matrix mout= idModule2Matrix( idCopy(rmat) );
1490  print_matrix(mout);
1491  }
1492  for ( i= 1; i <= numSet0; i++ )
1493  {
1494  Print(" row %d contains coeffs of f_%d\n",IMATELEM(*uRPos,i,1),linPolyS);
1495  }
1496  PrintS(" Sparse Matrix done\n");
1497 #endif
1498 
1499  return E->num;
1500 }
1501 
1502 // find a sufficiently generic and small vector
1503 void resMatrixSparse::randomVector( const int dim, mprfloat shift[] )
1504 {
1505  int i,j;
1506  i= 1;
1507 
1508  while ( i <= dim )
1509  {
1510  shift[i]= (mprfloat) (RVMULT*(siRand()%MAXRVVAL)/(mprfloat)MAXRVVAL);
1511  i++;
1512  for ( j= 1; j < i-1; j++ )
1513  {
1514  if ( (shift[j] < shift[i-1] + SIMPLEX_EPS) && (shift[j] > shift[i-1] - SIMPLEX_EPS) )
1515  {
1516  i--;
1517  break;
1518  }
1519  }
1520  }
1521 }
1522 
1524 {
1525  pointSet *vs;
1526  onePoint vert;
1527  int j,k,l;
1528 
1529  vert.point=(Coord_t*)omAlloc( ((currRing->N)+2) * sizeof(Coord_t) );
1530 
1531  vs= new pointSet( dim );
1532 
1533  for ( j= 1; j <= Q1->num; j++ )
1534  {
1535  for ( k= 1; k <= Q2->num; k++ )
1536  {
1537  for ( l= 1; l <= dim; l++ )
1538  {
1539  vert.point[l]= (*Q1)[j]->point[l] + (*Q2)[k]->point[l];
1540  }
1541  vs->mergeWithExp( &vert );
1542  //vs->addPoint( &vert );
1543  }
1544  }
1545 
1546  omFreeSize( (void *) vert.point, ((currRing->N)+2) * sizeof(Coord_t) );
1547 
1548  return vs;
1549 }
1550 
1552 {
1553  pointSet *vs,*vs_old;
1554  int j;
1555 
1556  vs= new pointSet( dim );
1557 
1558  for ( j= 1; j <= pQ[0]->num; j++ ) vs->addPoint( (*pQ[0])[j] );
1559 
1560  for ( j= 1; j < numq; j++ )
1561  {
1562  vs_old= vs;
1563  vs= minkSumTwo( vs_old, pQ[j], dim );
1564 
1565  delete vs_old;
1566  }
1567 
1568  return vs;
1569 }
1570 
1571 //----------------------------------------------------------------------------------------
1572 
1573 resMatrixSparse::resMatrixSparse( const ideal _gls, const int special )
1574  : resMatrixBase(), gls( _gls )
1575 {
1576  pointSet **Qi; // vertices sets of Conv(Supp(f_i)), i=0..idelem
1577  pointSet *E; // all integer lattice points of the minkowski sum of Q0...Qn
1578  int i,k;
1579  int pnt;
1580  int totverts; // total number of exponent vectors in ideal gls
1581  mprfloat shift[MAXVARS+2]; // shiftvector delta, index [1..dim]
1582 
1583  if ( (currRing->N) > MAXVARS )
1584  {
1585  WerrorS("resMatrixSparse::resMatrixSparse: Too many variables!");
1586  return;
1587  }
1588 
1589  rmat= NULL;
1590  numSet0= 0;
1591 
1592  if ( special == SNONE ) linPolyS= 0;
1593  else linPolyS= special;
1594 
1596 
1597  n= (currRing->N);
1598  idelem= IDELEMS(gls); // should be n+1
1599 
1600  // prepare matrix LP->LiPM for Linear Programming
1601  totverts = 0;
1602  for( i=0; i < idelem; i++) totverts += pLength( (gls->m)[i] );
1603 
1604  LP = new simplex( idelem+totverts*2+5, totverts+5 ); // rows, cols
1605 
1606  // get shift vector
1607 #ifdef mprTEST
1608  shift[0]=0.005; shift[1]=0.003; shift[2]=0.008; shift[3]=0.005; shift[4]=0.002;
1609  shift[5]=0.1; shift[6]=0.3; shift[7]=0.2; shift[8]=0.4; shift[9]=0.2;
1610 #else
1611  randomVector( idelem, shift );
1612 #endif
1613 #ifdef mprDEBUG_PROT
1614  PrintS(" shift vector: ");
1615  for ( i= 1; i <= idelem; i++ ) Print(" %.12f ",(double)shift[i]);
1616  PrintLn();
1617 #endif
1618 
1619  // evaluate convex hull for supports of gls
1620  convexHull chnp( LP );
1621  Qi= chnp.newtonPolytopesP( gls );
1622 
1623 #ifdef mprMINKSUM
1624  E= minkSumAll( Qi, n+1, n);
1625 #else
1626  // get inner points
1627  mayanPyramidAlg mpa( LP );
1628  E= mpa.getInnerPoints( Qi, shift );
1629 #endif
1630 
1631 #ifdef mprDEBUG_PROT
1632 #ifdef mprMINKSUM
1633  PrintS("(MinkSum)");
1634 #endif
1635  PrintS("\n E = (Q_0 + ... + Q_n) \\cap \\N :\n");
1636  for ( pnt= 1; pnt <= E->num; pnt++ )
1637  {
1638  Print("%d: <",pnt);print_exp( (*E)[pnt], E->dim );PrintS(">\n");
1639  }
1640  PrintLn();
1641 #endif
1642 
1643 #ifdef mprTEST
1644  int lift[5][5];
1645  lift[0][1]=3; lift[0][2]=4; lift[0][3]=8; lift[0][4]=2;
1646  lift[1][1]=6; lift[1][2]=1; lift[1][3]=7; lift[1][4]=4;
1647  lift[2][1]=2; lift[2][2]=5; lift[2][3]=9; lift[2][4]=6;
1648  lift[3][1]=2; lift[3][2]=1; lift[3][3]=9; lift[3][4]=5;
1649  lift[4][1]=3; lift[4][2]=7; lift[4][3]=1; lift[4][4]=5;
1650  // now lift everything
1651  for ( i= 0; i <= n; i++ ) Qi[i]->lift( lift[i] );
1652 #else
1653  // now lift everything
1654  for ( i= 0; i <= n; i++ ) Qi[i]->lift();
1655 #endif
1656  E->dim++;
1657 
1658  // run Row Content Function for every point in E
1659  for ( pnt= 1; pnt <= E->num; pnt++ )
1660  {
1661  RC( Qi, E, pnt, shift );
1662  }
1663 
1664  // remove points not in cells
1665  k= E->num;
1666  for ( pnt= k; pnt > 0; pnt-- )
1667  {
1668  if ( (*E)[pnt]->rcPnt == NULL )
1669  {
1670  E->removePoint(pnt);
1672  }
1673  }
1674  mprSTICKYPROT("\n");
1675 
1676 #ifdef mprDEBUG_PROT
1677  PrintS(" points which lie in a cell:\n");
1678  for ( pnt= 1; pnt <= E->num; pnt++ )
1679  {
1680  Print("%d: <",pnt);print_exp( (*E)[pnt], E->dim );PrintS(">\n");
1681  }
1682  PrintLn();
1683 #endif
1684 
1685  // unlift to old dimension, sort
1686  for ( i= 0; i <= n; i++ ) Qi[i]->unlift();
1687  E->unlift();
1688  E->sort();
1689 
1690 #ifdef mprDEBUG_PROT
1691  Print(" points with a[ij] (%d):\n",E->num);
1692  for ( pnt= 1; pnt <= E->num; pnt++ )
1693  {
1694  Print("Punkt p \\in E[%d]: <",pnt);print_exp( (*E)[pnt], E->dim );
1695  Print(">, RC(p) = (i:%d, j:%d), a[i,j] = <",(*E)[pnt]->rc.set,(*E)[pnt]->rc.pnt);
1696  //print_exp( (Qi[(*E)[pnt]->rc.set])[(*E)[pnt]->rc.pnt], E->dim );PrintS("> = <");
1697  print_exp( (*E)[pnt]->rcPnt, E->dim );PrintS(">\n");
1698  }
1699 #endif
1700 
1701  // now create matrix
1702  if (E->num <1)
1703  {
1704  WerrorS("could not handle a degenerate situation: no inner points found");
1705  goto theEnd;
1706  }
1707  if ( createMatrix( E ) != E->num )
1708  {
1709  // this can happen if the shiftvector shift is to large or not generic
1711  WerrorS("resMatrixSparse::resMatrixSparse: Error in resMatrixSparse::createMatrix!");
1712  goto theEnd;
1713  }
1714 
1715  theEnd:
1716  // clean up
1717  for ( i= 0; i < idelem; i++ )
1718  {
1719  delete Qi[i];
1720  }
1721  omFreeSize( (void *) Qi, idelem * sizeof(pointSet*) );
1722 
1723  delete E;
1724 
1725  delete LP;
1726 }
1727 
1728 //----------------------------------------------------------------------------------------
1729 
1731 {
1732  delete uRPos;
1733  idDelete( &rmat );
1734 }
1735 
1737 {
1738  int i,/*j,*/cp;
1739  poly pp,phelp,piter,pgls;
1740 
1741  // copy original sparse res matrix
1742  ideal rmat_out= idCopy(rmat);
1743 
1744  // now fill in coeffs of f0
1745  for ( i= 1; i <= numSet0; i++ )
1746  {
1747 
1748  pgls= (gls->m)[0]; // f0
1749 
1750  // get matrix row and delete it
1751  pp= (rmat_out->m)[IMATELEM(*uRPos,i,1)];
1752  pDelete( &pp );
1753  pp= NULL;
1754  phelp= pp;
1755  piter= NULL;
1756 
1757  // u_1,..,u_k
1758  cp=2;
1759  while ( pNext(pgls)!=NULL )
1760  {
1761  phelp= pOne();
1762  pSetCoeff( phelp, nCopy(pGetCoeff(pgls)) );
1763  pSetComp( phelp, IMATELEM(*uRPos,i,cp) );
1764  pSetmComp( phelp );
1765  if ( piter!=NULL )
1766  {
1767  pNext(piter)= phelp;
1768  piter= phelp;
1769  }
1770  else
1771  {
1772  pp= phelp;
1773  piter= phelp;
1774  }
1775  cp++;
1776  pIter( pgls );
1777  }
1778  // u0, now pgls points to last monom
1779  phelp= pOne();
1780  pSetCoeff( phelp, nCopy(pGetCoeff(pgls)) );
1781  //pSetComp( phelp, IMATELEM(*uRPos,i,idelem+1) );
1782  pSetComp( phelp, IMATELEM(*uRPos,i,pLength((gls->m)[0])+1) );
1783  pSetmComp( phelp );
1784  if (piter!=NULL) pNext(piter)= phelp;
1785  else pp= phelp;
1786  (rmat_out->m)[IMATELEM(*uRPos,i,1)]= pp;
1787  }
1788 
1789  return rmat_out;
1790 }
1791 
1792 // Fills in resMat[][] with evpoint[] and gets determinant
1793 // uRPos[i][1]: row of matrix
1794 // uRPos[i][idelem+1]: col of u(0)
1795 // uRPos[i][2..idelem]: col of u(1) .. u(n)
1796 // i= 1 .. numSet0
1797 number resMatrixSparse::getDetAt( const number* evpoint )
1798 {
1799  int i,cp;
1800  poly pp,phelp,piter;
1801 
1802  mprPROTnl("smCallDet");
1803 
1804  for ( i= 1; i <= numSet0; i++ )
1805  {
1806  pp= (rmat->m)[IMATELEM(*uRPos,i,1)];
1807  pDelete( &pp );
1808  pp= NULL;
1809  phelp= pp;
1810  piter= NULL;
1811  // u_1,..,u_n
1812  for ( cp= 2; cp <= idelem; cp++ )
1813  {
1814  if ( !nIsZero(evpoint[cp-1]) )
1815  {
1816  phelp= pOne();
1817  pSetCoeff( phelp, nCopy(evpoint[cp-1]) );
1818  pSetComp( phelp, IMATELEM(*uRPos,i,cp) );
1819  pSetmComp( phelp );
1820  if ( piter )
1821  {
1822  pNext(piter)= phelp;
1823  piter= phelp;
1824  }
1825  else
1826  {
1827  pp= phelp;
1828  piter= phelp;
1829  }
1830  }
1831  }
1832  // u0
1833  phelp= pOne();
1834  pSetCoeff( phelp, nCopy(evpoint[0]) );
1835  pSetComp( phelp, IMATELEM(*uRPos,i,idelem+1) );
1836  pSetmComp( phelp );
1837  pNext(piter)= phelp;
1838  (rmat->m)[IMATELEM(*uRPos,i,1)]= pp;
1839  }
1840 
1841  mprSTICKYPROT(ST__DET); // 1
1842 
1843  poly pres= sm_CallDet( rmat, currRing );
1844  number numres= nCopy( pGetCoeff( pres ) );
1845  pDelete( &pres );
1846 
1847  mprSTICKYPROT(ST__DET); // 2
1848 
1849  return ( numres );
1850 }
1851 
1852 // Fills in resMat[][] with evpoint[] and gets determinant
1853 // uRPos[i][1]: row of matrix
1854 // uRPos[i][idelem+1]: col of u(0)
1855 // uRPos[i][2..idelem]: col of u(1) .. u(n)
1856 // i= 1 .. numSet0
1857 poly resMatrixSparse::getUDet( const number* evpoint )
1858 {
1859  int i,cp;
1860  poly pp,phelp/*,piter*/;
1861 
1862  mprPROTnl("smCallDet");
1863 
1864  for ( i= 1; i <= numSet0; i++ )
1865  {
1866  pp= (rmat->m)[IMATELEM(*uRPos,i,1)];
1867  pDelete( &pp );
1868  phelp= NULL;
1869  // piter= NULL;
1870  for ( cp= 2; cp <= idelem; cp++ )
1871  { // u1 .. un
1872  if ( !nIsZero(evpoint[cp-1]) )
1873  {
1874  phelp= pOne();
1875  pSetCoeff( phelp, nCopy(evpoint[cp-1]) );
1876  pSetComp( phelp, IMATELEM(*uRPos,i,cp) );
1877  //pSetmComp( phelp );
1878  pSetm( phelp );
1879  //Print("comp %d\n",IMATELEM(*uRPos,i,cp));
1880  #if 0
1881  if ( piter!=NULL )
1882  {
1883  pNext(piter)= phelp;
1884  piter= phelp;
1885  }
1886  else
1887  {
1888  pp= phelp;
1889  piter= phelp;
1890  }
1891  #else
1892  pp=pAdd(pp,phelp);
1893  #endif
1894  }
1895  }
1896  // u0
1897  phelp= pOne();
1898  pSetExp(phelp,1,1);
1899  pSetComp( phelp, IMATELEM(*uRPos,i,idelem+1) );
1900  // Print("comp %d\n",IMATELEM(*uRPos,i,idelem+1));
1901  pSetm( phelp );
1902  #if 0
1903  pNext(piter)= phelp;
1904  #else
1905  pp=pAdd(pp,phelp);
1906  #endif
1907  pTest(pp);
1908  (rmat->m)[IMATELEM(*uRPos,i,1)]= pp;
1909  }
1910 
1911  mprSTICKYPROT(ST__DET); // 1
1912 
1913  poly pres= sm_CallDet( rmat, currRing );
1914 
1915  mprSTICKYPROT(ST__DET); // 2
1916 
1917  return ( pres );
1918 }
1919 //<-
1920 
1921 //-----------------------------------------------------------------------------
1922 //-------------- dense resultant matrix ---------------------------------------
1923 //-----------------------------------------------------------------------------
1924 
1925 //-> dense resultant matrix
1926 //
1927 struct resVector;
1928 
1929 /* dense resultant matrix */
1930 class resMatrixDense : virtual public resMatrixBase
1931 {
1932 public:
1933  /**
1934  * _gls: system of multivariate polynoms
1935  * special: -1 -> resMatrixDense is a symbolic matrix
1936  * 0,1, ... -> resMatrixDense ist eine u-Resultante, wobei special das
1937  * lineare u-Polynom angibt
1938  */
1939  resMatrixDense( const ideal _gls, const int special = SNONE );
1940  ~resMatrixDense();
1941 
1942  /** column vector of matrix, index von 0 ... numVectors-1 */
1943  resVector *getMVector( const int i );
1944 
1945  /** Returns the matrix M in an usable presentation */
1946  ideal getMatrix();
1947 
1948  /** Returns the submatrix M' of M in an usable presentation */
1949  ideal getSubMatrix();
1950 
1951  /** Evaluate the determinant of the matrix M at the point evpoint
1952  * where the ui's are replaced by the components of evpoint.
1953  * Uses singclap_det from factory.
1954  */
1955  number getDetAt( const number* evpoint );
1956 
1957  /** Evaluates the determinant of the submatrix M'.
1958  * Since the matrix is numerically, no evaluation point is needed.
1959  * Uses singclap_det from factory.
1960  */
1961  number getSubDet();
1962 
1963 private:
1964  /** deactivated copy constructor */
1965  resMatrixDense( const resMatrixDense & );
1966 
1967  /** Generate the "matrix" M. Each column is presented by a resVector
1968  * holding all entries for this column.
1969  */
1970  void generateBaseData();
1971 
1972  /** Generates needed set of monoms, split them into sets S0, ... Sn and
1973  * check if reduced/nonreduced and calculate size of submatrix.
1974  */
1975  void generateMonomData( int deg, intvec* polyDegs , intvec* iVO );
1976 
1977  /** Recursively generate all homogeneous monoms of
1978  * (currRing->N) variables of degree deg.
1979  */
1980  void generateMonoms( poly m, int var, int deg );
1981 
1982  /** Creates quadratic matrix M of size numVectors for later use.
1983  * u0, u1, ...,un are replaced by 0.
1984  * Entries equal to 0 are not initialized ( == NULL)
1985  */
1986  void createMatrix();
1987 
1988 private:
1990 
1994  int subSize;
1995 
1997 };
1998 //<-
1999 
2000 //-> struct resVector
2001 /* Holds a row vector of the dense resultant matrix */
2003 {
2004 public:
2005  void init()
2006  {
2007  isReduced = FALSE;
2008  elementOfS = SFREE;
2009  mon = NULL;
2010  }
2011  void init( const poly m )
2012  {
2013  isReduced = FALSE;
2014  elementOfS = SFREE;
2015  mon = m;
2016  }
2017 
2018  /** index von 0 ... numVectors-1 */
2019  poly getElem( const int i );
2020 
2021  /** index von 0 ... numVectors-1 */
2022  number getElemNum( const int i );
2023 
2024  // variables
2028 
2029  /** number of the set S mon is element of */
2031 
2032  /** holds the index of u0, u1, ..., un, if (elementOfS == linPolyS)
2033  * the size is given by (currRing->N)
2034  */
2036 
2037  /** holds the column vector if (elementOfS != linPolyS) */
2038  number *numColVector;
2039 
2040  /** size of numColVector */
2042 
2043  number *numColVecCopy;
2044 };
2045 //<-
2046 
2047 //-> resVector::*
2048 poly resVector::getElem( const int i ) // inline ???
2049 {
2050  assume( 0 < i || i > numColVectorSize );
2051  poly out= pOne();
2052  pSetCoeff( out, numColVector[i] );
2053  pTest( out );
2054  return( out );
2055 }
2056 
2057 number resVector::getElemNum( const int i ) // inline ??
2058 {
2059  assume( i >= 0 && i < numColVectorSize );
2060  return( numColVector[i] );
2061 }
2062 //<-
2063 
2064 //-> resMatrixDense::*
2065 resMatrixDense::resMatrixDense( const ideal _gls, const int special )
2066  : resMatrixBase()
2067 {
2068  int i;
2069 
2071  gls= idCopy( _gls );
2072  linPolyS= special;
2073  m=NULL;
2074 
2075  // init all
2076  generateBaseData();
2077 
2078  totDeg= 1;
2079  for ( i= 0; i < IDELEMS(gls); i++ )
2080  {
2081  totDeg*=pTotaldegree( (gls->m)[i] );
2082  }
2083 
2084  mprSTICKYPROT2(" resultant deg: %d\n",totDeg);
2085 
2087 }
2088 
2090 {
2091  int i,j;
2092  for (i=0; i < numVectors; i++)
2093  {
2094  pDelete( &resVectorList[i].mon );
2095  pDelete( &resVectorList[i].dividedBy );
2096  for ( j=0; j < resVectorList[i].numColVectorSize; j++ )
2097  {
2098  nDelete( resVectorList[i].numColVector+j );
2099  }
2100  // OB: ????? (solve_s.tst)
2101  if (resVectorList[i].numColVector!=NULL)
2102  omfreeSize( (void *)resVectorList[i].numColVector,
2103  numVectors * sizeof( number ) );
2104  if (resVectorList[i].numColParNr!=NULL)
2105  omfreeSize( (void *)resVectorList[i].numColParNr,
2106  ((currRing->N)+1) * sizeof(int) );
2107  }
2108 
2109  omFreeSize( (void *)resVectorList, veclistmax*sizeof( resVector ) );
2110 
2111  // free matrix m
2112  if ( m != NULL )
2113  {
2114  idDelete((ideal *)&m);
2115  }
2116 }
2117 
2118 // mprSTICKYPROT:
2119 // ST_DENSE_FR: found row S0
2120 // ST_DENSE_NR: normal row
2122 {
2123  int k,i,j;
2124  resVector *vecp;
2125 
2127 
2128  for ( i= 1; i <= MATROWS( m ); i++ )
2129  for ( j= 1; j <= MATCOLS( m ); j++ )
2130  {
2131  MATELEM(m,i,j)= pInit();
2132  pSetCoeff0( MATELEM(m,i,j), nInit(0) );
2133  }
2134 
2135 
2136  for ( k= 0; k <= numVectors - 1; k++ )
2137  {
2138  if ( linPolyS == getMVector(k)->elementOfS )
2139  {
2141  for ( i= 0; i < (currRing->N); i++ )
2142  {
2143  MATELEM(m,numVectors-k,numVectors-(getMVector(k)->numColParNr)[i])= pInit();
2144  }
2145  }
2146  else
2147  {
2149  vecp= getMVector(k);
2150  for ( i= 0; i < numVectors; i++)
2151  {
2152  if ( !nIsZero( vecp->getElemNum(i) ) )
2153  {
2154  MATELEM(m,numVectors - k,i + 1)= pInit();
2155  pSetCoeff0( MATELEM(m,numVectors - k,i + 1), nCopy(vecp->getElemNum(i)) );
2156  }
2157  }
2158  }
2159  } // for
2160  mprSTICKYPROT("\n");
2161 
2162 #ifdef mprDEBUG_ALL
2163  for ( k= numVectors - 1; k >= 0; k-- )
2164  {
2165  if ( linPolyS == getMVector(k)->elementOfS )
2166  {
2167  for ( i=0; i < (currRing->N); i++ )
2168  {
2169  Print(" %d ",(getMVector(k)->numColParNr)[i]);
2170  }
2171  PrintLn();
2172  }
2173  }
2174  for (i=1; i <= numVectors; i++)
2175  {
2176  for (j=1; j <= numVectors; j++ )
2177  {
2178  pWrite0(MATELEM(m,i,j));PrintS(" ");
2179  }
2180  PrintLn();
2181  }
2182 #endif
2183 }
2184 
2185 // mprSTICKYPROT:
2186 // ST_DENSE_MEM: more mem allocated
2187 // ST_DENSE_NMON: new monom added
2188 void resMatrixDense::generateMonoms( poly mm, int var, int deg )
2189 {
2190  if ( deg == 0 )
2191  {
2192  poly mon = pCopy( mm );
2193 
2194  if ( numVectors == veclistmax )
2195  {
2197  (veclistmax) * sizeof( resVector ),
2198  (veclistmax + veclistblock) * sizeof( resVector ) );
2199  int k;
2200  for ( k= veclistmax; k < (veclistmax + veclistblock); k++ )
2201  resVectorList[k].init();
2204 
2205  }
2206  resVectorList[numVectors].init( mon );
2207  numVectors++;
2209  return;
2210  }
2211  else
2212  {
2213  if ( var == (currRing->N)+1 ) return;
2214  poly newm = pCopy( mm );
2215  while ( deg >= 0 )
2216  {
2217  generateMonoms( newm, var+1, deg );
2218  pIncrExp( newm, var );
2219  pSetm( newm );
2220  deg--;
2221  }
2222  pDelete( & newm );
2223  }
2224 
2225  return;
2226 }
2227 
2228 void resMatrixDense::generateMonomData( int deg, intvec* polyDegs , intvec* iVO )
2229 {
2230  int i,j,k;
2231 
2232  // init monomData
2233  veclistblock= 512;
2236 
2237  // Init resVector()s
2238  for ( j= veclistmax - 1; j >= 0; j-- ) resVectorList[j].init();
2239  numVectors= 0;
2240 
2241  // Generate all monoms of degree deg
2242  poly start= pOne();
2243  generateMonoms( start, 1, deg );
2244  pDelete( & start );
2245 
2246  mprSTICKYPROT("\n");
2247 
2248  // Check for reduced monoms
2249  // First generate polyDegs.rows() monoms
2250  // x(k)^(polyDegs[k]), 0 <= k < polyDegs.rows()
2251  ideal pDegDiv= idInit( polyDegs->rows(), 1 );
2252  for ( k= 0; k < polyDegs->rows(); k++ )
2253  {
2254  poly p= pOne();
2255  pSetExp( p, k + 1, (*polyDegs)[k] );
2256  pSetm( p );
2257  (pDegDiv->m)[k]= p;
2258  }
2259 
2260  // Now check each monom if it is reduced.
2261  // A monom monom is called reduced if there exists
2262  // exactly one x(k)^(polyDegs[k]) that divides the monom.
2263  int divCount;
2264  for ( j= numVectors - 1; j >= 0; j-- )
2265  {
2266  divCount= 0;
2267  for ( k= 0; k < IDELEMS(pDegDiv); k++ )
2268  if ( pLmDivisibleByNoComp( (pDegDiv->m)[k], resVectorList[j].mon ) )
2269  divCount++;
2270  resVectorList[j].isReduced= (divCount == 1);
2271  }
2272 
2273  // create the sets S(k)s
2274  // a monom x(i)^deg, deg given, is element of the set S(i)
2275  // if all x(0)^(polyDegs[0]) ... x(i-1)^(polyDegs[i-1]) DONT divide
2276  // x(i)^deg and only x(i)^(polyDegs[i]) divides x(i)^deg
2277  bool doInsert;
2278  for ( k= 0; k < iVO->rows(); k++)
2279  {
2280  //mprPROTInl(" ------------ var:",(*iVO)[k]);
2281  for ( j= numVectors - 1; j >= 0; j-- )
2282  {
2283  //mprPROTPnl("testing monom",resVectorList[j].mon);
2284  if ( resVectorList[j].elementOfS == SFREE )
2285  {
2286  //mprPROTnl("\tfree");
2287  if ( pLmDivisibleByNoComp( (pDegDiv->m)[ (*iVO)[k] ], resVectorList[j].mon ) )
2288  {
2289  //mprPROTPnl("\tdivisible by ",(pDegDiv->m)[ (*iVO)[k] ]);
2290  doInsert=TRUE;
2291  for ( i= 0; i < k; i++ )
2292  {
2293  //mprPROTPnl("\tchecking db ",(pDegDiv->m)[ (*iVO)[i] ]);
2294  if ( pLmDivisibleByNoComp( (pDegDiv->m)[ (*iVO)[i] ], resVectorList[j].mon ) )
2295  {
2296  //mprPROTPnl("\t and divisible by",(pDegDiv->m)[ (*iVO)[i] ]);
2297  doInsert=FALSE;
2298  break;
2299  }
2300  }
2301  if ( doInsert )
2302  {
2303  //mprPROTInl("\t------------------> S ",(*iVO)[k]);
2304  resVectorList[j].elementOfS= (*iVO)[k];
2305  resVectorList[j].dividedBy= pCopy( (pDegDiv->m)[ (*iVO)[i] ] );
2306  }
2307  }
2308  }
2309  }
2310  }
2311 
2312  // size of submatrix M', equal to number of nonreduced monoms
2313  // (size of matrix M is equal to number of monoms=numVectors)
2314  subSize= 0;
2315  int sub;
2316  for ( i= 0; i < polyDegs->rows(); i++ )
2317  {
2318  sub= 1;
2319  for ( k= 0; k < polyDegs->rows(); k++ )
2320  if ( i != k ) sub*= (*polyDegs)[k];
2321  subSize+= sub;
2322  }
2324 
2325  // pDegDiv wieder freigeben!
2326  idDelete( &pDegDiv );
2327 
2328 #ifdef mprDEBUG_ALL
2329  // Print a list of monoms and their properties
2330  PrintS("// \n");
2331  for ( j= numVectors - 1; j >= 0; j-- )
2332  {
2333  Print("// %s, S(%d), db ",
2334  resVectorList[j].isReduced?"reduced":"nonreduced",
2335  resVectorList[j].elementOfS);
2336  pWrite0(resVectorList[j].dividedBy);
2337  PrintS(" monom ");
2338  pWrite(resVectorList[j].mon);
2339  }
2340  Print("// size: %d, subSize: %d\n",numVectors,subSize);
2341 #endif
2342 }
2343 
2345 {
2346  int k,j,i;
2347  number matEntry;
2348  poly pmatchPos;
2349  poly pi,factor,pmp;
2350 
2351  // holds the degrees of F0, F1, ..., Fn
2352  intvec polyDegs( IDELEMS(gls) );
2353  for ( k= 0; k < IDELEMS(gls); k++ )
2354  polyDegs[k]= pTotaldegree( (gls->m)[k] );
2355 
2356  // the internal Variable Ordering
2357  // make sure that the homogenization variable goes last!
2358  intvec iVO( (currRing->N) );
2359  if ( linPolyS != SNONE )
2360  {
2361  iVO[(currRing->N) - 1]= linPolyS;
2362  int p=0;
2363  for ( k= (currRing->N) - 1; k >= 0; k-- )
2364  {
2365  if ( k != linPolyS )
2366  {
2367  iVO[p]= k;
2368  p++;
2369  }
2370  }
2371  }
2372  else
2373  {
2374  linPolyS= 0;
2375  for ( k= 0; k < (currRing->N); k++ )
2376  iVO[k]= (currRing->N) - k - 1;
2377  }
2378 
2379  // the critical degree d= sum( deg(Fi) ) - n
2380  int sumDeg= 0;
2381  for ( k= 0; k < polyDegs.rows(); k++ )
2382  sumDeg+= polyDegs[k];
2383  sumDeg-= polyDegs.rows() - 1;
2384 
2385  // generate the base data
2386  generateMonomData( sumDeg, &polyDegs, &iVO );
2387 
2388  // generate "matrix"
2389  for ( k= numVectors - 1; k >= 0; k-- )
2390  {
2391  if ( resVectorList[k].elementOfS != linPolyS )
2392  {
2393  // column k is a normal column with numerical or symbolic entries
2394  // init stuff
2397  resVectorList[k].numColVector= (number *)omAlloc( numVectors*sizeof( number ) );
2398  for ( i= 0; i < numVectors; i++ ) resVectorList[k].numColVector[i]= nInit(0);
2399 
2400  // compute row poly
2401  poly pi= ppMult_qq( (gls->m)[ resVectorList[k].elementOfS ] , resVectorList[k].mon );
2402  pi= pDivideM( pCopy( pi ), pCopy( resVectorList[k].dividedBy ) );
2403 
2404  // fill in "matrix"
2405  while ( pi != NULL )
2406  {
2407  matEntry= nCopy(pGetCoeff(pi));
2408  pmatchPos= pLmInit( pi );
2409  pSetCoeff0( pmatchPos, nInit(1) );
2410 
2411  for ( i= 0; i < numVectors; i++)
2412  if ( pLmEqual( pmatchPos, resVectorList[i].mon ) )
2413  break;
2414 
2415  resVectorList[k].numColVector[numVectors - i - 1] = nCopy(matEntry);
2416 
2417  pDelete( &pmatchPos );
2418  nDelete( &matEntry );
2419 
2420  pIter( pi );
2421  }
2422  pDelete( &pi );
2423  }
2424  else
2425  {
2426  // column is a special column, i.e. is generated by S0 and F0
2427  // safe only the positions of the ui's in the column
2428  //mprPROTInl(" setup of numColParNr ",k);
2431  resVectorList[k].numColParNr= (int *)omAlloc0( ((currRing->N)+1) * sizeof(int) );
2432 
2433  pi= (gls->m)[ resVectorList[k].elementOfS ];
2434  factor= pDivideM( pCopy( resVectorList[k].mon ), pCopy( resVectorList[k].dividedBy ) );
2435 
2436  j=0;
2437  while ( pi != NULL )
2438  { // fill in "matrix"
2439  pmp= pMult( pCopy( factor ), pHead( pi ) );
2440  pTest( pmp );
2441 
2442  for ( i= 0; i < numVectors; i++)
2443  if ( pLmEqual( pmp, resVectorList[i].mon ) )
2444  break;
2445 
2447  pDelete( &pmp );
2448  pIter( pi );
2449  j++;
2450  }
2451  pDelete( &pi );
2452  pDelete( &factor );
2453  }
2454  } // for ( k= numVectors - 1; k >= 0; k-- )
2455 
2456  mprSTICKYPROT2(" size of matrix: %d\n",numVectors);
2457  mprSTICKYPROT2(" size of submatrix: %d\n",subSize);
2458 
2459  // create the matrix M
2460  createMatrix();
2461 
2462 }
2463 
2465 {
2466  assume( i >= 0 && i < numVectors );
2467  return &resVectorList[i];
2468 }
2469 
2471 {
2472  int i,j;
2473 
2474  // copy matrix
2475  matrix resmat= mpNew(numVectors,numVectors);
2476  poly p;
2477  for (i=1; i <= numVectors; i++)
2478  {
2479  for (j=1; j <= numVectors; j++ )
2480  {
2481  p=MATELEM(m,i,j);
2482  if (( p!=NULL)
2483  && (!nIsZero(pGetCoeff(p)))
2484  && (pGetCoeff(p)!=NULL)
2485  )
2486  {
2487  MATELEM(resmat,i,j)= pCopy( p );
2488  }
2489  }
2490  }
2491  for (i=0; i < numVectors; i++)
2492  {
2493  if ( resVectorList[i].elementOfS == linPolyS )
2494  {
2495  for (j=1; j <= (currRing->N); j++ )
2496  {
2497  if ( MATELEM(resmat,numVectors-i,
2498  numVectors-resVectorList[i].numColParNr[j-1])!=NULL )
2499  pDelete( &MATELEM(resmat,numVectors-i,numVectors-resVectorList[i].numColParNr[j-1]) );
2500  MATELEM(resmat,numVectors-i,numVectors-resVectorList[i].numColParNr[j-1])= pOne();
2501  // FIX ME
2502  if ( FALSE )
2503  {
2504  pSetCoeff( MATELEM(resmat,numVectors-i,numVectors-resVectorList[i].numColParNr[j-1]), n_Param(j,currRing) );
2505  }
2506  else
2507  {
2508  pSetExp( MATELEM(resmat,numVectors-i,numVectors-resVectorList[i].numColParNr[j-1]), j, 1 );
2509  pSetm(MATELEM(resmat,numVectors-i,numVectors-resVectorList[i].numColParNr[j-1]));
2510  }
2511  }
2512  }
2513  }
2514 
2515  // obachman: idMatrix2Module frees resmat !!
2516  ideal resmod= id_Matrix2Module(resmat,currRing);
2517  return resmod;
2518 }
2519 
2521 {
2522  int k,i,j,l;
2523  resVector *vecp;
2524 
2525  // generate quadratic matrix resmat of size subSize
2526  matrix resmat= mpNew( subSize, subSize );
2527 
2528  j=1;
2529  for ( k= numVectors - 1; k >= 0; k-- )
2530  {
2531  vecp= getMVector(k);
2532  if ( vecp->isReduced ) continue;
2533  l=1;
2534  for ( i= numVectors - 1; i >= 0; i-- )
2535  {
2536  if ( getMVector(i)->isReduced ) continue;
2537  if ( !nIsZero(vecp->getElemNum(numVectors - i - 1)) )
2538  {
2539  MATELEM(resmat,j,l)= pCopy( vecp->getElem(numVectors-i-1) );
2540  }
2541  l++;
2542  }
2543  j++;
2544  }
2545 
2546  // obachman: idMatrix2Module frees resmat !!
2547  ideal resmod= id_Matrix2Module(resmat,currRing);
2548  return resmod;
2549 }
2550 
2551 number resMatrixDense::getDetAt( const number* evpoint )
2552 {
2553  int k,i;
2554 
2555  // copy evaluation point into matrix
2556  // p0, p1, ..., pn replace u0, u1, ..., un
2557  for ( k= numVectors - 1; k >= 0; k-- )
2558  {
2559  if ( linPolyS == getMVector(k)->elementOfS )
2560  {
2561  for ( i= 0; i < (currRing->N); i++ )
2562  {
2563  pSetCoeff( 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
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
int status int void size_t count
Definition: si_signals.h:59
#define ST_SPARSE_MREC1
Definition: mpr_global.h:73
#define pSetmComp(p)
TODO:
Definition: polys.h:255
#define mprSTICKYPROT(msg)
Definition: mpr_global.h:54
complex root finder for univariate polynomials based on laguers algorithm
Definition: mpr_numeric.h:65
pointSet(const int _dim, const int _index=0, const int count=MAXINITELEMS)
Definition: mpr_base.cc:414
#define ppMult_qq(p, q)
Definition: polys.h:191
#define ST_SPARSE_MREC2
Definition: mpr_global.h:74
~pointSet()
Definition: mpr_base.cc:427
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:1239
long isReduced(const nmod_mat_t M)
Definition: facFqBivar.cc:1447
#define pSetm(p)
Definition: polys.h:253
void randomVector(const int dim, mprfloat shift[])
Definition: mpr_base.cc:1503
number * interpolateDense(const number *q)
Solves the Vandermode linear system {i=1}^{n} x_i^k-1 w_i = q_k, k=1,..,n.
Definition: mpr_numeric.cc:159
const poly a
Definition: syzextra.cc:212
void PrintLn()
Definition: reporter.cc:310
#define Print
Definition: emacs.cc:83
#define pAdd(p, q)
Definition: polys.h:186
Definition: mpr_base.cc:151
simplex * pLP
Definition: mpr_base.cc:327
vandermonde system solver for interpolating polynomials from their values
Definition: mpr_numeric.h:28
bool lifted
Definition: mpr_base.cc:166
int dim
Definition: mpr_base.cc:171
CanonicalForm cd(bCommonDen(FF))
Definition: cfModGcd.cc:4030
#define idDelete(H)
delete an ideal
Definition: ideals.h:29
#define nNormalize(n)
Definition: numbers.h:30
CanonicalForm num(const CanonicalForm &f)
#define pSetExp(p, i, v)
Definition: polys.h:42
#define FALSE
Definition: auxiliary.h:94
Linear Programming / Linear Optimization using Simplex - Algorithm.
Definition: mpr_numeric.h:194
Compatiblity layer for legacy polynomial operations (over currRing)
return P p
Definition: myNF.cc:203
#define ST_SPARSE_VREJ
Definition: mpr_global.h:71
#define nPower(a, b, res)
Definition: numbers.h:38
#define mprPROTL(msg, intval)
Definition: mpr_global.h:46
poly sm_CallDet(ideal I, const ring R)
Definition: sparsmat.cc:358
#define pTest(p)
Definition: polys.h:398
bool larger(int, int)
points[a] > points[b] ?
Definition: mpr_base.cc:630
#define mprSTICKYPROT2(msg, arg)
Definition: mpr_global.h:55
int rows() const
Definition: intvec.h:88
#define omFreeSize(addr, size)
Definition: omAllocDecl.h:260
bool checkMem()
Checks, if more mem is needed ( i.e.
Definition: mpr_base.cc:445
#define omfreeSize(addr, size)
Definition: omAllocDecl.h:236
#define SNONE
Definition: mpr_base.h:14
CFFListIterator iter
Definition: facAbsBiFact.cc:54
onePointP * points
Definition: mpr_base.cc:165
poly linearPoly(const resMatType rmt)
Definition: mpr_base.cc:2742
#define TRUE
Definition: auxiliary.h:98
virtual poly getUDet(const number *)
Definition: mpr_base.h:34
void getRowMP(const int indx, int *vert)
Definition: mpr_base.cc:601
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:998
mprfloat * shift
Definition: mpr_base.cc:321
poly mon
Definition: mpr_base.cc:2025
poly singclap_det(const matrix m, const ring s)
Definition: clapsing.cc:1579
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:911
void pWrite(poly p)
Definition: polys.h:290
void WerrorS(const char *s)
Definition: feFopen.cc:24
int k
Definition: cfEzgcd.cc:93
#define ST_BASE_EV
Definition: mpr_global.h:62
resMatrixBase * resMat
Definition: mpr_base.h:92
#define Q
Definition: sirandom.c:25
#define nEqual(n1, n2)
Definition: numbers.h:20
#define MINVDIST
Definition: mpr_base.cc:54
ideal loNewtonPolytope(const ideal id)
Definition: mpr_base.cc:3190
static number & pGetCoeff(poly p)
return an alias to the leading coefficient of p assumes that p != NULL NOTE: not copy ...
Definition: monomials.h:51
Base class for sparse and dense u-Resultant computation.
Definition: mpr_base.h:22
rootContainer ** specializeInU(BOOLEAN matchUp=false, const number subDetVal=NULL)
Definition: mpr_base.cc:3059
#define omAlloc(size)
Definition: omAllocDecl.h:210
void pWrite0(poly p)
Definition: polys.h:291
double mprfloat
Definition: mpr_global.h:17
int index
Definition: mpr_base.cc:172
#define mprPROTnl(msg)
Definition: mpr_global.h:42
poly pp
Definition: myNF.cc:296
#define ST_DENSE_NMON
Definition: mpr_global.h:67
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:2065
resVector * resVectorList
Definition: mpr_base.cc:1989
virtual number getDetAt(const number *)
Definition: mpr_base.h:36
bool found
Definition: facFactorize.cc:56
resMatType rmt
Definition: mpr_base.h:91
#define nPrint(a)
only for debug, over any initalized currRing
Definition: numbers.h:46
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..n_NumberOfParameters(...)
Definition: coeffs.h:817
poly getElem(const int i)
index von 0 ...
Definition: mpr_base.cc:2048
#define ST_DENSE_MEM
Definition: mpr_global.h:66
#define pIter(p)
Definition: monomials.h:44
ideal lift(const ideal J, const ring r, const ideal inI, const ring s)
Definition: lift.cc:26
poly res
Definition: myNF.cc:322
number num
Definition: mpr_base.cc:153
int getExpPos(const poly p)
Definition: mpr_base.cc:580
#define omReallocSize(addr, o_size, size)
Definition: omAllocDecl.h:220
ring currRing
Widely used global variable which specifies the current polynomial ring for Singular interpreter and ...
Definition: polys.cc:10
#define pGetExp(p, i)
Exponent.
Definition: polys.h:41
#define pIncrExp(p, i)
Definition: polys.h:43
#define ST_SPARSE_RC
Definition: mpr_global.h:75
number * numColVector
holds the column vector if (elementOfS != linPolyS)
Definition: mpr_base.cc:2038
const ring r
Definition: syzextra.cc:208
poly monomAt(poly p, int i)
Definition: mpr_base.cc:722
number getElemNum(const int i)
index von 0 ...
Definition: mpr_base.cc:2057
ideal newtonPolytopesI(const ideal gls)
Definition: mpr_base.cc:836
#define nTest(a)
Definition: numbers.h:35
Definition: intvec.h:14
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:551
const CanonicalForm CFMap CFMap & N
Definition: cfEzgcd.cc:49
pointSet * E
Definition: mpr_base.cc:320
for(int i=0;i< R->ExpL_Size;i++) Print("%09lx "
Definition: cfEzgcd.cc:66
int j
Definition: myNF.cc:70
static int max(int a, int b)
Definition: fast_mult.cc:264
pointSet * getInnerPoints(pointSet **_q_i, mprfloat _shift[])
Drive Mayan Pyramid Algorithm.
Definition: mpr_base.cc:893
#define pLmDivisibleByNoComp(a, b)
like pLmDivisibleBy, does not check components
Definition: polys.h:142
static long pTotaldegree(poly p)
Definition: polys.h:264
#define SCALEDOWN
Definition: mpr_base.cc:53
#define MAXRVVAL
Definition: mpr_base.cc:56
#define assume(x)
Definition: mod2.h:394
void unlift()
Definition: mpr_base.cc:231
#define pSetExpV(p, e)
Definition: polys.h:97
int nrows
Definition: cf_linsys.cc:32
#define pLmInit(p)
like pInit, except that expvector is initialized to that of p, p must be != NULL
Definition: polys.h:64
#define pDivideM(a, b)
Definition: polys.h:276
#define ST_DENSE_FR
Definition: mpr_global.h:64
void lift(int *l=NULL)
Lifts the point set using sufficiently generic linear lifting homogeneous forms l[1]..l[dim] in Z.
Definition: mpr_base.cc:672
poly dividedBy
Definition: mpr_base.cc:2026
void mergeWithPoly(const poly p)
Definition: mpr_base.cc:552
convexHull(simplex *_pLP)
Definition: mpr_base.cc:254
pointSet ** Q
Definition: mpr_base.cc:271
P bucket
Definition: myNF.cc:79
#define pSetComp(p, v)
Definition: polys.h:38
int m
Definition: cfEzgcd.cc:119
bool isReduced
Definition: mpr_base.cc:2027
void fillContainer(number *_coeffs, number *_ievpoint, const int _var, const int _tdg, const rootType _rt, const int _anz)
Definition: mpr_numeric.cc:312
#define MAXEVPOINT
Definition: mpr_base.cc:2652
int dim(ideal I, ring r)
int nextPrime(const int p)
Definition: mpr_base.cc:3172
int i
Definition: cfEzgcd.cc:123
poly getUDet(const number *evpoint)
Definition: mpr_base.cc:1857
void PrintS(const char *s)
Definition: reporter.cc:284
mayanPyramidAlg(simplex *_pLP)
Definition: mpr_base.cc:282
#define SFREE
Definition: mpr_base.h:15
int elementOfS
number of the set S mon is element of
Definition: mpr_base.cc:2030
int IsPrime(int p)
Definition: prime.cc:61
#define pOne()
Definition: polys.h:297
virtual ideal getSubMatrix()
Definition: mpr_base.h:32
CanonicalForm factor
Definition: facAbsFact.cc:101
static unsigned pLength(poly a)
Definition: p_polys.h:189
#define ST__DET
Definition: mpr_global.h:78
#define pHead(p)
returns newly allocated copy of Lm(p), coef is copied, next=NULL, p might be NULL ...
Definition: polys.h:67
#define IDELEMS(i)
Definition: simpleideals.h:24
Coord_t * point
Definition: mpr_base.cc:143
rootContainer ** interpolateDenseSP(BOOLEAN matchUp=false, const number subDetVal=NULL)
Definition: mpr_base.cc:2921
pointSet * minkSumTwo(pointSet *Q1, pointSet *Q2, int dim)
Definition: mpr_base.cc:1523
#define nDelete(n)
Definition: numbers.h:16
#define SIMPLEX_EPS
Definition: mpr_numeric.h:181
ideal idCopy(ideal A)
Definition: ideals.h:60
static int index(p_Length length, p_Ord ord)
Definition: p_Procs_Impl.h:592
int numColVectorSize
size of numColVector
Definition: mpr_base.cc:2041
simplex * pLP
Definition: mpr_base.cc:273
int size(const CanonicalForm &f, const Variable &v)
int size ( const CanonicalForm & f, const Variable & v )
Definition: cf_ops.cc:600
matrix mpNew(int r, int c)
create a r x c zero-matrix
Definition: matpol.cc:47
#define ST_DENSE_NR
Definition: mpr_global.h:65
#define pi
Definition: libparse.cc:1143
void init(const poly m)
Definition: mpr_base.cc:2011
ideal idInit(int idsize, int rank)
initialise an ideal / module
Definition: simpleideals.cc:38
struct onePoint * rcPnt
Definition: mpr_base.cc:145
int col
Definition: mpr_base.cc:154
ring sourceRing
Definition: mpr_base.h:48
#define nDiv(a, b)
Definition: numbers.h:32
onePointP operator[](const int index)
Definition: mpr_base.cc:439
#define MATCOLS(i)
Definition: matpol.h:28
#define nIsZero(n)
Definition: numbers.h:19
#define NULL
Definition: omList.c:10
#define RVMULT
Definition: mpr_base.cc:55
int num
Definition: mpr_base.cc:169
REvaluation E(1, terms.length(), IntRandom(25))
number * numColVecCopy
Definition: mpr_base.cc:2043
bool addPoint(const onePointP vert)
Adds a point to pointSet, copy vert[0,...,dim] ot point[num+1][0,...,dim].
Definition: mpr_base.cc:466
#define MAXVARS
Definition: mpr_base.cc:57
unsigned int Coord_t
Definition: mpr_base.cc:133
bool removePoint(const int indx)
Definition: mpr_base.cc:499
resMatrixSparse(const ideal _gls, const int special=SNONE)
Definition: mpr_base.cc:1573
void generateBaseData()
Generate the "matrix" M.
Definition: mpr_base.cc:2344
int linPolyS
Definition: mpr_base.h:47
pointSet * minkSumAll(pointSet **pQ, int numq, int dim)
Definition: mpr_base.cc:1551
ideal getSubMatrix()
Returns the submatrix M&#39; of M in an usable presentation.
Definition: mpr_base.cc:2520
#define pMult(p, q)
Definition: polys.h:190
IStateType istate
Definition: mpr_base.h:44
void sort()
sort lex
Definition: mpr_base.cc:649
bool smaller(int, int)
points[a] < points[b] ?
Definition: mpr_base.cc:611
void createMatrix()
Creates quadratic matrix M of size numVectors for later use.
Definition: mpr_base.cc:2121
#define ST_SPARSE_MPEND
Definition: mpr_global.h:72
int set
Definition: mpr_base.cc:137
bool remapXiToPoint(const int indx, pointSet **pQ, int *set, int *vtx)
Definition: mpr_base.cc:1220
#define pInit()
allocates a new monomial and initializes everything to 0
Definition: polys.h:61
int int ncols
Definition: cf_linsys.cc:32
poly interpolateDense(const number subDetVal=NULL)
Definition: mpr_base.cc:2769
#define pDelete(p_ptr)
Definition: polys.h:169
virtual long getDetDeg()
Definition: mpr_base.h:39
#define ST_SPARSE_RCRJ
Definition: mpr_global.h:76
#define pGetExpV(p, e)
Gets a copy of (resp. set) the exponent vector, where e is assumed to point to (r->N +1)*sizeof(long)...
Definition: polys.h:96
void init()
Definition: mpr_base.cc:2005
#define nCopy(n)
Definition: numbers.h:15
number getSubDet()
Evaluates the determinant of the submatrix M&#39;.
Definition: mpr_base.cc:2592
#define pNext(p)
Definition: monomials.h:43
bool mergeWithExp(const onePointP vert)
Adds point to pointSet, iff pointSet point = .
Definition: mpr_base.cc:514
#define mprPROTN(msg, nval)
Definition: mpr_global.h:48
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:1797
#define pSetCoeff0(p, n)
Definition: monomials.h:67
ideal getMatrix()
Returns the matrix M in an usable presentation.
Definition: mpr_base.cc:2470
int max
Definition: mpr_base.cc:170
#define phelp
Definition: libparse.cc:1240
void sort(CFArray &A, int l=0)
quick sort A
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:778
void generateMonomData(int deg, intvec *polyDegs, intvec *iVO)
Generates needed set of monoms, split them into sets S0, ...
Definition: mpr_base.cc:2228
p exp[i]
Definition: DebugPrint.cc:39
#define LIFT_COOR
Definition: mpr_base.cc:52
setID rc
Definition: mpr_base.cc:144
#define MATROWS(i)
Definition: matpol.h:27
struct _entry * next
Definition: mpr_base.cc:155
polyrec * poly
Definition: hilb.h:10
#define ST_SPARSE_MEM
Definition: mpr_global.h:69
int * numColParNr
holds the index of u0, u1, ..., un, if (elementOfS == linPolyS) the size is given by (currRing->N) ...
Definition: mpr_base.cc:2035
#define mprPROT(msg)
Definition: mpr_global.h:41
int createMatrix(pointSet *E)
create coeff matrix uRPos[i][1]: row of matrix uRPos[i][idelem+1]: col of u(0) uRPos[i][2..idelem]: col of u(1) .
Definition: mpr_base.cc:1411
#define nInit(i)
Definition: numbers.h:24
int pnt
Definition: mpr_base.cc:138
uResultant(const ideal _gls, const resMatType _rmt=sparseResMat, BOOLEAN extIdeal=true)
Definition: mpr_base.cc:2684
Rational pow(const Rational &a, int e)
Definition: GMPrat.cc:418
int BOOLEAN
Definition: auxiliary.h:85
#define IMATELEM(M, I, J)
Definition: intvec.h:77
const poly b
Definition: syzextra.cc:213
#define pSetCoeff(p, n)
deletes old coeff before setting the new one
Definition: polys.h:31
#define ST_SPARSE_VADD
Definition: mpr_global.h:70
#define MAXINITELEMS
Definition: mpr_base.cc:51
unsigned long over(const unsigned long n, const unsigned long d)
Definition: mpr_base.cc:2658
ideal gls
Definition: mpr_base.h:46
ideal id_Matrix2Module(matrix mat, const ring R)
resVector * getMVector(const int i)
column vector of matrix, index von 0 ...
Definition: mpr_base.cc:2464
void Werror(const char *fmt,...)
Definition: reporter.cc:189
#define mprPROTNnl(msg, nval)
Definition: mpr_global.h:49
ideal extendIdeal(const ideal gls, poly linPoly, const resMatType rmt)
Definition: mpr_base.cc:2714
virtual number getSubDet()
Definition: mpr_base.h:37
#define pLmEqual(p1, p2)
Definition: polys.h:111
void generateMonoms(poly m, int var, int deg)
Recursively generate all homogeneous monoms of (currRing->N) variables of degree deg.
Definition: mpr_base.cc:2188
int siRand()
Definition: sirandom.c:41
#define omAlloc0(size)
Definition: omAllocDecl.h:211
return result
Definition: facAbsBiFact.cc:76
int l
Definition: cfEzgcd.cc:94
simplex * LP
Definition: mpr_base.cc:126
intvec * uRPos
Definition: mpr_base.cc:122
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:732
number getDetAt(const number *evpoint)
Evaluate the determinant of the matrix M at the point evpoint where the ui&#39;s are replaced by the comp...
Definition: mpr_base.cc:2551
void runMayanPyramid(int dim)
Recursive Mayan Pyramid algorithm for directly computing MinkowskiSum lattice points for (n+1)-fold M...
Definition: mpr_base.cc:1164
bool storeMinkowskiSumPoint()
Stores point in E->points[pt], iff v-distance != 0 Returns true iff point was stored, else flase.
Definition: mpr_base.cc:1140
ideal getMatrix()
Definition: mpr_base.cc:1736
#define pCopy(p)
return a copy of the poly
Definition: polys.h:168
#define MATELEM(mat, i, j)
Definition: matpol.h:29
ideal gls
Definition: mpr_base.h:88
pointSet ** Qi
Definition: mpr_base.cc:319