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