Changeset 178bd3 in git for libpolys

Timestamp:

Dec 17, 2018, 7:56:25 AM (5 years ago)

Author:

Karim Abou Zeid <karim23697@…>

Branches:

(u'spielwiese', 'fe61d9c35bf7c61f2b6cbf1b56e25e2f08d536cc')

Children:

c7231560b85e2e1683beb9629db0560ee8ca4255

Parents:

04d263ab0fcd4cc5a8215c12fdf4423b782d3283ee3e7cd650e2cf6d6e770d6b57b1681594c7cda7

Message:

Merge branch 'spielwiese' into christmas-release

Location:

libpolys/polys

Files:

• matpol.cc (modified) (7 diffs)
• matpol.h (modified) (3 diffs)
• shiftop.cc (modified) (3 diffs)
• sparsmat.cc (modified) (1 diff)
• sparsmat.h (modified) (1 diff)

libpolys/polys/matpol.cc

@@ -23 +23 @@
 #include "matpol.h"
 #include "prCopy.h"
+#include "clapsing.h"
 
 #include "sparsmat.h"
@@ -2029 +2030 @@
 
 //  mu-Matrix
-static void mu(matrix A, matrix &X, const ring R)
+static matrix mu(matrix A, const ring R)
 {
   int n=MATROWS(A);
@@ -2044 +2045 @@
 
     // X als n*n Null-Matrix initalisieren
-    X=mpNew(n,n);
+    matrix X=mpNew(n,n);
 
     //  Diagonaleintraege von X berrechnen
@@ -2053 +2054 @@
         sum=p_Sub(sum,p_Copy(MATELEM0(A,i,i),R),R);
     }
+    p_Delete(&sum,R);
 
     //  Eintraege aus dem oberen Dreieck von A nach X uebertragen
-    for (int i = 0; i < n; i++)
+    for (int i = n-1; i >=0; i--)
     {
         for (int j = i+1; j < n; j++)
@@ -2062 +2064 @@
         }
     }
+    return X;
 }
 
@@ -2080 +2083 @@
 
     //speichere A ab:
-    matrix B=mp_Copy(A,R);
-    A=mp_Copy(A,R);
+    matrix workA=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);
+    for (int i = n-1; i >0; i--)
+    {
+        X=mu(workA,R);
+        id_Delete((ideal*)&workA,R);
+        workA=mp_Mult(X,A,R);
         id_Delete((ideal*)&X,R);
     }
@@ -2097 +2099 @@
     if (n%2 == 0)
     {
-        res=p_Neg(MATELEM0(A,0,0),R);
+        res=p_Neg(MATELEM0(workA,0,0),R);
     }
     else
     {
-        res=MATELEM0(A,0,0);
-    }
-    MATELEM0(A,0,0)=NULL;
-    id_Delete((ideal*)&A,R);
+        res=MATELEM0(workA,0,0);
+    }
+    MATELEM0(workA,0,0)=NULL;
+    id_Delete((ideal*)&workA,R);
     return res;
 }
 
+DetVariant mp_GetAlgorithmDet(matrix m, const ring r)
+{
+  if (MATROWS(m)+2*r->N>20+5*rField_is_Zp(r)) return DetMu;
+  if (MATROWS(m)<10+5*rField_is_Zp(r)) return DetSBareiss;
+  BOOLEAN isConst=TRUE;
+  int s=0;
+  for(int i=MATCOLS(m)*MATROWS(m)-1;i>=0;i--)
+  {
+    poly p=m->m[i];
+    if (p!=NULL)
+    {
+      if(!p_IsConstant(p,r)) isConst=FALSE;
+      s++;
+    }
+  }
+  if (isConst && rField_is_Q(r)) return DetFactory;
+  if (s*2<MATCOLS(m)*MATROWS(m)) // few entries
+    return DetSBareiss;
+  return DetMu;
+}
+DetVariant mp_GetAlgorithmDet(const char *s)
+{
+  if (strcmp(s,"Bareiss")==0) return DetBareiss;
+  if (strcmp(s,"SBareiss")==0) return DetSBareiss;
+  if (strcmp(s,"Mu")==0) return DetMu;
+  if (strcmp(s,"Factory")==0) return DetFactory;
+  WarnS("unknown method for det");
+  return DetDefault;
+}
+
+
+poly mp_Det(matrix a, const ring r, DetVariant d/*=DetDefault*/)
+{
+  if ((MATCOLS(a)==0)
+  && (MATROWS(a)==0))
+    return p_One(r);
+  if (d==DetDefault) d=mp_GetAlgorithmDet(a,r);
+  switch (d)
+  {
+    case DetBareiss: return mp_DetBareiss(a,r);
+    case DetMu: return mp_DetMu(a,r);
+    case DetFactory: return singclap_det(a,r);
+    case DetSBareiss:
+    {
+      ideal I=id_Matrix2Module(mp_Copy(a, r),r);
+      poly p=sm_CallDet(I, r);
+      id_Delete(&I, r);
+      return p;
+    }
+    default:
+      WerrorS("unknown algorith for det");
+  }
+}
+
+poly sm_Det(ideal a, const ring r, DetVariant d/*=DetDefault*/)
+{
+  if ((MATCOLS(a)==0)
+  && (MATROWS(a)==0))
+    return p_One(r);
+  if (d==DetDefault) d=mp_GetAlgorithmDet((matrix)a,r);
+  if (d==DetSBareiss) return sm_CallDet(a,r);
+  matrix m=id_Module2Matrix(id_Copy(a,r),r);
+  poly p=mp_Det(m,r,d);
+  id_Delete((ideal *)&m,r);
+  return p;
+}

