1 /****************************************
2 * Computer Algebra System SINGULAR *
3 ****************************************/
4
5 /*
6 * ABSTRACT:
7 */
8
9 #include "misc/auxiliary.h"
10
11 #include "misc/mylimits.h"
12
13 #include "misc/intvec.h"
14 #include "coeffs/numbers.h"
15
16 #include "reporter/reporter.h"
17
18
19 #include "monomials/ring.h"
20 #include "monomials/p_polys.h"
21
22 #include "simpleideals.h"
23 #include "matpol.h"
24 #include "prCopy.h"
25 #include "clapsing.h"
26
27 #include "sparsmat.h"
28
29 //omBin sip_sideal_bin = omGetSpecBin(sizeof(ip_smatrix));
30 /*0 implementation*/
31
32 static poly mp_Exdiv ( poly m, poly d, poly vars, const ring);
33 static poly mp_Select (poly fro, poly what, const ring);
34 static poly mp_SelectId (ideal I, poly what, const ring R);
35
36 /// create a r x c zero-matrix
mpNew(int r,int c)37 matrix mpNew(int r, int c)
38 {
39 int rr=r;
40 if (rr<=0) rr=1;
41 //if ( (((int)(MAX_INT_VAL/sizeof(poly))) / rr) <= c)
42 //{
43 // Werror("internal error: creating matrix[%d][%d]",r,c);
44 // return NULL;
45 //}
46 matrix rc = (matrix)omAllocBin(sip_sideal_bin);
47 rc->nrows = r;
48 rc->ncols = c;
49 rc->rank = r;
50 if ((c != 0)&&(r!=0))
51 {
52 size_t s=((size_t)r)*((size_t)c)*sizeof(poly);
53 rc->m = (poly*)omAlloc0(s);
54 //if (rc->m==NULL)
55 //{
56 // Werror("internal error: creating matrix[%d][%d]",r,c);
57 // return NULL;
58 //}
59 }
60 return rc;
61 }
62
63 /// copies matrix a (from ring r to r)
mp_Copy(matrix a,const ring r)64 matrix mp_Copy (matrix a, const ring r)
65 {
66 id_Test((ideal)a, r);
67 poly t;
68 int m=MATROWS(a), n=MATCOLS(a);
69 matrix b = mpNew(m, n);
70
71 for (long i=(long)m*(long)n-1; i>=0; i--)
72 {
73 t = a->m[i];
74 if (t!=NULL)
75 {
76 p_Normalize(t, r);
77 b->m[i] = p_Copy(t, r);
78 }
79 }
80 b->rank=a->rank;
81 return b;
82 }
83
84 /// copies matrix a from rSrc into rDst
mp_Copy(const matrix a,const ring rSrc,const ring rDst)85 matrix mp_Copy(const matrix a, const ring rSrc, const ring rDst)
86 {
87 id_Test((ideal)a, rSrc);
88
89 poly t;
90 int i, m=MATROWS(a), n=MATCOLS(a);
91
92 matrix b = mpNew(m, n);
93
94 for (i=m*n-1; i>=0; i--)
95 {
96 t = a->m[i];
97 if (t!=NULL)
98 {
99 b->m[i] = prCopyR_NoSort(t, rSrc, rDst);
100 p_Normalize(b->m[i], rDst);
101 }
102 }
103 b->rank=a->rank;
104
105 id_Test((ideal)b, rDst);
106
107 return b;
108 }
109
110
111
112 /// make it a p * unit matrix
mp_InitP(int r,int c,poly p,const ring R)113 matrix mp_InitP(int r, int c, poly p, const ring R)
114 {
115 matrix rc = mpNew(r,c);
116 int i=si_min(r,c), n = c*(i-1)+i-1, inc = c+1;
117
118 p_Normalize(p, R);
119 while (n>0)
120 {
121 rc->m[n] = p_Copy(p, R);
122 n -= inc;
123 }
124 rc->m[0]=p;
125 return rc;
126 }
127
128 /// make it a v * unit matrix
mp_InitI(int r,int c,int v,const ring R)129 matrix mp_InitI(int r, int c, int v, const ring R)
130 {
131 return mp_InitP(r, c, p_ISet(v, R), R);
132 }
133
134 /// c = f*a
mp_MultI(matrix a,int f,const ring R)135 matrix mp_MultI(matrix a, int f, const ring R)
136 {
137 int k, n = a->nrows, m = a->ncols;
138 poly p = p_ISet(f, R);
139 matrix c = mpNew(n,m);
140
141 for (k=m*n-1; k>0; k--)
142 c->m[k] = pp_Mult_qq(a->m[k], p, R);
143 c->m[0] = p_Mult_q(p_Copy(a->m[0], R), p, R);
144 return c;
145 }
146
147 /// multiply a matrix 'a' by a poly 'p', destroy the args
mp_MultP(matrix a,poly p,const ring R)148 matrix mp_MultP(matrix a, poly p, const ring R)
149 {
150 int k, n = a->nrows, m = a->ncols;
151
152 p_Normalize(p, R);
153 for (k=m*n-1; k>0; k--)
154 {
155 if (a->m[k]!=NULL)
156 a->m[k] = p_Mult_q(a->m[k], p_Copy(p, R), R);
157 }
158 a->m[0] = p_Mult_q(a->m[0], p, R);
159 return a;
160 }
161
162 /*2
163 * multiply a poly 'p' by a matrix 'a', destroy the args
164 */
pMultMp(poly p,matrix a,const ring R)165 matrix pMultMp(poly p, matrix a, const ring R)
166 {
167 int k, n = a->nrows, m = a->ncols;
168
169 p_Normalize(p, R);
170 for (k=m*n-1; k>0; k--)
171 {
172 if (a->m[k]!=NULL)
173 a->m[k] = p_Mult_q(p_Copy(p, R), a->m[k], R);
174 }
175 a->m[0] = p_Mult_q(p, a->m[0], R);
176 return a;
177 }
178
mp_Add(matrix a,matrix b,const ring R)179 matrix mp_Add(matrix a, matrix b, const ring R)
180 {
181 int k, n = a->nrows, m = a->ncols;
182 if ((n != b->nrows) || (m != b->ncols))
183 {
184 /*
185 * Werror("cannot add %dx%d matrix and %dx%d matrix",
186 * m,n,b->cols(),b->rows());
187 */
188 return NULL;
189 }
190 matrix c = mpNew(n,m);
191 for (k=m*n-1; k>=0; k--)
192 c->m[k] = p_Add_q(p_Copy(a->m[k], R), p_Copy(b->m[k], R), R);
193 return c;
194 }
195
mp_Sub(matrix a,matrix b,const ring R)196 matrix mp_Sub(matrix a, matrix b, const ring R)
197 {
198 int k, n = a->nrows, m = a->ncols;
199 if ((n != b->nrows) || (m != b->ncols))
200 {
201 /*
202 * Werror("cannot sub %dx%d matrix and %dx%d matrix",
203 * m,n,b->cols(),b->rows());
204 */
205 return NULL;
206 }
207 matrix c = mpNew(n,m);
208 for (k=m*n-1; k>=0; k--)
209 c->m[k] = p_Sub(p_Copy(a->m[k], R), p_Copy(b->m[k], R), R);
210 return c;
211 }
212
mp_Mult(matrix a,matrix b,const ring R)213 matrix mp_Mult(matrix a, matrix b, const ring R)
214 {
215 int i, j, k;
216 int m = MATROWS(a);
217 int p = MATCOLS(a);
218 int q = MATCOLS(b);
219
220 if (p!=MATROWS(b))
221 {
222 /*
223 * Werror("cannot multiply %dx%d matrix and %dx%d matrix",
224 * m,p,b->rows(),q);
225 */
226 return NULL;
227 }
228 matrix c = mpNew(m,q);
229
230 for (i=0; i<m; i++)
231 {
232 for (k=0; k<p; k++)
233 {
234 poly aik;
235 if ((aik=MATELEM0(a,i,k))!=NULL)
236 {
237 for (j=0; j<q; j++)
238 {
239 poly bkj;
240 if ((bkj=MATELEM0(b,k,j))!=NULL)
241 {
242 poly *cij=&(MATELEM0(c,i,j));
243 poly s = pp_Mult_qq(aik /*MATELEM0(a,i,k)*/, bkj/*MATELEM0(b,k,j)*/, R);
244 (*cij)/*MATELEM0(c,i,j)*/ = p_Add_q((*cij) /*MATELEM0(c,i,j)*/ ,s, R);
245 }
246 }
247 }
248 }
249 }
250 for(i=m*q-1;i>=0;i--) p_Normalize(c->m[i], R);
251 return c;
252 }
253
mp_Transp(matrix a,const ring R)254 matrix mp_Transp(matrix a, const ring R)
255 {
256 int i, j, r = MATROWS(a), c = MATCOLS(a);
257 poly *p;
258 matrix b = mpNew(c,r);
259
260 p = b->m;
261 for (i=0; i<c; i++)
262 {
263 for (j=0; j<r; j++)
264 {
265 if (a->m[j*c+i]!=NULL) *p = p_Copy(a->m[j*c+i], R);
266 p++;
267 }
268 }
269 return b;
270 }
271
272 /*2
273 *returns the trace of matrix a
274 */
mp_Trace(matrix a,const ring R)275 poly mp_Trace ( matrix a, const ring R)
276 {
277 int i;
278 int n = (MATCOLS(a)<MATROWS(a)) ? MATCOLS(a) : MATROWS(a);
279 poly t = NULL;
280
281 for (i=1; i<=n; i++)
282 t = p_Add_q(t, p_Copy(MATELEM(a,i,i), R), R);
283 return t;
284 }
285
286 /*2
287 *returns the trace of the product of a and b
288 */
TraceOfProd(matrix a,matrix b,int n,const ring R)289 poly TraceOfProd ( matrix a, matrix b, int n, const ring R)
290 {
291 int i, j;
292 poly p, t = NULL;
293
294 for (i=1; i<=n; i++)
295 {
296 for (j=1; j<=n; j++)
297 {
298 p = pp_Mult_qq(MATELEM(a,i,j), MATELEM(b,j,i), R);
299 t = p_Add_q(t, p, R);
300 }
301 }
302 return t;
303 }
304
305 // #ifndef SIZE_OF_SYSTEM_PAGE
306 // #define SIZE_OF_SYSTEM_PAGE 4096
307 // #endif
308
309 /*2
310 * corresponds to Maple's coeffs:
311 * var has to be the number of a variable
312 */
mp_Coeffs(ideal I,int var,const ring R)313 matrix mp_Coeffs (ideal I, int var, const ring R)
314 {
315 poly h,f;
316 int l, i, c, m=0;
317 /* look for maximal power m of x_var in I */
318 for (i=IDELEMS(I)-1; i>=0; i--)
319 {
320 f=I->m[i];
321 while (f!=NULL)
322 {
323 l=p_GetExp(f,var, R);
324 if (l>m) m=l;
325 pIter(f);
326 }
327 }
328 matrix co=mpNew((m+1)*I->rank,IDELEMS(I));
329 /* divide each monomial by a power of x_var,
330 * remember the power in l and the component in c*/
331 for (i=IDELEMS(I)-1; i>=0; i--)
332 {
333 f=I->m[i];
334 I->m[i]=NULL;
335 while (f!=NULL)
336 {
337 l=p_GetExp(f,var, R);
338 p_SetExp(f,var,0, R);
339 c=si_max((int)p_GetComp(f, R),1);
340 p_SetComp(f,0, R);
341 p_Setm(f, R);
342 /* now add the resulting monomial to co*/
343 h=pNext(f);
344 pNext(f)=NULL;
345 //MATELEM(co,c*(m+1)-l,i+1)
346 // =p_Add_q(MATELEM(co,c*(m+1)-l,i+1),f, R);
347 MATELEM(co,(c-1)*(m+1)+l+1,i+1)
348 =p_Add_q(MATELEM(co,(c-1)*(m+1)+l+1,i+1),f, R);
349 /* iterate f*/
350 f=h;
351 }
352 }
353 id_Delete(&I, R);
354 return co;
355 }
356
357 /*2
358 * given the result c of mpCoeffs(ideal/module i, var)
359 * i of rank r
360 * build the matrix of the corresponding monomials in m
361 */
mp_Monomials(matrix c,int r,int var,matrix m,const ring R)362 void mp_Monomials(matrix c, int r, int var, matrix m, const ring R)
363 {
364 /* clear contents of m*/
365 int k,l;
366 for (k=MATROWS(m);k>0;k--)
367 {
368 for(l=MATCOLS(m);l>0;l--)
369 {
370 p_Delete(&MATELEM(m,k,l), R);
371 }
372 }
373 omfreeSize((ADDRESS)m->m,MATROWS(m)*MATCOLS(m)*sizeof(poly));
374 /* allocate monoms in the right size r x MATROWS(c)*/
375 m->m=(poly*)omAlloc0(r*MATROWS(c)*sizeof(poly));
376 MATROWS(m)=r;
377 MATCOLS(m)=MATROWS(c);
378 m->rank=r;
379 /* the maximal power p of x_var: MATCOLS(m)=r*(p+1) */
380 int p=MATCOLS(m)/r-1;
381 /* fill in the powers of x_var=h*/
382 poly h=p_One(R);
383 for(k=r;k>0; k--)
384 {
385 MATELEM(m,k,k*(p+1))=p_One(R);
386 }
387 for(l=p;l>=0; l--)
388 {
389 p_SetExp(h,var,p-l, R);
390 p_Setm(h, R);
391 for(k=r;k>0; k--)
392 {
393 MATELEM(m,k,k*(p+1)-l)=p_Copy(h, R);
394 }
395 }
396 p_Delete(&h, R);
397 }
398
mp_CoeffProc(poly f,poly vars,const ring R)399 matrix mp_CoeffProc (poly f, poly vars, const ring R)
400 {
401 assume(vars!=NULL);
402 poly sel, h;
403 int l, i;
404 int pos_of_1 = -1;
405 matrix co;
406
407 if (f==NULL)
408 {
409 co = mpNew(2, 1);
410 MATELEM(co,1,1) = p_One(R);
411 //MATELEM(co,2,1) = NULL;
412 return co;
413 }
414 sel = mp_Select(f, vars, R);
415 l = pLength(sel);
416 co = mpNew(2, l);
417
418 if (rHasLocalOrMixedOrdering(R))
419 {
420 for (i=l; i>=1; i--)
421 {
422 h = sel;
423 pIter(sel);
424 pNext(h)=NULL;
425 MATELEM(co,1,i) = h;
426 //MATELEM(co,2,i) = NULL;
427 if (p_IsConstant(h, R)) pos_of_1 = i;
428 }
429 }
430 else
431 {
432 for (i=1; i<=l; i++)
433 {
434 h = sel;
435 pIter(sel);
436 pNext(h)=NULL;
437 MATELEM(co,1,i) = h;
438 //MATELEM(co,2,i) = NULL;
439 if (p_IsConstant(h, R)) pos_of_1 = i;
440 }
441 }
442 while (f!=NULL)
443 {
444 i = 1;
445 loop
446 {
447 if (i!=pos_of_1)
448 {
449 h = mp_Exdiv(f, MATELEM(co,1,i),vars, R);
450 if (h!=NULL)
451 {
452 MATELEM(co,2,i) = p_Add_q(MATELEM(co,2,i), h, R);
453 break;
454 }
455 }
456 if (i == l)
457 {
458 // check monom 1 last:
459 if (pos_of_1 != -1)
460 {
461 h = mp_Exdiv(f, MATELEM(co,1,pos_of_1),vars, R);
462 if (h!=NULL)
463 {
464 MATELEM(co,2,pos_of_1) = p_Add_q(MATELEM(co,2,pos_of_1), h, R);
465 }
466 }
467 break;
468 }
469 i ++;
470 }
471 pIter(f);
472 }
473 return co;
474 }
475
mp_CoeffProcId(ideal I,poly vars,const ring R)476 matrix mp_CoeffProcId (ideal I, poly vars, const ring R)
477 {
478 assume(vars!=NULL);
479 poly sel, h;
480 int l, i;
481 int pos_of_1 = -1;
482 matrix co;
483
484 if (idIs0(I))
485 {
486 co = mpNew(IDELEMS(I)+1,1);
487 MATELEM(co,1,1) = p_One(R);
488 return co;
489 }
490 sel = mp_SelectId(I, vars, R);
491 l = pLength(sel);
492 co = mpNew(IDELEMS(I)+1, l);
493
494 if (rHasLocalOrMixedOrdering(R))
495 {
496 for (i=l; i>=1; i--)
497 {
498 h = sel;
499 pIter(sel);
500 pNext(h)=NULL;
501 MATELEM(co,1,i) = h;
502 //MATELEM(co,2,i) = NULL;
503 if (p_IsConstant(h, R)) pos_of_1 = i;
504 }
505 }
506 else
507 {
508 for (i=1; i<=l; i++)
509 {
510 h = sel;
511 pIter(sel);
512 pNext(h)=NULL;
513 MATELEM(co,1,i) = h;
514 //MATELEM(co,2,i) = NULL;
515 if (p_IsConstant(h, R)) pos_of_1 = i;
516 }
517 }
518 for(int j=0;j<IDELEMS(I);j++)
519 {
520 poly f=I->m[j];
521 while (f!=NULL)
522 {
523 i = 1;
524 loop
525 {
526 if (i!=pos_of_1)
527 {
528 h = mp_Exdiv(f, MATELEM(co,1,i),vars, R);
529 if (h!=NULL)
530 {
531 MATELEM(co,j+2,i) = p_Add_q(MATELEM(co,j+2,i), h, R);
532 break;
533 }
534 }
535 if (i == l)
536 {
537 // check monom 1 last:
538 if (pos_of_1 != -1)
539 {
540 h = mp_Exdiv(f, MATELEM(co,1,pos_of_1),vars, R);
541 if (h!=NULL)
542 {
543 MATELEM(co,j+2,pos_of_1) = p_Add_q(MATELEM(co,j+2,pos_of_1), h, R);
544 }
545 }
546 break;
547 }
548 i ++;
549 }
550 pIter(f);
551 }
552 }
553 return co;
554 }
555
556 /*2
557 *exact divisor: let d == x^i*y^j, m is thought to have only one term;
558 * return m/d iff d divides m, and no x^k*y^l (k>i or l>j) divides m
559 * consider all variables in vars
560 */
mp_Exdiv(poly m,poly d,poly vars,const ring R)561 static poly mp_Exdiv ( poly m, poly d, poly vars, const ring R)
562 {
563 int i;
564 poly h = p_Head(m, R);
565 for (i=1; i<=rVar(R); i++)
566 {
567 if (p_GetExp(vars,i, R) > 0)
568 {
569 if (p_GetExp(d,i, R) != p_GetExp(h,i, R))
570 {
571 p_Delete(&h, R);
572 return NULL;
573 }
574 p_SetExp(h,i,0, R);
575 }
576 }
577 p_Setm(h, R);
578 return h;
579 }
580
mp_Coef2(poly v,poly mon,matrix * c,matrix * m,const ring R)581 void mp_Coef2(poly v, poly mon, matrix *c, matrix *m, const ring R)
582 {
583 poly* s;
584 poly p;
585 int sl,i,j;
586 int l=0;
587 poly sel=mp_Select(v,mon, R);
588
589 p_Vec2Polys(sel,&s,&sl, R);
590 for (i=0; i<sl; i++)
591 l=si_max(l,pLength(s[i]));
592 *c=mpNew(sl,l);
593 *m=mpNew(sl,l);
594 poly h;
595 int isConst;
596 for (j=1; j<=sl;j++)
597 {
598 p=s[j-1];
599 if (p_IsConstant(p, R)) /*p != NULL */
600 {
601 isConst=-1;
602 i=l;
603 }
604 else
605 {
606 isConst=1;
607 i=1;
608 }
609 while(p!=NULL)
610 {
611 h = p_Head(p, R);
612 MATELEM(*m,j,i) = h;
613 i+=isConst;
614 p = p->next;
615 }
616 }
617 while (v!=NULL)
618 {
619 i = 1;
620 j = __p_GetComp(v, R);
621 loop
622 {
623 poly mp=MATELEM(*m,j,i);
624 if (mp!=NULL)
625 {
626 h = mp_Exdiv(v, mp /*MATELEM(*m,j,i)*/, mp, R);
627 if (h!=NULL)
628 {
629 p_SetComp(h,0, R);
630 MATELEM(*c,j,i) = p_Add_q(MATELEM(*c,j,i), h, R);
631 break;
632 }
633 }
634 if (i < l)
635 i++;
636 else
637 break;
638 }
639 v = v->next;
640 }
641 }
642
mp_Compare(matrix a,matrix b,const ring R)643 int mp_Compare(matrix a, matrix b, const ring R)
644 {
645 if (MATCOLS(a)<MATCOLS(b)) return -1;
646 else if (MATCOLS(a)>MATCOLS(b)) return 1;
647 if (MATROWS(a)<MATROWS(b)) return -1;
648 else if (MATROWS(a)<MATROWS(b)) return 1;
649
650 unsigned ii=MATCOLS(a)*MATROWS(a)-1;
651 unsigned j=0;
652 int r=0;
653 while (j<=ii)
654 {
655 r=p_Compare(a->m[j],b->m[j],R);
656 if (r!=0) return r;
657 j++;
658 }
659 return r;
660 }
661
mp_Equal(matrix a,matrix b,const ring R)662 BOOLEAN mp_Equal(matrix a, matrix b, const ring R)
663 {
664 if ((MATCOLS(a)!=MATCOLS(b)) || (MATROWS(a)!=MATROWS(b)))
665 return FALSE;
666 int i=MATCOLS(a)*MATROWS(a)-1;
667 while (i>=0)
668 {
669 if (a->m[i]==NULL)
670 {
671 if (b->m[i]!=NULL) return FALSE;
672 }
673 else if (b->m[i]==NULL) return FALSE;
674 else if (p_Cmp(a->m[i],b->m[i], R)!=0) return FALSE;
675 i--;
676 }
677 i=MATCOLS(a)*MATROWS(a)-1;
678 while (i>=0)
679 {
680 if(!p_EqualPolys(a->m[i],b->m[i], R)) return FALSE;
681 i--;
682 }
683 return TRUE;
684 }
685
686 /*2
687 * insert a monomial into a list, avoid duplicates
688 * arguments are destroyed
689 */
p_Insert(poly p1,poly p2,const ring R)690 static poly p_Insert(poly p1, poly p2, const ring R)
691 {
692 poly a1, p, a2, a;
693 int c;
694
695 if (p1==NULL) return p2;
696 if (p2==NULL) return p1;
697 a1 = p1;
698 a2 = p2;
699 a = p = p_One(R);
700 loop
701 {
702 c = p_Cmp(a1, a2, R);
703 if (c == 1)
704 {
705 a = pNext(a) = a1;
706 pIter(a1);
707 if (a1==NULL)
708 {
709 pNext(a) = a2;
710 break;
711 }
712 }
713 else if (c == -1)
714 {
715 a = pNext(a) = a2;
716 pIter(a2);
717 if (a2==NULL)
718 {
719 pNext(a) = a1;
720 break;
721 }
722 }
723 else
724 {
725 p_LmDelete(&a2, R);
726 a = pNext(a) = a1;
727 pIter(a1);
728 if (a1==NULL)
729 {
730 pNext(a) = a2;
731 break;
732 }
733 else if (a2==NULL)
734 {
735 pNext(a) = a1;
736 break;
737 }
738 }
739 }
740 p_LmDelete(&p, R);
741 return p;
742 }
743
744 /*2
745 *if what == xy the result is the list of all different power products
746 * x^i*y^j (i, j >= 0) that appear in fro
747 */
mp_Select(poly fro,poly what,const ring R)748 static poly mp_Select (poly fro, poly what, const ring R)
749 {
750 int i;
751 poly h, res;
752 res = NULL;
753 while (fro!=NULL)
754 {
755 h = p_One(R);
756 for (i=1; i<=rVar(R); i++)
757 p_SetExp(h,i, p_GetExp(fro,i, R) * p_GetExp(what, i, R), R);
758 p_SetComp(h, p_GetComp(fro, R), R);
759 p_Setm(h, R);
760 res = p_Insert(h, res, R);
761 fro = fro->next;
762 }
763 return res;
764 }
765
mp_SelectId(ideal I,poly what,const ring R)766 static poly mp_SelectId (ideal I, poly what, const ring R)
767 {
768 int i;
769 poly h, res;
770 res = NULL;
771 for(int j=0;j<IDELEMS(I);j++)
772 {
773 poly fro=I->m[j];
774 while (fro!=NULL)
775 {
776 h = p_One(R);
777 for (i=1; i<=rVar(R); i++)
778 p_SetExp(h,i, p_GetExp(fro,i, R) * p_GetExp(what, i, R), R);
779 p_SetComp(h, p_GetComp(fro, R), R);
780 p_Setm(h, R);
781 res = p_Insert(h, res, R);
782 fro = fro->next;
783 }
784 }
785 return res;
786 }
787
788 /*
789 *static void ppp(matrix a)
790 *{
791 * int j,i,r=a->nrows,c=a->ncols;
792 * for(j=1;j<=r;j++)
793 * {
794 * for(i=1;i<=c;i++)
795 * {
796 * if(MATELEM(a,j,i)!=NULL) PrintS("X");
797 * else PrintS("0");
798 * }
799 * PrintLn();
800 * }
801 *}
802 */
803
mp_PartClean(matrix a,int lr,int lc,const ring R)804 static void mp_PartClean(matrix a, int lr, int lc, const ring R)
805 {
806 poly *q1;
807 int i,j;
808
809 for (i=lr-1;i>=0;i--)
810 {
811 q1 = &(a->m)[i*a->ncols];
812 for (j=lc-1;j>=0;j--) if(q1[j]) p_Delete(&q1[j], R);
813 }
814 }
815
mp_IsDiagUnit(matrix U,const ring R)816 BOOLEAN mp_IsDiagUnit(matrix U, const ring R)
817 {
818 if(MATROWS(U)!=MATCOLS(U))
819 return FALSE;
820 for(int i=MATCOLS(U);i>=1;i--)
821 {
822 for(int j=MATCOLS(U); j>=1; j--)
823 {
824 if (i==j)
825 {
826 if (!p_IsUnit(MATELEM(U,i,i), R)) return FALSE;
827 }
828 else if (MATELEM(U,i,j)!=NULL) return FALSE;
829 }
830 }
831 return TRUE;
832 }
833
iiWriteMatrix(matrix im,const char * n,int dim,const ring r,int spaces)834 void iiWriteMatrix(matrix im, const char *n, int dim, const ring r, int spaces)
835 {
836 int i,ii = MATROWS(im)-1;
837 int j,jj = MATCOLS(im)-1;
838 poly *pp = im->m;
839
840 for (i=0; i<=ii; i++)
841 {
842 for (j=0; j<=jj; j++)
843 {
844 if (spaces>0)
845 Print("%-*.*s",spaces,spaces," ");
846 if (dim == 2) Print("%s[%u,%u]=",n,i+1,j+1);
847 else if (dim == 1) Print("%s[%u]=",n,j+1);
848 else if (dim == 0) Print("%s=",n);
849 if ((i<ii)||(j<jj)) p_Write(*pp++, r);
850 else p_Write0(*pp, r);
851 }
852 }
853 }
854
iiStringMatrix(matrix im,int dim,const ring r,char ch)855 char * iiStringMatrix(matrix im, int dim, const ring r, char ch)
856 {
857 int i,ii = MATROWS(im);
858 int j,jj = MATCOLS(im);
859 poly *pp = im->m;
860 char ch_s[2];
861 ch_s[0]=ch;
862 ch_s[1]='\0';
863
864 StringSetS("");
865
866 for (i=0; i<ii; i++)
867 {
868 for (j=0; j<jj; j++)
869 {
870 p_String0(*pp++, r);
871 StringAppendS(ch_s);
872 if (dim > 1) StringAppendS("\n");
873 }
874 }
875 char *s=StringEndS();
876 s[strlen(s)- (dim > 1 ? 2 : 1)]='\0';
877 return s;
878 }
879
mp_Delete(matrix * a,const ring r)880 void mp_Delete(matrix* a, const ring r)
881 {
882 id_Delete((ideal *) a, r);
883 }
884
885 /*
886 * C++ classes for Bareiss algorithm
887 */
888 class row_col_weight
889 {
890 private:
891 int ym, yn;
892 public:
893 float *wrow, *wcol;
row_col_weight()894 row_col_weight() : ym(0) {}
895 row_col_weight(int, int);
896 ~row_col_weight();
897 };
898
row_col_weight(int i,int j)899 row_col_weight::row_col_weight(int i, int j)
900 {
901 ym = i;
902 yn = j;
903 wrow = (float *)omAlloc(i*sizeof(float));
904 wcol = (float *)omAlloc(j*sizeof(float));
905 }
~row_col_weight()906 row_col_weight::~row_col_weight()
907 {
908 if (ym!=0)
909 {
910 omFreeSize((ADDRESS)wcol, yn*sizeof(float));
911 omFreeSize((ADDRESS)wrow, ym*sizeof(float));
912 }
913 }
914
915 /*2
916 * a submatrix M of a matrix X[m,n]:
917 * 0 <= i < s_m <= a_m
918 * 0 <= j < s_n <= a_n
919 * M = ( Xarray[qrow[i],qcol[j]] )
920 * if a_m = a_n and s_m = s_n
921 * det(X) = sign*div^(s_m-1)*det(M)
922 * resticted pivot for elimination
923 * 0 <= j < piv_s
924 */
925 class mp_permmatrix
926 {
927 private:
928 int a_m, a_n, s_m, s_n, sign, piv_s;
929 int *qrow, *qcol;
930 poly *Xarray;
931 ring _R;
932 void mpInitMat();
mpRowAdr(int r)933 poly * mpRowAdr(int r)
934 { return &(Xarray[a_n*qrow[r]]); }
mpColAdr(int c)935 poly * mpColAdr(int c)
936 { return &(Xarray[qcol[c]]); }
937 void mpRowWeight(float *);
938 void mpColWeight(float *);
939 void mpRowSwap(int, int);
940 void mpColSwap(int, int);
941 public:
mp_permmatrix()942 mp_permmatrix() : a_m(0) {}
943 mp_permmatrix(matrix, ring);
944 mp_permmatrix(mp_permmatrix *);
945 ~mp_permmatrix();
946 int mpGetRow();
947 int mpGetCol();
mpGetRdim()948 int mpGetRdim() { return s_m; }
mpGetCdim()949 int mpGetCdim() { return s_n; }
mpGetSign()950 int mpGetSign() { return sign; }
951 void mpSetSearch(int s);
mpSaveArray()952 void mpSaveArray() { Xarray = NULL; }
953 poly mpGetElem(int, int);
954 void mpSetElem(poly, int, int);
955 void mpDelElem(int, int);
956 void mpElimBareiss(poly);
957 int mpPivotBareiss(row_col_weight *);
958 int mpPivotRow(row_col_weight *, int);
959 void mpToIntvec(intvec *);
960 void mpRowReorder();
961 void mpColReorder();
962 };
mp_permmatrix(matrix A,ring R)963 mp_permmatrix::mp_permmatrix(matrix A, ring R) : sign(1)
964 {
965 a_m = A->nrows;
966 a_n = A->ncols;
967 this->mpInitMat();
968 Xarray = A->m;
969 _R=R;
970 }
971
mp_permmatrix(mp_permmatrix * M)972 mp_permmatrix::mp_permmatrix(mp_permmatrix *M)
973 {
974 poly p, *athis, *aM;
975 int i, j;
976
977 _R=M->_R;
978 a_m = M->s_m;
979 a_n = M->s_n;
980 sign = M->sign;
981 this->mpInitMat();
982 Xarray = (poly *)omAlloc0(a_m*a_n*sizeof(poly));
983 for (i=a_m-1; i>=0; i--)
984 {
985 athis = this->mpRowAdr(i);
986 aM = M->mpRowAdr(i);
987 for (j=a_n-1; j>=0; j--)
988 {
989 p = aM[M->qcol[j]];
990 if (p)
991 {
992 athis[j] = p_Copy(p,_R);
993 }
994 }
995 }
996 }
997
~mp_permmatrix()998 mp_permmatrix::~mp_permmatrix()
999 {
1000 int k;
1001
1002 if (a_m != 0)
1003 {
1004 omFreeSize((ADDRESS)qrow,a_m*sizeof(int));
1005 omFreeSize((ADDRESS)qcol,a_n*sizeof(int));
1006 if (Xarray != NULL)
1007 {
1008 for (k=a_m*a_n-1; k>=0; k--)
1009 p_Delete(&Xarray[k],_R);
1010 omFreeSize((ADDRESS)Xarray,a_m*a_n*sizeof(poly));
1011 }
1012 }
1013 }
1014
1015
1016 static float mp_PolyWeight(poly p, const ring r);
mpColWeight(float * wcol)1017 void mp_permmatrix::mpColWeight(float *wcol)
1018 {
1019 poly p, *a;
1020 int i, j;
1021 float count;
1022
1023 for (j=s_n; j>=0; j--)
1024 {
1025 a = this->mpColAdr(j);
1026 count = 0.0;
1027 for(i=s_m; i>=0; i--)
1028 {
1029 p = a[a_n*qrow[i]];
1030 if (p)
1031 count += mp_PolyWeight(p,_R);
1032 }
1033 wcol[j] = count;
1034 }
1035 }
mpRowWeight(float * wrow)1036 void mp_permmatrix::mpRowWeight(float *wrow)
1037 {
1038 poly p, *a;
1039 int i, j;
1040 float count;
1041
1042 for (i=s_m; i>=0; i--)
1043 {
1044 a = this->mpRowAdr(i);
1045 count = 0.0;
1046 for(j=s_n; j>=0; j--)
1047 {
1048 p = a[qcol[j]];
1049 if (p)
1050 count += mp_PolyWeight(p,_R);
1051 }
1052 wrow[i] = count;
1053 }
1054 }
1055
mpRowSwap(int i1,int i2)1056 void mp_permmatrix::mpRowSwap(int i1, int i2)
1057 {
1058 poly p, *a1, *a2;
1059 int j;
1060
1061 a1 = &(Xarray[a_n*i1]);
1062 a2 = &(Xarray[a_n*i2]);
1063 for (j=a_n-1; j>= 0; j--)
1064 {
1065 p = a1[j];
1066 a1[j] = a2[j];
1067 a2[j] = p;
1068 }
1069 }
1070
mpColSwap(int j1,int j2)1071 void mp_permmatrix::mpColSwap(int j1, int j2)
1072 {
1073 poly p, *a1, *a2;
1074 int i, k = a_n*a_m;
1075
1076 a1 = &(Xarray[j1]);
1077 a2 = &(Xarray[j2]);
1078 for (i=0; i< k; i+=a_n)
1079 {
1080 p = a1[i];
1081 a1[i] = a2[i];
1082 a2[i] = p;
1083 }
1084 }
mpInitMat()1085 void mp_permmatrix::mpInitMat()
1086 {
1087 int k;
1088
1089 s_m = a_m;
1090 s_n = a_n;
1091 piv_s = 0;
1092 qrow = (int *)omAlloc(a_m*sizeof(int));
1093 qcol = (int *)omAlloc(a_n*sizeof(int));
1094 for (k=a_m-1; k>=0; k--) qrow[k] = k;
1095 for (k=a_n-1; k>=0; k--) qcol[k] = k;
1096 }
1097
mpColReorder()1098 void mp_permmatrix::mpColReorder()
1099 {
1100 int k, j, j1, j2;
1101
1102 if (a_n > a_m)
1103 k = a_n - a_m;
1104 else
1105 k = 0;
1106 for (j=a_n-1; j>=k; j--)
1107 {
1108 j1 = qcol[j];
1109 if (j1 != j)
1110 {
1111 this->mpColSwap(j1, j);
1112 j2 = 0;
1113 while (qcol[j2] != j) j2++;
1114 qcol[j2] = j1;
1115 }
1116 }
1117 }
1118
mpRowReorder()1119 void mp_permmatrix::mpRowReorder()
1120 {
1121 int k, i, i1, i2;
1122
1123 if (a_m > a_n)
1124 k = a_m - a_n;
1125 else
1126 k = 0;
1127 for (i=a_m-1; i>=k; i--)
1128 {
1129 i1 = qrow[i];
1130 if (i1 != i)
1131 {
1132 this->mpRowSwap(i1, i);
1133 i2 = 0;
1134 while (qrow[i2] != i) i2++;
1135 qrow[i2] = i1;
1136 }
1137 }
1138 }
1139
1140 /*
1141 * perform replacement for pivot strategy in Bareiss algorithm
1142 * change sign of determinant
1143 */
mpReplace(int j,int n,int & sign,int * perm)1144 static void mpReplace(int j, int n, int &sign, int *perm)
1145 {
1146 int k;
1147
1148 if (j != n)
1149 {
1150 k = perm[n];
1151 perm[n] = perm[j];
1152 perm[j] = k;
1153 sign = -sign;
1154 }
1155 }
1156 /*2
1157 * pivot strategy for Bareiss algorithm
1158 */
mpPivotBareiss(row_col_weight * C)1159 int mp_permmatrix::mpPivotBareiss(row_col_weight *C)
1160 {
1161 poly p, *a;
1162 int i, j, iopt, jopt;
1163 float sum, f1, f2, fo, r, ro, lp;
1164 float *dr = C->wrow, *dc = C->wcol;
1165
1166 fo = 1.0e20;
1167 ro = 0.0;
1168 iopt = jopt = -1;
1169
1170 s_n--;
1171 s_m--;
1172 if (s_m == 0)
1173 return 0;
1174 if (s_n == 0)
1175 {
1176 for(i=s_m; i>=0; i--)
1177 {
1178 p = this->mpRowAdr(i)[qcol[0]];
1179 if (p)
1180 {
1181 f1 = mp_PolyWeight(p,_R);
1182 if (f1 < fo)
1183 {
1184 fo = f1;
1185 if (iopt >= 0)
1186 p_Delete(&(this->mpRowAdr(iopt)[qcol[0]]),_R);
1187 iopt = i;
1188 }
1189 else
1190 p_Delete(&(this->mpRowAdr(i)[qcol[0]]),_R);
1191 }
1192 }
1193 if (iopt >= 0)
1194 mpReplace(iopt, s_m, sign, qrow);
1195 return 0;
1196 }
1197 this->mpRowWeight(dr);
1198 this->mpColWeight(dc);
1199 sum = 0.0;
1200 for(i=s_m; i>=0; i--)
1201 sum += dr[i];
1202 for(i=s_m; i>=0; i--)
1203 {
1204 r = dr[i];
1205 a = this->mpRowAdr(i);
1206 for(j=s_n; j>=0; j--)
1207 {
1208 p = a[qcol[j]];
1209 if (p)
1210 {
1211 lp = mp_PolyWeight(p,_R);
1212 ro = r - lp;
1213 f1 = ro * (dc[j]-lp);
1214 if (f1 != 0.0)
1215 {
1216 f2 = lp * (sum - ro - dc[j]);
1217 f2 += f1;
1218 }
1219 else
1220 f2 = lp-r-dc[j];
1221 if (f2 < fo)
1222 {
1223 fo = f2;
1224 iopt = i;
1225 jopt = j;
1226 }
1227 }
1228 }
1229 }
1230 if (iopt < 0)
1231 return 0;
1232 mpReplace(iopt, s_m, sign, qrow);
1233 mpReplace(jopt, s_n, sign, qcol);
1234 return 1;
1235 }
mpGetElem(int r,int c)1236 poly mp_permmatrix::mpGetElem(int r, int c)
1237 {
1238 return Xarray[a_n*qrow[r]+qcol[c]];
1239 }
1240
1241 /*
1242 * the Bareiss-type elimination with division by div (div != NULL)
1243 */
mpElimBareiss(poly div)1244 void mp_permmatrix::mpElimBareiss(poly div)
1245 {
1246 poly piv, elim, q1, q2, *ap, *a;
1247 int i, j, jj;
1248
1249 ap = this->mpRowAdr(s_m);
1250 piv = ap[qcol[s_n]];
1251 for(i=s_m-1; i>=0; i--)
1252 {
1253 a = this->mpRowAdr(i);
1254 elim = a[qcol[s_n]];
1255 if (elim != NULL)
1256 {
1257 elim = p_Neg(elim,_R);
1258 for (j=s_n-1; j>=0; j--)
1259 {
1260 q2 = NULL;
1261 jj = qcol[j];
1262 if (ap[jj] != NULL)
1263 {
1264 q2 = SM_MULT(ap[jj], elim, div,_R);
1265 if (a[jj] != NULL)
1266 {
1267 q1 = SM_MULT(a[jj], piv, div,_R);
1268 p_Delete(&a[jj],_R);
1269 q2 = p_Add_q(q2, q1, _R);
1270 }
1271 }
1272 else if (a[jj] != NULL)
1273 {
1274 q2 = SM_MULT(a[jj], piv, div, _R);
1275 }
1276 if ((q2!=NULL) && div)
1277 SM_DIV(q2, div, _R);
1278 a[jj] = q2;
1279 }
1280 p_Delete(&a[qcol[s_n]], _R);
1281 }
1282 else
1283 {
1284 for (j=s_n-1; j>=0; j--)
1285 {
1286 jj = qcol[j];
1287 if (a[jj] != NULL)
1288 {
1289 q2 = SM_MULT(a[jj], piv, div, _R);
1290 p_Delete(&a[jj], _R);
1291 if (div)
1292 SM_DIV(q2, div, _R);
1293 a[jj] = q2;
1294 }
1295 }
1296 }
1297 }
1298 }
1299 /*
1300 * weigth of a polynomial, for pivot strategy
1301 */
mp_PolyWeight(poly p,const ring r)1302 static float mp_PolyWeight(poly p, const ring r)
1303 {
1304 int i;
1305 float res;
1306
1307 if (pNext(p) == NULL)
1308 {
1309 res = (float)n_Size(pGetCoeff(p),r->cf);
1310 for (i=r->N;i>0;i--)
1311 {
1312 if(p_GetExp(p,i,r)!=0)
1313 {
1314 res += 2.0;
1315 break;
1316 }
1317 }
1318 }
1319 else
1320 {
1321 res = 0.0;
1322 do
1323 {
1324 res += (float)n_Size(pGetCoeff(p),r->cf)+2.0;
1325 pIter(p);
1326 }
1327 while (p);
1328 }
1329 return res;
1330 }
1331 /*
1332 * find best row
1333 */
mp_PivBar(matrix a,int lr,int lc,const ring r)1334 static int mp_PivBar(matrix a, int lr, int lc, const ring r)
1335 {
1336 float f1, f2;
1337 poly *q1;
1338 int i,j,io;
1339
1340 io = -1;
1341 f1 = 1.0e30;
1342 for (i=lr-1;i>=0;i--)
1343 {
1344 q1 = &(a->m)[i*a->ncols];
1345 f2 = 0.0;
1346 for (j=lc-1;j>=0;j--)
1347 {
1348 if (q1[j]!=NULL)
1349 f2 += mp_PolyWeight(q1[j],r);
1350 }
1351 if ((f2!=0.0) && (f2<f1))
1352 {
1353 f1 = f2;
1354 io = i;
1355 }
1356 }
1357 if (io<0) return 0;
1358 else return io+1;
1359 }
1360
mpSwapRow(matrix a,int pos,int lr,int lc)1361 static void mpSwapRow(matrix a, int pos, int lr, int lc)
1362 {
1363 poly sw;
1364 int j;
1365 poly* a2 = a->m;
1366 poly* a1 = &a2[a->ncols*(pos-1)];
1367
1368 a2 = &a2[a->ncols*(lr-1)];
1369 for (j=lc-1; j>=0; j--)
1370 {
1371 sw = a1[j];
1372 a1[j] = a2[j];
1373 a2[j] = sw;
1374 }
1375 }
1376
1377 /*2
1378 * prepare one step of 'Bareiss' algorithm
1379 * for application in minor
1380 */
mp_PrepareRow(matrix a,int lr,int lc,const ring R)1381 static int mp_PrepareRow (matrix a, int lr, int lc, const ring R)
1382 {
1383 int r;
1384
1385 r = mp_PivBar(a,lr,lc,R);
1386 if(r==0) return 0;
1387 if(r<lr) mpSwapRow(a, r, lr, lc);
1388 return 1;
1389 }
1390
1391 /*
1392 * find pivot in the last row
1393 */
mp_PivRow(matrix a,int lr,int lc,const ring r)1394 static int mp_PivRow(matrix a, int lr, int lc, const ring r)
1395 {
1396 float f1, f2;
1397 poly *q1;
1398 int j,jo;
1399
1400 jo = -1;
1401 f1 = 1.0e30;
1402 q1 = &(a->m)[(lr-1)*a->ncols];
1403 for (j=lc-1;j>=0;j--)
1404 {
1405 if (q1[j]!=NULL)
1406 {
1407 f2 = mp_PolyWeight(q1[j],r);
1408 if (f2<f1)
1409 {
1410 f1 = f2;
1411 jo = j;
1412 }
1413 }
1414 }
1415 if (jo<0) return 0;
1416 else return jo+1;
1417 }
1418
mpSwapCol(matrix a,int pos,int lr,int lc)1419 static void mpSwapCol(matrix a, int pos, int lr, int lc)
1420 {
1421 poly sw;
1422 int j;
1423 poly* a2 = a->m;
1424 poly* a1 = &a2[pos-1];
1425
1426 a2 = &a2[lc-1];
1427 for (j=a->ncols*(lr-1); j>=0; j-=a->ncols)
1428 {
1429 sw = a1[j];
1430 a1[j] = a2[j];
1431 a2[j] = sw;
1432 }
1433 }
1434
1435 /*2
1436 * prepare one step of 'Bareiss' algorithm
1437 * for application in minor
1438 */
mp_PreparePiv(matrix a,int lr,int lc,const ring r)1439 static int mp_PreparePiv (matrix a, int lr, int lc,const ring r)
1440 {
1441 int c;
1442
1443 c = mp_PivRow(a, lr, lc,r);
1444 if(c==0) return 0;
1445 if(c<lc) mpSwapCol(a, c, lr, lc);
1446 return 1;
1447 }
1448
mp_ElimBar(matrix a0,matrix re,poly div,int lr,int lc,const ring R)1449 static void mp_ElimBar(matrix a0, matrix re, poly div, int lr, int lc, const ring R)
1450 {
1451 int r=lr-1, c=lc-1;
1452 poly *b = a0->m, *x = re->m;
1453 poly piv, elim, q1, *ap, *a, *q;
1454 int i, j;
1455
1456 ap = &b[r*a0->ncols];
1457 piv = ap[c];
1458 for(j=c-1; j>=0; j--)
1459 if (ap[j] != NULL) ap[j] = p_Neg(ap[j],R);
1460 for(i=r-1; i>=0; i--)
1461 {
1462 a = &b[i*a0->ncols];
1463 q = &x[i*re->ncols];
1464 if (a[c] != NULL)
1465 {
1466 elim = a[c];
1467 for (j=c-1; j>=0; j--)
1468 {
1469 q1 = NULL;
1470 if (a[j] != NULL)
1471 {
1472 q1 = sm_MultDiv(a[j], piv, div,R);
1473 if (ap[j] != NULL)
1474 {
1475 poly q2 = sm_MultDiv(ap[j], elim, div, R);
1476 q1 = p_Add_q(q1,q2,R);
1477 }
1478 }
1479 else if (ap[j] != NULL)
1480 q1 = sm_MultDiv(ap[j], elim, div, R);
1481 if (q1 != NULL)
1482 {
1483 if (div)
1484 sm_SpecialPolyDiv(q1, div,R);
1485 q[j] = q1;
1486 }
1487 }
1488 }
1489 else
1490 {
1491 for (j=c-1; j>=0; j--)
1492 {
1493 if (a[j] != NULL)
1494 {
1495 q1 = sm_MultDiv(a[j], piv, div, R);
1496 if (div)
1497 sm_SpecialPolyDiv(q1, div, R);
1498 q[j] = q1;
1499 }
1500 }
1501 }
1502 }
1503 }
1504
1505 /*2*/
1506 /// entries of a are minors and go to result (only if not in R)
mp_MinorToResult(ideal result,int & elems,matrix a,int r,int c,ideal R,const ring)1507 void mp_MinorToResult(ideal result, int &elems, matrix a, int r, int c,
1508 ideal R, const ring)
1509 {
1510 poly *q1;
1511 int e=IDELEMS(result);
1512 int i,j;
1513
1514 if (R != NULL)
1515 {
1516 for (i=r-1;i>=0;i--)
1517 {
1518 q1 = &(a->m)[i*a->ncols];
1519 //for (j=c-1;j>=0;j--)
1520 //{
1521 // if (q1[j]!=NULL) q1[j] = kNF(R,currRing->qideal,q1[j]);
1522 //}
1523 }
1524 }
1525 for (i=r-1;i>=0;i--)
1526 {
1527 q1 = &(a->m)[i*a->ncols];
1528 for (j=c-1;j>=0;j--)
1529 {
1530 if (q1[j]!=NULL)
1531 {
1532 if (elems>=e)
1533 {
1534 pEnlargeSet(&(result->m),e,e);
1535 e += e;
1536 IDELEMS(result) =e;
1537 }
1538 result->m[elems] = q1[j];
1539 q1[j] = NULL;
1540 elems++;
1541 }
1542 }
1543 }
1544 }
1545 /*
1546 // from linalg_from_matpol.cc: TODO: compare with above & remove...
1547 void mp_MinorToResult(ideal result, int &elems, matrix a, int r, int c,
1548 ideal R, const ring R)
1549 {
1550 poly *q1;
1551 int e=IDELEMS(result);
1552 int i,j;
1553
1554 if (R != NULL)
1555 {
1556 for (i=r-1;i>=0;i--)
1557 {
1558 q1 = &(a->m)[i*a->ncols];
1559 for (j=c-1;j>=0;j--)
1560 {
1561 if (q1[j]!=NULL) q1[j] = kNF(R,currRing->qideal,q1[j]);
1562 }
1563 }
1564 }
1565 for (i=r-1;i>=0;i--)
1566 {
1567 q1 = &(a->m)[i*a->ncols];
1568 for (j=c-1;j>=0;j--)
1569 {
1570 if (q1[j]!=NULL)
1571 {
1572 if (elems>=e)
1573 {
1574 if(e<SIZE_OF_SYSTEM_PAGE)
1575 {
1576 pEnlargeSet(&(result->m),e,e);
1577 e += e;
1578 }
1579 else
1580 {
1581 pEnlargeSet(&(result->m),e,SIZE_OF_SYSTEM_PAGE);
1582 e += SIZE_OF_SYSTEM_PAGE;
1583 }
1584 IDELEMS(result) =e;
1585 }
1586 result->m[elems] = q1[j];
1587 q1[j] = NULL;
1588 elems++;
1589 }
1590 }
1591 }
1592 }
1593 */
1594
mpFinalClean(matrix a)1595 static void mpFinalClean(matrix a)
1596 {
1597 omFreeSize((ADDRESS)a->m,a->nrows*a->ncols*sizeof(poly));
1598 omFreeBin((ADDRESS)a, sip_sideal_bin);
1599 }
1600
1601 /*2*/
1602 /// produces recursively the ideal of all arxar-minors of a
mp_RecMin(int ar,ideal result,int & elems,matrix a,int lr,int lc,poly barDiv,ideal R,const ring r)1603 void mp_RecMin(int ar,ideal result,int &elems,matrix a,int lr,int lc,
1604 poly barDiv, ideal R, const ring r)
1605 {
1606 int k;
1607 int kr=lr-1,kc=lc-1;
1608 matrix nextLevel=mpNew(kr,kc);
1609
1610 loop
1611 {
1612 /*--- look for an optimal row and bring it to last position ------------*/
1613 if(mp_PrepareRow(a,lr,lc,r)==0) break;
1614 /*--- now take all pivots from the last row ------------*/
1615 k = lc;
1616 loop
1617 {
1618 if(mp_PreparePiv(a,lr,k,r)==0) break;
1619 mp_ElimBar(a,nextLevel,barDiv,lr,k,r);
1620 k--;
1621 if (ar>1)
1622 {
1623 mp_RecMin(ar-1,result,elems,nextLevel,kr,k,a->m[kr*a->ncols+k],R,r);
1624 mp_PartClean(nextLevel,kr,k, r);
1625 }
1626 else mp_MinorToResult(result,elems,nextLevel,kr,k,R,r);
1627 if (ar>k-1) break;
1628 }
1629 if (ar>=kr) break;
1630 /*--- now we have to take out the last row...------------*/
1631 lr = kr;
1632 kr--;
1633 }
1634 mpFinalClean(nextLevel);
1635 }
1636 /*
1637 // from linalg_from_matpol.cc: TODO: compare with above & remove...
1638 void mp_RecMin(int ar,ideal result,int &elems,matrix a,int lr,int lc,
1639 poly barDiv, ideal R, const ring R)
1640 {
1641 int k;
1642 int kr=lr-1,kc=lc-1;
1643 matrix nextLevel=mpNew(kr,kc);
1644
1645 loop
1646 {
1647 // --- look for an optimal row and bring it to last position ------------
1648 if(mpPrepareRow(a,lr,lc)==0) break;
1649 // --- now take all pivots from the last row ------------
1650 k = lc;
1651 loop
1652 {
1653 if(mpPreparePiv(a,lr,k)==0) break;
1654 mpElimBar(a,nextLevel,barDiv,lr,k);
1655 k--;
1656 if (ar>1)
1657 {
1658 mpRecMin(ar-1,result,elems,nextLevel,kr,k,a->m[kr*a->ncols+k],R);
1659 mpPartClean(nextLevel,kr,k);
1660 }
1661 else mpMinorToResult(result,elems,nextLevel,kr,k,R);
1662 if (ar>k-1) break;
1663 }
1664 if (ar>=kr) break;
1665 // --- now we have to take out the last row...------------
1666 lr = kr;
1667 kr--;
1668 }
1669 mpFinalClean(nextLevel);
1670 }
1671 */
1672
1673 /*2*/
1674 /// returns the determinant of the matrix m;
1675 /// uses Bareiss algorithm
mp_DetBareiss(matrix a,const ring r)1676 poly mp_DetBareiss (matrix a, const ring r)
1677 {
1678 int s;
1679 poly div, res;
1680 if (MATROWS(a) != MATCOLS(a))
1681 {
1682 Werror("det of %d x %d matrix",MATROWS(a),MATCOLS(a));
1683 return NULL;
1684 }
1685 matrix c = mp_Copy(a,r);
1686 mp_permmatrix *Bareiss = new mp_permmatrix(c,r);
1687 row_col_weight w(Bareiss->mpGetRdim(), Bareiss->mpGetCdim());
1688
1689 /* Bareiss */
1690 div = NULL;
1691 while(Bareiss->mpPivotBareiss(&w))
1692 {
1693 Bareiss->mpElimBareiss(div);
1694 div = Bareiss->mpGetElem(Bareiss->mpGetRdim(), Bareiss->mpGetCdim());
1695 }
1696 Bareiss->mpRowReorder();
1697 Bareiss->mpColReorder();
1698 Bareiss->mpSaveArray();
1699 s = Bareiss->mpGetSign();
1700 delete Bareiss;
1701
1702 /* result */
1703 res = MATELEM(c,1,1);
1704 MATELEM(c,1,1) = NULL;
1705 id_Delete((ideal *)&c,r);
1706 if (s < 0)
1707 res = p_Neg(res,r);
1708 return res;
1709 }
1710 /*
1711 // from linalg_from_matpol.cc: TODO: compare with above & remove...
1712 poly mp_DetBareiss (matrix a, const ring R)
1713 {
1714 int s;
1715 poly div, res;
1716 if (MATROWS(a) != MATCOLS(a))
1717 {
1718 Werror("det of %d x %d matrix",MATROWS(a),MATCOLS(a));
1719 return NULL;
1720 }
1721 matrix c = mp_Copy(a, R);
1722 mp_permmatrix *Bareiss = new mp_permmatrix(c, R);
1723 row_col_weight w(Bareiss->mpGetRdim(), Bareiss->mpGetCdim());
1724
1725 // Bareiss
1726 div = NULL;
1727 while(Bareiss->mpPivotBareiss(&w))
1728 {
1729 Bareiss->mpElimBareiss(div);
1730 div = Bareiss->mpGetElem(Bareiss->mpGetRdim(), Bareiss->mpGetCdim());
1731 }
1732 Bareiss->mpRowReorder();
1733 Bareiss->mpColReorder();
1734 Bareiss->mpSaveArray();
1735 s = Bareiss->mpGetSign();
1736 delete Bareiss;
1737
1738 // result
1739 res = MATELEM(c,1,1);
1740 MATELEM(c,1,1) = NULL;
1741 id_Delete((ideal *)&c, R);
1742 if (s < 0)
1743 res = p_Neg(res, R);
1744 return res;
1745 }
1746 */
1747
1748 /*2
1749 * compute all ar-minors of the matrix a
1750 */
mp_Wedge(matrix a,int ar,const ring R)1751 matrix mp_Wedge(matrix a, int ar, const ring R)
1752 {
1753 int i,j,k,l;
1754 int *rowchoise,*colchoise;
1755 BOOLEAN rowch,colch;
1756 matrix result;
1757 matrix tmp;
1758 poly p;
1759
1760 i = binom(a->nrows,ar);
1761 j = binom(a->ncols,ar);
1762
1763 rowchoise=(int *)omAlloc(ar*sizeof(int));
1764 colchoise=(int *)omAlloc(ar*sizeof(int));
1765 result = mpNew(i,j);
1766 tmp = mpNew(ar,ar);
1767 l = 1; /* k,l:the index in result*/
1768 idInitChoise(ar,1,a->nrows,&rowch,rowchoise);
1769 while (!rowch)
1770 {
1771 k=1;
1772 idInitChoise(ar,1,a->ncols,&colch,colchoise);
1773 while (!colch)
1774 {
1775 for (i=1; i<=ar; i++)
1776 {
1777 for (j=1; j<=ar; j++)
1778 {
1779 MATELEM(tmp,i,j) = MATELEM(a,rowchoise[i-1],colchoise[j-1]);
1780 }
1781 }
1782 p = mp_DetBareiss(tmp, R);
1783 if ((k+l) & 1) p=p_Neg(p, R);
1784 MATELEM(result,l,k) = p;
1785 k++;
1786 idGetNextChoise(ar,a->ncols,&colch,colchoise);
1787 }
1788 idGetNextChoise(ar,a->nrows,&rowch,rowchoise);
1789 l++;
1790 }
1791
1792 /*delete the matrix tmp*/
1793 for (i=1; i<=ar; i++)
1794 {
1795 for (j=1; j<=ar; j++) MATELEM(tmp,i,j) = NULL;
1796 }
1797 id_Delete((ideal *) &tmp, R);
1798 return (result);
1799 }
1800
1801 // helper for sm_Tensor
1802 // destroyes f, keeps B
sm_MultAndShift(poly f,ideal B,int s,const ring r)1803 static ideal sm_MultAndShift(poly f, ideal B, int s, const ring r)
1804 {
1805 assume(f!=NULL);
1806 ideal res=idInit(IDELEMS(B),B->rank+s);
1807 int q=IDELEMS(B); // p x q
1808 for(int j=0;j<q;j++)
1809 {
1810 poly h=pp_Mult_qq(f,B->m[j],r);
1811 if (h!=NULL)
1812 {
1813 if (s>0) p_Shift(&h,s,r);
1814 res->m[j]=h;
1815 }
1816 }
1817 p_Delete(&f,r);
1818 return res;
1819 }
1820 // helper for sm_Tensor
1821 // updates res, destroyes contents of sm
sm_AddSubMat(ideal res,ideal sm,int col,const ring r)1822 static void sm_AddSubMat(ideal res, ideal sm, int col, const ring r)
1823 {
1824 for(int i=0;i<IDELEMS(sm);i++)
1825 {
1826 res->m[col+i]=p_Add_q(res->m[col+i],sm->m[i],r);
1827 sm->m[i]=NULL;
1828 }
1829 }
1830
sm_Tensor(ideal A,ideal B,const ring r)1831 ideal sm_Tensor(ideal A, ideal B, const ring r)
1832 {
1833 // size of the result m*p x n*q
1834 int n=IDELEMS(A); // m x n
1835 int m=A->rank;
1836 int q=IDELEMS(B); // p x q
1837 int p=B->rank;
1838 ideal res=idInit(n*q,m*p);
1839 poly *a=(poly*)omAlloc(m*sizeof(poly));
1840 for(int i=0; i<n; i++)
1841 {
1842 memset(a,0,m*sizeof(poly));
1843 p_Vec2Array(A->m[i],a,m,r);
1844 for(int j=0;j<m;j++)
1845 {
1846 if (a[j]!=NULL)
1847 {
1848 ideal sm=sm_MultAndShift(a[j], // A_i_j
1849 B,
1850 j*p, // shift j*p down
1851 r);
1852 sm_AddSubMat(res,sm,i*q,r); // add this columns to col i*q ff
1853 id_Delete(&sm,r); // delete the now empty ideal
1854 }
1855 }
1856 }
1857 omFreeSize(a,m*sizeof(poly));
1858 return res;
1859 }
1860 // --------------------------------------------------------------------------
1861 /****************************************
1862 * Computer Algebra System SINGULAR *
1863 ****************************************/
1864
1865 /*
1866 * ABSTRACT: basic operation for sparse matrices:
1867 * type: ideal (of column vectors)
1868 * nrows: I->rank, ncols: IDELEMS(I)
1869 */
1870
sm_Add(ideal a,ideal b,const ring R)1871 ideal sm_Add(ideal a, ideal b, const ring R)
1872 {
1873 assume(IDELEMS(a)==IDELEMS(b));
1874 assume(a->rank==b->rank);
1875 ideal c=idInit(IDELEMS(a),a->rank);
1876 for (int k=IDELEMS(a)-1; k>=0; k--)
1877 c->m[k] = p_Add_q(p_Copy(a->m[k], R), p_Copy(b->m[k], R), R);
1878 return c;
1879 }
1880
sm_Sub(ideal a,ideal b,const ring R)1881 ideal sm_Sub(ideal a, ideal b, const ring R)
1882 {
1883 assume(IDELEMS(a)==IDELEMS(b));
1884 assume(a->rank==b->rank);
1885 ideal c=idInit(IDELEMS(a),a->rank);
1886 for (int k=IDELEMS(a)-1; k>=0; k--)
1887 c->m[k] = p_Sub(p_Copy(a->m[k], R), p_Copy(b->m[k], R), R);
1888 return c;
1889 }
1890
sm_Mult(ideal a,ideal b,const ring R)1891 ideal sm_Mult(ideal a, ideal b, const ring R)
1892 {
1893 int i, j, k;
1894 int m = a->rank;
1895 int p = IDELEMS(a);
1896 int q = IDELEMS(b);
1897
1898 assume (IDELEMS(a)==b->rank);
1899 ideal c = idInit(q,m);
1900
1901 for (i=0; i<m; i++)
1902 {
1903 for (k=0; k<p; k++)
1904 {
1905 poly aik;
1906 if ((aik=SMATELEM(a,i,k,R))!=NULL)
1907 {
1908 for (j=0; j<q; j++)
1909 {
1910 poly bkj=SMATELEM(b,k,j,R);
1911 if (bkj!=NULL)
1912 {
1913 poly s = p_Mult_q(p_Copy(aik,R) /*SMATELEM(a,i,k)*/, bkj/*SMATELEM(b,k,j)*/, R);
1914 if (s!=NULL) p_SetComp(s,i+1,R);
1915 c->m[j]=p_Add_q(c->m[j],s, R);
1916 }
1917 }
1918 p_Delete(&aik,R);
1919 }
1920 }
1921 }
1922 for(i=q-1;i>=0;i--) p_Normalize(c->m[i], R);
1923 return c;
1924 }
1925
sm_Flatten(ideal a,const ring R)1926 ideal sm_Flatten(ideal a, const ring R)
1927 {
1928 if (IDELEMS(a)==0) return id_Copy(a,R);
1929 ideal res=idInit(1,IDELEMS(a)*a->rank);
1930 for(int i=0;i<IDELEMS(a);i++)
1931 {
1932 if(a->m[i]!=NULL)
1933 {
1934 poly p=p_Copy(a->m[i],R);
1935 if (i==0) res->m[0]=p;
1936 else
1937 {
1938 p_Shift(&p,i*a->rank,R);
1939 res->m[0]=p_Add_q(res->m[0],p,R);
1940 }
1941 }
1942 }
1943 return res;
1944 }
1945
sm_UnFlatten(ideal a,int col,const ring R)1946 ideal sm_UnFlatten(ideal a, int col, const ring R)
1947 {
1948 if ((IDELEMS(a)!=1)
1949 ||((a->rank % col)!=0))
1950 {
1951 Werror("wrong format: %d x %d for unflatten",(int)a->rank,IDELEMS(a));
1952 return NULL;
1953 }
1954 int row=a->rank/col;
1955 ideal res=idInit(col,row);
1956 poly p=a->m[0];
1957 while(p!=NULL)
1958 {
1959 poly h=p_Head(p,R);
1960 int comp=p_GetComp(h,R);
1961 int c=(comp-1)/row;
1962 int r=comp%row; if (r==0) r=row;
1963 p_SetComp(h,r,R); p_SetmComp(h,R);
1964 res->m[c]=p_Add_q(res->m[c],h,R);
1965 pIter(p);
1966 }
1967 return res;
1968 }
1969
1970 /*2
1971 *returns the trace of matrix a
1972 */
sm_Trace(ideal a,const ring R)1973 poly sm_Trace ( ideal a, const ring R)
1974 {
1975 int i;
1976 int n = (IDELEMS(a)<a->rank) ? IDELEMS(a) : a->rank;
1977 poly t = NULL;
1978
1979 for (i=0; i<=n; i++)
1980 t = p_Add_q(t, p_Copy(SMATELEM(a,i,i,R), R), R);
1981 return t;
1982 }
1983
sm_Compare(ideal a,ideal b,const ring R)1984 int sm_Compare(ideal a, ideal b, const ring R)
1985 {
1986 if (IDELEMS(a)<IDELEMS(b)) return -1;
1987 else if (IDELEMS(a)>IDELEMS(b)) return 1;
1988 if ((a->rank)<(b->rank)) return -1;
1989 else if ((a->rank)<(b->rank)) return 1;
1990
1991 unsigned ii=IDELEMS(a)-1;
1992 unsigned j=0;
1993 int r=0;
1994 while (j<=ii)
1995 {
1996 r=p_Compare(a->m[j],b->m[j],R);
1997 if (r!=0) return r;
1998 j++;
1999 }
2000 return r;
2001 }
2002
sm_Equal(ideal a,ideal b,const ring R)2003 BOOLEAN sm_Equal(ideal a, ideal b, const ring R)
2004 {
2005 if ((a->rank!=b->rank) || (IDELEMS(a)!=IDELEMS(b)))
2006 return FALSE;
2007 int i=IDELEMS(a)-1;
2008 while (i>=0)
2009 {
2010 if (a->m[i]==NULL)
2011 {
2012 if (b->m[i]!=NULL) return FALSE;
2013 }
2014 else if (b->m[i]==NULL) return FALSE;
2015 else if (p_Cmp(a->m[i],b->m[i], R)!=0) return FALSE;
2016 i--;
2017 }
2018 i=IDELEMS(a)-1;
2019 while (i>=0)
2020 {
2021 if(!p_EqualPolys(a->m[i],b->m[i], R)) return FALSE;
2022 i--;
2023 }
2024 return TRUE;
2025 }
2026
2027 /*
2028 * mu-Algorithmus:
2029 */
2030
2031 // mu-Matrix
mu(matrix A,const ring R)2032 static matrix mu(matrix A, const ring R)
2033 {
2034 int n=MATROWS(A);
2035 assume(MATCOLS(A)==n);
2036 /* Die Funktion erstellt die Matrix mu
2037 *
2038 * Input:
2039 * int n: Dimension der Matrix
2040 * int A: Matrix der Groesse n*n
2041 * int X: Speicherplatz fuer Output
2042 *
2043 * In der Matrix X speichert man die Matrix mu
2044 */
2045
2046 // X als n*n Null-Matrix initalisieren
2047 matrix X=mpNew(n,n);
2048
2049 // Diagonaleintraege von X berrechnen
2050 poly sum = NULL;
2051 for (int i = n-1; i >= 0; i--)
2052 {
2053 MATELEM0(X,i,i) = p_Copy(sum,R);
2054 sum=p_Sub(sum,p_Copy(MATELEM0(A,i,i),R),R);
2055 }
2056 p_Delete(&sum,R);
2057
2058 // Eintraege aus dem oberen Dreieck von A nach X uebertragen
2059 for (int i = n-1; i >=0; i--)
2060 {
2061 for (int j = i+1; j < n; j++)
2062 {
2063 MATELEM0(X,i,j)=p_Copy(MATELEM0(A,i,j),R);
2064 }
2065 }
2066 return X;
2067 }
2068
2069 // Funktion muDet
mp_DetMu(matrix A,const ring R)2070 poly mp_DetMu(matrix A, const ring R)
2071 {
2072 int n=MATROWS(A);
2073 assume(MATCOLS(A)==n);
2074 /*
2075 * Intput:
2076 * int n: Dimension der Matrix
2077 * int A: n*n Matrix
2078 *
2079 * Berechnet n-1 mal: X = mu(X)*A
2080 *
2081 * Output: det(A)
2082 */
2083
2084 //speichere A ab:
2085 matrix workA=mp_Copy(A,R);
2086
2087 // berechen X = mu(X)*A
2088 matrix X;
2089 for (int i = n-1; i >0; i--)
2090 {
2091 X=mu(workA,R);
2092 id_Delete((ideal*)&workA,R);
2093 workA=mp_Mult(X,A,R);
2094 id_Delete((ideal*)&X,R);
2095 }
2096
2097 // berrechne det(A)
2098 poly res;
2099 if (n%2 == 0)
2100 {
2101 res=p_Neg(MATELEM0(workA,0,0),R);
2102 }
2103 else
2104 {
2105 res=MATELEM0(workA,0,0);
2106 }
2107 MATELEM0(workA,0,0)=NULL;
2108 id_Delete((ideal*)&workA,R);
2109 return res;
2110 }
2111
mp_GetAlgorithmDet(matrix m,const ring r)2112 DetVariant mp_GetAlgorithmDet(matrix m, const ring r)
2113 {
2114 if (MATROWS(m)+2*r->N>20+5*rField_is_Zp(r)) return DetMu;
2115 if (MATROWS(m)<10+5*rField_is_Zp(r)) return DetSBareiss;
2116 BOOLEAN isConst=TRUE;
2117 int s=0;
2118 for(int i=MATCOLS(m)*MATROWS(m)-1;i>=0;i--)
2119 {
2120 poly p=m->m[i];
2121 if (p!=NULL)
2122 {
2123 if(!p_IsConstant(p,r)) isConst=FALSE;
2124 s++;
2125 }
2126 }
2127 if (isConst && rField_is_Q(r)) return DetFactory;
2128 if (s*2<MATCOLS(m)*MATROWS(m)) // few entries
2129 return DetSBareiss;
2130 return DetMu;
2131 }
mp_GetAlgorithmDet(const char * s)2132 DetVariant mp_GetAlgorithmDet(const char *s)
2133 {
2134 if (strcmp(s,"Bareiss")==0) return DetBareiss;
2135 if (strcmp(s,"SBareiss")==0) return DetSBareiss;
2136 if (strcmp(s,"Mu")==0) return DetMu;
2137 if (strcmp(s,"Factory")==0) return DetFactory;
2138 WarnS("unknown method for det");
2139 return DetDefault;
2140 }
2141
2142
mp_Det(matrix a,const ring r,DetVariant d)2143 poly mp_Det(matrix a, const ring r, DetVariant d/*=DetDefault*/)
2144 {
2145 if ((MATCOLS(a)==0)
2146 && (MATROWS(a)==0))
2147 return p_One(r);
2148 if (d==DetDefault) d=mp_GetAlgorithmDet(a,r);
2149 switch (d)
2150 {
2151 case DetBareiss: return mp_DetBareiss(a,r);
2152 case DetMu: return mp_DetMu(a,r);
2153 case DetFactory: return singclap_det(a,r);
2154 case DetSBareiss:
2155 {
2156 ideal I=id_Matrix2Module(mp_Copy(a, r),r);
2157 poly p=sm_CallDet(I, r);
2158 id_Delete(&I, r);
2159 return p;
2160 }
2161 default:
2162 WerrorS("unknown algorithm for det");
2163 return NULL;
2164 }
2165 }
2166
sm_Det(ideal a,const ring r,DetVariant d)2167 poly sm_Det(ideal a, const ring r, DetVariant d/*=DetDefault*/)
2168 {
2169 if ((MATCOLS(a)==0)
2170 && (MATROWS(a)==0))
2171 return p_One(r);
2172 if (d==DetDefault) d=mp_GetAlgorithmDet((matrix)a,r);
2173 if (d==DetSBareiss) return sm_CallDet(a,r);
2174 matrix m=id_Module2Matrix(id_Copy(a,r),r);
2175 poly p=mp_Det(m,r,d);
2176 id_Delete((ideal *)&m,r);
2177 return p;
2178 }
2179