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