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