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