Changeset d2da42 in git

Timestamp:

Dec 6, 2018, 3:31:44 PM (5 years ago)

Author:

Hans Schoenemann <hannes@…>

Branches:

(u'spielwiese', 'fe61d9c35bf7c61f2b6cbf1b56e25e2f08d536cc')

Children:

757eac36df59510348bf034d6ebdc03389e843f8

Parents:

7d0c5307af92569d0d63f695cef1d0431f7d1b9f

Message:

add: mp_DetMu: determinant via mu algorithm

Location:

libpolys/polys

Files:

• matpol.cc (modified) (2 diffs)
• matpol.h (modified) (2 diffs)

libpolys/polys/matpol.cc

@@ -241 +241 @@
             poly *cij=&(MATELEM(c,i,j));
             poly s = pp_Mult_qq(aik /*MATELEM(a,i,k)*/, bkj/*MATELEM(b,k,j)*/, R);
-            if (/*MATELEM(c,i,j)*/ (*cij)==NULL) (*cij)=s;
-            else (*cij) = p_Add_q((*cij) /*MATELEM(c,i,j)*/ ,s, R);
+            (*cij)/*MATELEM(c,i,j)*/ = p_Add_q((*cij) /*MATELEM(c,i,j)*/ ,s, R);
           }
         }
@@ -2027 +2026 @@
 }
 
+/*
+*   mu-Algorithmus:
+*/
+
+//  mu-Matrix
+static void mu(matrix A, matrix &X, const ring R)
+{
+  int n=MATROWS(A);
+  assume(MATCOLS(A)==n);
+    /*  Die Funktion erstellt die Matrix mu
+    *
+    *   Input:
+    *   int n: Dimension der Matrix
+    *   int A: Matrix der Größe n*n
+    *   int X: Speicherplatz fÃŒr Output
+    *
+    *   In der Matrix X speichert man die Matrix mu
+    */
+
+    // X als n*n Null-Matrix initalisieren
+    X=mpNew(n,n);
+
+    //  DiagonaleintrÀge von X berrechnen
+    poly sum = NULL;
+    for (int i = n-1; i >= 0; i--)
+    {
+        MATELEM0(X,i,i) = p_Copy(sum,R);
+        sum=p_Sub(sum,p_Copy(MATELEM0(A,i,i),R),R);
+    }
+
+    //  EintrÀge aus dem oberen Dreieck von A nach X Ìbertragen
+    for (int i = 0; i < n; i++)
+    {
+        for (int j = i+1; j < n; j++)
+        {
+            MATELEM0(X,i,j)=p_Copy(MATELEM0(A,i,j),R);
+        }
+    }
+}
+
+// Funktion muDet
+poly mp_DetMu(matrix A, const ring R)
+{
+  int n=MATROWS(A);
+  assume(MATCOLS(A)==n);
+    /*
+    *   Intput:
+    *   int n: Dimension der Matrix
+    *   int A: n*n Matrix
+    *
+    *   Berechnet n-1 mal: X = mu(X)*A
+    *
+    *   Output: det(A)
+    */
+
+    //speichere A ab:
+    matrix B=mp_Copy(A,R);
+    A=mp_Copy(A,R);
+
+    // berechen X = mu(X)*A
+    matrix X;
+    for (int i = 0; i < n-1; i++)
+    {
+        mu(A,X,R);
+        id_Delete((ideal*)&A,R);
+        A=mp_Mult(X,B,R);
+        id_Delete((ideal*)&X,R);
+    }
+
+    // berrechne det(A)
+    poly res;
+    if (n%2 == 0)
+    {
+        res=p_Neg(MATELEM0(A,0,0),R);
+    }
+    else
+    {
+        res=MATELEM0(A,0,0);
+    }
+    MATELEM0(A,0,0)=NULL;
+    id_Delete((ideal*)&A,R);
+    return res;
+}
+

libpolys/polys/matpol.h

@@ -26 +26 @@
   #define MATROWS(i) ((i)->nrows)
   #define MATCOLS(i) ((i)->ncols)
+  /// 1-based access to matrix
   #define MATELEM(mat,i,j) ((mat)->m)[MATCOLS((mat)) * ((i)-1) + (j)-1]
+  /// 0-based access to matrix
+  #define MATELEM0(mat,i,j) ((mat)->m)[MATCOLS((mat)) * (i) + (j)]
 };
 
@@ -57 +60 @@
 
 poly mp_DetBareiss (matrix a, const ring r);
+poly mp_DetMu(matrix A, const ring R);
+
 
 //matrix mp_Homogen(matrix a, int v, const ring r);