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