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}
635
636int mp_Compare(matrix a, matrix b, const ring R)
637{
638 if (MATCOLS(a)<MATCOLS(b)) return -1;
639 else if (MATCOLS(a)>MATCOLS(b)) return 1;
640 if (MATROWS(a)<MATROWS(b)) return -1;
641 else if (MATROWS(a)<MATROWS(b)) return 1;
642
643 unsigned ii=MATCOLS(a)*MATROWS(a)-1;
644 unsigned j=0;
645 int r=0;
646 while (j<=ii)
647 {
648 r=p_Compare(a->m[j],b->m[j],R);
649 if (r!=0) return r;
650 j++;
651 }
652 return r;
653}
654
656{
657 if ((MATCOLS(a)!=MATCOLS(b)) || (MATROWS(a)!=MATROWS(b)))
658 return FALSE;
659 int i=MATCOLS(a)*MATROWS(a)-1;
660 while (i>=0)
661 {
662 if (a->m[i]==NULL)
663 {
664 if (b->m[i]!=NULL) return FALSE;
665 }
666 else if (b->m[i]==NULL) return FALSE;
667 else if (p_Cmp(a->m[i],b->m[i], R)!=0) return FALSE;
668 i--;
669 }
670 i=MATCOLS(a)*MATROWS(a)-1;
671 while (i>=0)
672 {
673 if(!p_EqualPolys(a->m[i],b->m[i], R)) return FALSE;
674 i--;
675 }
676 return TRUE;
677}
678
679/*2
680* insert a monomial into a list, avoid duplicates
681* arguments are destroyed
682*/
683static poly p_Insert(poly p1, poly p2, const ring R)
684{
685 poly a1, p, a2, a;
686 int c;
687
688 if (p1==NULL) return p2;
689 if (p2==NULL) return p1;
690 a1 = p1;
691 a2 = p2;
692 a = p = p_One(R);
693 loop
694 {
695 c = p_Cmp(a1, a2, R);
696 if (c == 1)
697 {
698 a = pNext(a) = a1;
699 pIter(a1);
700 if (a1==NULL)
701 {
702 pNext(a) = a2;
703 break;
704 }
705 }
706 else if (c == -1)
707 {
708 a = pNext(a) = a2;
709 pIter(a2);
710 if (a2==NULL)
711 {
712 pNext(a) = a1;
713 break;
714 }
715 }
716 else
717 {
718 p_LmDelete(&a2, R);
719 a = pNext(a) = a1;
720 pIter(a1);
721 if (a1==NULL)
722 {
723 pNext(a) = a2;
724 break;
725 }
726 else if (a2==NULL)
727 {
728 pNext(a) = a1;
729 break;
730 }
731 }
732 }
733 p_LmDelete(&p, R);
734 return p;
735}
736
737/*2
738*if what == xy the result is the list of all different power products
739* x^i*y^j (i, j >= 0) that appear in fro
740*/
741static poly mp_Select (poly fro, poly what, const ring R)
742{
743 int i;
744 poly h, res;
745 res = NULL;
746 while (fro!=NULL)
747 {
748 h = p_One(R);
749 for (i=1; i<=rVar(R); i++)
750 p_SetExp(h,i, p_GetExp(fro,i, R) * p_GetExp(what, i, R), R);
751 p_SetComp(h, p_GetComp(fro, R), R);
752 p_Setm(h, R);
753 res = p_Insert(h, res, R);
754 fro = fro->next;
755 }
756 return res;
757}
758
759static poly mp_SelectId (ideal I, poly what, const ring R)
760{
761 int i;
762 poly h, res;
763 res = NULL;
764 for(int j=0;j<IDELEMS(I);j++)
765 {
766 poly fro=I->m[j];
767 while (fro!=NULL)
768 {
769 h = p_One(R);
770 for (i=1; i<=rVar(R); i++)
771 p_SetExp(h,i, p_GetExp(fro,i, R) * p_GetExp(what, i, R), R);
772 p_SetComp(h, p_GetComp(fro, R), R);
773 p_Setm(h, R);
774 res = p_Insert(h, res, R);
775 fro = fro->next;
776 }
777 }
778 return res;
779}
780
781/*
782*static void ppp(matrix a)
783*{
784* int j,i,r=a->nrows,c=a->ncols;
785* for(j=1;j<=r;j++)
786* {
787* for(i=1;i<=c;i++)
788* {
789* if(MATELEM(a,j,i)!=NULL) PrintS("X");
790* else PrintS("0");
791* }
792* PrintLn();
793* }
794*}
795*/
796
797static void mp_PartClean(matrix a, int lr, int lc, const ring R)
798{
799 poly *q1;
800 int i,j;
801
802 for (i=lr-1;i>=0;i--)
803 {
804 q1 = &(a->m)[i*a->ncols];
805 for (j=lc-1;j>=0;j--) if(q1[j]) p_Delete(&q1[j], R);
806 }
807}
808
810{
811 if(MATROWS(U)!=MATCOLS(U))
812 return FALSE;
813 for(int i=MATCOLS(U);i>=1;i--)
814 {
815 for(int j=MATCOLS(U); j>=1; j--)
816 {
817 if (i==j)
818 {
819 if (!p_IsUnit(MATELEM(U,i,i), R)) return FALSE;
820 }
821 else if (MATELEM(U,i,j)!=NULL) return FALSE;
822 }
823 }
824 return TRUE;
825}
826
827void iiWriteMatrix(matrix im, const char *n, int dim, const ring r, int spaces)
828{
829 int i,ii = MATROWS(im)-1;
830 int j,jj = MATCOLS(im)-1;
831 poly *pp = im->m;
832
833 for (i=0; i<=ii; i++)
834 {
835 for (j=0; j<=jj; j++)
836 {
837 if (spaces>0)
838 Print("%-*.*s",spaces,spaces," ");
839 if (dim == 2) Print("%s[%u,%u]=",n,i+1,j+1);
840 else if (dim == 1) Print("%s[%u]=",n,j+1);
841 else if (dim == 0) Print("%s=",n);
842 if ((i<ii)||(j<jj)) p_Write(*pp++, r);
843 else p_Write0(*pp, r);
844 }
845 }
846}
847
848char * iiStringMatrix(matrix im, int dim, const ring r, char ch)
849{
850 int i,ii = MATROWS(im);
851 int j,jj = MATCOLS(im);
852 poly *pp = im->m;
853 char ch_s[2];
854 ch_s[0]=ch;
855 ch_s[1]='\0';
856
857 StringSetS("");
858
859 for (i=0; i<ii; i++)
860 {
861 for (j=0; j<jj; j++)
862 {
863 p_String0(*pp++, r);
864 StringAppendS(ch_s);
865 if (dim > 1) StringAppendS("\n");
866 }
867 }
868 char *s=StringEndS();
869 s[strlen(s)- (dim > 1 ? 2 : 1)]='\0';
870 return s;
871}
872
873void mp_Delete(matrix* a, const ring r)
874{
875 id_Delete((ideal *) a, r);
876}
877
878/*
879* C++ classes for Bareiss algorithm
880*/
882{
883 private:
884 int ym, yn;
885 public:
886 float *wrow, *wcol;
888 row_col_weight(int, int);
890};
891
893{
894 ym = i;
895 yn = j;
896 wrow = (float *)omAlloc(i*sizeof(float));
897 wcol = (float *)omAlloc(j*sizeof(float));
898}
900{
901 if (ym!=0)
902 {
903 omFreeSize((ADDRESS)wcol, yn*sizeof(float));
904 omFreeSize((ADDRESS)wrow, ym*sizeof(float));
905 }
906}
907
908/*2
909* a submatrix M of a matrix X[m,n]:
910* 0 <= i < s_m <= a_m
911* 0 <= j < s_n <= a_n
912* M = ( Xarray[qrow[i],qcol[j]] )
913* if a_m = a_n and s_m = s_n
914* det(X) = sign*div^(s_m-1)*det(M)
915* resticted pivot for elimination
916* 0 <= j < piv_s
917*/
919{
920 private:
922 int *qrow, *qcol;
923 poly *Xarray;
924 ring _R;
925 void mpInitMat();
926 poly * mpRowAdr(int r)
927 { return &(Xarray[a_n*qrow[r]]); }
928 poly * mpColAdr(int c)
929 { return &(Xarray[qcol[c]]); }
930 void mpRowWeight(float *);
931 void mpColWeight(float *);
932 void mpRowSwap(int, int);
933 void mpColSwap(int, int);
934 public:
936 mp_permmatrix(matrix, ring);
939 int mpGetRow();
940 int mpGetCol();
941 int mpGetRdim() { return s_m; }
942 int mpGetCdim() { return s_n; }
943 int mpGetSign() { return sign; }
944 void mpSetSearch(int s);
945 void mpSaveArray() { Xarray = NULL; }
946 poly mpGetElem(int, int);
947 void mpSetElem(poly, int, int);
948 void mpDelElem(int, int);
949 void mpElimBareiss(poly);
953 void mpRowReorder();
954 void mpColReorder();
955};
957{
958 a_m = A->nrows;
959 a_n = A->ncols;
960 this->mpInitMat();
961 Xarray = A->m;
962 _R=R;
963}
964
966{
967 poly p, *athis, *aM;
968 int i, j;
969
970 _R=M->_R;
971 a_m = M->s_m;
972 a_n = M->s_n;
973 sign = M->sign;
974 this->mpInitMat();
975 Xarray = (poly *)omAlloc0(a_m*a_n*sizeof(poly));
976 for (i=a_m-1; i>=0; i--)
977 {
978 athis = this->mpRowAdr(i);
979 aM = M->mpRowAdr(i);
980 for (j=a_n-1; j>=0; j--)
981 {
982 p = aM[M->qcol[j]];
983 if (p)
984 {
985 athis[j] = p_Copy(p,_R);
986 }
987 }
988 }
989}
990
992{
993 int k;
994
995 if (a_m != 0)
996 {
997 omFreeSize((ADDRESS)qrow,a_m*sizeof(int));
998 omFreeSize((ADDRESS)qcol,a_n*sizeof(int));
999 if (Xarray != NULL)
1000 {
1001 for (k=a_m*a_n-1; k>=0; k--)
1002 p_Delete(&Xarray[k],_R);
1003 omFreeSize((ADDRESS)Xarray,a_m*a_n*sizeof(poly));
1004 }
1005 }
1006}
1007
1008
1009static float mp_PolyWeight(poly p, const ring r);
1011{
1012 poly p, *a;
1013 int i, j;
1014 float count;
1015
1016 for (j=s_n; j>=0; j--)
1017 {
1018 a = this->mpColAdr(j);
1019 count = 0.0;
1020 for(i=s_m; i>=0; i--)
1021 {
1022 p = a[a_n*qrow[i]];
1023 if (p)
1024 count += mp_PolyWeight(p,_R);
1025 }
1026 wcol[j] = count;
1027 }
1028}
1030{
1031 poly p, *a;
1032 int i, j;
1033 float count;
1034
1035 for (i=s_m; i>=0; i--)
1036 {
1037 a = this->mpRowAdr(i);
1038 count = 0.0;
1039 for(j=s_n; j>=0; j--)
1040 {
1041 p = a[qcol[j]];
1042 if (p)
1043 count += mp_PolyWeight(p,_R);
1044 }
1045 wrow[i] = count;
1046 }
1047}
1048
1049void mp_permmatrix::mpRowSwap(int i1, int i2)
1050{
1051 poly p, *a1, *a2;
1052 int j;
1053
1054 a1 = &(Xarray[a_n*i1]);
1055 a2 = &(Xarray[a_n*i2]);
1056 for (j=a_n-1; j>= 0; j--)
1057 {
1058 p = a1[j];
1059 a1[j] = a2[j];
1060 a2[j] = p;
1061 }
1062}
1063
1064void mp_permmatrix::mpColSwap(int j1, int j2)
1065{
1066 poly p, *a1, *a2;
1067 int i, k = a_n*a_m;
1068
1069 a1 = &(Xarray[j1]);
1070 a2 = &(Xarray[j2]);
1071 for (i=0; i< k; i+=a_n)
1072 {
1073 p = a1[i];
1074 a1[i] = a2[i];
1075 a2[i] = p;
1076 }
1077}
1079{
1080 int k;
1081
1082 s_m = a_m;
1083 s_n = a_n;
1084 piv_s = 0;
1085 qrow = (int *)omAlloc(a_m*sizeof(int));
1086 qcol = (int *)omAlloc(a_n*sizeof(int));
1087 for (k=a_m-1; k>=0; k--) qrow[k] = k;
1088 for (k=a_n-1; k>=0; k--) qcol[k] = k;
1089}
1090
1092{
1093 int k, j, j1, j2;
1094
1095 if (a_n > a_m)
1096 k = a_n - a_m;
1097 else
1098 k = 0;
1099 for (j=a_n-1; j>=k; j--)
1100 {
1101 j1 = qcol[j];
1102 if (j1 != j)
1103 {
1104 this->mpColSwap(j1, j);
1105 j2 = 0;
1106 while (qcol[j2] != j) j2++;
1107 qcol[j2] = j1;
1108 }
1109 }
1110}
1111
1113{
1114 int k, i, i1, i2;
1115
1116 if (a_m > a_n)
1117 k = a_m - a_n;
1118 else
1119 k = 0;
1120 for (i=a_m-1; i>=k; i--)
1121 {
1122 i1 = qrow[i];
1123 if (i1 != i)
1124 {
1125 this->mpRowSwap(i1, i);
1126 i2 = 0;
1127 while (qrow[i2] != i) i2++;
1128 qrow[i2] = i1;
1129 }
1130 }
1131}
1132
1133/*
1134* perform replacement for pivot strategy in Bareiss algorithm
1135* change sign of determinant
1136*/
1137static void mpReplace(int j, int n, int &sign, int *perm)
1138{
1139 int k;
1140
1141 if (j != n)
1142 {
1143 k = perm[n];
1144 perm[n] = perm[j];
1145 perm[j] = k;
1146 sign = -sign;
1147 }
1148}
1149/*2
1150* pivot strategy for Bareiss algorithm
1151*/
1153{
1154 poly p, *a;
1155 int i, j, iopt, jopt;
1156 float sum, f1, f2, fo, r, ro, lp;
1157 float *dr = C->wrow, *dc = C->wcol;
1158
1159 fo = 1.0e20;
1160 ro = 0.0;
1161 iopt = jopt = -1;
1162
1163 s_n--;
1164 s_m--;
1165 if (s_m == 0)
1166 return 0;
1167 if (s_n == 0)
1168 {
1169 for(i=s_m; i>=0; i--)
1170 {
1171 p = this->mpRowAdr(i)[qcol[0]];
1172 if (p)
1173 {
1174 f1 = mp_PolyWeight(p,_R);
1175 if (f1 < fo)
1176 {
1177 fo = f1;
1178 if (iopt >= 0)
1179 p_Delete(&(this->mpRowAdr(iopt)[qcol[0]]),_R);
1180 iopt = i;
1181 }
1182 else
1183 p_Delete(&(this->mpRowAdr(i)[qcol[0]]),_R);
1184 }
1185 }
1186 if (iopt >= 0)
1187 mpReplace(iopt, s_m, sign, qrow);
1188 return 0;
1189 }
1190 this->mpRowWeight(dr);
1191 this->mpColWeight(dc);
1192 sum = 0.0;
1193 for(i=s_m; i>=0; i--)
1194 sum += dr[i];
1195 for(i=s_m; i>=0; i--)
1196 {
1197 r = dr[i];
1198 a = this->mpRowAdr(i);
1199 for(j=s_n; j>=0; j--)
1200 {
1201 p = a[qcol[j]];
1202 if (p)
1203 {
1204 lp = mp_PolyWeight(p,_R);
1205 ro = r - lp;
1206 f1 = ro * (dc[j]-lp);
1207 if (f1 != 0.0)
1208 {
1209 f2 = lp * (sum - ro - dc[j]);
1210 f2 += f1;
1211 }
1212 else
1213 f2 = lp-r-dc[j];
1214 if (f2 < fo)
1215 {
1216 fo = f2;
1217 iopt = i;
1218 jopt = j;
1219 }
1220 }
1221 }
1222 }
1223 if (iopt < 0)
1224 return 0;
1225 mpReplace(iopt, s_m, sign, qrow);
1226 mpReplace(jopt, s_n, sign, qcol);
1227 return 1;
1228}
1230{
1231 return Xarray[a_n*qrow[r]+qcol[c]];
1232}
1233
1234/*
1235* the Bareiss-type elimination with division by div (div != NULL)
1236*/
1238{
1239 poly piv, elim, q1, q2, *ap, *a;
1240 int i, j, jj;
1241
1242 ap = this->mpRowAdr(s_m);
1243 piv = ap[qcol[s_n]];
1244 for(i=s_m-1; i>=0; i--)
1245 {
1246 a = this->mpRowAdr(i);
1247 elim = a[qcol[s_n]];
1248 if (elim != NULL)
1249 {
1250 elim = p_Neg(elim,_R);
1251 for (j=s_n-1; j>=0; j--)
1252 {
1253 q2 = NULL;
1254 jj = qcol[j];
1255 if (ap[jj] != NULL)
1256 {
1257 q2 = SM_MULT(ap[jj], elim, div,_R);
1258 if (a[jj] != NULL)
1259 {
1260 q1 = SM_MULT(a[jj], piv, div,_R);
1261 p_Delete(&a[jj],_R);
1262 q2 = p_Add_q(q2, q1, _R);
1263 }
1264 }
1265 else if (a[jj] != NULL)
1266 {
1267 q2 = SM_MULT(a[jj], piv, div, _R);
1268 }
1269 if ((q2!=NULL) && div)
1270 SM_DIV(q2, div, _R);
1271 a[jj] = q2;
1272 }
1273 p_Delete(&a[qcol[s_n]], _R);
1274 }
1275 else
1276 {
1277 for (j=s_n-1; j>=0; j--)
1278 {
1279 jj = qcol[j];
1280 if (a[jj] != NULL)
1281 {
1282 q2 = SM_MULT(a[jj], piv, div, _R);
1283 p_Delete(&a[jj], _R);
1284 if (div)
1285 SM_DIV(q2, div, _R);
1286 a[jj] = q2;
1287 }
1288 }
1289 }
1290 }
1291}
1292/*
1293* weigth of a polynomial, for pivot strategy
1294*/
1295static float mp_PolyWeight(poly p, const ring r)
1296{
1297 int i;
1298 float res;
1299
1300 if (pNext(p) == NULL)
1301 {
1302 res = (float)n_Size(pGetCoeff(p),r->cf);
1303 for (i=r->N;i>0;i--)
1304 {
1305 if(p_GetExp(p,i,r)!=0)
1306 {
1307 res += 2.0;
1308 break;
1309 }
1310 }
1311 }
1312 else
1313 {
1314 res = 0.0;
1315 do
1316 {
1317 res += (float)n_Size(pGetCoeff(p),r->cf)+2.0;
1318 pIter(p);
1319 }
1320 while (p);
1321 }
1322 return res;
1323}
1324/*
1325* find best row
1326*/
1327static int mp_PivBar(matrix a, int lr, int lc, const ring r)
1328{
1329 float f1, f2;
1330 poly *q1;
1331 int i,j,io;
1332
1333 io = -1;
1334 f1 = 1.0e30;
1335 for (i=lr-1;i>=0;i--)
1336 {
1337 q1 = &(a->m)[i*a->ncols];
1338 f2 = 0.0;
1339 for (j=lc-1;j>=0;j--)
1340 {
1341 if (q1[j]!=NULL)
1342 f2 += mp_PolyWeight(q1[j],r);
1343 }
1344 if ((f2!=0.0) && (f2<f1))
1345 {
1346 f1 = f2;
1347 io = i;
1348 }
1349 }
1350 if (io<0) return 0;
1351 else return io+1;
1352}
1353
1354static void mpSwapRow(matrix a, int pos, int lr, int lc)
1355{
1356 poly sw;
1357 int j;
1358 poly* a2 = a->m;
1359 poly* a1 = &a2[a->ncols*(pos-1)];
1360
1361 a2 = &a2[a->ncols*(lr-1)];
1362 for (j=lc-1; j>=0; j--)
1363 {
1364 sw = a1[j];
1365 a1[j] = a2[j];
1366 a2[j] = sw;
1367 }
1368}
1369
1370/*2
1371* prepare one step of 'Bareiss' algorithm
1372* for application in minor
1373*/
1374static int mp_PrepareRow (matrix a, int lr, int lc, const ring R)
1375{
1376 int r;
1377
1378 r = mp_PivBar(a,lr,lc,R);
1379 if(r==0) return 0;
1380 if(r<lr) mpSwapRow(a, r, lr, lc);
1381 return 1;
1382}
1383
1384/*
1385* find pivot in the last row
1386*/
1387static int mp_PivRow(matrix a, int lr, int lc, const ring r)
1388{
1389 float f1, f2;
1390 poly *q1;
1391 int j,jo;
1392
1393 jo = -1;
1394 f1 = 1.0e30;
1395 q1 = &(a->m)[(lr-1)*a->ncols];
1396 for (j=lc-1;j>=0;j--)
1397 {
1398 if (q1[j]!=NULL)
1399 {
1400 f2 = mp_PolyWeight(q1[j],r);
1401 if (f2<f1)
1402 {
1403 f1 = f2;
1404 jo = j;
1405 }
1406 }
1407 }
1408 if (jo<0) return 0;
1409 else return jo+1;
1410}
1411
1412static void mpSwapCol(matrix a, int pos, int lr, int lc)
1413{
1414 poly sw;
1415 int j;
1416 poly* a2 = a->m;
1417 poly* a1 = &a2[pos-1];
1418
1419 a2 = &a2[lc-1];
1420 for (j=a->ncols*(lr-1); j>=0; j-=a->ncols)
1421 {
1422 sw = a1[j];
1423 a1[j] = a2[j];
1424 a2[j] = sw;
1425 }
1426}
1427
1428/*2
1429* prepare one step of 'Bareiss' algorithm
1430* for application in minor
1431*/
1432static int mp_PreparePiv (matrix a, int lr, int lc,const ring r)
1433{
1434 int c;
1435
1436 c = mp_PivRow(a, lr, lc,r);
1437 if(c==0) return 0;
1438 if(c<lc) mpSwapCol(a, c, lr, lc);
1439 return 1;
1440}
1441
1442static void mp_ElimBar(matrix a0, matrix re, poly div, int lr, int lc, const ring R)
1443{
1444 int r=lr-1, c=lc-1;
1445 poly *b = a0->m, *x = re->m;
1446 poly piv, elim, q1, *ap, *a, *q;
1447 int i, j;
1448
1449 ap = &b[r*a0->ncols];
1450 piv = ap[c];
1451 for(j=c-1; j>=0; j--)
1452 if (ap[j] != NULL) ap[j] = p_Neg(ap[j],R);
1453 for(i=r-1; i>=0; i--)
1454 {
1455 a = &b[i*a0->ncols];
1456 q = &x[i*re->ncols];
1457 if (a[c] != NULL)
1458 {
1459 elim = a[c];
1460 for (j=c-1; j>=0; j--)
1461 {
1462 q1 = NULL;
1463 if (a[j] != NULL)
1464 {
1465 q1 = sm_MultDiv(a[j], piv, div,R);
1466 if (ap[j] != NULL)
1467 {
1468 poly q2 = sm_MultDiv(ap[j], elim, div, R);
1469 q1 = p_Add_q(q1,q2,R);
1470 }
1471 }
1472 else if (ap[j] != NULL)
1473 q1 = sm_MultDiv(ap[j], elim, div, R);
1474 if (q1 != NULL)
1475 {
1476 if (div)
1477 sm_SpecialPolyDiv(q1, div,R);
1478 q[j] = q1;
1479 }
1480 }
1481 }
1482 else
1483 {
1484 for (j=c-1; j>=0; j--)
1485 {
1486 if (a[j] != NULL)
1487 {
1488 q1 = sm_MultDiv(a[j], piv, div, R);
1489 if (div)
1490 sm_SpecialPolyDiv(q1, div, R);
1491 q[j] = q1;
1492 }
1493 }
1494 }
1495 }
1496}
1497
1498/*2*/
1499/// entries of a are minors and go to result (only if not in R)
1500void mp_MinorToResult(ideal result, int &elems, matrix a, int r, int c,
1501 ideal R, const ring)
1502{
1503 poly *q1;
1504 int e=IDELEMS(result);
1505 int i,j;
1506
1507 if (R != NULL)
1508 {
1509 for (i=r-1;i>=0;i--)
1510 {
1511 q1 = &(a->m)[i*a->ncols];
1512 //for (j=c-1;j>=0;j--)
1513 //{
1514 // if (q1[j]!=NULL) q1[j] = kNF(R,currRing->qideal,q1[j]);
1515 //}
1516 }
1517 }
1518 for (i=r-1;i>=0;i--)
1519 {
1520 q1 = &(a->m)[i*a->ncols];
1521 for (j=c-1;j>=0;j--)
1522 {
1523 if (q1[j]!=NULL)
1524 {
1525 if (elems>=e)
1526 {
1527 pEnlargeSet(&(result->m),e,e);
1528 e += e;
1529 IDELEMS(result) =e;
1530 }
1531 result->m[elems] = q1[j];
1532 q1[j] = NULL;
1533 elems++;
1534 }
1535 }
1536 }
1537}
1538/*
1539// from linalg_from_matpol.cc: TODO: compare with above & remove...
1540void mp_MinorToResult(ideal result, int &elems, matrix a, int r, int c,
1541 ideal R, const ring R)
1542{
1543 poly *q1;
1544 int e=IDELEMS(result);
1545 int i,j;
1546
1547 if (R != NULL)
1548 {
1549 for (i=r-1;i>=0;i--)
1550 {
1551 q1 = &(a->m)[i*a->ncols];
1552 for (j=c-1;j>=0;j--)
1553 {
1554 if (q1[j]!=NULL) q1[j] = kNF(R,currRing->qideal,q1[j]);
1555 }
1556 }
1557 }
1558 for (i=r-1;i>=0;i--)
1559 {
1560 q1 = &(a->m)[i*a->ncols];
1561 for (j=c-1;j>=0;j--)
1562 {
1563 if (q1[j]!=NULL)
1564 {
1565 if (elems>=e)
1566 {
1567 if(e<SIZE_OF_SYSTEM_PAGE)
1568 {
1569 pEnlargeSet(&(result->m),e,e);
1570 e += e;
1571 }
1572 else
1573 {
1574 pEnlargeSet(&(result->m),e,SIZE_OF_SYSTEM_PAGE);
1575 e += SIZE_OF_SYSTEM_PAGE;
1576 }
1577 IDELEMS(result) =e;
1578 }
1579 result->m[elems] = q1[j];
1580 q1[j] = NULL;
1581 elems++;
1582 }
1583 }
1584 }
1585}
1586*/
1587
1588static void mpFinalClean(matrix a)
1589{
1590 omFreeSize((ADDRESS)a->m,a->nrows*a->ncols*sizeof(poly));
1592}
1593
1594/*2*/
1595/// produces recursively the ideal of all arxar-minors of a
1596void mp_RecMin(int ar,ideal result,int &elems,matrix a,int lr,int lc,
1597 poly barDiv, ideal R, const ring r)
1598{
1599 int k;
1600 int kr=lr-1,kc=lc-1;
1601 matrix nextLevel=mpNew(kr,kc);
1602
1603 loop
1604 {
1605/*--- look for an optimal row and bring it to last position ------------*/
1606 if(mp_PrepareRow(a,lr,lc,r)==0) break;
1607/*--- now take all pivots from the last row ------------*/
1608 k = lc;
1609 loop
1610 {
1611 if(mp_PreparePiv(a,lr,k,r)==0) break;
1612 mp_ElimBar(a,nextLevel,barDiv,lr,k,r);
1613 k--;
1614 if (ar>1)
1615 {
1616 mp_RecMin(ar-1,result,elems,nextLevel,kr,k,a->m[kr*a->ncols+k],R,r);
1617 mp_PartClean(nextLevel,kr,k, r);
1618 }
1619 else mp_MinorToResult(result,elems,nextLevel,kr,k,R,r);
1620 if (ar>k-1) break;
1621 }
1622 if (ar>=kr) break;
1623/*--- now we have to take out the last row...------------*/
1624 lr = kr;
1625 kr--;
1626 }
1627 mpFinalClean(nextLevel);
1628}
1629/*
1630// from linalg_from_matpol.cc: TODO: compare with above & remove...
1631void mp_RecMin(int ar,ideal result,int &elems,matrix a,int lr,int lc,
1632 poly barDiv, ideal R, const ring R)
1633{
1634 int k;
1635 int kr=lr-1,kc=lc-1;
1636 matrix nextLevel=mpNew(kr,kc);
1637
1638 loop
1639 {
1640// --- look for an optimal row and bring it to last position ------------
1641 if(mpPrepareRow(a,lr,lc)==0) break;
1642// --- now take all pivots from the last row ------------
1643 k = lc;
1644 loop
1645 {
1646 if(mpPreparePiv(a,lr,k)==0) break;
1647 mpElimBar(a,nextLevel,barDiv,lr,k);
1648 k--;
1649 if (ar>1)
1650 {
1651 mpRecMin(ar-1,result,elems,nextLevel,kr,k,a->m[kr*a->ncols+k],R);
1652 mpPartClean(nextLevel,kr,k);
1653 }
1654 else mpMinorToResult(result,elems,nextLevel,kr,k,R);
1655 if (ar>k-1) break;
1656 }
1657 if (ar>=kr) break;
1658// --- now we have to take out the last row...------------
1659 lr = kr;
1660 kr--;
1661 }
1662 mpFinalClean(nextLevel);
1663}
1664*/
1665
1666/*2*/
1667/// returns the determinant of the matrix m;
1668/// uses Bareiss algorithm
1669poly mp_DetBareiss (matrix a, const ring r)
1670{
1671 int s;
1672 poly div, res;
1673 if (MATROWS(a) != MATCOLS(a))
1674 {
1675 Werror("det of %d x %d matrix",MATROWS(a),MATCOLS(a));
1676 return NULL;
1677 }
1678 matrix c = mp_Copy(a,r);
1679 mp_permmatrix *Bareiss = new mp_permmatrix(c,r);
1680 row_col_weight w(Bareiss->mpGetRdim(), Bareiss->mpGetCdim());
1681
1682 /* Bareiss */
1683 div = NULL;
1684 while(Bareiss->mpPivotBareiss(&w))
1685 {
1686 Bareiss->mpElimBareiss(div);
1687 div = Bareiss->mpGetElem(Bareiss->mpGetRdim(), Bareiss->mpGetCdim());
1688 }
1689 Bareiss->mpRowReorder();
1690 Bareiss->mpColReorder();
1691 Bareiss->mpSaveArray();
1692 s = Bareiss->mpGetSign();
1693 delete Bareiss;
1694
1695 /* result */
1696 res = MATELEM(c,1,1);
1697 MATELEM(c,1,1) = NULL;
1698 id_Delete((ideal *)&c,r);
1699 if (s < 0)
1700 res = p_Neg(res,r);
1701 return res;
1702}
1703/*
1704// from linalg_from_matpol.cc: TODO: compare with above & remove...
1705poly mp_DetBareiss (matrix a, const ring R)
1706{
1707 int s;
1708 poly div, res;
1709 if (MATROWS(a) != MATCOLS(a))
1710 {
1711 Werror("det of %d x %d matrix",MATROWS(a),MATCOLS(a));
1712 return NULL;
1713 }
1714 matrix c = mp_Copy(a, R);
1715 mp_permmatrix *Bareiss = new mp_permmatrix(c, R);
1716 row_col_weight w(Bareiss->mpGetRdim(), Bareiss->mpGetCdim());
1717
1718 // Bareiss
1719 div = NULL;
1720 while(Bareiss->mpPivotBareiss(&w))
1721 {
1722 Bareiss->mpElimBareiss(div);
1723 div = Bareiss->mpGetElem(Bareiss->mpGetRdim(), Bareiss->mpGetCdim());
1724 }
1725 Bareiss->mpRowReorder();
1726 Bareiss->mpColReorder();
1727 Bareiss->mpSaveArray();
1728 s = Bareiss->mpGetSign();
1729 delete Bareiss;
1730
1731 // result
1732 res = MATELEM(c,1,1);
1733 MATELEM(c,1,1) = NULL;
1734 id_Delete((ideal *)&c, R);
1735 if (s < 0)
1736 res = p_Neg(res, R);
1737 return res;
1738}
1739*/
1740
1741/*2
1742* compute all ar-minors of the matrix a
1743*/
1744matrix mp_Wedge(matrix a, int ar, const ring R)
1745{
1746 int i,j,k,l;
1747 int *rowchoise,*colchoise;
1748 BOOLEAN rowch,colch;
1749 matrix result;
1750 matrix tmp;
1751 poly p;
1752
1753 i = binom(a->nrows,ar);
1754 j = binom(a->ncols,ar);
1755
1756 rowchoise=(int *)omAlloc(ar*sizeof(int));
1757 colchoise=(int *)omAlloc(ar*sizeof(int));
1758 result = mpNew(i,j);
1759 tmp = mpNew(ar,ar);
1760 l = 1; /* k,l:the index in result*/
1761 idInitChoise(ar,1,a->nrows,&rowch,rowchoise);
1762 while (!rowch)
1763 {
1764 k=1;
1765 idInitChoise(ar,1,a->ncols,&colch,colchoise);
1766 while (!colch)
1767 {
1768 for (i=1; i<=ar; i++)
1769 {
1770 for (j=1; j<=ar; j++)
1771 {
1772 MATELEM(tmp,i,j) = MATELEM(a,rowchoise[i-1],colchoise[j-1]);
1773 }
1774 }
1775 p = mp_DetBareiss(tmp, R);
1776 if ((k+l) & 1) p=p_Neg(p, R);
1777 MATELEM(result,l,k) = p;
1778 k++;
1779 idGetNextChoise(ar,a->ncols,&colch,colchoise);
1780 }
1781 idGetNextChoise(ar,a->nrows,&rowch,rowchoise);
1782 l++;
1783 }
1784
1785 /*delete the matrix tmp*/
1786 for (i=1; i<=ar; i++)
1787 {
1788 for (j=1; j<=ar; j++) MATELEM(tmp,i,j) = NULL;
1789 }
1790 id_Delete((ideal *) &tmp, R);
1791 return (result);
1792}
1793
1794// helper for sm_Tensor
1795// destroyes f, keeps B
1796static ideal sm_MultAndShift(poly f, ideal B, int s, const ring r)
1797{
1798 assume(f!=NULL);
1799 ideal res=idInit(IDELEMS(B),B->rank+s);
1800 int q=IDELEMS(B); // p x q
1801 for(int j=0;j<q;j++)
1802 {
1803 poly h=pp_Mult_qq(f,B->m[j],r);
1804 if (h!=NULL)
1805 {
1806 if (s>0) p_Shift(&h,s,r);
1807 res->m[j]=h;
1808 }
1809 }
1810 p_Delete(&f,r);
1811 return res;
1812}
1813// helper for sm_Tensor
1814// updates res, destroyes contents of sm
1815static void sm_AddSubMat(ideal res, ideal sm, int col, const ring r)
1816{
1817 for(int i=0;i<IDELEMS(sm);i++)
1818 {
1819 res->m[col+i]=p_Add_q(res->m[col+i],sm->m[i],r);
1820 sm->m[i]=NULL;
1821 }
1822}
1823
1824ideal sm_Tensor(ideal A, ideal B, const ring r)
1825{
1826 // size of the result m*p x n*q
1827 int n=IDELEMS(A); // m x n
1828 int m=A->rank;
1829 int q=IDELEMS(B); // p x q
1830 int p=B->rank;
1831 ideal res=idInit(n*q,m*p);
1832 poly *a=(poly*)omAlloc(m*sizeof(poly));
1833 for(int i=0; i<n; i++)
1834 {
1835 memset(a,0,m*sizeof(poly));
1836 p_Vec2Array(A->m[i],a,m,r);
1837 for(int j=0;j<m;j++)
1838 {
1839 if (a[j]!=NULL)
1840 {
1841 ideal sm=sm_MultAndShift(a[j], // A_i_j
1842 B,
1843 j*p, // shift j*p down
1844 r);
1845 sm_AddSubMat(res,sm,i*q,r); // add this columns to col i*q ff
1846 id_Delete(&sm,r); // delete the now empty ideal
1847 }
1848 }
1849 }
1850 omFreeSize(a,m*sizeof(poly));
1851 return res;
1852}
1853// --------------------------------------------------------------------------
1854/****************************************
1855* Computer Algebra System SINGULAR *
1856****************************************/
1857
1858/*
1859* ABSTRACT: basic operation for sparse matrices:
1860* type: ideal (of column vectors)
1861* nrows: I->rank, ncols: IDELEMS(I)
1862*/
1863
1864ideal sm_Add(ideal a, ideal b, const ring R)
1865{
1866 assume(IDELEMS(a)==IDELEMS(b));
1867 assume(a->rank==b->rank);
1868 ideal c=idInit(IDELEMS(a),a->rank);
1869 for (int k=IDELEMS(a)-1; k>=0; k--)
1870 c->m[k] = p_Add_q(p_Copy(a->m[k], R), p_Copy(b->m[k], R), R);
1871 return c;
1872}
1873
1874ideal sm_Sub(ideal a, ideal b, const ring R)
1875{
1876 assume(IDELEMS(a)==IDELEMS(b));
1877 assume(a->rank==b->rank);
1878 ideal c=idInit(IDELEMS(a),a->rank);
1879 for (int k=IDELEMS(a)-1; k>=0; k--)
1880 c->m[k] = p_Sub(p_Copy(a->m[k], R), p_Copy(b->m[k], R), R);
1881 return c;
1882}
1883
1884ideal sm_Mult(ideal a, ideal b, const ring R)
1885{
1886 int i, j, k;
1887 int m = a->rank;
1888 int p = IDELEMS(a);
1889 int q = IDELEMS(b);
1890
1891 assume (IDELEMS(a)==b->rank);
1892 ideal c = idInit(q,m);
1893
1894 for (i=0; i<m; i++)
1895 {
1896 for (k=0; k<p; k++)
1897 {
1898 poly aik;
1899 if ((aik=SMATELEM(a,i,k,R))!=NULL)
1900 {
1901 for (j=0; j<q; j++)
1902 {
1903 poly bkj=SMATELEM(b,k,j,R);
1904 if (bkj!=NULL)
1905 {
1906 poly s = p_Mult_q(p_Copy(aik,R) /*SMATELEM(a,i,k)*/, bkj/*SMATELEM(b,k,j)*/, R);
1907 if (s!=NULL) p_SetComp(s,i+1,R);
1908 c->m[j]=p_Add_q(c->m[j],s, R);
1909 }
1910 }
1911 p_Delete(&aik,R);
1912 }
1913 }
1914 }
1915 for(i=q-1;i>=0;i--) p_Normalize(c->m[i], R);
1916 return c;
1917}
1918
1919ideal sm_Flatten(ideal a, const ring R)
1920{
1921 if (IDELEMS(a)==0) return id_Copy(a,R);
1922 ideal res=idInit(1,IDELEMS(a)*a->rank);
1923 for(int i=0;i<IDELEMS(a);i++)
1924 {
1925 if(a->m[i]!=NULL)
1926 {
1927 poly p=p_Copy(a->m[i],R);
1928 if (i==0) res->m[0]=p;
1929 else
1930 {
1931 p_Shift(&p,i*a->rank,R);
1932 res->m[0]=p_Add_q(res->m[0],p,R);
1933 }
1934 }
1935 }
1936 return res;
1937}
1938
1939ideal sm_UnFlatten(ideal a, int col, const ring R)
1940{
1941 if ((IDELEMS(a)!=1)
1942 ||((a->rank % col)!=0))
1943 {
1944 Werror("wrong format: %d x %d for unflatten",(int)a->rank,IDELEMS(a));
1945 return NULL;
1946 }
1947 int row=a->rank/col;
1948 ideal res=idInit(col,row);
1949 poly p=a->m[0];
1950 while(p!=NULL)
1951 {
1952 poly h=p_Head(p,R);
1953 int comp=p_GetComp(h,R);
1954 int c=(comp-1)/row;
1955 int r=comp%row; if (r==0) r=row;
1956 p_SetComp(h,r,R); p_SetmComp(h,R);
1957 res->m[c]=p_Add_q(res->m[c],h,R);
1958 pIter(p);
1959 }
1960 return res;
1961}
1962
1963/*2
1964*returns the trace of matrix a
1965*/
1966poly sm_Trace ( ideal a, const ring R)
1967{
1968 int i;
1969 int n = (IDELEMS(a)<a->rank) ? IDELEMS(a) : a->rank;
1970 poly t = NULL;
1971
1972 for (i=0; i<=n; i++)
1973 t = p_Add_q(t, p_Copy(SMATELEM(a,i,i,R), R), R);
1974 return t;
1975}
1976
1977int sm_Compare(ideal a, ideal b, const ring R)
1978{
1979 if (IDELEMS(a)<IDELEMS(b)) return -1;
1980 else if (IDELEMS(a)>IDELEMS(b)) return 1;
1981 if ((a->rank)<(b->rank)) return -1;
1982 else if ((a->rank)<(b->rank)) return 1;
1983
1984 unsigned ii=IDELEMS(a)-1;
1985 unsigned j=0;
1986 int r=0;
1987 while (j<=ii)
1988 {
1989 r=p_Compare(a->m[j],b->m[j],R);
1990 if (r!=0) return r;
1991 j++;
1992 }
1993 return r;
1994}
1995
1996BOOLEAN sm_Equal(ideal a, ideal b, const ring R)
1997{
1998 if ((a->rank!=b->rank) || (IDELEMS(a)!=IDELEMS(b)))
1999 return FALSE;
2000 int i=IDELEMS(a)-1;
2001 while (i>=0)
2002 {
2003 if (a->m[i]==NULL)
2004 {
2005 if (b->m[i]!=NULL) return FALSE;
2006 }
2007 else if (b->m[i]==NULL) return FALSE;
2008 else if (p_Cmp(a->m[i],b->m[i], R)!=0) return FALSE;
2009 i--;
2010 }
2011 i=IDELEMS(a)-1;
2012 while (i>=0)
2013 {
2014 if(!p_EqualPolys(a->m[i],b->m[i], R)) return FALSE;
2015 i--;
2016 }
2017 return TRUE;
2018}
2019
2020/*
2021* mu-Algorithmus:
2022*/
2023
2024// mu-Matrix
2025static matrix mu(matrix A, const ring R)
2026{
2027 int n=MATROWS(A);
2028 assume(MATCOLS(A)==n);
2029 /* Die Funktion erstellt die Matrix mu
2030 *
2031 * Input:
2032 * int n: Dimension der Matrix
2033 * int A: Matrix der Groesse n*n
2034 * int X: Speicherplatz fuer Output
2035 *
2036 * In der Matrix X speichert man die Matrix mu
2037 */
2038
2039 // X als n*n Null-Matrix initalisieren
2040 matrix X=mpNew(n,n);
2041
2042 // Diagonaleintraege von X berrechnen
2043 poly sum = NULL;
2044 for (int i = n-1; i >= 0; i--)
2045 {
2046 MATELEM0(X,i,i) = p_Copy(sum,R);
2047 sum=p_Sub(sum,p_Copy(MATELEM0(A,i,i),R),R);
2048 }
2049 p_Delete(&sum,R);
2050
2051 // Eintraege aus dem oberen Dreieck von A nach X uebertragen
2052 for (int i = n-1; i >=0; i--)
2053 {
2054 for (int j = i+1; j < n; j++)
2055 {
2056 MATELEM0(X,i,j)=p_Copy(MATELEM0(A,i,j),R);
2057 }
2058 }
2059 return X;
2060}
2061
2062// Funktion muDet
2063poly mp_DetMu(matrix A, const ring R)
2064{
2065 int n=MATROWS(A);
2066 assume(MATCOLS(A)==n);
2067 /*
2068 * Intput:
2069 * int n: Dimension der Matrix
2070 * int A: n*n Matrix
2071 *
2072 * Berechnet n-1 mal: X = mu(X)*A
2073 *
2074 * Output: det(A)
2075 */
2076
2077 //speichere A ab:
2078 matrix workA=mp_Copy(A,R);
2079
2080 // berechen X = mu(X)*A
2081 matrix X;
2082 for (int i = n-1; i >0; i--)
2083 {
2084 X=mu(workA,R);
2085 id_Delete((ideal*)&workA,R);
2086 workA=mp_Mult(X,A,R);
2087 id_Delete((ideal*)&X,R);
2088 }
2089
2090 // berrechne det(A)
2091 poly res;
2092 if (n%2 == 0)
2093 {
2094 res=p_Neg(MATELEM0(workA,0,0),R);
2095 }
2096 else
2097 {
2098 res=MATELEM0(workA,0,0);
2099 }
2100 MATELEM0(workA,0,0)=NULL;
2101 id_Delete((ideal*)&workA,R);
2102 return res;
2103}
2104
2106{
2107 if (MATROWS(m)+2*r->N>20+5*rField_is_Zp(r)) return DetMu;
2108 if (MATROWS(m)<10+5*rField_is_Zp(r)) return DetSBareiss;
2109 BOOLEAN isConst=TRUE;
2110 int s=0;
2111 for(int i=MATCOLS(m)*MATROWS(m)-1;i>=0;i--)
2112 {
2113 poly p=m->m[i];
2114 if (p!=NULL)
2115 {
2116 if(!p_IsConstant(p,r)) isConst=FALSE;
2117 s++;
2118 }
2119 }
2120 if (isConst && rField_is_Q(r)) return DetFactory;
2121 if (s*2<MATCOLS(m)*MATROWS(m)) // few entries
2122 return DetSBareiss;
2123 return DetMu;
2124}
2126{
2127 if (strcmp(s,"Bareiss")==0) return DetBareiss;
2128 if (strcmp(s,"SBareiss")==0) return DetSBareiss;
2129 if (strcmp(s,"Mu")==0) return DetMu;
2130 if (strcmp(s,"Factory")==0) return DetFactory;
2131 WarnS("unknown method for det");
2132 return DetDefault;
2133}
2134
2135
2136poly mp_Det(matrix a, const ring r, DetVariant d/*=DetDefault*/)
2137{
2138 if ((MATCOLS(a)==0)
2139 && (MATROWS(a)==0))
2140 return p_One(r);
2141 if (d==DetDefault) d=mp_GetAlgorithmDet(a,r);
2142 switch (d)
2143 {
2144 case DetBareiss: return mp_DetBareiss(a,r);
2145 case DetMu: return mp_DetMu(a,r);
2146 case DetFactory: return singclap_det(a,r);
2147 case DetSBareiss:
2148 {
2149 ideal I=id_Matrix2Module(mp_Copy(a, r),r);
2150 poly p=sm_CallDet(I, r);
2151 id_Delete(&I, r);
2152 return p;
2153 }
2154 default:
2155 WerrorS("unknown algorithm for det");
2156 return NULL;
2157 }
2158}
2159
2160poly sm_Det(ideal a, const ring r, DetVariant d/*=DetDefault*/)
2161{
2162 if ((MATCOLS(a)==0)
2163 && (MATROWS(a)==0))
2164 return p_One(r);
2165 if (d==DetDefault) d=mp_GetAlgorithmDet((matrix)a,r);
2166 if (d==DetSBareiss) return sm_CallDet(a,r);
2168 poly p=mp_Det(m,r,d);
2169 id_Delete((ideal *)&m,r);
2170 return p;
2171}
All the auxiliary stuff.
static int si_max(const int a, const int b)
Definition: auxiliary.h:124
int BOOLEAN
Definition: auxiliary.h:87
#define TRUE
Definition: auxiliary.h:100
#define FALSE
Definition: auxiliary.h:96
void * ADDRESS
Definition: auxiliary.h:119
static int si_min(const int a, const int b)
Definition: auxiliary.h:125
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:4082
int p
Definition: cfModGcd.cc:4078
CanonicalForm b
Definition: cfModGcd.cc:4103
FILE * f
Definition: checklibs.c:9
poly singclap_det(const matrix m, const ring s)
Definition: clapsing.cc:1757
Variable next() const
Definition: variable.h:52
Definition: intvec.h:23
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:1064
poly * mpRowAdr(int r)
Definition: matpol.cc:926
int mpGetCdim()
Definition: matpol.cc:942
void mpRowReorder()
Definition: matpol.cc:1112
poly mpGetElem(int, int)
Definition: matpol.cc:1229
void mpColReorder()
Definition: matpol.cc:1091
int mpPivotRow(row_col_weight *, int)
void mpSaveArray()
Definition: matpol.cc:945
void mpElimBareiss(poly)
Definition: matpol.cc:1237
poly * Xarray
Definition: matpol.cc:923
void mpDelElem(int, int)
int * qrow
Definition: matpol.cc:922
void mpSetElem(poly, int, int)
int mpPivotBareiss(row_col_weight *)
Definition: matpol.cc:1152
void mpColWeight(float *)
Definition: matpol.cc:1010
void mpInitMat()
Definition: matpol.cc:1078
void mpSetSearch(int s)
int * qcol
Definition: matpol.cc:922
int mpGetRdim()
Definition: matpol.cc:941
poly * mpColAdr(int c)
Definition: matpol.cc:928
void mpRowSwap(int, int)
Definition: matpol.cc:1049
int mpGetSign()
Definition: matpol.cc:943
void mpToIntvec(intvec *)
~mp_permmatrix()
Definition: matpol.cc:991
void mpRowWeight(float *)
Definition: matpol.cc:1029
float * wcol
Definition: matpol.cc:886
float * wrow
Definition: matpol.cc:886
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:567
#define Print
Definition: emacs.cc:80
#define WarnS
Definition: emacs.cc:78
return result
Definition: facAbsBiFact.cc:75
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
for(j=0;j< factors.length();j++)
Definition: facHensel.cc:129
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:809
matrix mp_Wedge(matrix a, int ar, const ring R)
Definition: matpol.cc:1744
static poly mp_SelectId(ideal I, poly what, const ring R)
Definition: matpol.cc:759
static void mp_PartClean(matrix a, int lr, int lc, const ring R)
Definition: matpol.cc:797
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:1824
static ideal sm_MultAndShift(poly f, ideal B, int s, const ring r)
Definition: matpol.cc:1796
ideal sm_UnFlatten(ideal a, int col, const ring R)
Definition: matpol.cc:1939
int sm_Compare(ideal a, ideal b, const ring R)
Definition: matpol.cc:1977
ideal sm_Add(ideal a, ideal b, const ring R)
Definition: matpol.cc:1864
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:1966
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:873
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:2025
DetVariant mp_GetAlgorithmDet(matrix m, const ring r)
Definition: matpol.cc:2105
static float mp_PolyWeight(poly p, const ring r)
Definition: matpol.cc:1295
static int mp_PivRow(matrix a, int lr, int lc, const ring r)
Definition: matpol.cc:1387
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:1412
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:2160
static poly p_Insert(poly p1, poly p2, const ring R)
Definition: matpol.cc:683
static int mp_PrepareRow(matrix a, int lr, int lc, const ring R)
Definition: matpol.cc:1374
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:1874
ideal sm_Mult(ideal a, ideal b, const ring R)
Definition: matpol.cc:1884
static void mpFinalClean(matrix a)
Definition: matpol.cc:1588
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:2136
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:636
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:1996
static void mp_ElimBar(matrix a0, matrix re, poly div, int lr, int lc, const ring R)
Definition: matpol.cc:1442
static void mpReplace(int j, int n, int &sign, int *perm)
Definition: matpol.cc:1137
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:655
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:848
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:741
static int mp_PreparePiv(matrix a, int lr, int lc, const ring r)
Definition: matpol.cc:1432
void iiWriteMatrix(matrix im, const char *n, int dim, const ring r, int spaces)
set spaces to zero by default
Definition: matpol.cc:827
poly mp_DetMu(matrix A, const ring R)
Definition: matpol.cc:2063
static void sm_AddSubMat(ideal res, ideal sm, int col, const ring r)
Definition: matpol.cc:1815
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:1500
ideal sm_Flatten(ideal a, const ring R)
Definition: matpol.cc:1919
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:1596
static int mp_PivBar(matrix a, int lr, int lc, const ring r)
Definition: matpol.cc:1327
static void mpSwapRow(matrix a, int pos, int lr, int lc)
Definition: matpol.cc:1354
poly mp_DetBareiss(matrix a, const ring r)
returns the determinant of the matrix m; uses Bareiss algorithm
Definition: matpol.cc:1669
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)
Definition: omAllocDecl.h:260
#define omAlloc(size)
Definition: omAllocDecl.h:210
#define omAllocBin(bin)
Definition: omAllocDecl.h:205
#define omAlloc0(size)
Definition: omAllocDecl.h:211
#define omFreeBin(addr, bin)
Definition: omAllocDecl.h:259
#define omfreeSize(addr, size)
Definition: omAllocDecl.h:236
#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:1297
void p_Vec2Polys(poly v, poly **p, int *len, const ring r)
Definition: p_polys.cc:3625
void p_Shift(poly *p, int i, const ring r)
shifts components of the vector p by i
Definition: p_polys.cc:4706
void p_Normalize(poly p, const ring r)
Definition: p_polys.cc:3813
int p_Compare(const poly a, const poly b, const ring R)
Definition: p_polys.cc:4849
poly p_One(const ring r)
Definition: p_polys.cc:1313
poly p_Sub(poly p1, poly p2, const ring r)
Definition: p_polys.cc:1990
void pEnlargeSet(poly **p, int l, int increment)
Definition: p_polys.cc:3696
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:3595
BOOLEAN p_EqualPolys(poly p1, poly p2, const ring r)
Definition: p_polys.cc:4512
static poly p_Neg(poly p, const ring r)
Definition: p_polys.h:1105
static int pLength(poly a)
Definition: p_polys.h:188
static poly p_Add_q(poly p, poly q, const ring r)
Definition: p_polys.h:934
static void p_LmDelete(poly p, const ring r)
Definition: p_polys.h:721
static poly p_Mult_q(poly p, poly q, const ring r)
Definition: p_polys.h:1112
static int p_Cmp(poly p1, poly p2, ring r)
Definition: p_polys.h:1725
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:486
static unsigned long p_SetComp(poly p, unsigned long c, ring r)
Definition: p_polys.h:245
static void p_Setm(poly p, const ring r)
Definition: p_polys.h:231
#define p_SetmComp
Definition: p_polys.h:242
static poly p_Head(const poly p, const ring r)
copy the (leading) term of p
Definition: p_polys.h:858
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:467
static BOOLEAN p_IsConstant(const poly p, const ring r)
Definition: p_polys.h:1962
static void p_Delete(poly *p, const ring r)
Definition: p_polys.h:899
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:1149
static BOOLEAN p_IsUnit(const poly p, const ring r)
Definition: p_polys.h:1989
static poly p_Copy(poly p, const ring r)
returns a copy of p
Definition: p_polys.h:844
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:3427
static BOOLEAN rField_is_Zp(const ring r)
Definition: ring.h:500
static BOOLEAN rField_is_Q(const ring r)
Definition: ring.h:506
static short rVar(const ring r)
#define rVar(r) (r->N)
Definition: ring.h:592
BOOLEAN rHasLocalOrMixedOrdering(const ring r)
Definition: ring.h:760
int status int void size_t count
Definition: si_signals.h:59
ideal idInit(int idsize, int rank)
initialise an ideal / module
Definition: simpleideals.cc:35
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
Definition: simpleideals.cc:27
#define IDELEMS(i)
Definition: simpleideals.h:23
#define id_Test(A, lR)
Definition: simpleideals.h:87
#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:75
int dim(ideal I, ring r)