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