My Project
Loading...
Searching...
No Matches
matpol.cc File Reference
#include "misc/auxiliary.h"
#include "misc/mylimits.h"
#include "misc/intvec.h"
#include "coeffs/numbers.h"
#include "reporter/reporter.h"
#include "monomials/ring.h"
#include "monomials/p_polys.h"
#include "simpleideals.h"
#include "matpol.h"
#include "prCopy.h"
#include "clapsing.h"
#include "sparsmat.h"

Go to the source code of this file.

Data Structures

class  row_col_weight
 
class  mp_permmatrix
 

Functions

static poly mp_Exdiv (poly m, poly d, poly vars, const ring)
 
static poly mp_Select (poly fro, poly what, const ring)
 
static poly mp_SelectId (ideal I, poly what, const ring R)
 
matrix mpNew (int r, int c)
 create a r x c zero-matrix
 
matrix mp_Copy (matrix a, const ring r)
 copies matrix a (from ring r to r)
 
matrix mp_Copy (const matrix a, const ring rSrc, const ring rDst)
 copies matrix a from rSrc into rDst
 
matrix mp_InitP (int r, int c, poly p, const ring R)
 make it a p * unit matrix
 
matrix mp_InitI (int r, int c, int v, const ring R)
 make it a v * unit matrix
 
matrix mp_MultI (matrix a, long f, const ring R)
 c = f*a
 
matrix mp_MultP (matrix a, poly p, const ring R)
 multiply a matrix 'a' by a poly 'p', destroy the args
 
matrix pMultMp (poly p, matrix a, const ring R)
 
matrix mp_Add (matrix a, matrix b, const ring R)
 
matrix mp_Sub (matrix a, matrix b, const ring R)
 
matrix mp_Mult (matrix a, matrix b, const ring R)
 
matrix mp_Transp (matrix a, const ring R)
 
poly mp_Trace (matrix a, const ring R)
 
poly TraceOfProd (matrix a, matrix b, int n, const ring R)
 
matrix mp_Coeffs (ideal I, int var, const ring R)
 corresponds to Maple's coeffs: var has to be the number of a variable
 
void mp_Monomials (matrix c, int r, int var, matrix m, const ring R)
 
matrix mp_CoeffProc (poly f, poly vars, const ring R)
 
matrix mp_CoeffProcId (ideal I, poly vars, const ring R)
 
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, eg 'xy'; then the poly f is searched for monomials in x and y, these monomials are written to the first row of the matrix co. the second row of co contains the respective factors in f. Thus f = sum co[1,i]*co[2,i], i = 1..cols, rows equals 2.
 
int mp_Compare (matrix a, matrix b, const ring R)
 
BOOLEAN mp_Equal (matrix a, matrix b, const ring R)
 
static poly p_Insert (poly p1, poly p2, const ring R)
 
static void mp_PartClean (matrix a, int lr, int lc, const ring R)
 
BOOLEAN mp_IsDiagUnit (matrix U, const ring R)
 
void iiWriteMatrix (matrix im, const char *n, int dim, const ring r, int spaces)
 set spaces to zero by default
 
char * iiStringMatrix (matrix im, int dim, const ring r, char ch)
 
void mp_Delete (matrix *a, const ring r)
 
static float mp_PolyWeight (poly p, const ring r)
 
static void mpReplace (int j, int n, int &sign, int *perm)
 
static int mp_PivBar (matrix a, int lr, int lc, const ring r)
 
static void mpSwapRow (matrix a, int pos, int lr, int lc)
 
static int mp_PrepareRow (matrix a, int lr, int lc, const ring R)
 
static int mp_PivRow (matrix a, int lr, int lc, const ring r)
 
static void mpSwapCol (matrix a, int pos, int lr, int lc)
 
static int mp_PreparePiv (matrix a, int lr, int lc, const ring r)
 
static void mp_ElimBar (matrix a0, matrix re, poly div, int lr, int lc, const ring R)
 
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)
 
static void mpFinalClean (matrix a)
 
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
 
poly mp_DetBareiss (matrix a, const ring r)
 returns the determinant of the matrix m; uses Bareiss algorithm
 
matrix mp_Wedge (matrix a, int ar, const ring R)
 
static ideal sm_MultAndShift (poly f, ideal B, int s, const ring r)
 
static void sm_AddSubMat (ideal res, ideal sm, int col, const ring r)
 
