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