source: git/kernel/numeric/mpr_base.cc @ 1b2b4f

spielwiese
Last change on this file since 1b2b4f was 1b2b4f, checked in by Hans Schoenemann <hannes@…>, 19 months ago
fix: mpresmat
  • Property mode set to 100644
File size: 74.0 KB
Line 
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 */
61class pointSet;
62
63
64
65/* sparse resultant matrix class */
66class resMatrixSparse : virtual public resMatrixBase
67{
68public:
69  resMatrixSparse( const ideal _gls, const int special = SNONE );
70  ~resMatrixSparse();
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
85private:
86  resMatrixSparse( const resMatrixSparse & );
87
88  void randomVector( const int dim, mprfloat shift[] );
89
90  /** Row Content Function
91   * Finds the largest i such that F[i] is a point, F[i]= a[ij] in A[i] for some j.
92   * Returns -1 iff the point vert does not lie in a cell
93   */
94  int RC( pointSet **pQ, pointSet *E, int vert, mprfloat shift[] );
95
96  /* Remaps a result of LP to the according point set Qi.
97   * Returns false iff remaping was not possible, otherwise true.
98   */
99  bool remapXiToPoint( const int indx, pointSet **pQ, int *set, int *vtx );
100
101  /** create coeff matrix
102   * uRPos[i][1]: row of matrix
103   * uRPos[i][idelem+1]: col of u(0)
104   *  uRPos[i][2..idelem]: col of u(1) .. u(n)
105   *  i= 1 .. numSet0
106   * Returns the dimension of the matrix or -1 in case of an error
107   */
108  int createMatrix( pointSet *E );
109
110  pointSet * minkSumAll( pointSet **pQ, int numq, int dim );
111  pointSet * minkSumTwo( pointSet *Q1, pointSet *Q2, int dim );
112
113private:
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
120  intvec *uRPos;
121
122  ideal rmat;        // sparse matrix representation
123
124  simplex * LP;      // linear programming stuff
125};
126//<-
127
128//-> typedefs and structs
129poly monomAt( poly p, int i );
130
131typedef unsigned int Coord_t;
132
133struct setID
134{
135  int set;
136  int pnt;
137};
138
139struct 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
146typedef struct onePoint * onePointP;
147
148/* sparse matrix entry */
149struct _entry
150{
151  number num;
152  int col;
153  struct _entry * next;
154};
155
156typedef struct _entry * entry;
157//<-
158
159//-> class pointSet
160class pointSet
161{
162private:
163  onePointP *points;     // set of onePoint's, index [1..num], supports of monoms
164  bool lifted;
165
166public:
167  int num;               // number of elements in points
168  int max;               // maximal entries in points, i.e. allocated mem
169  int dim;               // dimension, i.e. valid coord entries in point
170  int index;             // should hold unique identifier of point set
171
172  pointSet( const int _dim, const int _index= 0, const int count= MAXINITELEMS );
173   ~pointSet();
174
175  // pointSet.points[i] equals pointSet[i]
176  inline onePointP operator[] ( const int index );
177
178  /** Adds a point to pointSet, copy vert[0,...,dim] ot point[num+1][0,...,dim].
179   * Returns false, iff additional memory was allocated ( i.e. num >= max )
180   * else returns true
181   */
182  bool addPoint( const onePointP vert );
183
184  /** Adds a point to pointSet, copy vert[0,...,dim] ot point[num+1][0,...,dim].
185   * Returns false, iff additional memory was allocated ( i.e. num >= max )
186   * else returns true
187   */
188  bool addPoint( const int * vert );
189
190  /** Adds a point to pointSet, copy vert[0,...,dim] ot point[num+1][0,...,dim].
191   * Returns false, iff additional memory was allocated ( i.e. num >= max )
192   * else returns true
193   */
194  bool addPoint( const Coord_t * vert );
195
196  /* Removes the point at intex indx */
197  bool removePoint( const int indx );
198
199  /** Adds point to pointSet, iff pointSet \cap point = \emptyset.
200   * Returns true, iff added, else false.
201   */
202  bool mergeWithExp( const onePointP vert );
203
204  /** Adds point to pointSet, iff pointSet \cap point = \emptyset.
205   * Returns true, iff added, else false.
206   */
207  bool mergeWithExp( const int * vert );
208
209  /* Adds support of poly p to pointSet, iff pointSet \cap point = \emptyset. */
210  void mergeWithPoly( const poly p );
211
212  /* Returns the row polynom multiplicator in vert[] */
213  void getRowMP( const int indx, int * vert );
214
215  /* Returns index of supp(LT(p)) in pointSet. */
216  int getExpPos( const poly p );
217
218  /** sort lex
219   */
220  void sort();
221
222  /** Lifts the point set using sufficiently generic linear lifting
223   * homogeneous forms l[1]..l[dim] in Z. Every l[i] is of the form
224   * L1x1+...+Lnxn, for generic L1..Ln in Z.
225   *
226   * Lifting raises dimension by one!
227   */
228  void lift( int *l= NULL );     // !! increments dim by 1
229  void unlift() { dim--; lifted= false; }
230
231private:
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 */
249class convexHull
250{
251public:
252  convexHull( simplex * _pLP ) : pLP(_pLP) {}
253  ~convexHull() {}
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
262private:
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
268private:
269  pointSet **Q;
270  int n;
271  simplex * pLP;
272};
273//<-
274
275//-> class mayanPyramidAlg
276/* Compute all lattice points in a given convex hull */
277class mayanPyramidAlg
278{
279public:
280  mayanPyramidAlg( simplex * _pLP ) : n((currRing->N)), pLP(_pLP) {}
281  ~mayanPyramidAlg() {}
282
283  /** Drive Mayan Pyramid Algorithm.
284   * The Alg computes conv(Qi[]+shift[]).
285   */
286  pointSet * getInnerPoints( pointSet **_q_i, mprfloat _shift[] );
287
288private:
289
290  /** Recursive Mayan Pyramid algorithm for directly computing MinkowskiSum
291   * lattice points for (n+1)-fold MinkowskiSum of given point sets Qi[].
292   * Recursively for range of dim: dim in [0..n); acoords[0..var) fixed.
293   * Stores only MinkowskiSum points of udist > 0: done by storeMinkowskiSumPoints.
294   */
295  void runMayanPyramid( int dim );
296
297  /**  Compute v-distance via Linear Programing
298   * Linear Program finds the v-distance of the point in accords[].
299   * The v-distance is the distance along the direction v to boundary of
300   * Minkowski Sum of Qi (here vector v is represented by shift[]).
301   * Returns the v-distance or -1.0 if an error occurred.
302   */
303  mprfloat vDistance( Coord_t * acoords, int dim );
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
316private:
317  pointSet **Qi;
318  pointSet *E;
319  mprfloat *shift;
320
321  int n,idelem;
322
323  Coord_t acoords[MAXVARS+2];
324
325  simplex * pLP;
326};
327//<-
328
329//-> debug output stuff
330#if defined(mprDEBUG_PROT) || defined(mprDEBUG_ALL)
331void 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}
342void 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
359void 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}
370void 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::*
412pointSet::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
425pointSet::~pointSet()
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
437inline onePointP pointSet::operator[] ( const int index_i )
438{
439  assume( index_i > 0 && index_i <= num );
440  return points[index_i];
441}
442
443inline 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;
458    mprSTICKYPROT(ST_SPARSE_MEM);
459    return false;
460  }
461  return true;
462}
463
464bool 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
475bool 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
486bool 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
497bool 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
512bool 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
531bool 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
550void 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
578int 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
599void 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
609inline 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
628inline 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
647void pointSet::sort()
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
670void 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
720poly 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::*
730bool 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)
776pointSet ** convexHull::newtonPolytopesP( const ideal gls )
777{
778  int i, j, k;
779  int m;  // Anzahl der Exponentvektoren im i-ten Polynom (gls->m)[i] des Ideals gls
780  int idelem= IDELEMS(gls);
781  int * vert;
782
783  n= (currRing->N);
784  vert= (int *)omAlloc( (idelem+1) * sizeof(int) );
785
786  Q = (pointSet **)omAlloc( idelem * sizeof(pointSet*) );        // support hulls
787  for ( i= 0; i < idelem; i++ )
788    Q[i] = new pointSet( (currRing->N), i+1, pLength((gls->m)[i])+1 );
789
790  for( i= 0; i < idelem; i++ )
791  {
792    k=1;
793    m = pLength( (gls->m)[i] );
794
795    poly p= (gls->m)[i];
796    for( j= 1; j <= m; j++) {  // für jeden Exponentvektor
797      if( !inHull( (gls->m)[i], p, m, j ) )
798      {
799        p_GetExpV( p, vert, currRing );
800        Q[i]->addPoint( vert );
801        k++;
802        mprSTICKYPROT(ST_SPARSE_VADD);
803      }
804      else
805      {
806        mprSTICKYPROT(ST_SPARSE_VREJ);
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)
834ideal convexHull::newtonPolytopesI( const ideal gls )
835{
836  int i, j;
837  int m;  // Anzahl der Exponentvektoren im i-ten Polynom (gls->m)[i] des Ideals gls
838  int idelem= IDELEMS(gls);
839  ideal id;
840  poly p,pid;
841  int * vert;
842
843  n= (currRing->N);
844  vert= (int *)omAlloc( (idelem+1) * sizeof(int) );
845  id= idInit( idelem, 1 );
846
847  for( i= 0; i < idelem; i++ )
848  {
849    m = pLength( (gls->m)[i] );
850
851    p= (gls->m)[i];
852    for( j= 1; j <= m; j++) {  // für jeden Exponentvektor
853      if( !inHull( (gls->m)[i], p, m, j ) )
854      {
855        if ( (id->m)[i] == NULL )
856        {
857          (id->m)[i]= pHead(p);
858          pid=(id->m)[i];
859        }
860        else
861        {
862          pNext(pid)= pHead(p);
863          pIter(pid);
864          pNext(pid)= NULL;
865        }
866        mprSTICKYPROT(ST_SPARSE_VADD);
867      }
868      else
869      {
870        mprSTICKYPROT(ST_SPARSE_VREJ);
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::*
891pointSet * mayanPyramidAlg::getInnerPoints( pointSet **_q_i, mprfloat _shift[] )
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
909mprfloat mayanPyramidAlg::vDistance( Coord_t * acoords_a, int dim )
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
996void  mayanPyramidAlg::mn_mx_MinkowskiSum( int dim, Coord_t *minR, Coord_t *maxR )
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
1138bool mayanPyramidAlg::storeMinkowskiSumPoint()
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  {
1148    mprSTICKYPROT(ST_SPARSE_VREJ);
1149    return false;
1150  }
1151
1152  E->addPoint( &(acoords[0]) );
1153  mprSTICKYPROT(ST_SPARSE_VADD);
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
1162void mayanPyramidAlg::runMayanPyramid( int dim )
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    }
1188    mprSTICKYPROT(ST_SPARSE_MPEND);
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  ??
1198      mprSTICKYPROT(ST_SPARSE_MREC1);
1199      runMayanPyramid( dim + 1 );         // recurse with higer dimension
1200    }
1201    else
1202    {
1203      // get v-distance of pt
1204      dist= vDistance( &(acoords[0]), dim + 1 );// dim+1 == known coordinates
1205
1206      if( dist >= SIMPLEX_EPS )
1207      {
1208        mprSTICKYPROT(ST_SPARSE_MREC2);
1209        runMayanPyramid( dim + 1 );       // recurse with higer dimension
1210      }
1211    }
1212    acoords[dim]++;
1213  } // while
1214}
1215//<-
1216
1217//-> resMatrixSparse::*
1218bool 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
1237int resMatrixSparse::RC( pointSet **pQ, pointSet *E, int vert, mprfloat shift[] )
1238{
1239  int i, j, k,c ;
1240  int size;
1241  bool found= true;
1242  mprfloat cd;
1243  int onum;
1244  int bucket[MAXVARS+2];
1245  setID *optSum;
1246
1247  LP->n = 1;
1248  LP->m = n + n + 1;   // number of constrains
1249
1250  // fill in LP matrix
1251  for ( i= 0; i <= n; i++ )
1252  {
1253    size= pQ[i]->num;
1254    for ( k= 1; k <= size; k++ )
1255    {
1256      LP->n++;
1257
1258      // objective funtion, minimize
1259      LP->LiPM[1][LP->n] = - ( (mprfloat) (*pQ[i])[k]->point[pQ[i]->dim] / SCALEDOWN );
1260
1261      // lambdas sum up to 1
1262      for ( j = 0; j <= n; j++ )
1263      {
1264        if ( i==j )
1265          LP->LiPM[j+2][LP->n] = -1.0;
1266        else
1267          LP->LiPM[j+2][LP->n] = 0.0;
1268      }
1269
1270      // the points
1271      for ( j = 1; j <= n; j++ )
1272      {
1273        LP->LiPM[j+n+2][LP->n] =  - ( (mprfloat) (*pQ[i])[k]->point[j] );
1274      }
1275    }
1276  }
1277
1278  for ( j = 0; j <= n; j++ ) LP->LiPM[j+2][1] = 1.0;
1279  for ( j= 1; j <= n; j++ )
1280  {
1281    LP->LiPM[j+n+2][1]= (mprfloat)(*E)[vert]->point[j] - shift[j];
1282  }
1283  LP->n--;
1284
1285  LP->LiPM[1][1] = 0.0;
1286
1287#ifdef mprDEBUG_ALL
1288  PrintLn();
1289  Print(" n= %d, LP->m=M= %d, LP->n=N= %d\n",n,LP->m,LP->n);
1290  print_mat(LP->LiPM, LP->m+1, LP->n+1);
1291#endif
1292
1293  LP->m3= LP->m;
1294
1295  LP->compute();
1296
1297  if ( LP->icase < 0 )
1298  {
1299    // infeasibility: the point does not lie in a cell -> remove it
1300    return -1;
1301  }
1302
1303  // store result
1304  (*E)[vert]->point[E->dim]= (int)(-LP->LiPM[1][1] * SCALEDOWN);
1305
1306#ifdef mprDEBUG_ALL
1307  Print(" simplx returned %d, Objective value = %f\n", LP->icase, LP->LiPM[1][1]);
1308  //print_bmat(LP->LiPM, NumCons + 1, LP->n+1-NumCons, LP->n+1, LP->iposv); // ( rows= M+1, cols= N+1-m3 )
1309  //print_mat(LP->LiPM, NumCons+1, LP->n);
1310#endif
1311
1312#if 1
1313  // sort LP results
1314  while (found)
1315  {
1316    found=false;
1317    for ( i= 1; i < LP->m; i++ )
1318    {
1319      if ( LP->iposv[i] > LP->iposv[i+1] )
1320      {
1321
1322        c= LP->iposv[i];
1323        LP->iposv[i]=LP->iposv[i+1];
1324        LP->iposv[i+1]=c;
1325
1326        cd=LP->LiPM[i+1][1];
1327        LP->LiPM[i+1][1]=LP->LiPM[i+2][1];
1328        LP->LiPM[i+2][1]=cd;
1329
1330        found= true;
1331      }
1332    }
1333  }
1334#endif
1335
1336#ifdef mprDEBUG_ALL
1337  print_bmat(LP->LiPM, LP->m + 1, LP->n+1-LP->m, LP->n+1, LP->iposv);
1338  PrintS(" now split into sets\n");
1339#endif
1340
1341
1342  // init bucket
1343  for ( i= 0; i <= E->dim; i++ ) bucket[i]= 0;
1344  // remap results of LP to sets Qi
1345  c=0;
1346  optSum= (setID*)omAlloc( (LP->m) * sizeof(struct setID) );
1347  for ( i= 0; i < LP->m; i++ )
1348  {
1349    //Print("% .15f\n",LP->LiPM[i+2][1]);
1350    if ( LP->LiPM[i+2][1] > 1e-12 )
1351    {
1352      if ( !remapXiToPoint( LP->iposv[i+1], pQ, &(optSum[c].set), &(optSum[c].pnt) ) )
1353      {
1354        Werror(" resMatrixSparse::RC: Found bad solution in LP: %d!",LP->iposv[i+1]);
1355        WerrorS(" resMatrixSparse::RC: remapXiToPoint failed!");
1356        return -1;
1357      }
1358      bucket[optSum[c].set]++;
1359      c++;
1360    }
1361  }
1362
1363  onum= c;
1364  // find last min in bucket[]: maximum i such that Fi is a point
1365  c= 0;
1366  for ( i= 1; i < E->dim; i++ )
1367  {
1368    if ( bucket[c] >= bucket[i] )
1369    {
1370      c= i;
1371    }
1372  }
1373  // find matching point set
1374  for ( i= onum - 1; i >= 0; i-- )
1375  {
1376    if ( optSum[i].set == c )
1377      break;
1378  }
1379  // store
1380  (*E)[vert]->rc.set= c;
1381  (*E)[vert]->rc.pnt= optSum[i].pnt;
1382  (*E)[vert]->rcPnt= (*pQ[c])[optSum[i].pnt];
1383  // count
1384  if ( (*E)[vert]->rc.set == linPolyS ) numSet0++;
1385
1386#ifdef mprDEBUG_PROT
1387  Print("\n Point E[%d] was <",vert);print_exp((*E)[vert],E->dim-1);Print(">, bucket={");
1388  for ( j= 0; j < E->dim; j++ )
1389  {
1390    Print(" %d",bucket[j]);
1391  }
1392  PrintS(" }\n optimal Sum: Qi ");
1393  for ( j= 0; j < LP->m; j++ )
1394  {
1395    Print(" [ %d, %d ]",optSum[j].set,optSum[j].pnt);
1396  }
1397  Print(" -> i= %d, j = %d\n",(*E)[vert]->rc.set,optSum[i].pnt);
1398#endif
1399
1400  // clean up
1401  omFreeSize( (void *) optSum, (LP->m) * sizeof(struct setID) );
1402
1403  mprSTICKYPROT(ST_SPARSE_RC);
1404
1405  return (int)(-LP->LiPM[1][1] * SCALEDOWN);
1406}
1407
1408// create coeff matrix
1409int resMatrixSparse::createMatrix( pointSet *E )
1410{
1411  // sparse matrix
1412  //    uRPos[i][1]: row of matrix
1413  //    uRPos[i][idelem+1]: col of u(0)
1414  //    uRPos[i][2..idelem]: col of u(1) .. u(n)
1415  //    i= 1 .. numSet0
1416  int i,epos;
1417  int rp,cp;
1418  poly rowp,epp;
1419  poly iterp;
1420  int *epp_mon, *eexp;
1421
1422  epp_mon= (int *)omAlloc( (n+2) * sizeof(int) );
1423  eexp= (int *)omAlloc0(((currRing->N)+1)*sizeof(int));
1424
1425  totDeg= numSet0;
1426
1427  mprSTICKYPROT2(" size of matrix: %d\n", E->num);
1428  mprSTICKYPROT2("  resultant deg: %d\n", numSet0);
1429
1430  uRPos= new intvec( numSet0, pLength((gls->m)[0])+1, 0 );
1431
1432  // sparse Matrix represented as a module where
1433  // each poly is column vector ( pSetComp(p,k) gives the row )
1434  rmat= idInit( E->num, E->num );    // cols, rank= number of rows
1435  msize= E->num;
1436
1437  rp= 1;
1438  rowp= NULL;
1439  epp= pOne();
1440  for ( i= 1; i <= E->num; i++ )
1441  {       // for every row
1442    E->getRowMP( i, epp_mon );           // compute (p-a[ij]), (i,j) = RC(p)
1443    pSetExpV( epp, epp_mon );
1444
1445    //
1446    rowp= ppMult_qq( epp, (gls->m)[(*E)[i]->rc.set] );  // x^(p-a[ij]) * f(i)
1447
1448    cp= 2;
1449    // get column for every monomial in rowp and store it
1450    iterp= rowp;
1451    while ( iterp!=NULL )
1452    {
1453      epos= E->getExpPos( iterp );
1454      if ( epos == 0 )
1455      {
1456        // this can happen, if the shift vektor or the lift funktions
1457        // are not generically chosen.
1458        Werror("resMatrixSparse::createMatrix: Found exponent not in E, id %d, set [%d, %d]!",
1459               i,(*E)[i]->rc.set,(*E)[i]->rc.pnt);
1460        return i;
1461      }
1462      pSetExpV(iterp,eexp);
1463      pSetComp(iterp, epos );
1464      pSetm(iterp);
1465      if ( (*E)[i]->rc.set == linPolyS )
1466      { // store coeff positions
1467        IMATELEM(*uRPos,rp,cp)= epos;
1468        cp++;
1469      }
1470      pIter( iterp );
1471    } // while
1472    if ( (*E)[i]->rc.set == linPolyS )
1473    {   // store row
1474      IMATELEM(*uRPos,rp,1)= i-1;
1475      rp++;
1476    }
1477    (rmat->m)[i-1]= rowp;
1478  } // for
1479
1480  pDelete( &epp );
1481  omFreeSize( (void *) epp_mon, (n+2) * sizeof(int) );
1482  omFreeSize( (void *) eexp, ((currRing->N)+1)*sizeof(int));
1483
1484#ifdef mprDEBUG_ALL
1485  if ( E->num <= 40 )
1486  {
1487    matrix mout= idModule2Matrix( idCopy(rmat) );
1488    print_matrix(mout);
1489  }
1490  for ( i= 1; i <= numSet0; i++ )
1491  {
1492    Print(" row  %d contains coeffs of f_%d\n",IMATELEM(*uRPos,i,1),linPolyS);
1493  }
1494  PrintS(" Sparse Matrix done\n");
1495#endif
1496
1497  return E->num;
1498}
1499
1500// find a sufficiently generic and small vector
1501void 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
1521pointSet * resMatrixSparse::minkSumTwo( pointSet *Q1, pointSet *Q2, int dim )
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
1549pointSet * resMatrixSparse::minkSumAll( pointSet **pQ, int numq, int dim )
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
1571resMatrixSparse::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
1593  istate= resMatrixBase::ready;
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);
1669      mprSTICKYPROT(ST_SPARSE_RCRJ);
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
1708    istate= resMatrixBase::fatalError;
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
1728resMatrixSparse::~resMatrixSparse()
1729{
1730  delete uRPos;
1731  idDelete( &rmat );
1732}
1733
1734ideal resMatrixSparse::getMatrix()
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
1796number 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
1856poly 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//
1926struct resVector;
1927
1928/* dense resultant matrix */
1929class resMatrixDense : virtual public resMatrixBase
1930{
1931public:
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
1962private:
1963  /** deactivated copy constructor */
1964  resMatrixDense( const resMatrixDense & );
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
1987private:
1988  resVector *resVectorList;
1989
1990  int veclistmax;
1991  int veclistblock;
1992  int numVectors;
1993  int subSize;
1994
1995  matrix m;
1996};
1997//<-
1998
1999//-> struct resVector
2000/* Holds a row vector of the dense resultant matrix */
2001struct resVector
2002{
2003public:
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;
2025  poly dividedBy;
2026  bool isReduced;
2027
2028  /** number of the set S mon is element of */
2029  int elementOfS;
2030
2031  /** holds the index of u0, u1, ..., un, if (elementOfS == linPolyS)
2032   *  the size is given by (currRing->N)
2033   */
2034  int *numColParNr;
2035
2036  /** holds the column vector if (elementOfS != linPolyS) */
2037  number *numColVector;
2038
2039  /** size of numColVector */
2040  int numColVectorSize;
2041
2042  number *numColVecCopy;
2043};
2044//<-
2045
2046//-> resVector::*
2047poly 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
2056number resVector::getElemNum( const int i ) // inline ??
2057{
2058  assume( i >= 0 && i < numColVectorSize );
2059  return( numColVector[i] );
2060}
2061//<-
2062
2063//-> resMatrixDense::*
2064resMatrixDense::resMatrixDense( const ideal _gls, const int special )
2065  : resMatrixBase()
2066{
2067  int i;
2068
2069  sourceRing=currRing;
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
2085  istate= resMatrixBase::ready;
2086}
2087
2088resMatrixDense::~resMatrixDense()
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
2120void resMatrixDense::createMatrix()
2121{
2122  int k,i,j;
2123  resVector *vecp;
2124
2125  m= mpNew( numVectors, numVectors );
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    {
2139      mprSTICKYPROT(ST_DENSE_FR);
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    {
2147      mprSTICKYPROT(ST_DENSE_NR);
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
2187void 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    {
2195      resVectorList= (resVector * )omReallocSize( resVectorList,
2196                                            (veclistmax) * sizeof( resVector ),
2197                                            (veclistmax + veclistblock) * sizeof( resVector ) );
2198      int k;
2199      for ( k= veclistmax; k < (veclistmax + veclistblock); k++ )
2200        resVectorList[k].init();
2201      veclistmax+= veclistblock;
2202      mprSTICKYPROT(ST_DENSE_MEM);
2203
2204    }
2205    resVectorList[numVectors].init( mon );
2206    numVectors++;
2207    mprSTICKYPROT(ST_DENSE_NMON);
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
2227void resMatrixDense::generateMonomData( int deg, intvec* polyDegs , intvec* iVO )
2228{
2229  int i,j,k;
2230
2231  // init monomData
2232  veclistblock= 512;
2233  veclistmax= veclistblock;
2234  resVectorList= (resVector *)omAlloc( veclistmax*sizeof( resVector ) );
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  }
2322  subSize= numVectors - subSize;
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
2343void resMatrixDense::generateBaseData()
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
2394      resVectorList[k].numColParNr= NULL;
2395      resVectorList[k].numColVectorSize= numVectors;
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);
2428      resVectorList[k].numColVectorSize= 0;
2429      resVectorList[k].numColVector= NULL;
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
2445        resVectorList[k].numColParNr[j]= i;
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
2463resVector *resMatrixDense::getMVector(const int i)
2464{
2465  assume( i >= 0 && i < numVectors );
2466  return &resVectorList[i];
2467}
2468
2469ideal resMatrixDense::getMatrix()
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
2519ideal resMatrixDense::getSubMatrix()
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
2550number 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
2570  mprSTICKYPROT(ST__DET);
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
2588  mprSTICKYPROT(ST__DET);
2589
2590  return( numres );
2591}
2592
2593number resMatrixDense::getSubDet()
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//
2659unsigned 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::*
2685uResultant::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
2710uResultant::~uResultant( )
2711{
2712  delete resMat;
2713}
2714
2715ideal 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
2743poly uResultant::linearPoly( const resMatType rrmt )
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
2770poly 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
2827    mprSTICKYPROT(ST_BASE_EV);
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
2922rootContainer ** 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
3014      mprSTICKYPROT(ST_BASE_EV);
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,
3046                                (matchUp?rootContainer::cspecialmu:rootContainer::cspecial),
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
3060rootContainer ** 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,
3160                                (matchUp?rootContainer::cspecialmu:rootContainer::cspecial),
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
3173int 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(...)
3191ideal 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
Note: See TracBrowser for help on using the repository browser.