ideal sm_Tensor (ideal A, ideal B, const ring r)
 
ideal sm_Add (ideal a, ideal b, const ring R)
 
ideal sm_Sub (ideal a, ideal b, const ring R)
 
ideal sm_Mult (ideal a, ideal b, const ring R)
 
ideal sm_Flatten (ideal a, const ring R)
 
ideal sm_UnFlatten (ideal a, int col, const ring R)
 
poly sm_Trace (ideal a, const ring R)
 
int sm_Compare (ideal a, ideal b, const ring R)
 
BOOLEAN sm_Equal (ideal a, ideal b, const ring R)
 
static matrix mu (matrix A, const ring R)
 
poly mp_DetMu (matrix A, const ring R)
 
DetVariant mp_GetAlgorithmDet (matrix m, const ring r)
 
DetVariant mp_GetAlgorithmDet (const char *s)
 
poly mp_Det (matrix a, const ring r, DetVariant d)
 
poly sm_Det (ideal a, const ring r, DetVariant d)
 

Function Documentation

◆ iiStringMatrix()

char * iiStringMatrix ( matrix im,
int dim,
const ring r,
char ch )

Definition at line 849 of file matpol.cc.

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}
CanonicalForm FACTORY_PUBLIC pp(const CanonicalForm &)
CanonicalForm pp ( const CanonicalForm & f )
Definition cf_gcd.cc:676
int i
Definition cfEzgcd.cc:132
poly * m
Definition matpol.h:18
const CanonicalForm int s
Definition facAbsFact.cc:51
int j
Definition facHensel.cc:110
#define MATROWS(i)
Definition matpol.h:26
#define MATCOLS(i)
Definition matpol.h:27
void p_String0(poly p, ring lmRing, ring tailRing)
print p according to ShortOut in lmRing & tailRing
Definition polys0.cc:223
void StringSetS(const char *st)
Definition reporter.cc:128
void StringAppendS(const char *st)
Definition reporter.cc:107
char * StringEndS()
Definition reporter.cc:151
int dim(ideal I, ring r)

◆ iiWriteMatrix()

void iiWriteMatrix ( matrix im,
const char * n,
int dim,
const ring r,
int spaces )

set spaces to zero by default

Definition at line 828 of file matpol.cc.

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}
#define Print
Definition emacs.cc:80
void p_Write(poly p, ring lmRing, ring tailRing)
Definition polys0.cc:342
void p_Write0(poly p, ring lmRing, ring tailRing)
Definition polys0.cc:332

◆ mp_Add()

matrix mp_Add ( matrix a,
matrix b,
const ring R )

Definition at line 172 of file matpol.cc.

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}
int m
Definition cfEzgcd.cc:128
int k
Definition cfEzgcd.cc:99
CanonicalForm b
Definition cfModGcd.cc:4111
int nrows
Definition matpol.h:20
int ncols
Definition matpol.h:21
matrix mpNew(int r, int c)
create a r x c zero-matrix
Definition matpol.cc:37
ip_smatrix * matrix
Definition matpol.h:43
#define NULL
Definition omList.c:12
static poly p_Add_q(poly p, poly q, const ring r)
Definition p_polys.h:938
static poly p_Copy(poly p, const ring r)
returns a copy of p
Definition p_polys.h:848
#define R
Definition sirandom.c:27

◆ mp_Coef2()

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, eg 'xy'; then the poly f is searched for monomials in x and y, these monomials are written to the first row of the matrix co. the second row of co contains the respective factors in f. Thus f = sum co[1,i]*co[2,i], i = 1..cols, rows equals 2.

Definition at line 574 of file matpol.cc.

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}
static int si_max(const int a, const int b)
Definition auxiliary.h:125
int l
Definition cfEzgcd.cc:100
int p
Definition cfModGcd.cc:4086
const Variable & v
< [in] a sqrfree bivariate poly
Definition facBivar.h:39
STATIC_VAR Poly * h
Definition janet.cc:971
static poly mp_Exdiv(poly m, poly d, poly vars, const ring)
Definition matpol.cc:554
static poly mp_Select(poly fro, poly what, const ring)
Definition matpol.cc:742
#define MATELEM(mat, i, j)
1-based access to matrix
Definition matpol.h:29
#define __p_GetComp(p, r)
Definition monomials.h:63
#define omFree(addr)
void p_Vec2Polys(poly v, poly **p, int *len, const ring r)
Definition p_polys.cc:3705
static int pLength(poly a)
Definition p_polys.h:190
static unsigned long p_SetComp(poly p, unsigned long c, ring r)
Definition p_polys.h:249
static poly p_Head(const poly p, const ring r)
copy the (leading) term of p
Definition p_polys.h:862
static BOOLEAN p_IsConstant(const poly p, const ring r)
Definition p_polys.h:1980
#define loop
Definition structs.h:71

