source: git/numeric/mpr_base.cc @ cd4f24

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