1 | /**************************************** |
---|
2 | * Computer Algebra System SINGULAR * |
---|
3 | ****************************************/ |
---|
4 | /* |
---|
5 | * ABSTRACT: gauss implementation for F4 |
---|
6 | */ |
---|
7 | #ifdef HAVE_CONFIG_H |
---|
8 | #include "singularconfig.h" |
---|
9 | #endif /* HAVE_CONFIG_H */ |
---|
10 | #include <kernel/mod2.h> |
---|
11 | #include <misc/options.h> |
---|
12 | #include <kernel/tgbgauss.h> |
---|
13 | #include <omalloc/omalloc.h> |
---|
14 | #include <stdlib.h> |
---|
15 | #include <kernel/kutil.h> |
---|
16 | #include <kernel/febase.h> |
---|
17 | #include <kernel/polys.h> |
---|
18 | static const int bundle_size=100; |
---|
19 | |
---|
20 | mac_poly mac_p_add_ff_qq(mac_poly a, number f,mac_poly b) |
---|
21 | { |
---|
22 | mac_poly erg; |
---|
23 | mac_poly* set_this; |
---|
24 | set_this=&erg; |
---|
25 | while((a!=NULL) &&(b!=NULL)) |
---|
26 | { |
---|
27 | if (a->exp<b->exp) |
---|
28 | { |
---|
29 | (*set_this)=a; |
---|
30 | a=a->next; |
---|
31 | set_this= &((*set_this)->next); |
---|
32 | } |
---|
33 | else |
---|
34 | { |
---|
35 | if (a->exp>b->exp) |
---|
36 | { |
---|
37 | mac_poly in =new mac_poly_r(); |
---|
38 | in->exp=b->exp; |
---|
39 | in->coef=nMult(b->coef,f); |
---|
40 | (*set_this)=in; |
---|
41 | b=b->next; |
---|
42 | set_this= &((*set_this)->next); |
---|
43 | } |
---|
44 | else |
---|
45 | { |
---|
46 | //a->exp==b->ecp |
---|
47 | number n=nMult(b->coef,f); |
---|
48 | number n2=nAdd(a->coef,n); |
---|
49 | nDelete(&n); |
---|
50 | nDelete(&(a->coef)); |
---|
51 | if (nIsZero(n2)) |
---|
52 | { |
---|
53 | nDelete(&n2); |
---|
54 | mac_poly ao=a; |
---|
55 | a=a->next; |
---|
56 | delete ao; |
---|
57 | b=b->next; |
---|
58 | } |
---|
59 | else |
---|
60 | { |
---|
61 | a->coef=n2; |
---|
62 | b=b->next; |
---|
63 | (*set_this)=a; |
---|
64 | a=a->next; |
---|
65 | set_this= &((*set_this)->next); |
---|
66 | } |
---|
67 | } |
---|
68 | } |
---|
69 | } |
---|
70 | if((a==NULL)&&(b==NULL)) |
---|
71 | { |
---|
72 | (*set_this)=NULL; |
---|
73 | return erg; |
---|
74 | } |
---|
75 | if (b==NULL) |
---|
76 | { |
---|
77 | (*set_this=a); |
---|
78 | return erg; |
---|
79 | } |
---|
80 | |
---|
81 | //a==NULL |
---|
82 | while(b!=NULL) |
---|
83 | { |
---|
84 | mac_poly mp= new mac_poly_r(); |
---|
85 | mp->exp=b->exp; |
---|
86 | mp->coef=nMult(f,b->coef); |
---|
87 | (*set_this)=mp; |
---|
88 | set_this=&(mp->next); |
---|
89 | b=b->next; |
---|
90 | } |
---|
91 | (*set_this)=NULL; |
---|
92 | return erg; |
---|
93 | } |
---|
94 | |
---|
95 | void mac_mult_cons(mac_poly p,number c) |
---|
96 | { |
---|
97 | while(p) |
---|
98 | { |
---|
99 | number m=nMult(p->coef,c); |
---|
100 | nDelete(&(p->coef)); |
---|
101 | p->coef=m; |
---|
102 | p=p->next; |
---|
103 | } |
---|
104 | } |
---|
105 | |
---|
106 | int mac_length(mac_poly p) |
---|
107 | { |
---|
108 | int l=0; |
---|
109 | while(p){ |
---|
110 | l++; |
---|
111 | p=p->next; |
---|
112 | } |
---|
113 | return l; |
---|
114 | } |
---|
115 | |
---|
116 | //contrary to delete on the mac_poly_r, the coefficients are also destroyed here |
---|
117 | void mac_destroy(mac_poly p) |
---|
118 | { |
---|
119 | mac_poly iter=p; |
---|
120 | while(iter) |
---|
121 | { |
---|
122 | mac_poly next=iter->next; |
---|
123 | nDelete(&iter->coef); |
---|
124 | delete iter; |
---|
125 | iter=next; |
---|
126 | } |
---|
127 | } |
---|
128 | |
---|
129 | void simple_gauss(tgb_sparse_matrix* mat, slimgb_alg* /*c*/) |
---|
130 | { |
---|
131 | int col, row; |
---|
132 | int* row_cache=(int*) omAlloc(mat->get_rows()*sizeof(int)); |
---|
133 | col=0; |
---|
134 | row=0; |
---|
135 | int i; |
---|
136 | int pn=mat->get_rows(); |
---|
137 | int matcol=mat->get_columns(); |
---|
138 | int* area=(int*) omAlloc(sizeof(int)*((matcol-1)/bundle_size+1)); |
---|
139 | const int max_area_index=(matcol-1)/bundle_size; |
---|
140 | //rows are divided in areas |
---|
141 | //if row begins with columns col, it is located in [area[col/bundle_size],area[col/bundle_size+1]-1] |
---|
142 | assume(pn>0); |
---|
143 | //first clear zeroes |
---|
144 | for(i=0;i<pn;i++) |
---|
145 | { |
---|
146 | if(mat->zero_row(i)) |
---|
147 | { |
---|
148 | mat->perm_rows(i,pn-1); |
---|
149 | pn--; |
---|
150 | if(i!=pn){i--;} |
---|
151 | } |
---|
152 | |
---|
153 | } |
---|
154 | mat->sort_rows(); |
---|
155 | for(i=0;i<pn;i++) |
---|
156 | { |
---|
157 | row_cache[i]=mat->min_col_not_zero_in_row(i); |
---|
158 | // Print("row_cache:%d\n",row_cache[i]); |
---|
159 | } |
---|
160 | int last_area=-1; |
---|
161 | for(i=0;i<pn;i++) |
---|
162 | { |
---|
163 | int this_area=row_cache[i]/bundle_size; |
---|
164 | assume(this_area>=last_area); |
---|
165 | if(this_area>last_area) |
---|
166 | { |
---|
167 | int j; |
---|
168 | for(j=last_area+1;j<=this_area;j++) |
---|
169 | area[j]=i; |
---|
170 | last_area=this_area; |
---|
171 | } |
---|
172 | } |
---|
173 | for(i=last_area+1;i<=max_area_index;i++) |
---|
174 | { |
---|
175 | area[i]=pn; |
---|
176 | } |
---|
177 | while(row<pn-1) |
---|
178 | { |
---|
179 | //row is the row where pivot should be |
---|
180 | // row== pn-1 means we have only to act on one row so no red nec. |
---|
181 | //we assume further all rows till the pn-1 row are non-zero |
---|
182 | |
---|
183 | //select column |
---|
184 | |
---|
185 | //col=mat->min_col_not_zero_in_row(row); |
---|
186 | int max_in_area; |
---|
187 | { |
---|
188 | int tai=row_cache[row]/bundle_size; |
---|
189 | assume(tai<=max_area_index); |
---|
190 | if(tai==max_area_index) |
---|
191 | max_in_area=pn-1; |
---|
192 | else |
---|
193 | max_in_area=area[tai+1]-1; |
---|
194 | } |
---|
195 | assume(row_cache[row]==mat->min_col_not_zero_in_row(row)); |
---|
196 | col=row_cache[row]; |
---|
197 | |
---|
198 | assume(col!=matcol); |
---|
199 | int found_in_row; |
---|
200 | |
---|
201 | found_in_row=row; |
---|
202 | BOOLEAN must_reduce=FALSE; |
---|
203 | assume(pn<=mat->get_rows()); |
---|
204 | for(i=row+1;i<=max_in_area;i++) |
---|
205 | { |
---|
206 | int first;//=mat->min_col_not_zero_in_row(i); |
---|
207 | assume(row_cache[i]==mat->min_col_not_zero_in_row(i)); |
---|
208 | first=row_cache[i]; |
---|
209 | assume(first!=matcol); |
---|
210 | if(first<col) |
---|
211 | { |
---|
212 | col=first; |
---|
213 | found_in_row=i; |
---|
214 | must_reduce=FALSE; |
---|
215 | } |
---|
216 | else |
---|
217 | { |
---|
218 | if(first==col) |
---|
219 | must_reduce=TRUE; |
---|
220 | } |
---|
221 | } |
---|
222 | //select pivot |
---|
223 | int act_l=nSize(mat->get(found_in_row,col))*mat->non_zero_entries(found_in_row); |
---|
224 | if(must_reduce) |
---|
225 | { |
---|
226 | for(i=found_in_row+1;i<=max_in_area;i++) |
---|
227 | { |
---|
228 | assume(mat->min_col_not_zero_in_row(i)>=col); |
---|
229 | assume(row_cache[i]==mat->min_col_not_zero_in_row(i)); |
---|
230 | #ifndef NDEBUG |
---|
231 | int first=row_cache[i]; |
---|
232 | assume(first!=matcol); |
---|
233 | #endif |
---|
234 | // if((!(mat->is_zero_entry(i,col)))&&(mat->non_zero_entries(i)<act_l)) |
---|
235 | int nz; |
---|
236 | if((row_cache[i]==col)&&((nz=nSize(mat->get(i,col))*mat->non_zero_entries(i))<act_l)) |
---|
237 | { |
---|
238 | found_in_row=i; |
---|
239 | act_l=nz; |
---|
240 | } |
---|
241 | |
---|
242 | } |
---|
243 | } |
---|
244 | mat->perm_rows(row,found_in_row); |
---|
245 | int h=row_cache[row]; |
---|
246 | row_cache[row]=row_cache[found_in_row]; |
---|
247 | row_cache[found_in_row]=h; |
---|
248 | |
---|
249 | if(!must_reduce) |
---|
250 | { |
---|
251 | row++; |
---|
252 | continue; |
---|
253 | } |
---|
254 | //reduction |
---|
255 | //must extract content and normalize here |
---|
256 | mat->row_content(row); |
---|
257 | mat->row_normalize(row); |
---|
258 | |
---|
259 | //for(i=row+1;i<pn;i++){ |
---|
260 | for(i=max_in_area;i>row;i--) |
---|
261 | { |
---|
262 | int col_area_index=col/bundle_size; |
---|
263 | assume(col_area_index<=max_area_index); |
---|
264 | assume(mat->min_col_not_zero_in_row(i)>=col); |
---|
265 | assume(row_cache[i]==mat->min_col_not_zero_in_row(i)); |
---|
266 | #ifndef NDEBUG |
---|
267 | int first=row_cache[i]; |
---|
268 | assume(first!=matcol); |
---|
269 | #endif |
---|
270 | if(row_cache[i]==col) |
---|
271 | { |
---|
272 | |
---|
273 | number c1=mat->get(i,col); |
---|
274 | number c2=mat->get(row,col); |
---|
275 | number n1=c1; |
---|
276 | number n2=c2; |
---|
277 | |
---|
278 | ksCheckCoeff(&n1,&n2,currRing->cf); |
---|
279 | //nDelete(&c1); |
---|
280 | n1=nNeg(n1); |
---|
281 | mat->mult_row(i,n2); |
---|
282 | mat->add_lambda_times_row(i,row,n1); |
---|
283 | nDelete(&n1); |
---|
284 | nDelete(&n2); |
---|
285 | assume(mat->is_zero_entry(i,col)); |
---|
286 | row_cache[i]=mat->min_col_not_zero_in_row(i); |
---|
287 | assume(mat->min_col_not_zero_in_row(i)>col); |
---|
288 | if(row_cache[i]==matcol) |
---|
289 | { |
---|
290 | int index; |
---|
291 | index=i; |
---|
292 | int last_in_area; |
---|
293 | int this_cai=col_area_index; |
---|
294 | while(this_cai<max_area_index) |
---|
295 | { |
---|
296 | last_in_area=area[this_cai+1]-1; |
---|
297 | int h_c=row_cache[last_in_area]; |
---|
298 | row_cache[last_in_area]=row_cache[index]; |
---|
299 | row_cache[index]=h_c; |
---|
300 | mat->perm_rows(index,last_in_area); |
---|
301 | index=last_in_area; |
---|
302 | this_cai++; |
---|
303 | area[this_cai]--; |
---|
304 | } |
---|
305 | mat->perm_rows(index,pn-1); |
---|
306 | row_cache[index]=row_cache[pn-1]; |
---|
307 | row_cache[pn-1]=matcol; |
---|
308 | pn--; |
---|
309 | } |
---|
310 | else |
---|
311 | { |
---|
312 | int index; |
---|
313 | index=i; |
---|
314 | int last_in_area; |
---|
315 | int this_cai=col_area_index; |
---|
316 | int final_cai=row_cache[index]/bundle_size; |
---|
317 | assume(final_cai<=max_area_index); |
---|
318 | while(this_cai<final_cai) |
---|
319 | { |
---|
320 | last_in_area=area[this_cai+1]-1; |
---|
321 | int h_c=row_cache[last_in_area]; |
---|
322 | row_cache[last_in_area]=row_cache[index]; |
---|
323 | row_cache[index]=h_c; |
---|
324 | mat->perm_rows(index,last_in_area); |
---|
325 | index=last_in_area; |
---|
326 | this_cai++; |
---|
327 | area[this_cai]--; |
---|
328 | } |
---|
329 | } |
---|
330 | } |
---|
331 | else |
---|
332 | assume(mat->min_col_not_zero_in_row(i)>col); |
---|
333 | } |
---|
334 | // for(i=row+1;i<pn;i++) |
---|
335 | // { |
---|
336 | // assume(mat->min_col_not_zero_in_row(i)==row_cache[i]); |
---|
337 | // // if(mat->zero_row(i)) |
---|
338 | // assume(matcol==mat->get_columns()); |
---|
339 | // if(row_cache[i]==matcol) |
---|
340 | // { |
---|
341 | // assume(mat->zero_row(i)); |
---|
342 | // mat->perm_rows(i,pn-1); |
---|
343 | // row_cache[i]=row_cache[pn-1]; |
---|
344 | // row_cache[pn-1]=matcol; |
---|
345 | // pn--; |
---|
346 | // if(i!=pn){i--;} |
---|
347 | // } |
---|
348 | // } |
---|
349 | #ifdef TGB_DEBUG |
---|
350 | { |
---|
351 | int last=-1; |
---|
352 | for(i=0;i<pn;i++) |
---|
353 | { |
---|
354 | int act=mat->min_col_not_zero_in_row(i); |
---|
355 | assume(act>last); |
---|
356 | } |
---|
357 | for(i=pn;i<mat->get_rows();i++) |
---|
358 | { |
---|
359 | assume(mat->zero_row(i)); |
---|
360 | } |
---|
361 | } |
---|
362 | #endif |
---|
363 | row++; |
---|
364 | } |
---|
365 | omFree(area); |
---|
366 | omFree(row_cache); |
---|
367 | } |
---|
368 | |
---|
369 | void simple_gauss2(tgb_matrix* mat) |
---|
370 | { |
---|
371 | int col, row; |
---|
372 | col=0; |
---|
373 | row=0; |
---|
374 | int i; |
---|
375 | int pn=mat->get_rows(); |
---|
376 | assume(pn>0); |
---|
377 | //first clear zeroes |
---|
378 | // for(i=0;i<pn;i++) |
---|
379 | // { |
---|
380 | // if(mat->zero_row(i)) |
---|
381 | // { |
---|
382 | // mat->perm_rows(i,pn-1); |
---|
383 | // pn--; |
---|
384 | // if(i!=pn){i--;} |
---|
385 | // } |
---|
386 | // } |
---|
387 | while((row<pn-1)&&(col<mat->get_columns())){ |
---|
388 | //row is the row where pivot should be |
---|
389 | // row== pn-1 means we have only to act on one row so no red nec. |
---|
390 | //we assume further all rows till the pn-1 row are non-zero |
---|
391 | |
---|
392 | //select column |
---|
393 | |
---|
394 | // col=mat->min_col_not_zero_in_row(row); |
---|
395 | assume(col!=mat->get_columns()); |
---|
396 | int found_in_row=-1; |
---|
397 | |
---|
398 | // found_in_row=row; |
---|
399 | assume(pn<=mat->get_rows()); |
---|
400 | for(i=row;i<pn;i++) |
---|
401 | { |
---|
402 | // int first=mat->min_col_not_zero_in_row(i); |
---|
403 | // if(first<col) |
---|
404 | if(!(mat->is_zero_entry(i,col))) |
---|
405 | { |
---|
406 | found_in_row=i; |
---|
407 | break; |
---|
408 | } |
---|
409 | } |
---|
410 | if(found_in_row!=-1) |
---|
411 | { |
---|
412 | //select pivot |
---|
413 | int act_l=mat->non_zero_entries(found_in_row); |
---|
414 | for(i=i+1;i<pn;i++) |
---|
415 | { |
---|
416 | int vgl; |
---|
417 | assume(mat->min_col_not_zero_in_row(i)>=col); |
---|
418 | if((!(mat->is_zero_entry(i,col))) |
---|
419 | &&((vgl=mat->non_zero_entries(i))<act_l)) |
---|
420 | { |
---|
421 | found_in_row=i; |
---|
422 | act_l=vgl; |
---|
423 | } |
---|
424 | |
---|
425 | } |
---|
426 | mat->perm_rows(row,found_in_row); |
---|
427 | |
---|
428 | //reduction |
---|
429 | for(i=row+1;i<pn;i++){ |
---|
430 | assume(mat->min_col_not_zero_in_row(i)>=col); |
---|
431 | if(!(mat->is_zero_entry(i,col))) |
---|
432 | { |
---|
433 | number c1=nNeg(nCopy(mat->get(i,col))); |
---|
434 | number c2=mat->get(row,col); |
---|
435 | number n1=c1; |
---|
436 | number n2=c2; |
---|
437 | |
---|
438 | ksCheckCoeff(&n1,&n2,currRing->cf); |
---|
439 | nDelete(&c1); |
---|
440 | mat->mult_row(i,n2); |
---|
441 | mat->add_lambda_times_row(i,row,n1); |
---|
442 | assume(mat->is_zero_entry(i,col)); |
---|
443 | } |
---|
444 | assume(mat->min_col_not_zero_in_row(i)>col); |
---|
445 | } |
---|
446 | row++; |
---|
447 | } |
---|
448 | col++; |
---|
449 | // for(i=row+1;i<pn;i++) |
---|
450 | // { |
---|
451 | // if(mat->zero_row(i)) |
---|
452 | // { |
---|
453 | // mat->perm_rows(i,pn-1); |
---|
454 | // pn--; |
---|
455 | // if(i!=pn){i--;} |
---|
456 | // } |
---|
457 | // } |
---|
458 | } |
---|
459 | } |
---|
460 | |
---|
461 | |
---|
462 | tgb_matrix::tgb_matrix(int i, int j) |
---|
463 | { |
---|
464 | n=(number**) omAlloc(i*sizeof (number*));; |
---|
465 | int z; |
---|
466 | int z2; |
---|
467 | for(z=0;z<i;z++) |
---|
468 | { |
---|
469 | n[z]=(number*)omAlloc(j*sizeof(number)); |
---|
470 | for(z2=0;z2<j;z2++) |
---|
471 | { |
---|
472 | n[z][z2]=nInit(0); |
---|
473 | } |
---|
474 | } |
---|
475 | columns=j; |
---|
476 | rows=i; |
---|
477 | free_numbers=FALSE; |
---|
478 | } |
---|
479 | |
---|
480 | tgb_matrix::~tgb_matrix() |
---|
481 | { |
---|
482 | int z; |
---|
483 | for(z=0;z<rows;z++) |
---|
484 | { |
---|
485 | if(n[z]) |
---|
486 | { |
---|
487 | if(free_numbers) |
---|
488 | { |
---|
489 | int z2; |
---|
490 | for(z2=0;z2<columns;z2++) |
---|
491 | { |
---|
492 | nDelete(&(n[z][z2])); |
---|
493 | } |
---|
494 | } |
---|
495 | omFree(n[z]); |
---|
496 | } |
---|
497 | } |
---|
498 | omfree(n); |
---|
499 | } |
---|
500 | |
---|
501 | void tgb_matrix::print() |
---|
502 | { |
---|
503 | int i; |
---|
504 | int j; |
---|
505 | PrintLn(); |
---|
506 | for(i=0;i<rows;i++) |
---|
507 | { |
---|
508 | PrintS("("); |
---|
509 | for(j=0;j<columns;j++) |
---|
510 | { |
---|
511 | StringSetS(""); |
---|
512 | n_Write(n[i][j],currRing); |
---|
513 | char *s=StringEndS(); |
---|
514 | PrintS(s); |
---|
515 | omFree(s); |
---|
516 | PrintS("\t"); |
---|
517 | } |
---|
518 | PrintS(")\n"); |
---|
519 | } |
---|
520 | } |
---|
521 | |
---|
522 | //transfers ownership of n to the matrix |
---|
523 | void tgb_matrix::set(int i, int j, number nn) |
---|
524 | { |
---|
525 | assume(i<rows); |
---|
526 | assume(j<columns); |
---|
527 | n[i][j]=nn; |
---|
528 | } |
---|
529 | |
---|
530 | int tgb_matrix::get_rows() |
---|
531 | { |
---|
532 | return rows; |
---|
533 | } |
---|
534 | |
---|
535 | int tgb_matrix::get_columns() |
---|
536 | { |
---|
537 | return columns; |
---|
538 | } |
---|
539 | |
---|
540 | number tgb_matrix::get(int i, int j) |
---|
541 | { |
---|
542 | assume(i<rows); |
---|
543 | assume(j<columns); |
---|
544 | return n[i][j]; |
---|
545 | } |
---|
546 | |
---|
547 | BOOLEAN tgb_matrix::is_zero_entry(int i, int j) |
---|
548 | { |
---|
549 | return (nIsZero(n[i][j])); |
---|
550 | } |
---|
551 | |
---|
552 | void tgb_matrix::perm_rows(int i, int j) |
---|
553 | { |
---|
554 | number* h; |
---|
555 | h=n[i]; |
---|
556 | n[i]=n[j]; |
---|
557 | n[j]=h; |
---|
558 | } |
---|
559 | |
---|
560 | int tgb_matrix::min_col_not_zero_in_row(int row) |
---|
561 | { |
---|
562 | int i; |
---|
563 | for(i=0;i<columns;i++) |
---|
564 | { |
---|
565 | if(!(nIsZero(n[row][i]))) |
---|
566 | return i; |
---|
567 | } |
---|
568 | return columns;//error code |
---|
569 | } |
---|
570 | |
---|
571 | int tgb_matrix::next_col_not_zero(int row,int pre) |
---|
572 | { |
---|
573 | int i; |
---|
574 | for(i=pre+1;i<columns;i++) |
---|
575 | { |
---|
576 | if(!(nIsZero(n[row][i]))) |
---|
577 | return i; |
---|
578 | } |
---|
579 | return columns;//error code |
---|
580 | } |
---|
581 | |
---|
582 | BOOLEAN tgb_matrix::zero_row(int row) |
---|
583 | { |
---|
584 | int i; |
---|
585 | for(i=0;i<columns;i++) |
---|
586 | { |
---|
587 | if(!(nIsZero(n[row][i]))) |
---|
588 | return FALSE; |
---|
589 | } |
---|
590 | return TRUE; |
---|
591 | } |
---|
592 | |
---|
593 | int tgb_matrix::non_zero_entries(int row) |
---|
594 | { |
---|
595 | int i; |
---|
596 | int z=0; |
---|
597 | for(i=0;i<columns;i++) |
---|
598 | { |
---|
599 | if(!(nIsZero(n[row][i]))) |
---|
600 | z++; |
---|
601 | } |
---|
602 | return z; |
---|
603 | } |
---|
604 | |
---|
605 | //row add_to=row add_to +row summand*factor |
---|
606 | void tgb_matrix::add_lambda_times_row(int add_to,int summand,number factor) |
---|
607 | { |
---|
608 | int i; |
---|
609 | for(i=0;i<columns;i++) |
---|
610 | { |
---|
611 | if(!(nIsZero(n[summand][i]))) |
---|
612 | { |
---|
613 | number n1=n[add_to][i]; |
---|
614 | number n2=nMult(factor,n[summand][i]); |
---|
615 | n[add_to][i]=nAdd(n1,n2); |
---|
616 | nDelete(&n1); |
---|
617 | nDelete(&n2); |
---|
618 | } |
---|
619 | } |
---|
620 | } |
---|
621 | |
---|
622 | void tgb_matrix::mult_row(int row,number factor) |
---|
623 | { |
---|
624 | if (nIsOne(factor)) |
---|
625 | return; |
---|
626 | int i; |
---|
627 | for(i=0;i<columns;i++) |
---|
628 | { |
---|
629 | if(!(nIsZero(n[row][i]))) |
---|
630 | { |
---|
631 | number n1=n[row][i]; |
---|
632 | n[row][i]=nMult(n1,factor); |
---|
633 | nDelete(&n1); |
---|
634 | } |
---|
635 | } |
---|
636 | } |
---|
637 | |
---|
638 | void tgb_matrix::free_row(int row, BOOLEAN free_non_zeros) |
---|
639 | { |
---|
640 | int i; |
---|
641 | for(i=0;i<columns;i++) |
---|
642 | if((free_non_zeros)||(!(nIsZero(n[row][i])))) |
---|
643 | nDelete(&(n[row][i])); |
---|
644 | omFree(n[row]); |
---|
645 | n[row]=NULL; |
---|
646 | } |
---|
647 | |
---|
648 | tgb_sparse_matrix::tgb_sparse_matrix(int i, int j, ring rarg) |
---|
649 | { |
---|
650 | mp=(mac_poly*) omAlloc(i*sizeof (mac_poly));; |
---|
651 | int z; |
---|
652 | for(z=0;z<i;z++) |
---|
653 | { |
---|
654 | mp[z]=NULL; |
---|
655 | } |
---|
656 | columns=j; |
---|
657 | rows=i; |
---|
658 | free_numbers=FALSE; |
---|
659 | r=rarg; |
---|
660 | } |
---|
661 | |
---|
662 | tgb_sparse_matrix::~tgb_sparse_matrix() |
---|
663 | { |
---|
664 | int z; |
---|
665 | for(z=0;z<rows;z++) |
---|
666 | { |
---|
667 | if(mp[z]!=NULL) |
---|
668 | { |
---|
669 | if(free_numbers) |
---|
670 | { |
---|
671 | mac_destroy(mp[z]); |
---|
672 | } |
---|
673 | else { |
---|
674 | while(mp[z]!=NULL) |
---|
675 | { |
---|
676 | mac_poly next=mp[z]->next; |
---|
677 | delete mp[z]; |
---|
678 | mp[z]=next; |
---|
679 | } |
---|
680 | } |
---|
681 | } |
---|
682 | } |
---|
683 | omfree(mp); |
---|
684 | } |
---|
685 | |
---|
686 | static int row_cmp_gen(const void* a, const void* b) |
---|
687 | { |
---|
688 | const mac_poly ap= *((mac_poly*) a); |
---|
689 | const mac_poly bp=*((mac_poly*) b); |
---|
690 | if (ap==NULL) return 1; |
---|
691 | if (bp==NULL) return -1; |
---|
692 | if (ap->exp<bp->exp) return -1; |
---|
693 | return 1; |
---|
694 | } |
---|
695 | |
---|
696 | void tgb_sparse_matrix::sort_rows() |
---|
697 | { |
---|
698 | qsort(mp,rows,sizeof(mac_poly),row_cmp_gen); |
---|
699 | } |
---|
700 | |
---|
701 | void tgb_sparse_matrix::print() |
---|
702 | { |
---|
703 | int i; |
---|
704 | int j; |
---|
705 | PrintLn(); |
---|
706 | for(i=0;i<rows;i++) |
---|
707 | { |
---|
708 | PrintS("("); |
---|
709 | for(j=0;j<columns;j++) |
---|
710 | { |
---|
711 | StringSetS(""); |
---|
712 | number n=get(i,j); |
---|
713 | n_Write(n,currRing); |
---|
714 | char *s=StringEndS(); |
---|
715 | PrintS(s); |
---|
716 | omFree(s); |
---|
717 | PrintS("\t"); |
---|
718 | } |
---|
719 | PrintS(")\n"); |
---|
720 | } |
---|
721 | } |
---|
722 | |
---|
723 | //transfers ownership of n to the matrix |
---|
724 | void tgb_sparse_matrix::set(int i, int j, number n) |
---|
725 | { |
---|
726 | assume(i<rows); |
---|
727 | assume(j<columns); |
---|
728 | mac_poly* set_this=&mp[i]; |
---|
729 | // while(((*set_this)!=NULL)&&((*set_this)>exp<j)) |
---|
730 | while(((*set_this)!=NULL) && ((*set_this)->exp<j)) |
---|
731 | set_this=&((*set_this)->next); |
---|
732 | |
---|
733 | if (((*set_this)==NULL)||((*set_this)->exp>j)) |
---|
734 | { |
---|
735 | if (nIsZero(n)) return; |
---|
736 | mac_poly old=(*set_this); |
---|
737 | (*set_this)=new mac_poly_r(); |
---|
738 | (*set_this)->exp=j; |
---|
739 | (*set_this)->coef=n; |
---|
740 | (*set_this)->next=old; |
---|
741 | return; |
---|
742 | } |
---|
743 | assume((*set_this)->exp==j); |
---|
744 | if(!nIsZero(n)) |
---|
745 | { |
---|
746 | nDelete(&(*set_this)->coef); |
---|
747 | (*set_this)->coef=n; |
---|
748 | } |
---|
749 | else |
---|
750 | { |
---|
751 | nDelete(&(*set_this)->coef); |
---|
752 | mac_poly dt=(*set_this); |
---|
753 | (*set_this)=dt->next; |
---|
754 | delete dt; |
---|
755 | } |
---|
756 | return; |
---|
757 | } |
---|
758 | |
---|
759 | int tgb_sparse_matrix::get_rows() |
---|
760 | { |
---|
761 | return rows; |
---|
762 | } |
---|
763 | |
---|
764 | int tgb_sparse_matrix::get_columns() |
---|
765 | { |
---|
766 | return columns; |
---|
767 | } |
---|
768 | |
---|
769 | number tgb_sparse_matrix::get(int i, int j) |
---|
770 | { |
---|
771 | assume(i<rows); |
---|
772 | assume(j<columns); |
---|
773 | mac_poly rr=mp[i]; |
---|
774 | while((rr!=NULL)&&(rr->exp<j)) |
---|
775 | rr=rr->next; |
---|
776 | if ((rr==NULL)||(rr->exp>j)) |
---|
777 | { |
---|
778 | number n=nInit(0); |
---|
779 | return n; |
---|
780 | } |
---|
781 | assume(rr->exp==j); |
---|
782 | return rr->coef; |
---|
783 | } |
---|
784 | |
---|
785 | BOOLEAN tgb_sparse_matrix::is_zero_entry(int i, int j) |
---|
786 | { |
---|
787 | assume(i<rows); |
---|
788 | assume(j<columns); |
---|
789 | mac_poly rr=mp[i]; |
---|
790 | while((rr!=NULL)&&(rr->exp<j)) |
---|
791 | rr=rr->next; |
---|
792 | if ((rr==NULL)||(rr->exp>j)) |
---|
793 | { |
---|
794 | return TRUE; |
---|
795 | } |
---|
796 | assume(!nIsZero(rr->coef)); |
---|
797 | assume(rr->exp==j); |
---|
798 | return FALSE; |
---|
799 | } |
---|
800 | |
---|
801 | int tgb_sparse_matrix::min_col_not_zero_in_row(int row) |
---|
802 | { |
---|
803 | if(mp[row]!=NULL) |
---|
804 | { |
---|
805 | assume(!nIsZero(mp[row]->coef)); |
---|
806 | return mp[row]->exp; |
---|
807 | } |
---|
808 | else |
---|
809 | return columns;//error code |
---|
810 | } |
---|
811 | |
---|
812 | int tgb_sparse_matrix::next_col_not_zero(int row,int pre) |
---|
813 | { |
---|
814 | mac_poly rr=mp[row]; |
---|
815 | while((rr!=NULL)&&(rr->exp<=pre)) |
---|
816 | rr=rr->next; |
---|
817 | if(rr!=NULL) |
---|
818 | { |
---|
819 | assume(!nIsZero(rr->coef)); |
---|
820 | return rr->exp; |
---|
821 | } |
---|
822 | return columns;//error code |
---|
823 | } |
---|
824 | |
---|
825 | BOOLEAN tgb_sparse_matrix::zero_row(int row) |
---|
826 | { |
---|
827 | assume((mp[row]==NULL)||(!nIsZero(mp[row]->coef))); |
---|
828 | if (mp[row]==NULL) |
---|
829 | return TRUE; |
---|
830 | else |
---|
831 | return FALSE; |
---|
832 | } |
---|
833 | |
---|
834 | void tgb_sparse_matrix::row_normalize(int row) |
---|
835 | { |
---|
836 | if (!rField_has_simple_inverse(r)) /* Z/p, GF(p,n), R, long R/C */ |
---|
837 | { |
---|
838 | mac_poly m=mp[row]; |
---|
839 | while (m!=NULL) |
---|
840 | { |
---|
841 | #ifndef NDEBUG |
---|
842 | if (currRing==r) {nTest(m->coef);} |
---|
843 | #endif |
---|
844 | n_Normalize(m->coef,r); |
---|
845 | m=m->next; |
---|
846 | } |
---|
847 | } |
---|
848 | } |
---|
849 | |
---|
850 | void tgb_sparse_matrix::row_content(int row) |
---|
851 | { |
---|
852 | mac_poly ph=mp[row]; |
---|
853 | number h,d; |
---|
854 | mac_poly p; |
---|
855 | |
---|
856 | if(TEST_OPT_CONTENTSB) return; |
---|
857 | if(ph->next==NULL) |
---|
858 | { |
---|
859 | nDelete(&ph->coef); |
---|
860 | ph->coef=nInit(1); |
---|
861 | } |
---|
862 | else |
---|
863 | { |
---|
864 | nNormalize(ph->coef); |
---|
865 | if(!nGreaterZero(ph->coef)) |
---|
866 | { |
---|
867 | //ph = pNeg(ph); |
---|
868 | p=ph; |
---|
869 | while(p!=NULL) |
---|
870 | { |
---|
871 | p->coef=nNeg(p->coef); |
---|
872 | p=p->next; |
---|
873 | } |
---|
874 | } |
---|
875 | |
---|
876 | h=nCopy(ph->coef); |
---|
877 | p = ph->next; |
---|
878 | |
---|
879 | while (p!=NULL) |
---|
880 | { |
---|
881 | nNormalize(p->coef); |
---|
882 | d=nGcd(h,p->coef,currRing); |
---|
883 | nDelete(&h); |
---|
884 | h = d; |
---|
885 | if(nIsOne(h)) |
---|
886 | { |
---|
887 | break; |
---|
888 | } |
---|
889 | p=p->next; |
---|
890 | } |
---|
891 | p = ph; |
---|
892 | //number tmp; |
---|
893 | if(!nIsOne(h)) |
---|
894 | { |
---|
895 | while (p!=NULL) |
---|
896 | { |
---|
897 | d = nIntDiv(p->coef,h); |
---|
898 | nDelete(&p->coef); |
---|
899 | p->coef=d; |
---|
900 | p=p->next; |
---|
901 | } |
---|
902 | } |
---|
903 | nDelete(&h); |
---|
904 | } |
---|
905 | } |
---|
906 | int tgb_sparse_matrix::non_zero_entries(int row) |
---|
907 | { |
---|
908 | return mac_length(mp[row]); |
---|
909 | } |
---|
910 | |
---|
911 | //row add_to=row add_to +row summand*factor |
---|
912 | void tgb_sparse_matrix::add_lambda_times_row(int add_to,int summand,number factor) |
---|
913 | { |
---|
914 | mp[add_to]= mac_p_add_ff_qq(mp[add_to], factor,mp[summand]); |
---|
915 | } |
---|
916 | |
---|
917 | void tgb_sparse_matrix::mult_row(int row,number factor) |
---|
918 | { |
---|
919 | if (nIsZero(factor)) |
---|
920 | { |
---|
921 | mac_destroy(mp[row]); |
---|
922 | mp[row]=NULL; |
---|
923 | |
---|
924 | return; |
---|
925 | } |
---|
926 | if(nIsOne(factor)) |
---|
927 | return; |
---|
928 | mac_mult_cons(mp[row],factor); |
---|
929 | } |
---|
930 | |
---|
931 | void tgb_sparse_matrix::free_row(int row, BOOLEAN free_non_zeros) |
---|
932 | { |
---|
933 | if(free_non_zeros) |
---|
934 | mac_destroy(mp[row]); |
---|
935 | else |
---|
936 | { |
---|
937 | while(mp[row]!=NULL) |
---|
938 | { |
---|
939 | mac_poly next=mp[row]->next; |
---|
940 | delete mp[row]; |
---|
941 | mp[row]=next; |
---|
942 | } |
---|
943 | } |
---|
944 | mp[row]=NULL; |
---|
945 | } |
---|