◆ mp_CoeffProc()

matrix mp_CoeffProc ( poly f,
poly vars,
const ring R )

Definition at line 392 of file matpol.cc.

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}
FILE * f
Definition checklibs.c:9
#define assume(x)
Definition mod2.h:389
#define pIter(p)
Definition monomials.h:37
#define pNext(p)
Definition monomials.h:36
poly p_One(const ring r)
Definition p_polys.cc:1314
static BOOLEAN rHasLocalOrMixedOrdering(const ring r)
Definition ring.h:769

◆ mp_CoeffProcId()

matrix mp_CoeffProcId ( ideal I,
poly vars,
const ring R )

Definition at line 469 of file matpol.cc.

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}
BOOLEAN idIs0(ideal h)
returns true if h is the zero ideal
static poly mp_SelectId(ideal I, poly what, const ring R)
Definition matpol.cc:760
#define IDELEMS(i)

◆ mp_Coeffs()

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 at line 306 of file matpol.cc.

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}
#define p_GetComp(p, r)
Definition monomials.h:64
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 void p_Setm(poly p, const ring r)
Definition p_polys.h:235
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
void id_Delete(ideal *h, ring r)
deletes an ideal/module/matrix

◆ mp_Compare()

int mp_Compare ( matrix a,
matrix b,
const ring R )

Definition at line 637 of file matpol.cc.

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}
int p_Compare(const poly a, const poly b, const ring R)
Definition p_polys.cc:5005

◆ mp_Copy() [1/2]

matrix mp_Copy ( const matrix a,
const ring rSrc,
const ring rDst )

copies matrix a from rSrc into rDst

Definition at line 78 of file matpol.cc.

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}
long rank
Definition matpol.h:19
void p_Normalize(poly p, const ring r)
Definition p_polys.cc:3894
poly prCopyR_NoSort(poly p, ring src_r, ring dest_r)
Definition prCopy.cc:77
#define id_Test(A, lR)

◆ mp_Copy() [2/2]

matrix mp_Copy ( matrix a,
const ring r )

copies matrix a (from ring r to r)

Definition at line 57 of file matpol.cc.

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}

◆ mp_Delete()

void mp_Delete ( matrix * a,
const ring r )

Definition at line 874 of file matpol.cc.

875{
876 id_Delete((ideal *) a, r);
877}

◆ mp_Det()

poly mp_Det ( matrix a,
const ring r,
DetVariant d )

Definition at line 2139 of file matpol.cc.

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}
poly singclap_det(const matrix m, const ring s)
Definition clapsing.cc:1760
void WerrorS(const char *s)
Definition feFopen.cc:24
DetVariant mp_GetAlgorithmDet(matrix m, const ring r)
Definition matpol.cc:2108
matrix mp_Copy(matrix a, const ring r)
copies matrix a (from ring r to r)
Definition matpol.cc:57
poly mp_DetMu(matrix A, const ring R)
Definition matpol.cc:2066
poly mp_DetBareiss(matrix a, const ring r)
returns the determinant of the matrix m; uses Bareiss algorithm
Definition matpol.cc:1670
@ 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
ideal id_Matrix2Module(matrix mat, const ring R)
converts mat to module, destroys mat
poly sm_CallDet(ideal I, const ring R)
Definition sparsmat.cc:302

◆ mp_DetBareiss()

poly mp_DetBareiss ( matrix a,
const ring r )

returns the determinant of the matrix m; uses Bareiss algorithm