libpolys/polys/matpol.h

@@ -32 +32 @@
 };
 
+enum DetVariant
+{
+  DetDefault=0,
+  DetBareiss,
+  DetSBareiss,
+  DetMu,
+  DetFactory
+};
+
 typedef ip_smatrix *       matrix;
 
@@ -59 +68 @@
 // BOOLEAN mpKoszul(leftv res,leftv b/*in*/, leftv c/*ip*/, leftv id=NULL);
 
+poly mp_Det(matrix a, const ring r, DetVariant d=DetDefault);
 poly mp_DetBareiss (matrix a, const ring r);
 poly mp_DetMu(matrix A, const ring R);
@@ -106 +116 @@
 BOOLEAN sm_Equal(ideal a, ideal b, const ring R);
 ideal sm_Tensor(ideal A, ideal B, const ring r);
+poly sm_Det(ideal I, const ring, DetVariant d=DetDefault);
+DetVariant mp_GetAlgorithmDet(matrix m, const ring r);
+DetVariant mp_GetAlgorithmDet(const char *s);
 
 #define SMATELEM(A,i,j,R) p_Vec2Poly(A->m[j],i+1,R)

libpolys/polys/shiftop.cc

@@ -368 +368 @@
   p_GetExpV(m,e,ri);
 
+  if (p_mLastVblock(m, e, ri) + sh > ri->N/lV)
+  {
+    Werror("letterplace degree bound is %d, but at least %d is needed for this shift", ri->N/lV, p_mLastVblock(m, e, ri) + sh);
+  }
   for (int i = ri->N - sh*lV; i > 0; i--)
   {
@@ -506 +510 @@
   PrintLn(); WriteLPExpV(m2ExpV, ri);
 #endif
-  if (m1Length + m2Length > ri->N)
-  {
-    WarnS("letterplace degree bound too low for this multiplication");
-  }
-  for (int i = 1 + m1Length; i < 1 + m1Length + m2Length; ++i)
+  int last = m1Length + m2Length;
+  if (last > ri->N)
+  {
+    Werror("letterplace degree bound is %d, but at least %d is needed for this multiplication", ri->N/ri->isLPring, last/ri->isLPring);
+    last = ri->N;
+  }
+  for (int i = 1 + m1Length; i < 1 + last; ++i)
   {
     assume(m2ExpV[i - m1Length] <= 1);
@@ -531 +537 @@
   PrintLn(); WriteLPExpV(m2ExpV, ri);
 #endif
-  if (m1Length + m2Length > ri->N)
-  {
-    WarnS("letterplace degree bound too low for this multiplication");
+  int last = m1Length + m2Length;
+  if (last > ri->N)
+  {
+    Werror("letterplace degree bound is %d, but at least %d is needed for this multiplication", ri->N/ri->isLPring, last/ri->isLPring);
+    last = ri->N;
   }
 
   // shift m1 by m2Length
-  for (int i = m2Length + m1Length; i >= 1 + m2Length; --i)
+  for (int i = last; i >= 1 + m2Length; --i)
   {
     m1ExpV[i] = m1ExpV[i - m2Length];

libpolys/polys/sparsmat.cc

@@ -295 +295 @@
 }
 
-/*2
-* Bareiss or Chinese remainder ?
-* I is dXd
-* sw = TRUE  -> I Matrix
-*      FALSE -> I Module
-* return True  -> change Type
-*        FALSE -> same Type
-*/
-BOOLEAN sm_CheckDet(ideal I, int d, BOOLEAN sw, const ring r)
-{
-  int s,t,i;
-  poly p;
-
-  if (d>100)
-    return TRUE;
-  if (!rField_is_Q(r))
-    return TRUE;
-  s = t = 0;
-  // now: field is Q, n<=100
-  for(i=IDELEMS(I)-1;i>=0;i--)
-  {
-    p=I->m[i];
-    if (p!=NULL)
-    {
-      if(!p_IsConstant(p,r))
-        return TRUE;
-      s++;
-      t+=n_Size(pGetCoeff(p),r->cf);
-    }
-  }
-  s*=15;
-  if (t>s)
-    return FALSE;// few large constanst entries
-  else
-    return TRUE; //many small entries
-}
-
 /* ----------------- basics (used from 'C') ------------------ */
 /*2

libpolys/polys/sparsmat.h

@@ -32 +32 @@
 void sm_KillModifiedRing(ring r);
 long sm_ExpBound(ideal, int, int, int, const ring);
-BOOLEAN sm_CheckDet(ideal, int, BOOLEAN, const ring);
 #endif