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 ////////////////////////////////// |
---|
11 | typedef Integer* IntegerP; |
---|
12 | typedef BigInt* BigIntP; |
---|
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 | |
---|
32 | coefficients=new IntegerP[rows]; |
---|
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 | |
---|
59 | coefficients=new IntegerP[rows]; |
---|
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 | |
---|
102 | coefficients=new IntegerP[rows]; |
---|
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 | |
---|
142 | coefficients=new IntegerP[rows]; |
---|
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 | |
---|
174 | coefficients=new IntegerP[rows]; |
---|
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 | { |
---|
183 | H=new BigIntP[_kernel_dimension]; |
---|
184 | for(short k=0;k<_kernel_dimension;k++) |
---|
185 | H[k]=new BigInt[columns]; |
---|
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!) |
---|
264 | BigInt** b=new BigIntP[columns]; |
---|
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++) |
---|
324 | if(H[i][j]==(const BigInt&)0) |
---|
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++) |
---|
352 | if(M[i]>(const BigInt&)0) |
---|
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++) |
---|
389 | if(H[0][j]==(const BigInt&)0) |
---|
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++) |
---|
409 | if(H[0][j]==(const BigInt&)0) |
---|
410 | { |
---|
411 | remaining_zero_components++; |
---|
412 | for(short i=current_position;i<_kernel_dimension;i++) |
---|
413 | if(H[i][j]==(const BigInt&)0) |
---|
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 | |
---|
425 | if(min==(const BigInt&)remaining_zero_components) |
---|
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++) |
---|
446 | if(M[i]>(const BigInt&)0) |
---|
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++) |
---|
481 | if(H[0][j]!=(const BigInt&)0) |
---|
482 | if(H[0][j]+(const BigInt&)mult*H[current_position][j] |
---|
483 | ==(const BigInt&)0) |
---|
484 | found=FALSE; |
---|
485 | |
---|
486 | if(found==TRUE) |
---|
487 | for(short j=0;j<columns;j++) |
---|
488 | H[0][j]+=(const BigInt&)mult*H[current_position][j]; |
---|
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++) |
---|
498 | if(H[0][j]!=(const BigInt&)0) |
---|
499 | if(H[0][j]-(const BigInt&)mult*H[current_position][j] |
---|
500 | ==(const BigInt&)0) |
---|
501 | found=FALSE; |
---|
502 | |
---|
503 | if(found==TRUE) |
---|
504 | for(short j=0;j<columns;j++) |
---|
505 | H[0][j]-=(const BigInt&)mult*H[current_position][j]; |
---|
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++) |
---|
536 | if(H[0][j]<(const BigInt&)0) |
---|
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]; |
---|
550 | memset(F,0,r*sizeof(short)); |
---|
551 | short counter=0; |
---|
552 | for(short j=0;j<columns;j++) |
---|
553 | if(H[0][j]>(const BigInt&)0) |
---|
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]; |
---|
564 | memset(F,0,r*sizeof(short)); |
---|
565 | short counter=0; |
---|
566 | for(short j=0;j<columns;j++) |
---|
567 | if(H[0][j]<(const BigInt&)0) |
---|
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]; |
---|
608 | memset(sat_var,0,sizeof(short)*(columns/2)); |
---|
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 | { |
---|
627 | if(H[k][j]>(const BigInt&)0) |
---|
628 | pos_sat_var++; |
---|
629 | else |
---|
630 | if(H[k][j]<(const BigInt&)0) |
---|
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) |
---|
641 | if(H[k][j]>(const BigInt&)0) |
---|
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 |
---|
650 | if(H[k][j]<(const BigInt&)0) |
---|
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) |
---|
659 | if(H[k][j]<(const BigInt&)0) |
---|
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 |
---|
668 | if(H[k][j]>(const BigInt&)0) |
---|
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 |
---|