Definition at line 1670 of file matpol.cc.

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}
CF_NO_INLINE FACTORY_PUBLIC CanonicalForm div(const CanonicalForm &, const CanonicalForm &)
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
void mpSaveArray()
Definition matpol.cc:946
void mpElimBareiss(poly)
Definition matpol.cc:1238
int mpPivotBareiss(row_col_weight *)
Definition matpol.cc:1153
int mpGetRdim()
Definition matpol.cc:942
int mpGetSign()
Definition matpol.cc:944
CanonicalForm res
Definition facAbsFact.cc:60
const CanonicalForm & w
Definition facAbsFact.cc:51
static poly p_Neg(poly p, const ring r)
Definition p_polys.h:1109
void Werror(const char *fmt,...)
Definition reporter.cc:189

◆ mp_DetMu()

poly mp_DetMu ( matrix A,
const ring R )

Definition at line 2066 of file matpol.cc.

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}
static matrix mu(matrix A, const ring R)
Definition matpol.cc:2028
matrix mp_Mult(matrix a, matrix b, const ring R)
Definition matpol.cc:206
#define MATELEM0(mat, i, j)
0-based access to matrix
Definition matpol.h:31
#define A
Definition sirandom.c:24

◆ mp_ElimBar()

static void mp_ElimBar ( matrix a0,
matrix re,
poly div,
int lr,
int lc,
const ring R )
static

Definition at line 1443 of file matpol.cc.

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}
CanonicalForm lc(const CanonicalForm &f)
Variable x
Definition cfModGcd.cc:4090
Definition ap.h:40
poly sm_MultDiv(poly a, poly b, const poly c, const ring R)
Definition sparsmat.cc:1759
void sm_SpecialPolyDiv(poly a, poly b, const ring R)
Definition sparsmat.cc:1840

◆ mp_Equal()

BOOLEAN mp_Equal ( matrix a,
matrix b,
const ring R )

Definition at line 656 of file matpol.cc.

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}
#define TRUE
Definition auxiliary.h:101
#define FALSE
Definition auxiliary.h:97
BOOLEAN p_EqualPolys(poly p1, poly p2, const ring r)
Definition p_polys.cc:4621
static int p_Cmp(poly p1, poly p2, ring r)
Definition p_polys.h:1743

◆ mp_Exdiv()

static poly mp_Exdiv ( poly m,
poly d,
poly vars,
const ring R )
static

Definition at line 554 of file matpol.cc.

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}
static void p_Delete(poly *p, const ring r)
Definition p_polys.h:903
static short rVar(const ring r)
#define rVar(r) (r->N)
Definition ring.h:598

◆ mp_GetAlgorithmDet() [1/2]

DetVariant mp_GetAlgorithmDet ( const char * s)

Definition at line 2128 of file matpol.cc.

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}
#define WarnS
Definition emacs.cc:78

◆ mp_GetAlgorithmDet() [2/2]

DetVariant mp_GetAlgorithmDet ( matrix m,
const ring r )

Definition at line 2108 of file matpol.cc.

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}
int BOOLEAN
Definition auxiliary.h:88
static BOOLEAN rField_is_Zp(const ring r)
Definition ring.h:506
static BOOLEAN rField_is_Q(const ring r)
Definition ring.h:512

◆ mp_InitI()

matrix mp_InitI ( int r,
int c,
int v,
const ring R )

make it a v * unit matrix

Definition at line 122 of file matpol.cc.

123{
124 return mp_InitP(r, c, p_ISet(v, R), R);
125}
matrix mp_InitP(int r, int c, poly p, const ring R)
make it a p * unit matrix
Definition matpol.cc:106
poly p_ISet(long i, const ring r)
returns the poly representing the integer i
Definition p_polys.cc:1298

◆ mp_InitP()

matrix mp_InitP ( int r,
int c,
poly p,
const ring R )

make it a p * unit matrix

Definition at line 106 of file matpol.cc.

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}
static int si_min(const int a, const int b)
Definition auxiliary.h:126

◆ mp_IsDiagUnit()

BOOLEAN mp_IsDiagUnit ( matrix U,
const ring R )

Definition at line 810 of file matpol.cc.

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}
static BOOLEAN p_IsUnit(const poly p, const ring r)
Definition p_polys.h:2007

◆ mp_MinorToResult()

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 at line 1501 of file matpol.cc.

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}
return result
void pEnlargeSet(poly **p, int l, int increment)
Definition p_polys.cc:3776

◆ mp_Monomials()

void mp_Monomials ( matrix c,
int r,
int var,
matrix m,
const ring R )

Definition at line 355 of file matpol.cc.

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}
void * ADDRESS
Definition auxiliary.h:120
#define omAlloc0(size)
#define omfreeSize(addr, size)

