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