[6ba162] | 1 | // matrix.cc |
---|
| 2 | |
---|
| 3 | // implementation of class matrix |
---|
| 4 | |
---|
| 5 | #ifndef MATRIX_CC |
---|
| 6 | #define MATRIX_CC |
---|
| 7 | |
---|
| 8 | #include "matrix.h" |
---|
| 9 | |
---|
| 10 | ////////////// constructors and destructor ////////////////////////////////// |
---|
[36f8d9] | 11 | typedef Integer* IntegerP; |
---|
| 12 | typedef BigInt* BigIntP; |
---|
[6ba162] | 13 | |
---|
| 14 | matrix::matrix(const short& row_number, const short& column_number) |
---|
| 15 | :rows(row_number),columns(column_number) |
---|
| 16 | { |
---|
| 17 | _kernel_dimension=-2; |
---|
| 18 | // LLL-algorithm not yet performed |
---|
| 19 | |
---|
| 20 | // argument check |
---|
| 21 | if((rows<=0)||(columns<=0)) |
---|
| 22 | // bad input, set "error flag" |
---|
| 23 | { |
---|
| 24 | cerr<<"\nWARNING: matrix::matrix(const short&, const short&):\n" |
---|
| 25 | "argument out of range"<<endl; |
---|
| 26 | columns=-1; |
---|
| 27 | return; |
---|
| 28 | } |
---|
| 29 | |
---|
| 30 | // memory allocation and initialization |
---|
| 31 | |
---|
[36f8d9] | 32 | coefficients=new IntegerP[rows]; |
---|
[6ba162] | 33 | for(short i=0;i<rows;i++) |
---|
| 34 | coefficients[i]=new Integer[columns]; |
---|
| 35 | for(short i=0;i<rows;i++) |
---|
| 36 | for(short j=0;j<columns;j++) |
---|
| 37 | coefficients[i][j]=0; |
---|
| 38 | } |
---|
| 39 | |
---|
| 40 | matrix::matrix(const short& row_number, const short& column_number, |
---|
| 41 | Integer** entries) |
---|
| 42 | :rows(row_number),columns(column_number) |
---|
| 43 | { |
---|
| 44 | _kernel_dimension=-2; |
---|
| 45 | // LLL-algorithm not yet performed |
---|
| 46 | |
---|
| 47 | // argument check |
---|
| 48 | if((rows<=0)||(columns<=0)) |
---|
| 49 | // bad input, set "error flag" |
---|
| 50 | { |
---|
| 51 | cerr<<"\nWARNING: matrix::matrix(const short&, const short&, Integr**):\n" |
---|
| 52 | "argument out of range"<<endl; |
---|
| 53 | columns=-1; |
---|
| 54 | return; |
---|
| 55 | } |
---|
| 56 | |
---|
| 57 | // memory allocation and initialization |
---|
| 58 | |
---|
[36f8d9] | 59 | coefficients=new IntegerP[rows]; |
---|
[6ba162] | 60 | for(short i=0;i<rows;i++) |
---|
| 61 | coefficients[i]=new Integer[columns]; |
---|
| 62 | for(short i=0;i<rows;i++) |
---|
| 63 | for(short j=0;j<columns;j++) |
---|
| 64 | coefficients[i][j]=entries[i][j]; |
---|
| 65 | // coefficients[i] is the i-th row vector |
---|
| 66 | } |
---|
| 67 | |
---|
| 68 | |
---|
| 69 | |
---|
| 70 | |
---|
| 71 | matrix::matrix(ifstream& input) |
---|
| 72 | { |
---|
| 73 | _kernel_dimension=-2; |
---|
| 74 | // LLL-algorithm not yet performed |
---|
| 75 | |
---|
| 76 | input>>rows; |
---|
| 77 | if(!input) |
---|
| 78 | // input failure, set "error flag" |
---|
| 79 | { |
---|
| 80 | cerr<<"\nWARNING: matrix::matrix(ifstream&): input failure"<<endl; |
---|
| 81 | columns=-2; |
---|
| 82 | return; |
---|
| 83 | } |
---|
| 84 | |
---|
| 85 | input>>columns; |
---|
| 86 | if(!input) |
---|
| 87 | // input failure, set "error flag" |
---|
| 88 | { |
---|
| 89 | cerr<<"\nWARNING: matrix::matrix(ifstream&): input failure"<<endl; |
---|
| 90 | columns=-2; |
---|
| 91 | return; |
---|
| 92 | } |
---|
| 93 | |
---|
| 94 | if((rows<=0)||(columns<=0)) |
---|
| 95 | // bad input, set "error flag" |
---|
| 96 | { |
---|
| 97 | cerr<<"\nWARNING: matrix::matrix(ifstream&): bad input"<<endl; |
---|
| 98 | columns=-1; |
---|
| 99 | return; |
---|
| 100 | } |
---|
| 101 | |
---|
[36f8d9] | 102 | coefficients=new IntegerP[rows]; |
---|
[6ba162] | 103 | for(short i=0;i<rows;i++) |
---|
| 104 | coefficients[i]=new Integer[columns]; |
---|
| 105 | for(short i=0;i<rows;i++) |
---|
| 106 | for(short j=0;j<columns;j++) |
---|
| 107 | { |
---|
| 108 | input>>coefficients[i][j]; |
---|
| 109 | if(!input) |
---|
| 110 | // bad input, set error flag |
---|
| 111 | { |
---|
| 112 | cerr<<"\nWARNING: matrix::matrix(ifstream&): input failure"<<endl; |
---|
| 113 | columns=-2; |
---|
| 114 | return; |
---|
| 115 | } |
---|
| 116 | } |
---|
| 117 | } |
---|
| 118 | |
---|
| 119 | |
---|
| 120 | |
---|
| 121 | |
---|
| 122 | matrix::matrix(const short& m, const short& n, ifstream& input) |
---|
| 123 | { |
---|
| 124 | _kernel_dimension=-2; |
---|
| 125 | // LLL-algorithm not yet performed |
---|
| 126 | |
---|
| 127 | // argument check |
---|
| 128 | if((m<=0) || (n<=0)) |
---|
| 129 | // bad input, set "error flag" |
---|
| 130 | { |
---|
| 131 | cerr<<"\nWARNING: matrix::matrix(const short&, const short&, ifstream&):\n" |
---|
| 132 | "argument out of range"<<endl; |
---|
| 133 | columns=-1; |
---|
| 134 | return; |
---|
| 135 | } |
---|
| 136 | |
---|
| 137 | rows=m; |
---|
| 138 | columns=n; |
---|
| 139 | |
---|
| 140 | // memory allocation and initialization |
---|
| 141 | |
---|
[36f8d9] | 142 | coefficients=new IntegerP[rows]; |
---|
[6ba162] | 143 | for(short i=0;i<rows;i++) |
---|
| 144 | coefficients[i]=new Integer[columns]; |
---|
| 145 | for(short i=0;i<rows;i++) |
---|
| 146 | for(short j=0;j<columns;j++) |
---|
| 147 | { |
---|
| 148 | input>>coefficients[i][j]; |
---|
| 149 | if(!input) |
---|
| 150 | // bad input, set error flag |
---|
| 151 | { |
---|
| 152 | columns=-2; |
---|
| 153 | return; |
---|
| 154 | } |
---|
| 155 | } |
---|
| 156 | } |
---|
| 157 | |
---|
| 158 | |
---|
| 159 | |
---|
| 160 | |
---|
| 161 | matrix::matrix(const matrix& A) |
---|
| 162 | :rows(A.rows),columns(A.columns),_kernel_dimension(A._kernel_dimension) |
---|
| 163 | { |
---|
| 164 | |
---|
| 165 | if(columns<0) |
---|
| 166 | { |
---|
| 167 | cerr<<"\nWARNING: matrix::matrix(const matrix&):\n" |
---|
| 168 | "Building a matrix from a corrupt one"<<endl; |
---|
| 169 | return; |
---|
| 170 | } |
---|
| 171 | |
---|
| 172 | // memory allocation and initialization (also for H) |
---|
| 173 | |
---|
[36f8d9] | 174 | coefficients=new IntegerP[rows]; |
---|
[6ba162] | 175 | for(short i=0;i<rows;i++) |
---|
| 176 | coefficients[i]=new Integer[columns]; |
---|
| 177 | for(short i=0;i<rows;i++) |
---|
| 178 | for(short j=0;j<columns;j++) |
---|
| 179 | coefficients[i][j]=A.coefficients[i][j]; |
---|
| 180 | |
---|
| 181 | if(_kernel_dimension>0) |
---|
| 182 | { |
---|
[36f8d9] | 183 | H=new BigIntP[_kernel_dimension]; |
---|
[6ba162] | 184 | for(short k=0;k<_kernel_dimension;k++) |
---|
[ae9ccf] | 185 | H[k]=new BigInt[columns]; |
---|
[6ba162] | 186 | for(short k=0;k<_kernel_dimension;k++) |
---|
| 187 | for(short j=0;j<columns;j++) |
---|
| 188 | H[k][j]=(A.H)[k][j]; |
---|
| 189 | } |
---|
| 190 | } |
---|
| 191 | |
---|
| 192 | |
---|
| 193 | |
---|
| 194 | |
---|
| 195 | matrix::~matrix() |
---|
| 196 | { |
---|
| 197 | for(short i=0;i<rows;i++) |
---|
| 198 | delete[] coefficients[i]; |
---|
| 199 | delete[] coefficients; |
---|
| 200 | |
---|
| 201 | if(_kernel_dimension>0) |
---|
| 202 | // LLL-algorithm performed |
---|
| 203 | { |
---|
| 204 | for(short i=0;i<_kernel_dimension;i++) |
---|
| 205 | delete[] H[i]; |
---|
| 206 | delete[] H; |
---|
| 207 | } |
---|
| 208 | } |
---|
| 209 | |
---|
| 210 | |
---|
| 211 | |
---|
| 212 | |
---|
| 213 | //////////////////// object properties ////////////////////////////////////// |
---|
| 214 | |
---|
| 215 | |
---|
| 216 | |
---|
| 217 | |
---|
| 218 | BOOLEAN matrix::is_nonnegative() const |
---|
| 219 | { |
---|
| 220 | for(short i=0;i<rows;i++) |
---|
| 221 | for(short j=0;j<columns;j++) |
---|
| 222 | if(coefficients[i][j]<0) |
---|
| 223 | return FALSE; |
---|
| 224 | return TRUE; |
---|
| 225 | } |
---|
| 226 | |
---|
| 227 | |
---|
| 228 | |
---|
| 229 | short matrix::error_status() const |
---|
| 230 | { |
---|
| 231 | if(columns<0) |
---|
| 232 | return columns; |
---|
| 233 | else |
---|
| 234 | return 0; |
---|
| 235 | } |
---|
| 236 | |
---|
| 237 | |
---|
| 238 | |
---|
| 239 | short matrix::row_number() const |
---|
| 240 | { |
---|
| 241 | return rows; |
---|
| 242 | } |
---|
| 243 | |
---|
| 244 | |
---|
| 245 | |
---|
| 246 | short matrix::column_number() const |
---|
| 247 | { |
---|
| 248 | return columns; |
---|
| 249 | } |
---|
| 250 | |
---|
| 251 | |
---|
| 252 | |
---|
| 253 | |
---|
| 254 | ////////// special routines for the IP-algorithms ///////////////////////// |
---|
| 255 | |
---|
| 256 | |
---|
| 257 | |
---|
| 258 | |
---|
| 259 | short matrix::LLL_kernel_basis() |
---|
| 260 | { |
---|
| 261 | |
---|
| 262 | // copy the column vectors of the actual matrix |
---|
| 263 | // (They are modified by the LLL-algorithm!) |
---|
[36f8d9] | 264 | BigInt** b=new BigIntP[columns]; |
---|
[6ba162] | 265 | for(short n=0;n<columns;n++) |
---|
| 266 | b[n]=new BigInt[rows]; |
---|
| 267 | for(short n=0;n<columns;n++) |
---|
| 268 | for(short m=0;m<rows;m++) |
---|
| 269 | b[n][m]=coefficients[m][n]; |
---|
| 270 | |
---|
| 271 | // compute a LLL-reduced basis of the relations of b[0],...,b[columns-1] |
---|
| 272 | _kernel_dimension=relations(b,columns,rows,H); |
---|
| 273 | |
---|
| 274 | // The kernel lattice basis is now stored in the member H (vectors |
---|
| 275 | // H[0],...,H[_kernel_dimension-1]). |
---|
| 276 | |
---|
| 277 | // delete auxiliary vectors |
---|
| 278 | for(short n=0;n<columns;n++) |
---|
| 279 | delete[] b[n]; |
---|
| 280 | delete[] b; |
---|
| 281 | |
---|
| 282 | return _kernel_dimension; |
---|
| 283 | } |
---|
| 284 | |
---|
| 285 | |
---|
| 286 | |
---|
| 287 | |
---|
| 288 | short matrix::compute_nonzero_kernel_vector() |
---|
| 289 | { |
---|
| 290 | |
---|
| 291 | if(_kernel_dimension==-2) |
---|
| 292 | // lattice basis not yet computed |
---|
| 293 | LLL_kernel_basis(); |
---|
| 294 | |
---|
| 295 | if(_kernel_dimension==-1) |
---|
| 296 | { |
---|
| 297 | cerr<<"\nWARNING: short matrix::compute_non_zero_kernel_vector(BigInt*&):" |
---|
| 298 | "\nerror in kernel basis, cannot compute the desired vector"<<endl; |
---|
| 299 | return 0; |
---|
| 300 | } |
---|
| 301 | |
---|
| 302 | if(_kernel_dimension==0) |
---|
| 303 | { |
---|
| 304 | cerr<<"\nWARNING: short matrix::compute_non_zero_kernel_vector(BigInt*&): " |
---|
| 305 | "\nkernel dimension is zero"<<endl; |
---|
| 306 | return 0; |
---|
| 307 | } |
---|
| 308 | |
---|
| 309 | // Now, the kernel dimension is positive. |
---|
| 310 | |
---|
| 311 | BigInt *M=new BigInt[_kernel_dimension]; |
---|
| 312 | // M stores a number by which the algorithm decides which vector to |
---|
| 313 | // take next. |
---|
| 314 | |
---|
| 315 | |
---|
| 316 | // STEP 1: Determine the vector with the least zero components (if it is not |
---|
| 317 | // unique, choose the smallest). |
---|
| 318 | |
---|
| 319 | // determine number of zero components |
---|
| 320 | for(short i=0;i<_kernel_dimension;i++) |
---|
| 321 | { |
---|
| 322 | M[i]=0; |
---|
| 323 | for(short j=0;j<columns;j++) |
---|
[4dcfb4f] | 324 | if(H[i][j]==BigInt(0)) |
---|
[6ba162] | 325 | M[i]++; |
---|
| 326 | } |
---|
| 327 | |
---|
| 328 | // determine minimal number of zero components |
---|
| 329 | BigInt min=columns; |
---|
| 330 | // columns is an upper bound (not reached because the kernel basis cannot |
---|
| 331 | // contain the zero vector) |
---|
| 332 | for(short i=0;i<_kernel_dimension;i++) |
---|
| 333 | if(M[i]<min) |
---|
| 334 | min=M[i]; |
---|
| 335 | |
---|
| 336 | // add the square of the norm to the vectors with the least zero components |
---|
| 337 | // and discard the others (the norm computation is why we have chosen the |
---|
| 338 | // M[i] to be BigInts) |
---|
| 339 | for(short i=0;i<_kernel_dimension;i++) |
---|
| 340 | if(M[i]!=min) |
---|
| 341 | M[i]=-1; |
---|
| 342 | else |
---|
| 343 | for(short j=0;j<columns;j++) |
---|
| 344 | M[i]+=H[i][j]*H[i][j]; |
---|
| 345 | // As the lattice basis does not contain the zero vector, at least one M[i] |
---|
| 346 | // is positive! |
---|
| 347 | |
---|
| 348 | // determine the start vector, i.e. the one with least zero components, but |
---|
| 349 | // smallest possible (euclidian) norm |
---|
| 350 | short min_index=-1; |
---|
| 351 | for(short i=0;i<_kernel_dimension;i++) |
---|
[4dcfb4f] | 352 | if(M[i]>BigInt(0)) |
---|
[6ba162] | 353 | if(min_index==-1) |
---|
| 354 | min_index=i; |
---|
| 355 | else |
---|
| 356 | if(M[i]<M[min_index]) |
---|
| 357 | min_index=i; |
---|
| 358 | |
---|
| 359 | // Now, H[min_index] is the vector to be transformed into a nonnegative one. |
---|
| 360 | // For a better overview, it is swapped with the first vector |
---|
| 361 | // (only pointers). |
---|
| 362 | |
---|
| 363 | if(min_index!=0) |
---|
| 364 | { |
---|
| 365 | BigInt* swap=H[min_index]; |
---|
| 366 | H[min_index]=H[0]; |
---|
| 367 | H[0]=swap; |
---|
| 368 | } |
---|
| 369 | |
---|
| 370 | |
---|
| 371 | // Now construct the desired vector. |
---|
| 372 | // This is done by adding a linear combination of |
---|
| 373 | // H[1],...,H[_kernel_dimension-1] to H[0]. It is important that the final |
---|
| 374 | // result, written as a linear combination of |
---|
| 375 | // H[0],...,H[_kernel_dimension-1], has coefficient 1 or -1 at H[0] |
---|
| 376 | // (to make sure that it is together with H[1],...,H[_kernel_dimension] |
---|
| 377 | // still a l a t t i c e basis). |
---|
| 378 | |
---|
| 379 | for(short current_position=1;current_position<columns;current_position++) |
---|
| 380 | // in fact, this loop will terminate before the condition in the |
---|
| 381 | // for-statement is satisfied... |
---|
| 382 | { |
---|
| 383 | |
---|
| 384 | |
---|
| 385 | // STEP 2: Nonnegative vector already found? |
---|
| 386 | |
---|
| 387 | BOOLEAN found=TRUE; |
---|
| 388 | for(short j=0;j<columns;j++) |
---|
[4dcfb4f] | 389 | if(H[0][j]==BigInt(0)) |
---|
[6ba162] | 390 | found=FALSE; |
---|
| 391 | |
---|
| 392 | if(found==TRUE) |
---|
| 393 | // H[0] has only positive entries, |
---|
| 394 | return 1; |
---|
| 395 | // else there are further zero components |
---|
| 396 | |
---|
| 397 | |
---|
| 398 | // STEP 3: Can a furhter zero component be "eliminated"? |
---|
| 399 | // If this is the case, find a basis vector that can do this. |
---|
| 400 | |
---|
| 401 | // determine number of components in each remaining vector that are zero |
---|
| 402 | // in the vector itself as well as in the already constructed vector |
---|
| 403 | for(short i=current_position;i<_kernel_dimension;i++) |
---|
| 404 | M[i]=0; |
---|
| 405 | |
---|
| 406 | short remaining_zero_components=0; |
---|
| 407 | |
---|
| 408 | for(short j=0;j<columns;j++) |
---|
[4dcfb4f] | 409 | if(H[0][j]==BigInt(0)) |
---|
[6ba162] | 410 | { |
---|
| 411 | remaining_zero_components++; |
---|
| 412 | for(short i=current_position;i<_kernel_dimension;i++) |
---|
[4dcfb4f] | 413 | if(H[i][j]==BigInt(0)) |
---|
[6ba162] | 414 | M[i]++; |
---|
| 415 | } |
---|
| 416 | |
---|
| 417 | // determine minimal number of such components |
---|
| 418 | min=remaining_zero_components; |
---|
| 419 | // this is the number of zero components in H[0] and an upper bound |
---|
| 420 | // for the M[i] |
---|
| 421 | for(short i=current_position;i<_kernel_dimension;i++) |
---|
| 422 | if(M[i]<min) |
---|
| 423 | min=M[i]; |
---|
| 424 | |
---|
[d96b79] | 425 | if(min==(const BigInt&)remaining_zero_components) |
---|
[6ba162] | 426 | // all zero components in H[0] are zero in each remaining vector |
---|
| 427 | // => desired vector does not exist |
---|
| 428 | return 0; |
---|
| 429 | |
---|
| 430 | // add the square of the norm to the vectors with the least common zero |
---|
| 431 | // components |
---|
| 432 | // discard the others |
---|
| 433 | for(short i=current_position;i<_kernel_dimension;i++) |
---|
| 434 | if(M[i]!=min) |
---|
| 435 | M[i]=-1; |
---|
| 436 | else |
---|
| 437 | for(short j=0;j<columns;j++) |
---|
| 438 | M[i]+=H[i][j]*H[i][j]; |
---|
| 439 | // Again, at least one M[i] is positive! |
---|
| 440 | |
---|
| 441 | // determine vector to proceed with |
---|
| 442 | // This is the vector with the least common zero components with respect |
---|
| 443 | // to H[0], but the smallest possible norm. |
---|
| 444 | short min_index=0; |
---|
| 445 | for(short i=current_position;i<_kernel_dimension;i++) |
---|
[4dcfb4f] | 446 | if(M[i]>BigInt(0)) |
---|
[6ba162] | 447 | if(min_index==0) |
---|
| 448 | min_index=i; |
---|
| 449 | else |
---|
| 450 | if(M[i]<M[min_index]) |
---|
| 451 | min_index=i; |
---|
| 452 | |
---|
| 453 | // Now, a multiple of H[min_index] will be added to the already constructed |
---|
| 454 | // vector H[0]. |
---|
| 455 | // For a better handling, it is swapped with the vector at current_position |
---|
| 456 | // (only pointers). |
---|
| 457 | |
---|
| 458 | if(min_index!=current_position) |
---|
| 459 | { |
---|
| 460 | BigInt* swap=H[min_index]; |
---|
| 461 | H[min_index]=H[current_position]; |
---|
| 462 | H[current_position]=swap; |
---|
| 463 | } |
---|
| 464 | |
---|
| 465 | |
---|
| 466 | // STEP 4: Choose a convenient multiple of H[current_position] to add to H[0]. |
---|
| 467 | // The number of factors "mult" that have to be tested is bounded by the |
---|
| 468 | // number of nonzero components in H[0] (for each such components, there is at |
---|
| 469 | // most one such factor that will eliminate it in the linear combination |
---|
| 470 | // H[0] + mult*H[current_position]. |
---|
| 471 | |
---|
| 472 | found=FALSE; |
---|
| 473 | |
---|
| 474 | for(short mult=1;found==FALSE;mult++) |
---|
| 475 | { |
---|
| 476 | found=TRUE; |
---|
| 477 | |
---|
| 478 | // check if any component !=0 of H[0] becomes zero by adding |
---|
| 479 | // mult*H[current_position] |
---|
| 480 | for(short j=0;j<columns;j++) |
---|
[4dcfb4f] | 481 | if(H[0][j]!=BigInt(0)) |
---|
[d96b79] | 482 | if(H[0][j]+(const BigInt&)mult*H[current_position][j] |
---|
[4dcfb4f] | 483 | ==BigInt(0)) |
---|
[6ba162] | 484 | found=FALSE; |
---|
| 485 | |
---|
| 486 | if(found==TRUE) |
---|
| 487 | for(short j=0;j<columns;j++) |
---|
[d96b79] | 488 | H[0][j]+=(const BigInt&)mult*H[current_position][j]; |
---|
[6ba162] | 489 | else |
---|
| 490 | // try -mult |
---|
| 491 | { |
---|
| 492 | |
---|
| 493 | found=TRUE; |
---|
| 494 | |
---|
| 495 | // check if any component !=0 of H[0] becomes zero by subtracting |
---|
| 496 | // mult*H[current_position] |
---|
| 497 | for(short j=0;j<columns;j++) |
---|
[4dcfb4f] | 498 | if(H[0][j]!=BigInt(0)) |
---|
[d96b79] | 499 | if(H[0][j]-(const BigInt&)mult*H[current_position][j] |
---|
[4dcfb4f] | 500 | ==BigInt(0)) |
---|
[6ba162] | 501 | found=FALSE; |
---|
| 502 | |
---|
| 503 | if(found==TRUE) |
---|
| 504 | for(short j=0;j<columns;j++) |
---|
[d96b79] | 505 | H[0][j]-=(const BigInt&)mult*H[current_position][j]; |
---|
[6ba162] | 506 | } |
---|
| 507 | } |
---|
| 508 | } |
---|
| 509 | |
---|
| 510 | // When reaching this line, an error must have occurred. |
---|
| 511 | cerr<<"FATAL ERROR in short matrix::compute_nonnegative_vector()"<<endl; |
---|
| 512 | abort(); |
---|
| 513 | } |
---|
| 514 | |
---|
| 515 | short matrix::compute_flip_variables(short*& F) |
---|
| 516 | { |
---|
| 517 | // first compute nonzero vector |
---|
| 518 | int okay=compute_nonzero_kernel_vector(); |
---|
| 519 | |
---|
| 520 | if(!okay) |
---|
| 521 | { |
---|
| 522 | cout<<"\nWARNING: short matrix::compute_flip_variables(short*&):\n" |
---|
| 523 | "kernel of the matrix contains no vector with nonzero components,\n" |
---|
| 524 | "no flip variables computed"<<endl; |
---|
| 525 | return -1; |
---|
| 526 | } |
---|
| 527 | |
---|
| 528 | // compute variables to flip; these might either be those corresponding |
---|
| 529 | // to the positive components of the kernel vector without zero components |
---|
| 530 | // or those corresponding to the negative ones |
---|
| 531 | |
---|
| 532 | short r=0; |
---|
| 533 | // number of flip variables |
---|
| 534 | |
---|
| 535 | for(short j=0;j<columns;j++) |
---|
[4dcfb4f] | 536 | if(H[0][j]<BigInt(0)) |
---|
[6ba162] | 537 | r++; |
---|
| 538 | // remember that all components of H[0] are !=0 |
---|
| 539 | |
---|
| 540 | if(r==0) |
---|
| 541 | // no flip variables |
---|
| 542 | return 0; |
---|
| 543 | |
---|
| 544 | if(2*r>columns) |
---|
| 545 | // more negative than positive components in H[0] |
---|
| 546 | // all variables corresponding to positive components will be flipped |
---|
| 547 | { |
---|
| 548 | r=columns-r; |
---|
| 549 | F=new short[r]; |
---|
[056f98] | 550 | memset(F,0,r*sizeof(short)); |
---|
[6ba162] | 551 | short counter=0; |
---|
| 552 | for(short j=0;j<columns;j++) |
---|
[4dcfb4f] | 553 | if(H[0][j]> BigInt(0)) |
---|
[6ba162] | 554 | { |
---|
| 555 | F[counter]=j; |
---|
| 556 | counter++; |
---|
| 557 | } |
---|
| 558 | } |
---|
| 559 | else |
---|
| 560 | // more (or as many) positive than negative components in v |
---|
| 561 | // all variables corresponding to negative components will be flipped |
---|
| 562 | { |
---|
| 563 | F=new short[r]; |
---|
[056f98] | 564 | memset(F,0,r*sizeof(short)); |
---|
[6ba162] | 565 | short counter=0; |
---|
| 566 | for(short j=0;j<columns;j++) |
---|
[4dcfb4f] | 567 | if(H[0][j]< BigInt(0)) |
---|
[6ba162] | 568 | { |
---|
| 569 | F[counter]=j; |
---|
| 570 | counter++; |
---|
| 571 | } |
---|
| 572 | } |
---|
| 573 | |
---|
| 574 | return r; |
---|
| 575 | } |
---|
| 576 | |
---|
| 577 | |
---|
| 578 | |
---|
| 579 | |
---|
| 580 | short matrix::hosten_shapiro(short*& sat_var) |
---|
| 581 | { |
---|
| 582 | |
---|
| 583 | if(_kernel_dimension==-2) |
---|
| 584 | // lattice basis not yet computed |
---|
| 585 | LLL_kernel_basis(); |
---|
| 586 | |
---|
| 587 | if(_kernel_dimension==-1) |
---|
| 588 | { |
---|
| 589 | cerr<<"\nWARNING: short matrix::hosten_shapiro(short*&):\n" |
---|
| 590 | "error in kernel basis, cannot compute the saturation variables"<<endl; |
---|
| 591 | return 0; |
---|
| 592 | } |
---|
| 593 | |
---|
| 594 | if(_kernel_dimension==0) |
---|
| 595 | // the toric ideal corresponding to the kernel lattice is the zero ideal, |
---|
| 596 | // no saturation variables necessary |
---|
| 597 | return 0; |
---|
| 598 | |
---|
| 599 | // Now, the kernel dimension is positive. |
---|
| 600 | |
---|
| 601 | if(columns==1) |
---|
| 602 | // matrix consists of one zero column, kernel is generated by the vector |
---|
| 603 | // (1) corresponding to the toric ideal <x-1> which is already staurated |
---|
| 604 | return 0; |
---|
| 605 | |
---|
| 606 | short number_of_sat_var=0; |
---|
| 607 | sat_var=new short[columns/2]; |
---|
[056f98] | 608 | memset(sat_var,0,sizeof(short)*(columns/2)); |
---|
[6ba162] | 609 | |
---|
| 610 | BOOLEAN* ideal_saturated_by_var=new BOOLEAN[columns]; |
---|
| 611 | // auxiliary array used to remember by which variables the ideal has still to |
---|
| 612 | // be saturated |
---|
| 613 | for(short j=0;j<columns;j++) |
---|
| 614 | ideal_saturated_by_var[j]=FALSE; |
---|
| 615 | |
---|
| 616 | for(short k=0;k<_kernel_dimension;k++) |
---|
| 617 | { |
---|
| 618 | // determine number of positive and negative components in H[k] |
---|
| 619 | // corresponding to variables by which the ideal has still to be saturated |
---|
| 620 | short pos_sat_var=0; |
---|
| 621 | short neg_sat_var=0; |
---|
| 622 | |
---|
| 623 | for(short j=0;j<columns;j++) |
---|
| 624 | { |
---|
| 625 | if(ideal_saturated_by_var[j]==FALSE) |
---|
| 626 | { |
---|
[4dcfb4f] | 627 | if(H[k][j]> BigInt(0)) |
---|
[6ba162] | 628 | pos_sat_var++; |
---|
| 629 | else |
---|
[4dcfb4f] | 630 | if(H[k][j]< BigInt(0)) |
---|
[6ba162] | 631 | neg_sat_var++; |
---|
| 632 | } |
---|
| 633 | } |
---|
| 634 | |
---|
| 635 | |
---|
| 636 | // now add the smaller set to the saturation variables |
---|
| 637 | if(pos_sat_var<=neg_sat_var) |
---|
| 638 | { |
---|
| 639 | for(short j=0;j<columns;j++) |
---|
| 640 | if(ideal_saturated_by_var[j]==FALSE) |
---|
[4dcfb4f] | 641 | if(H[k][j]> BigInt(0)) |
---|
[6ba162] | 642 | // ideal has to be saturated by the variables corresponding |
---|
| 643 | // to positive components |
---|
| 644 | { |
---|
| 645 | sat_var[number_of_sat_var]=j; |
---|
| 646 | ideal_saturated_by_var[j]=TRUE; |
---|
| 647 | number_of_sat_var++; |
---|
| 648 | } |
---|
| 649 | else |
---|
[4dcfb4f] | 650 | if(H[k][j]< BigInt(0)) |
---|
[6ba162] | 651 | // then the ideal is automatically saturated by the variables |
---|
| 652 | // corresponding to negative components |
---|
| 653 | ideal_saturated_by_var[j]=TRUE; |
---|
| 654 | } |
---|
| 655 | else |
---|
| 656 | { |
---|
| 657 | for(short j=0;j<columns;j++) |
---|
| 658 | if(ideal_saturated_by_var[j]==FALSE) |
---|
[4dcfb4f] | 659 | if(H[k][j]< BigInt(0)) |
---|
[6ba162] | 660 | // ideal has to be saturated by the variables corresponding |
---|
| 661 | // to negative components |
---|
| 662 | { |
---|
| 663 | sat_var[number_of_sat_var]=j; |
---|
| 664 | ideal_saturated_by_var[j]=TRUE; |
---|
| 665 | number_of_sat_var++; |
---|
| 666 | } |
---|
| 667 | else |
---|
[4dcfb4f] | 668 | if(H[k][j]> BigInt(0)) |
---|
[6ba162] | 669 | // then the ideal is automatically saturated by the variables |
---|
| 670 | // corresponding to positive components |
---|
| 671 | ideal_saturated_by_var[j]=TRUE; |
---|
| 672 | } |
---|
| 673 | } |
---|
| 674 | |
---|
| 675 | // clean up memory |
---|
| 676 | delete[] ideal_saturated_by_var; |
---|
| 677 | |
---|
| 678 | return number_of_sat_var; |
---|
| 679 | } |
---|
| 680 | |
---|
| 681 | |
---|
| 682 | |
---|
| 683 | |
---|
| 684 | //////////////////// output /////////////////////////////////////////////// |
---|
| 685 | |
---|
| 686 | |
---|
| 687 | |
---|
| 688 | |
---|
| 689 | void matrix::print() const |
---|
| 690 | { |
---|
| 691 | printf("\n%3d x %3d\n",rows,columns); |
---|
| 692 | |
---|
| 693 | for(short i=0;i<rows;i++) |
---|
| 694 | { |
---|
| 695 | for(short j=0;j<columns;j++) |
---|
| 696 | printf("%6d",coefficients[i][j]); |
---|
| 697 | printf("\n"); |
---|
| 698 | } |
---|
| 699 | } |
---|
| 700 | |
---|
| 701 | |
---|
| 702 | void matrix::print(FILE *output) const |
---|
| 703 | { |
---|
| 704 | fprintf(output,"\n%3d x %3d\n",rows,columns); |
---|
| 705 | |
---|
| 706 | for(short i=0;i<rows;i++) |
---|
| 707 | { |
---|
| 708 | for(short j=0;j<columns;j++) |
---|
| 709 | fprintf(output,"%6d",coefficients[i][j]); |
---|
| 710 | fprintf(output,"\n"); |
---|
| 711 | } |
---|
| 712 | } |
---|
| 713 | |
---|
| 714 | |
---|
| 715 | void matrix::print(ofstream& output) const |
---|
| 716 | { |
---|
| 717 | output<<endl<<setw(3)<<rows<<" x "<<setw(3)<<columns<<endl; |
---|
| 718 | |
---|
| 719 | for(short i=0;i<rows;i++) |
---|
| 720 | { |
---|
| 721 | for(short j=0;j<columns;j++) |
---|
| 722 | output<<setw(6)<<coefficients[i][j]; |
---|
| 723 | output<<endl; |
---|
| 724 | } |
---|
| 725 | } |
---|
| 726 | #endif // matrix.cc |
---|