◆ mp_Mult()

matrix mp_Mult ( matrix a,
matrix b,
const ring R )

Definition at line 206 of file matpol.cc.

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}
static poly pp_Mult_qq(poly p, poly q, const ring r)
Definition p_polys.h:1162

◆ mp_MultI()

matrix mp_MultI ( matrix a,
long f,
const ring R )

c = f*a

Definition at line 128 of file matpol.cc.

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}
static poly p_Mult_q(poly p, poly q, const ring r)
Definition p_polys.h:1120

◆ mp_MultP()

matrix mp_MultP ( matrix a,
poly p,
const ring R )

multiply a matrix 'a' by a poly 'p', destroy the args

Definition at line 141 of file matpol.cc.

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}

◆ mp_PartClean()

static void mp_PartClean ( matrix a,
int lr,
int lc,
const ring R )
static

Definition at line 798 of file matpol.cc.

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}

◆ mp_PivBar()

static int mp_PivBar ( matrix a,
int lr,
int lc,
const ring r )
static

Definition at line 1328 of file matpol.cc.

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}
static float mp_PolyWeight(poly p, const ring r)
Definition matpol.cc:1296

◆ mp_PivRow()

static int mp_PivRow ( matrix a,
int lr,
int lc,
const ring r )
static

Definition at line 1388 of file matpol.cc.

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}

◆ mp_PolyWeight()

static float mp_PolyWeight ( poly p,
const ring r )
static

Definition at line 1296 of file matpol.cc.

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}
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
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

◆ mp_PreparePiv()

static int mp_PreparePiv ( matrix a,
int lr,
int lc,
const ring r )
static

Definition at line 1433 of file matpol.cc.

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}
static int mp_PivRow(matrix a, int lr, int lc, const ring r)
Definition matpol.cc:1388
static void mpSwapCol(matrix a, int pos, int lr, int lc)
Definition matpol.cc:1413

◆ mp_PrepareRow()

static int mp_PrepareRow ( matrix a,
int lr,
int lc,
const ring R )
static

Definition at line 1375 of file matpol.cc.

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}
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

◆ mp_RecMin()

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

for minors with Bareiss

Definition at line 1597 of file matpol.cc.

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}
static void mp_PartClean(matrix a, int lr, int lc, const ring R)
Definition matpol.cc:798
static int mp_PrepareRow(matrix a, int lr, int lc, const ring R)
Definition matpol.cc:1375
static void mpFinalClean(matrix a)
Definition matpol.cc:1589
static void mp_ElimBar(matrix a0, matrix re, poly div, int lr, int lc, const ring R)
Definition matpol.cc:1443
static int mp_PreparePiv(matrix a, int lr, int lc, const ring r)
Definition matpol.cc:1433
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
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

◆ mp_Select()

static poly mp_Select ( poly fro,
poly what,
const ring R )
static

Definition at line 742 of file matpol.cc.

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}
static poly p_Insert(poly p1, poly p2, const ring R)
Definition matpol.cc:684

◆ mp_SelectId()

static poly mp_SelectId ( ideal I,
poly what,
const ring R )
static

Definition at line 760 of file matpol.cc.

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}

◆ mp_Sub()

matrix mp_Sub ( matrix a,
matrix b,
const ring R )

Definition at line 189 of file matpol.cc.

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}
poly p_Sub(poly p1, poly p2, const ring r)
Definition p_polys.cc:1994

◆ mp_Trace()

poly mp_Trace ( matrix a,
const ring R )

Definition at line 268 of file matpol.cc.

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}

◆ mp_Transp()

matrix mp_Transp ( matrix a,
const ring R )

Definition at line 247 of file matpol.cc.

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}

◆ mp_Wedge()

matrix mp_Wedge ( matrix a,
int ar,
const ring R )

Definition at line 1745 of file matpol.cc.

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}
int binom(int n, int r)
void idGetNextChoise(int r, int end, BOOLEAN *endch, int *choise)
void idInitChoise(int r, int beg, int end, BOOLEAN *endch, int *choise)
#define omAlloc(size)

◆ mpFinalClean()

static void mpFinalClean ( matrix a)
static

Definition at line 1589 of file matpol.cc.

1590{
1591 omFreeSize((ADDRESS)a->m,a->nrows*a->ncols*sizeof(poly));
1593}
#define omFreeSize(addr, size)
#define omFreeBin(addr, bin)
VAR omBin sip_sideal_bin

