source: git/kernel/mpr_base.cc @ d0f98e

spielwiese
Last change on this file since d0f98e was e2ad5d, checked in by Hans Schönemann <hannes@…>, 17 years ago
*hannes: format git-svn-id: file:///usr/local/Singular/svn/trunk@10095 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 73.9 KB
Line 
1/****************************************
2*  Computer Algebra System SINGULAR     *
3****************************************/
4/* $Id: mpr_base.cc,v 1.6 2007-05-24 17:46:04 Singular Exp $ */
5
6/*
7 * ABSTRACT - multipolynomial resultants - resultant matrices
8 *            ( sparse, dense, u-resultant solver )
9 */
10
11#include <math.h>
12
13#include "mod2.h"
14
15#include <mylimits.h>
16#include "omalloc.h"
17
18//-> includes
19#include "structs.h"
20#include "polys.h"
21#include "ideals.h"
22#include "febase.h"
23#include "intvec.h"
24#include "matpol.h"
25#include "numbers.h"
26#ifdef HAVE_FACTORY
27#include "clapsing.h"
28#endif
29#include "sparsmat.h"
30
31#include "mpr_global.h"
32#include "mpr_base.h"
33#include "mpr_numeric.h"
34//<-
35
36extern void nPrint(number n);  // for debugging output
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  const 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  const number getDetAt( const number* evpoint );
82
83  const 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 const 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(pVariables), pLP(_pLP) {}
281  ~mayanPyramidAlg() {}
282
283  /** Drive Mayan Pyramid Algorithm.
284   * The Alg computes conv(Qi[]+shift[]).
285   */
286  pointSet * getInnerPoints( pointSet **_Qi, 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 occured.
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= nInt(pGetCoeff( MATELEM( omat, i, j) ));
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( (ADDRESS) points[i]->point, fdim * sizeof(Coord_t) );
432    omFreeSize( (ADDRESS) points[i], sizeof(onePoint) );
433  }
434  omFreeSize( (ADDRESS) points, (max+1) * sizeof(onePointP) );
435}
436
437inline const onePointP pointSet::operator[] ( const int index )
438{
439  assume( index > 0 && index <= num );
440  return points[index];
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    pGetExpV( piter, vert );
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( (ADDRESS) 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  pGetExpV( p, vert );
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( (ADDRESS) 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,j;
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( (ADDRESS) 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; j++ ) iter= 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= pVariables;
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( pVariables, 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        pGetExpV( p, vert );
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( (ADDRESS) 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] , pVariables );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,pd;
841  int * vert;
842
843  n= pVariables;
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( (ADDRESS) 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 **_Qi, mprfloat _shift[] )
892{
893  int i;
894
895  Qi= _Qi;
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, 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[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[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,n= pVariables;
1221  int loffset= 0;
1222  for ( i= 0; i <= n; 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 set,pnt,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 faild!");
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( (ADDRESS) 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,j,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((pVariables+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 choosen.
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( (ADDRESS) epp_mon, (n+2) * sizeof(int) );
1482  omFreeSize( (ADDRESS) eexp, (pVariables+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  time_t *tp = NULL;
1506
1507  while ( i <= dim )
1508  {
1509    shift[i]= (mprfloat) (RVMULT*(siRand()%MAXRVVAL)/(mprfloat)MAXRVVAL);
1510    i++;
1511    for ( j= 1; j < i-1; j++ )
1512    {
1513      if ( (shift[j] < shift[i-1] + SIMPLEX_EPS) && (shift[j] > shift[i-1] - SIMPLEX_EPS) )
1514      {
1515        i--;
1516        break;
1517      }
1518    }
1519  }
1520}
1521
1522pointSet * resMatrixSparse::minkSumTwo( pointSet *Q1, pointSet *Q2, int dim )
1523{
1524  pointSet *vs;
1525  onePoint vert;
1526  int j,k,l;
1527
1528  vert.point=(Coord_t*)omAlloc( (pVariables+2) * sizeof(Coord_t) );
1529
1530  vs= new pointSet( dim );
1531
1532  for ( j= 1; j <= Q1->num; j++ )
1533  {
1534    for ( k= 1; k <= Q2->num; k++ )
1535    {
1536      for ( l= 1; l <= dim; l++ )
1537      {
1538        vert.point[l]= (*Q1)[j]->point[l] + (*Q2)[k]->point[l];
1539      }
1540      vs->mergeWithExp( &vert );
1541      //vs->addPoint( &vert );
1542    }
1543  }
1544
1545  omFreeSize( (ADDRESS) vert.point, (pVariables+2) * sizeof(Coord_t) );
1546
1547  return vs;
1548}
1549
1550pointSet * resMatrixSparse::minkSumAll( pointSet **pQ, int numq, int dim )
1551{
1552  pointSet *vs,*vs_old;
1553  int j;
1554
1555  vs= new pointSet( dim );
1556
1557  for ( j= 1; j <= pQ[0]->num; j++ ) vs->addPoint( (*pQ[0])[j] );
1558
1559  for ( j= 1; j < numq; j++ )
1560  {
1561    vs_old= vs;
1562    vs= minkSumTwo( vs_old, pQ[j], dim );
1563
1564    delete vs_old;
1565  }
1566
1567  return vs;
1568}
1569
1570//----------------------------------------------------------------------------------------
1571
1572resMatrixSparse::resMatrixSparse( const ideal _gls, const int special )
1573  : resMatrixBase(), gls( _gls )
1574{
1575  pointSet **Qi; // vertices sets of Conv(Supp(f_i)), i=0..idelem
1576  pointSet *E;   // all integer lattice points of the minkowski sum of Q0...Qn
1577  int i,j,k;
1578  int pnt;
1579  int totverts;                // total number of exponent vectors in ideal gls
1580  mprfloat shift[MAXVARS+2];   // shiftvector delta, index [1..dim]
1581
1582  if ( pVariables > MAXVARS )
1583  {
1584    WerrorS("resMatrixSparse::resMatrixSparse: Too many variables!");
1585    return;
1586  }
1587
1588  rmat= NULL;
1589  numSet0= 0;
1590
1591  if ( special == SNONE ) linPolyS= 0;
1592  else linPolyS= special;
1593
1594  istate= resMatrixBase::ready;
1595
1596  n= pVariables;
1597  idelem= IDELEMS(gls);  // should be n+1
1598
1599  // prepare matrix LP->LiPM for Linear Programming
1600  totverts = 0;
1601  for( i=0; i < idelem; i++) totverts += pLength( (gls->m)[i] );
1602
1603  LP = new simplex( idelem+totverts*2+5, totverts+5 ); // rows, cols
1604
1605  // get shift vector
1606#ifdef mprTEST
1607  shift[0]=0.005; shift[1]=0.003; shift[2]=0.008; shift[3]=0.005; shift[4]=0.002;
1608  shift[5]=0.1; shift[6]=0.3; shift[7]=0.2; shift[8]=0.4; shift[9]=0.2;
1609#else
1610  randomVector( idelem, shift );
1611#endif
1612#ifdef mprDEBUG_PROT
1613  PrintS(" shift vector: ");
1614  for ( i= 1; i <= idelem; i++ ) Print(" %.12f ",(double)shift[i]);
1615  PrintLn();
1616#endif
1617
1618  // evaluate convex hull for supports of gls
1619  convexHull chnp( LP );
1620  Qi= chnp.newtonPolytopesP( gls );
1621
1622#ifdef mprMINKSUM
1623  E= minkSumAll( Qi, n+1, n);
1624#else
1625  // get inner points
1626  mayanPyramidAlg mpa( LP );
1627  E= mpa.getInnerPoints( Qi, shift );
1628#endif
1629
1630#ifdef mprDEBUG_PROT
1631#ifdef mprMINKSUM
1632  PrintS("(MinkSum)");
1633#endif
1634  PrintS("\n E = (Q_0 + ... + Q_n) \\cap \\N :\n");
1635  for ( pnt= 1; pnt <= E->num; pnt++ )
1636  {
1637    Print("%d: <",pnt);print_exp( (*E)[pnt], E->dim );PrintS(">\n");
1638  }
1639  PrintLn();
1640#endif
1641
1642#ifdef mprTEST
1643  int lift[5][5];
1644  lift[0][1]=3; lift[0][2]=4; lift[0][3]=8;  lift[0][4]=2;
1645  lift[1][1]=6; lift[1][2]=1; lift[1][3]=7;  lift[1][4]=4;
1646  lift[2][1]=2; lift[2][2]=5; lift[2][3]=9;  lift[2][4]=6;
1647  lift[3][1]=2; lift[3][2]=1; lift[3][3]=9;  lift[3][4]=5;
1648  lift[4][1]=3; lift[4][2]=7; lift[4][3]=1;  lift[4][4]=5;
1649  // now lift everything
1650  for ( i= 0; i <= n; i++ ) Qi[i]->lift( lift[i] );
1651#else
1652  // now lift everything
1653  for ( i= 0; i <= n; i++ ) Qi[i]->lift();
1654#endif
1655  E->dim++;
1656
1657  // run Row Content Function for every point in E
1658  for ( pnt= 1; pnt <= E->num; pnt++ )
1659  {
1660    RC( Qi, E, pnt, shift );
1661  }
1662
1663  // remove points not in cells
1664  k= E->num;
1665  for ( pnt= k; pnt > 0; pnt-- )
1666  {
1667    if ( (*E)[pnt]->rcPnt == NULL )
1668    {
1669      E->removePoint(pnt);
1670      mprSTICKYPROT(ST_SPARSE_RCRJ);
1671    }
1672  }
1673  mprSTICKYPROT("\n");
1674
1675#ifdef mprDEBUG_PROT
1676  PrintS(" points which lie in a cell:\n");
1677  for ( pnt= 1; pnt <= E->num; pnt++ )
1678  {
1679    Print("%d: <",pnt);print_exp( (*E)[pnt], E->dim );PrintS(">\n");
1680  }
1681  PrintLn();
1682#endif
1683
1684  // unlift to old dimension, sort
1685  for ( i= 0; i <= n; i++ ) Qi[i]->unlift();
1686  E->unlift();
1687  E->sort();
1688
1689#ifdef mprDEBUG_PROT
1690  Print(" points with a[ij] (%d):\n",E->num);
1691  for ( pnt= 1; pnt <= E->num; pnt++ )
1692  {
1693    Print("Punkt p \\in E[%d]: <",pnt);print_exp( (*E)[pnt], E->dim );
1694    Print(">, RC(p) = (i:%d, j:%d), a[i,j] = <",(*E)[pnt]->rc.set,(*E)[pnt]->rc.pnt);
1695    //print_exp( (Qi[(*E)[pnt]->rc.set])[(*E)[pnt]->rc.pnt], E->dim );PrintS("> = <");
1696    print_exp( (*E)[pnt]->rcPnt, E->dim );PrintS(">\n");
1697  }
1698#endif
1699
1700  // now create matrix
1701  if (E->num <1)
1702  {
1703    WerrorS("could not handle a degenerate situation: no inner points found");
1704    goto theEnd;
1705  }
1706  if ( createMatrix( E ) != E->num )
1707  {
1708    // this can happen if the shiftvector shift is to large or not generic
1709    istate= resMatrixBase::fatalError;
1710    WerrorS("resMatrixSparse::resMatrixSparse: Error in resMatrixSparse::createMatrix!");
1711    goto theEnd;
1712  }
1713
1714 theEnd:
1715  // clean up
1716  for ( i= 0; i < idelem; i++ )
1717  {
1718    delete Qi[i];
1719  }
1720  omFreeSize( (ADDRESS) Qi, idelem * sizeof(pointSet*) );
1721
1722  delete E;
1723
1724  delete LP;
1725}
1726
1727//----------------------------------------------------------------------------------------
1728
1729resMatrixSparse::~resMatrixSparse()
1730{
1731  delete uRPos;
1732  idDelete( &rmat );
1733}
1734
1735const ideal resMatrixSparse::getMatrix()
1736{
1737  int i,j,rp,cp;
1738  poly pp,phelp,piter,pgls;
1739
1740  // copy original sparse res matrix
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
1796const number resMatrixSparse::getDetAt( const number* evpoint )
1797{
1798  int i,cp,rp;
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= smCallDet( rmat );
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
1856const poly resMatrixSparse::getUDet( const number* evpoint )
1857{
1858  int i,cp,rp;
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= smCallDet( rmat );
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//
1926class 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  const ideal getMatrix();
1946
1947  /** Returns the submatrix M' of M in an usable presentation */
1948  const 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  const 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  const 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   * pVariables 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 pVariables
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    omfreeSize( (ADDRESS)resVectorList[i].numColVector,
2101                numVectors * sizeof( number ) );
2102    omfreeSize( (ADDRESS)resVectorList[i].numColParNr,
2103                (pVariables+1) * sizeof(int) );
2104  }
2105
2106  omFreeSize( (ADDRESS)resVectorList, veclistmax*sizeof( resVector ) );
2107
2108  // free matrix m
2109  if ( m != NULL )
2110  {
2111    idDelete((ideal *)&m);
2112  }
2113}
2114
2115// mprSTICKYPROT:
2116// ST_DENSE_FR: found row S0
2117// ST_DENSE_NR: normal row
2118void resMatrixDense::createMatrix()
2119{
2120  int k,i,j;
2121  resVector *vecp;
2122
2123  m= mpNew( numVectors, numVectors );
2124
2125  for ( i= 1; i <= MATROWS( m ); i++ )
2126    for ( j= 1; j <= MATCOLS( m ); j++ )
2127    {
2128      MATELEM(m,i,j)= pInit();
2129      pSetCoeff0( MATELEM(m,i,j), nInit(0) );
2130    }
2131
2132
2133  for ( k= 0; k <= numVectors - 1; k++ )
2134  {
2135    if ( linPolyS == getMVector(k)->elementOfS )
2136    {
2137      mprSTICKYPROT(ST_DENSE_FR);
2138      for ( i= 0; i < pVariables; i++ )
2139      {
2140        MATELEM(m,numVectors-k,numVectors-(getMVector(k)->numColParNr)[i])= pInit();
2141      }
2142    }
2143    else
2144    {
2145      mprSTICKYPROT(ST_DENSE_NR);
2146      vecp= getMVector(k);
2147      for ( i= 0; i < numVectors; i++)
2148      {
2149        if ( !nIsZero( vecp->getElemNum(i) ) )
2150        {
2151          MATELEM(m,numVectors - k,i + 1)= pInit();
2152          pSetCoeff0( MATELEM(m,numVectors - k,i + 1), nCopy(vecp->getElemNum(i)) );
2153        }
2154      }
2155    }
2156  } // for
2157  mprSTICKYPROT("\n");
2158
2159#ifdef mprDEBUG_ALL
2160  for ( k= numVectors - 1; k >= 0; k-- )
2161  {
2162    if ( linPolyS == getMVector(k)->elementOfS )
2163    {
2164      for ( i=0; i < pVariables; i++ )
2165      {
2166        Print(" %d ",(getMVector(k)->numColParNr)[i]);
2167      }
2168      PrintLn();
2169    }
2170  }
2171  for (i=1; i <= numVectors; i++)
2172  {
2173    for (j=1; j <= numVectors; j++ )
2174    {
2175      pWrite0(MATELEM(m,i,j));PrintS("  ");
2176    }
2177    PrintLn();
2178  }
2179#endif
2180}
2181
2182// mprSTICKYPROT:
2183// ST_DENSE_MEM: more mem allocated
2184// ST_DENSE_NMON: new monom added
2185void resMatrixDense::generateMonoms( poly m, int var, int deg )
2186{
2187  if ( deg == 0 )
2188  {
2189    poly mon = pCopy( m );
2190
2191    if ( numVectors == veclistmax )
2192    {
2193      resVectorList= (resVector * )omReallocSize( resVectorList,
2194                                            (veclistmax) * sizeof( resVector ),
2195                                            (veclistmax + veclistblock) * sizeof( resVector ) );
2196      int k;
2197      for ( k= veclistmax; k < (veclistmax + veclistblock); k++ )
2198        resVectorList[k].init();
2199      veclistmax+= veclistblock;
2200      mprSTICKYPROT(ST_DENSE_MEM);
2201
2202    }
2203    resVectorList[numVectors].init( mon );
2204    numVectors++;
2205    mprSTICKYPROT(ST_DENSE_NMON);
2206    return;
2207  }
2208  else
2209  {
2210    if ( var == pVariables+1 ) return;
2211    poly newm = pCopy( m );
2212    while ( deg >= 0 )
2213    {
2214      generateMonoms( newm, var+1, deg );
2215      pIncrExp( newm, var );
2216      pSetm( newm );
2217      deg--;
2218    }
2219    pDelete( & newm );
2220  }
2221
2222  return;
2223}
2224
2225void resMatrixDense::generateMonomData( int deg, intvec* polyDegs , intvec* iVO )
2226{
2227  int i,j,k;
2228
2229  // init monomData
2230  veclistblock= 512;
2231  veclistmax= veclistblock;
2232  resVectorList= (resVector *)omAlloc( veclistmax*sizeof( resVector ) );
2233
2234  // Init resVector()s
2235  for ( j= veclistmax - 1; j >= 0; j-- ) resVectorList[j].init();
2236  numVectors= 0;
2237
2238  // Generate all monoms of degree deg
2239  poly start= pOne();
2240  generateMonoms( start, 1, deg );
2241  pDelete( & start );
2242
2243  mprSTICKYPROT("\n");
2244
2245  // Check for reduced monoms
2246  // First generate polyDegs.rows() monoms
2247  //  x(k)^(polyDegs[k]),  0 <= k < polyDegs.rows()
2248  ideal pDegDiv= idInit( polyDegs->rows(), 1 );
2249  for ( k= 0; k < polyDegs->rows(); k++ )
2250  {
2251    poly p= pOne();
2252    pSetExp( p, k + 1, (*polyDegs)[k] );
2253    pSetm( p );
2254    (pDegDiv->m)[k]= p;
2255  }
2256
2257  // Now check each monom if it is reduced.
2258  // A monom monom is called reduced if there exists
2259  // exactly one x(k)^(polyDegs[k]) that divides the monom.
2260  int divCount;
2261  for ( j= numVectors - 1; j >= 0; j-- )
2262  {
2263    divCount= 0;
2264    for ( k= 0; k < IDELEMS(pDegDiv); k++ )
2265      if ( pLmDivisibleByNoComp( (pDegDiv->m)[k], resVectorList[j].mon ) )
2266        divCount++;
2267    resVectorList[j].isReduced= (divCount == 1);
2268  }
2269
2270  // create the sets S(k)s
2271  // a monom x(i)^deg, deg given, is element of the set S(i)
2272  // if all x(0)^(polyDegs[0]) ... x(i-1)^(polyDegs[i-1]) DONT divide
2273  // x(i)^deg and only x(i)^(polyDegs[i]) divides x(i)^deg
2274  bool doInsert;
2275  for ( k= 0; k < iVO->rows(); k++)
2276  {
2277    //mprPROTInl(" ------------ var:",(*iVO)[k]);
2278    for ( j= numVectors - 1; j >= 0; j-- )
2279    {
2280      //mprPROTPnl("testing monom",resVectorList[j].mon);
2281      if ( resVectorList[j].elementOfS == SFREE )
2282      {
2283        //mprPROTnl("\tfree");
2284        if ( pLmDivisibleByNoComp( (pDegDiv->m)[ (*iVO)[k] ], resVectorList[j].mon ) )
2285        {
2286          //mprPROTPnl("\tdivisible by ",(pDegDiv->m)[ (*iVO)[k] ]);
2287          doInsert=TRUE;
2288          for ( i= 0; i < k; i++ )
2289          {
2290            //mprPROTPnl("\tchecking db ",(pDegDiv->m)[ (*iVO)[i] ]);
2291            if ( pLmDivisibleByNoComp( (pDegDiv->m)[ (*iVO)[i] ], resVectorList[j].mon ) )
2292            {
2293              //mprPROTPnl("\t and divisible by",(pDegDiv->m)[ (*iVO)[i] ]);
2294              doInsert=FALSE;
2295              break;
2296            }
2297          }
2298          if ( doInsert )
2299          {
2300            //mprPROTInl("\t------------------> S ",(*iVO)[k]);
2301            resVectorList[j].elementOfS= (*iVO)[k];
2302            resVectorList[j].dividedBy= pCopy( (pDegDiv->m)[ (*iVO)[i] ] );
2303          }
2304        }
2305      }
2306    }
2307  }
2308
2309  // size of submatrix M', equal to number of nonreduced monoms
2310  // (size of matrix M is equal to number of monoms=numVectors)
2311  subSize= 0;
2312  int sub;
2313  for ( i= 0; i < polyDegs->rows(); i++ )
2314  {
2315    sub= 1;
2316    for ( k= 0; k < polyDegs->rows(); k++ )
2317      if ( i != k ) sub*= (*polyDegs)[k];
2318    subSize+= sub;
2319  }
2320  subSize= numVectors - subSize;
2321
2322  // pDegDiv wieder freigeben!
2323  idDelete( &pDegDiv );
2324
2325#ifdef mprDEBUG_ALL
2326  // Print a list of monoms and their properties
2327  PrintS("// \n");
2328  for ( j= numVectors - 1; j >= 0; j-- )
2329  {
2330    Print("// %s, S(%d),  db ",
2331          resVectorList[j].isReduced?"reduced":"nonreduced",
2332          resVectorList[j].elementOfS);
2333    pWrite0(resVectorList[j].dividedBy);
2334    PrintS("  monom ");
2335    pWrite(resVectorList[j].mon);
2336  }
2337  Print("// size: %d, subSize: %d\n",numVectors,subSize);
2338#endif
2339}
2340
2341void resMatrixDense::generateBaseData()
2342{
2343  int k,j,i;
2344  number matEntry;
2345  poly pmatchPos;
2346  poly pi,factor,pmp;
2347
2348  // holds the degrees of F0, F1, ..., Fn
2349  intvec polyDegs( IDELEMS(gls) );
2350  for ( k= 0; k < IDELEMS(gls); k++ )
2351    polyDegs[k]= pTotaldegree( (gls->m)[k] );
2352
2353  // the internal Variable Ordering
2354  // make sure that the homogenization variable goes last!
2355  intvec iVO( pVariables );
2356  if ( linPolyS != SNONE )
2357  {
2358    iVO[pVariables - 1]= linPolyS;
2359    int p=0;
2360    for ( k= pVariables - 1; k >= 0; k-- )
2361    {
2362      if ( k != linPolyS )
2363      {
2364        iVO[p]= k;
2365        p++;
2366      }
2367    }
2368  }
2369  else
2370  {
2371    linPolyS= 0;
2372    for ( k= 0; k < pVariables; k++ )
2373      iVO[k]= pVariables - k - 1;
2374  }
2375
2376  // the critical degree d= sum( deg(Fi) ) - n
2377  int sumDeg= 0;
2378  for ( k= 0; k < polyDegs.rows(); k++ )
2379    sumDeg+= polyDegs[k];
2380  sumDeg-= polyDegs.rows() - 1;
2381
2382  // generate the base data
2383  generateMonomData( sumDeg, &polyDegs, &iVO );
2384
2385  // generate "matrix"
2386  for ( k= numVectors - 1; k >= 0; k-- )
2387  {
2388    if ( resVectorList[k].elementOfS != linPolyS )
2389    {
2390      // column k is a normal column with numerical or symbolic entries
2391      // init stuff
2392      resVectorList[k].numColParNr= NULL;
2393      resVectorList[k].numColVectorSize= numVectors;
2394      resVectorList[k].numColVector= (number *)omAlloc( numVectors*sizeof( number ) );
2395      for ( i= 0; i < numVectors; i++ ) resVectorList[k].numColVector[i]= nInit(0);
2396
2397      // compute row poly
2398      poly pi= ppMult_qq( (gls->m)[ resVectorList[k].elementOfS ] , resVectorList[k].mon );
2399      pi= pDivideM( pCopy( pi ), pCopy( resVectorList[k].dividedBy ) );
2400
2401      // fill in "matrix"
2402      while ( pi != NULL )
2403      {
2404        matEntry= nCopy(pGetCoeff(pi));
2405        pmatchPos= pLmInit( pi );
2406        pSetCoeff0( pmatchPos, nInit(1) );
2407
2408        for ( i= 0; i < numVectors; i++)
2409          if ( pLmEqual( pmatchPos, resVectorList[i].mon ) )
2410            break;
2411
2412        resVectorList[k].numColVector[numVectors - i - 1] = nCopy(matEntry);
2413
2414        pDelete( &pmatchPos );
2415        nDelete( &matEntry );
2416
2417        pIter( pi );
2418      }
2419      pDelete( &pi );
2420    }
2421    else
2422    {
2423      // column is a special column, i.e. is generated by S0 and F0
2424      // safe only the positions of the ui's in the column
2425      //mprPROTInl(" setup of numColParNr ",k);
2426      resVectorList[k].numColVectorSize= 0;
2427      resVectorList[k].numColVector= NULL;
2428      resVectorList[k].numColParNr= (int *)omAlloc0( (pVariables+1) * sizeof(int) );
2429
2430      pi= (gls->m)[ resVectorList[k].elementOfS ];
2431      factor= pDivideM( pCopy( resVectorList[k].mon ), pCopy( resVectorList[k].dividedBy ) );
2432
2433      j=0;
2434      while ( pi  != NULL )
2435      { // fill in "matrix"
2436        pmp= pMult( pCopy( factor ), pHead( pi ) );
2437        pTest( pmp );
2438
2439        for ( i= 0; i < numVectors; i++)
2440          if ( pLmEqual( pmp, resVectorList[i].mon ) )
2441            break;
2442
2443        resVectorList[k].numColParNr[j]= i;
2444        pDelete( &pmp );
2445        pIter( pi );
2446        j++;
2447      }
2448      pDelete( &pi );
2449      pDelete( &factor );
2450    }
2451  } // for ( k= numVectors - 1; k >= 0; k-- )
2452
2453  mprSTICKYPROT2(" size of matrix:    %d\n",numVectors);
2454  mprSTICKYPROT2(" size of submatrix: %d\n",subSize);
2455
2456  // create the matrix M
2457  createMatrix();
2458
2459}
2460
2461resVector *resMatrixDense::getMVector(const int i)
2462{
2463  assume( i >= 0 && i < numVectors );
2464  return &resVectorList[i];
2465}
2466
2467const ideal resMatrixDense::getMatrix()
2468{
2469  int k,i,j;
2470
2471  // copy matrix
2472  matrix resmat= mpNew(numVectors,numVectors);
2473  poly p;
2474  for (i=1; i <= numVectors; i++)
2475  {
2476    for (j=1; j <= numVectors; j++ )
2477    {
2478      p=MATELEM(m,i,j);
2479      if (( p!=NULL)
2480      && (!nIsZero(pGetCoeff(p)))
2481      && (pGetCoeff(p)!=NULL)
2482      )
2483      {
2484        MATELEM(resmat,i,j)= pCopy( p );
2485      }
2486    }
2487  }
2488  for (i=0; i < numVectors; i++)
2489  {
2490    if ( resVectorList[i].elementOfS == linPolyS )
2491    {
2492      for (j=1; j <= pVariables; j++ )
2493      {
2494        if ( MATELEM(resmat,numVectors-i,
2495                     numVectors-resVectorList[i].numColParNr[j-1])!=NULL )
2496          pDelete( &MATELEM(resmat,numVectors-i,numVectors-resVectorList[i].numColParNr[j-1]) );
2497        MATELEM(resmat,numVectors-i,numVectors-resVectorList[i].numColParNr[j-1])= pOne();
2498        // FIX ME
2499        if ( FALSE )
2500        {
2501          pSetCoeff( MATELEM(resmat,numVectors-i,numVectors-resVectorList[i].numColParNr[j-1]), nPar(j) );
2502        }
2503        else
2504        {
2505          pSetExp( MATELEM(resmat,numVectors-i,numVectors-resVectorList[i].numColParNr[j-1]), j, 1 );
2506          pSetm(MATELEM(resmat,numVectors-i,numVectors-resVectorList[i].numColParNr[j-1]));
2507        }
2508      }
2509    }
2510  }
2511
2512  // obachman: idMatrix2Module frees resmat !!
2513  ideal resmod= idMatrix2Module(resmat);
2514  return resmod;
2515}
2516
2517const ideal resMatrixDense::getSubMatrix()
2518{
2519  int k,i,j,l;
2520  resVector *vecp;
2521
2522  // generate quadratic matrix resmat of size subSize
2523  matrix resmat= mpNew( subSize, subSize );
2524
2525  j=1;
2526  for ( k= numVectors - 1; k >= 0; k-- )
2527  {
2528    vecp= getMVector(k);
2529    if ( vecp->isReduced ) continue;
2530    l=1;
2531    for ( i= numVectors - 1; i >= 0; i-- )
2532    {
2533      if ( getMVector(i)->isReduced ) continue;
2534      if ( !nIsZero(vecp->getElemNum(numVectors - i - 1)) )
2535      {
2536        MATELEM(resmat,j,l)= pCopy( vecp->getElem(numVectors-i-1) );
2537      }
2538      l++;
2539    }
2540    j++;
2541  }
2542
2543  // obachman: idMatrix2Module frees resmat !!
2544  ideal resmod= idMatrix2Module(resmat);
2545  return resmod;
2546}
2547
2548const number resMatrixDense::getDetAt( const number* evpoint )
2549{
2550  int k,i,j;
2551
2552  // copy evaluation point into matrix
2553  // p0, p1, ..., pn replace u0, u1, ..., un
2554  for ( k= numVectors - 1; k >= 0; k-- )
2555  {
2556    if ( linPolyS == getMVector(k)->elementOfS )
2557    {
2558      for ( i= 0; i < pVariables; i++ )
2559      {
2560        pSetCoeff( MATELEM(m,numVectors-k,numVectors-(getMVector(k)->numColParNr)[i]),
2561                   nCopy(evpoint[i]) );
2562      }
2563    }
2564  }
2565
2566  mprSTICKYPROT(ST__DET);
2567
2568  // evaluate determinant of matrix m using factory singclap_det
2569#ifdef HAVE_FACTORY
2570  poly res= singclap_det( m );
2571#else
2572  poly res= NULL;
2573#endif
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
2593const number 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#ifdef HAVE_FACTORY
2634  poly res= singclap_det( mat );
2635#else
2636  poly res= NULL;
2637#endif
2638
2639  number numres;
2640  if ((res != NULL) && (!nIsZero(pGetCoeff( res ))) )
2641  {
2642    numres= nCopy(pGetCoeff( res ));
2643  }
2644  else
2645  {
2646    numres= nInit(0);
2647  }
2648  pDelete( &res );
2649  return numres;
2650}
2651//<--
2652
2653//-----------------------------------------------------------------------------
2654//-------------- uResultant ---------------------------------------------------
2655//-----------------------------------------------------------------------------
2656
2657#define MAXEVPOINT 1000000 // 0x7fffffff
2658//#define MPR_MASI
2659
2660//-> unsigned long over(unsigned long n,unsigned long d)
2661// Calculates (n+d \over d) using gmp functionality
2662//
2663unsigned long over( const unsigned long n , const unsigned long d )
2664{ // (d+n)! / ( d! n! )
2665  mpz_t res;
2666  mpz_init(res);
2667  mpz_t m,md,mn;
2668  mpz_init(m);mpz_set_ui(m,1);
2669  mpz_init(md);mpz_set_ui(md,1);
2670  mpz_init(mn);mpz_set_ui(mn,1);
2671
2672  mpz_fac_ui(m,n+d);
2673  mpz_fac_ui(md,d);
2674  mpz_fac_ui(mn,n);
2675
2676  mpz_mul(res,md,mn);
2677  mpz_tdiv_q(res,m,res);
2678
2679  mpz_clear(m);mpz_clear(md);mpz_clear(mn);
2680
2681  unsigned long result = mpz_get_ui(res);
2682  mpz_clear(res);
2683
2684  return result;
2685}
2686//<-
2687
2688//-> uResultant::*
2689uResultant::uResultant( const ideal _gls, const resMatType _rmt, BOOLEAN extIdeal )
2690  : rmt( _rmt )
2691{
2692  if ( extIdeal )
2693  {
2694    // extend given ideal by linear poly F0=u0x0 + u1x1 +...+ unxn
2695    gls= extendIdeal( _gls, linearPoly( rmt ), rmt );
2696    n= IDELEMS( gls );
2697  }
2698  else
2699    gls= idCopy( _gls );
2700
2701  switch ( rmt )
2702  {
2703  case sparseResMat:
2704    resMat= new resMatrixSparse( gls );
2705    break;
2706  case denseResMat:
2707#ifdef HAVE_FACTORY
2708    resMat= new resMatrixDense( gls );
2709    break;
2710#endif
2711  default:
2712    WerrorS("uResultant::uResultant: Unknown resultant matrix type choosen!");
2713  }
2714}
2715
2716uResultant::~uResultant( )
2717{
2718  delete resMat;
2719}
2720
2721ideal uResultant::extendIdeal( const ideal gls, poly linPoly, const resMatType rmt )
2722{
2723  ideal newGls= idCopy( gls );
2724  newGls->m= (poly *)omReallocSize( newGls->m,
2725                              IDELEMS(gls) * sizeof(poly),
2726                              (IDELEMS(gls) + 1) * sizeof(poly) );
2727  IDELEMS(newGls)++;
2728
2729  switch ( rmt )
2730  {
2731  case sparseResMat:
2732  case denseResMat:
2733    {
2734      int i;
2735      for ( i= IDELEMS(newGls)-1; i > 0; i-- )
2736      {
2737        newGls->m[i]= newGls->m[i-1];
2738      }
2739      newGls->m[0]= linPoly;
2740    }
2741    break;
2742  default:
2743    WerrorS("uResultant::extendIdeal: Unknown resultant matrix type choosen!");
2744  }
2745
2746  return( newGls );
2747}
2748
2749poly uResultant::linearPoly( const resMatType rmt )
2750{
2751  int i,j;
2752
2753  poly newlp= pOne();
2754  poly actlp, rootlp= newlp;
2755
2756  for ( i= 1; i <= pVariables; i++ )
2757  {
2758    actlp= newlp;
2759    pSetExp( actlp, i, 1 );
2760    pSetm( actlp );
2761    newlp= pOne();
2762    actlp->next= newlp;
2763  }
2764  actlp->next= NULL;
2765  pDelete( &newlp );
2766
2767  if ( rmt == sparseResMat )
2768  {
2769    newlp= pOne();
2770    actlp->next= newlp;
2771    newlp->next= NULL;
2772  }
2773  return ( rootlp );
2774}
2775
2776poly uResultant::interpolateDense( const number subDetVal )
2777{
2778  int i,j,p;
2779  long tdg;
2780
2781  // D is a Polynom homogeneous in the coeffs of F0 and of degree tdg = d0*d1*...*dn
2782  tdg= resMat->getDetDeg();
2783
2784  // maximum number of terms in polynom D (homogeneous, of degree tdg)
2785  // long mdg= (facul(tdg+n-1) / facul( tdg )) / facul( n - 1 );
2786  long mdg= over( n-1, tdg );
2787
2788  // maximal number of terms in a polynom of degree tdg
2789  long l=(long)pow( (double)(tdg+1), n );
2790
2791#ifdef mprDEBUG_PROT
2792  Print("// total deg of D: tdg %d\n",tdg);
2793  Print("// maximum number of terms in D: mdg: %d\n",mdg);
2794  Print("// maximum number of terms in polynom of deg tdg: l %d\n",l);
2795#endif
2796
2797  // we need mdg results of D(p0,p1,...,pn)
2798  number *presults;
2799  presults= (number *)omAlloc( mdg * sizeof( number ) );
2800  for (i=0; i < mdg; i++) presults[i]= nInit(0);
2801
2802  number *pevpoint= (number *)omAlloc( n * sizeof( number ) );
2803  number *pev= (number *)omAlloc( n * sizeof( number ) );
2804  for (i=0; i < n; i++) pev[i]= nInit(0);
2805
2806  mprPROTnl("// initial evaluation point: ");
2807  // initial evaluatoin point
2808  p=1;
2809  for (i=0; i < n; i++)
2810  {
2811    // init pevpoint with primes 3,5,7,11, ...
2812    p= nextPrime( p );
2813    pevpoint[i]=nInit( p );
2814    nTest(pevpoint[i]);
2815    mprPROTNnl(" ",pevpoint[i]);
2816  }
2817
2818  // evaluate the determinant in the points pev^0, pev^1, ..., pev^mdg
2819  mprPROTnl("// evaluating:");
2820  for ( i=0; i < mdg; i++ )
2821  {
2822    for (j=0; j < n; j++)
2823    {
2824      nDelete( &pev[j] );
2825      nPower(pevpoint[j],i,&pev[j]);
2826      mprPROTN(" ",pev[j]);
2827    }
2828    mprPROTnl("");
2829
2830    nDelete( &presults[i] );
2831    presults[i]=resMat->getDetAt( pev );
2832
2833    mprSTICKYPROT(ST_BASE_EV);
2834  }
2835  mprSTICKYPROT("\n");
2836
2837  // now interpolate using vandermode interpolation
2838  mprPROTnl("// interpolating:");
2839  number *ncpoly;
2840  {
2841    vandermonde vm( mdg, n, tdg, pevpoint );
2842    ncpoly= vm.interpolateDense( presults );
2843  }
2844
2845  if ( subDetVal != NULL )
2846  {   // divide by common factor
2847    number detdiv;
2848    for ( i= 0; i <= mdg; i++ )
2849    {
2850      detdiv= nDiv( ncpoly[i], subDetVal );
2851      nNormalize( detdiv );
2852      nDelete( &ncpoly[i] );
2853      ncpoly[i]= detdiv;
2854    }
2855  }
2856
2857#ifdef mprDEBUG_ALL
2858  PrintLn();
2859  for ( i=0; i < mdg; i++ )
2860  {
2861    nPrint(ncpoly[i]); PrintS(" --- ");
2862  }
2863  PrintLn();
2864#endif
2865
2866  // prepare ncpoly for later use
2867  number nn=nInit(0);
2868  for ( i=0; i < mdg; i++ )
2869  {
2870    if ( nEqual(ncpoly[i],nn) )
2871    {
2872      nDelete( &ncpoly[i] );
2873      ncpoly[i]=NULL;
2874    }
2875  }
2876  nDelete( &nn );
2877
2878  // create poly presenting the determinat of the uResultant
2879  intvec exp( n );
2880  for ( i= 0; i < n; i++ ) exp[i]=0;
2881
2882  poly result= NULL;
2883
2884  long sum=0;
2885  long c=0;
2886
2887  for ( i=0; i < l; i++ )
2888  {
2889    if ( sum == tdg )
2890    {
2891      if ( !nIsZero(ncpoly[c]) )
2892      {
2893        poly p= pOne();
2894        if ( rmt == denseResMat )
2895        {
2896          for ( j= 0; j < n; j++ ) pSetExp( p, j+1, exp[j] );
2897        }
2898        else if ( rmt == sparseResMat )
2899        {
2900          for ( j= 1; j < n; j++ ) pSetExp( p, j, exp[j] );
2901        }
2902        pSetCoeff( p, ncpoly[c] );
2903        pSetm( p );
2904        if (result!=NULL) result= pAdd( result, p );
2905        else result=  p;
2906      }
2907      c++;
2908    }
2909    sum=0;
2910    exp[0]++;
2911    for ( j= 0; j < n - 1; j++ )
2912    {
2913      if ( exp[j] > tdg )
2914      {
2915        exp[j]= 0;
2916        exp[j + 1]++;
2917      }
2918      sum+=exp[j];
2919    }
2920    sum+=exp[n-1];
2921  }
2922
2923  pTest( result );
2924
2925  return result;
2926}
2927
2928rootContainer ** uResultant::interpolateDenseSP( BOOLEAN matchUp, const number subDetVal )
2929{
2930  int i,j,p,uvar;
2931  long tdg;
2932  int loops= (matchUp?n-2:n-1);
2933
2934  mprPROTnl("uResultant::interpolateDenseSP");
2935
2936  tdg= resMat->getDetDeg();
2937
2938  // evaluate D in tdg+1 distinct points, so
2939  // we need tdg+1 results of D(p0,1,0,...,0) =
2940  //              c(0)*u0^tdg + c(1)*u0^tdg-1 + ... + c(tdg-1)*u0 + c(tdg)
2941  number *presults;
2942  presults= (number *)omAlloc( (tdg + 1) * sizeof( number ) );
2943  for ( i=0; i <= tdg; i++ ) presults[i]= nInit(0);
2944
2945  rootContainer ** roots;
2946  roots= (rootContainer **) omAlloc( loops * sizeof(rootContainer*) );
2947  for ( i=0; i < loops; i++ ) roots[i]= new rootContainer(); // 0..n-2
2948
2949  number *pevpoint= (number *)omAlloc( n * sizeof( number ) );
2950  for (i=0; i < n; i++) pevpoint[i]= nInit(0);
2951
2952  number *pev= (number *)omAlloc( n * sizeof( number ) );
2953  for (i=0; i < n; i++) pev[i]= nInit(0);
2954
2955  // now we evaluate D(u0,-1,0,...0), D(u0,0,-1,0,...,0), ..., D(u0,0,..,0,-1)
2956  // or D(u0,k1,k2,0,...,0), D(u0,k1,k2,k3,0,...,0), ..., D(u0,k1,k2,k3,...,kn)
2957  // this gives us n-1 evaluations
2958  p=3;
2959  for ( uvar= 0; uvar < loops; uvar++ )
2960  {
2961    // generate initial evaluation point
2962    if ( matchUp )
2963    {
2964      for (i=0; i < n; i++)
2965      {
2966        // prime(random number) between 1 and MAXEVPOINT
2967        nDelete( &pevpoint[i] );
2968        if ( i == 0 )
2969        {
2970          //p= nextPrime( p );
2971          pevpoint[i]= nInit( p );
2972        }
2973        else if ( i <= uvar + 2 )
2974        {
2975          pevpoint[i]=nInit(1+siRand()%MAXEVPOINT);
2976          //pevpoint[i]=nInit(383);
2977        }
2978        else
2979          pevpoint[i]=nInit(0);
2980        mprPROTNnl(" ",pevpoint[i]);
2981      }
2982    }
2983    else
2984    {
2985      for (i=0; i < n; i++)
2986      {
2987        // init pevpoint with  prime,0,...0,1,0,...,0
2988        nDelete( &pevpoint[i] );
2989        if ( i == 0 )
2990        {
2991          //p=nextPrime( p );
2992          pevpoint[i]=nInit( p );
2993        }
2994        else
2995        {
2996          if ( i == (uvar + 1) ) pevpoint[i]= nInit(-1);
2997          else pevpoint[i]= nInit(0);
2998        }
2999        mprPROTNnl(" ",pevpoint[i]);
3000      }
3001    }
3002
3003    // prepare aktual evaluation point
3004    for (i=0; i < n; i++)
3005    {
3006      nDelete( &pev[i] );
3007      pev[i]= nCopy( pevpoint[i] );
3008    }
3009    // evaluate the determinant in the points pev^0, pev^1, ..., pev^tdg
3010    for ( i=0; i <= tdg; i++ )
3011    {
3012      nDelete( &pev[0] );
3013      nPower(pevpoint[0],i,&pev[0]);          // new evpoint
3014
3015      nDelete( &presults[i] );
3016      presults[i]=resMat->getDetAt( pev );   // evaluate det at point evpoint
3017
3018      mprPROTNnl("",presults[i]);
3019
3020      mprSTICKYPROT(ST_BASE_EV);
3021      mprPROTI("",tdg-i);
3022    }
3023    mprSTICKYPROT("\n");
3024
3025    // now interpolate
3026    vandermonde vm( tdg + 1, 1, tdg, pevpoint, FALSE );
3027    number *ncpoly= vm.interpolateDense( presults );
3028
3029    if ( subDetVal != NULL )
3030    {  // divide by common factor
3031      number detdiv;
3032      for ( i= 0; i <= tdg; i++ )
3033      {
3034        detdiv= nDiv( ncpoly[i], subDetVal );
3035        nNormalize( detdiv );
3036        nDelete( &ncpoly[i] );
3037        ncpoly[i]= detdiv;
3038      }
3039    }
3040
3041#ifdef mprDEBUG_ALL
3042    PrintLn();
3043    for ( i=0; i <= tdg; i++ )
3044    {
3045      nPrint(ncpoly[i]); PrintS(" --- ");
3046    }
3047    PrintLn();
3048#endif
3049
3050    // save results
3051    roots[uvar]->fillContainer( ncpoly, pevpoint, uvar+1, tdg,
3052                                (matchUp?rootContainer::cspecialmu:rootContainer::cspecial),
3053                                loops );
3054  }
3055
3056  // free some stuff: pev, presult
3057  for ( i=0; i < n; i++ ) nDelete( pev + i );
3058  omFreeSize( (ADDRESS)pev, n * sizeof( number ) );
3059
3060  for ( i=0; i <= tdg; i++ ) nDelete( presults+i );
3061  omFreeSize( (ADDRESS)presults, (tdg + 1) * sizeof( number ) );
3062
3063  return roots;
3064}
3065
3066rootContainer ** uResultant::specializeInU( BOOLEAN matchUp, const number subDetVal )
3067{
3068  int i,j,p,uvar;
3069  long tdg;
3070  poly pures,piter;
3071  int loops=(matchUp?n-2:n-1);
3072  int nn=n;
3073  if (loops==0) { loops=1;nn++;}
3074
3075  mprPROTnl("uResultant::specializeInU");
3076
3077  tdg= resMat->getDetDeg();
3078
3079  rootContainer ** roots;
3080  roots= (rootContainer **) omAlloc( loops * sizeof(rootContainer*) );
3081  for ( i=0; i < loops; i++ ) roots[i]= new rootContainer(); // 0..n-2
3082
3083  number *pevpoint= (number *)omAlloc( nn * sizeof( number ) );
3084  for (i=0; i < nn; i++) pevpoint[i]= nInit(0);
3085
3086  // now we evaluate D(u0,-1,0,...0), D(u0,0,-1,0,...,0), ..., D(u0,0,..,0,-1)
3087  // or D(u0,k1,k2,0,...,0), D(u0,k1,k2,k3,0,...,0), ..., D(u0,k1,k2,k3,...,kn)
3088  p=3;
3089  for ( uvar= 0; uvar < loops; uvar++ )
3090  {
3091    // generate initial evaluation point
3092    if ( matchUp )
3093    {
3094      for (i=0; i < n; i++)
3095      {
3096        // prime(random number) between 1 and MAXEVPOINT
3097        nDelete( &pevpoint[i] );
3098        if ( i <= uvar + 2 )
3099        {
3100          pevpoint[i]=nInit(1+siRand()%MAXEVPOINT);
3101          //pevpoint[i]=nInit(383);
3102        }
3103        else pevpoint[i]=nInit(0);
3104        mprPROTNnl(" ",pevpoint[i]);
3105      }
3106    }
3107    else
3108    {
3109      for (i=0; i < n; i++)
3110      {
3111        // init pevpoint with  prime,0,...0,-1,0,...,0
3112        nDelete( &(pevpoint[i]) );
3113        if ( i == (uvar + 1) ) pevpoint[i]= nInit(-1);
3114        else pevpoint[i]= nInit(0);
3115        mprPROTNnl(" ",pevpoint[i]);
3116      }
3117    }
3118
3119    pures= resMat->getUDet( pevpoint );
3120
3121    number *ncpoly= (number *)omAlloc( (tdg+1) * sizeof(number) );
3122
3123#ifdef MPR_MASI
3124    BOOLEAN masi=true;
3125#endif
3126
3127    piter= pures;
3128    for ( i= tdg; i >= 0; i-- )
3129    {
3130      //if ( piter ) Print("deg %d, pDeg(piter) %d\n",i,pTotaldegree(piter));
3131      if ( piter && pTotaldegree(piter) == i )
3132      {
3133        ncpoly[i]= nCopy( pGetCoeff( piter ) );
3134        pIter( piter );
3135#ifdef MPR_MASI
3136        masi=false;
3137#endif
3138      }
3139      else
3140      {
3141        ncpoly[i]= nInit(0);
3142      }
3143      mprPROTNnl("", ncpoly[i] );
3144    }
3145#ifdef MPR_MASI
3146    if ( masi ) mprSTICKYPROT("MASI");
3147#endif
3148
3149    mprSTICKYPROT(ST_BASE_EV); // .
3150
3151    if ( subDetVal != NULL )  // divide by common factor
3152    {
3153      number detdiv;
3154      for ( i= 0; i <= tdg; i++ )
3155      {
3156        detdiv= nDiv( ncpoly[i], subDetVal );
3157        nNormalize( detdiv );
3158        nDelete( &ncpoly[i] );
3159        ncpoly[i]= detdiv;
3160      }
3161    }
3162
3163    pDelete( &pures );
3164
3165    // save results
3166    roots[uvar]->fillContainer( ncpoly, pevpoint, uvar+1, tdg,
3167                                (matchUp?rootContainer::cspecialmu:rootContainer::cspecial),
3168                                loops );
3169  }
3170
3171  mprSTICKYPROT("\n");
3172
3173  // free some stuff: pev, presult
3174  for ( i=0; i < n; i++ ) nDelete( pevpoint + i );
3175  omFreeSize( (ADDRESS)pevpoint, n * sizeof( number ) );
3176
3177  return roots;
3178}
3179
3180int uResultant::nextPrime( const int i )
3181{
3182  int init=i;
3183  int ii=i+2;
3184  int j= IsPrime( ii );
3185  while ( j <= init )
3186  {
3187    ii+=2;
3188    j= IsPrime( ii );
3189  }
3190  return j;
3191}
3192//<-
3193
3194//-----------------------------------------------------------------------------
3195
3196//-> loNewtonPolytope(...)
3197ideal loNewtonPolytope( const ideal id )
3198{
3199  simplex * LP;
3200  int i;
3201  int n,totverts,idelem;
3202  ideal idr;
3203
3204  n= pVariables;
3205  idelem= IDELEMS(id);  // should be n+1
3206
3207  totverts = 0;
3208  for( i=0; i < idelem; i++) totverts += pLength( (id->m)[i] );
3209
3210  LP = new simplex( idelem+totverts*2+5, totverts+5 ); // rows, cols
3211
3212  // evaluate convex hull for supports of id
3213  convexHull chnp( LP );
3214  idr = chnp.newtonPolytopesI( id );
3215
3216  delete LP;
3217
3218  return idr;
3219}
3220//<-
3221
3222//%e
3223
3224//-----------------------------------------------------------------------------
3225
3226// local Variables: ***
3227// folded-file: t ***
3228// compile-command-1: "make installg" ***
3229// compile-command-2: "make install" ***
3230// End: ***
3231
3232// in folding: C-c x
3233// leave fold: C-c y
3234//   foldmode: F10
Note: See TracBrowser for help on using the repository browser.