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