◆ mpNew()

matrix mpNew ( int r,
int c )

create a r x c zero-matrix

Definition at line 37 of file matpol.cc.

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}
#define omAllocBin(bin)

◆ mpReplace()

static void mpReplace ( int j,
int n,
int & sign,
int * perm )
static

Definition at line 1138 of file matpol.cc.

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}
static int sign(int x)
Definition ring.cc:3503

◆ mpSwapCol()

static void mpSwapCol ( matrix a,
int pos,
int lr,
int lc )
static

Definition at line 1413 of file matpol.cc.

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}

◆ mpSwapRow()

static void mpSwapRow ( matrix a,
int pos,
int lr,
int lc )
static

Definition at line 1355 of file matpol.cc.

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}

◆ mu()

static matrix mu ( matrix A,
const ring R )
static

Definition at line 2028 of file matpol.cc.

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}

◆ p_Insert()

static poly p_Insert ( poly p1,
poly p2,
const ring R )
static

Definition at line 684 of file matpol.cc.

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}
static void p_LmDelete(poly p, const ring r)
Definition p_polys.h:725

◆ pMultMp()

matrix pMultMp ( poly p,
matrix a,
const ring R )

Definition at line 158 of file matpol.cc.

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}

◆ sm_Add()

ideal sm_Add ( ideal a,
ideal b,
const ring R )

Definition at line 1867 of file matpol.cc.

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}
ideal idInit(int idsize, int rank)
initialise an ideal / module

◆ sm_AddSubMat()

static void sm_AddSubMat ( ideal res,
ideal sm,
int col,
const ring r )
static

Definition at line 1818 of file matpol.cc.

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}

◆ sm_Compare()

int sm_Compare ( ideal a,
ideal b,
const ring R )

Definition at line 1980 of file matpol.cc.

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}

◆ sm_Det()

poly sm_Det ( ideal a,
const ring r,
DetVariant d )

Definition at line 2163 of file matpol.cc.

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}
ideal id_Copy(ideal h1, const ring r)
copy an ideal
poly mp_Det(matrix a, const ring r, DetVariant d)
Definition matpol.cc:2139
matrix id_Module2Matrix(ideal mod, const ring R)

◆ sm_Equal()

BOOLEAN sm_Equal ( ideal a,
ideal b,
const ring R )

Definition at line 1999 of file matpol.cc.

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}

◆ sm_Flatten()

ideal sm_Flatten ( ideal a,
const ring R )

Definition at line 1922 of file matpol.cc.

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}
void p_Shift(poly *p, int i, const ring r)
shifts components of the vector p by i
Definition p_polys.cc:4815

◆ sm_Mult()

ideal sm_Mult ( ideal a,
ideal b,
const ring R )

Definition at line 1887 of file matpol.cc.

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}
#define SMATELEM(A, i, j, R)
Definition matpol.h:123

◆ sm_MultAndShift()

static ideal sm_MultAndShift ( poly f,
ideal B,
int s,
const ring r )
static

Definition at line 1799 of file matpol.cc.

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}
b *CanonicalForm B
Definition facBivar.cc:52

◆ sm_Sub()

ideal sm_Sub ( ideal a,
ideal b,
const ring R )

Definition at line 1877 of file matpol.cc.

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}

◆ sm_Tensor()

ideal sm_Tensor ( ideal A,
ideal B,
const ring r )

Definition at line 1827 of file matpol.cc.

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}
static ideal sm_MultAndShift(poly f, ideal B, int s, const ring r)
Definition matpol.cc:1799
static void sm_AddSubMat(ideal res, ideal sm, int col, const ring r)
Definition matpol.cc:1818
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

◆ sm_Trace()

poly sm_Trace ( ideal a,
const ring R )

Definition at line 1969 of file matpol.cc.

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}

◆ sm_UnFlatten()

ideal sm_UnFlatten ( ideal a,
int col,
const ring R )

Definition at line 1942 of file matpol.cc.

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}
int comp(const CanonicalForm &A, const CanonicalForm &B)
compare polynomials
#define p_SetmComp
Definition p_polys.h:246

◆ TraceOfProd()

poly TraceOfProd ( matrix a,
matrix b,
int n,
const ring R )

Definition at line 282 of file matpol.cc.

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}