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