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