Changeset e68029 in git

Timestamp:

Jan 19, 2022, 9:40:11 AM (17 months ago)

Author:

Hans Schoenemann <hannes@…>

Branches:

(u'jengelh-datetime', 'ceac47cbc86fe4a15902392bdbb9bd2ae0ea02c6')(u'spielwiese', '23b0a9c07070b684aa12355009488dc00d9501e3')

Children:

3882094d003ecd5487dfad0aa941dee9d650a30e

Parents:

844f3006400f1bbcfdb84aa86a873413a954da98

Message:

rref for smatrix and ntl

Files:

• Singular/extra.cc (modified) (1 diff)
• libpolys/polys/clapsing.cc (modified) (2 diffs)
• libpolys/polys/clapsing.h (modified) (1 diff)
• libpolys/polys/flintconv.cc (modified) (2 diffs)
• libpolys/polys/flintconv.h (modified) (1 diff)

Singular/extra.cc

@@ -1239 +1239 @@
   #endif
 /* ====== rref ============================*/
-  #ifdef HAVE_FLINT
+  #if defined(HAVE_FLINT) || defined(HAVE_NTL)
   if(strcmp(sys_cmd,"rref")==0)
   {
     const short t1[]={1,MATRIX_CMD};
-    if (iiCheckTypes(h,t1,1))
+    const short t2[]={1,SMATRIX_CMD};
+    if (iiCheckTypes(h,t1,0))
     {
       matrix M=(matrix)h->Data();
+      #if defined(HAVE_FLINT)
       res->data=(void*)singflint_rref(M,currRing);
+      #elif defined(HAVE_NTL)
+      res->data=(void*)singntl_rref(M,currRing);
+      #endif
       res->rtyp=MATRIX_CMD;
       return FALSE;
     }
-    else
-    {
-      WerrorS("expected system(\"rref\",<matrix>)");
+    else if (iiCheckTypes(h,t2,1))
+    {
+      ideal M=(ideal)h->Data();
+      #if defined(HAVE_FLINT)
+      res->data=(void*)singflint_rref(M,currRing);
+      #elif defined(HAVE_NTL)
+      res->data=(void*)singntl_rref(M,currRing);
+      #endif
+      res->rtyp=SMATRIX_CMD;
+      return FALSE;
+    }
+    else
+    {
+      WerrorS("expected system(\"rref\",<matrix>/<smatrix>)");
       return TRUE;
     }

libpolys/polys/clapsing.cc

@@ -24 +24 @@
 #include "polys/flint_mpoly.h"
 
+#ifdef HAVE_NTL
+#include <NTL/config.h>
+#ifdef NTL_STD_CXX
+#ifdef NOSTREAMIO
+#  ifdef HAVE_IOSTREAM
+#    include <iostream>
+#    define OSTREAM std::ostream
+#    define ISTREAM std::istream
+#  elif defined(HAVE_IOSTREAM_H)
+#    include <iostream.h>
+#    define OSTREAM ostream
+#    define ISTREAM istream
+#  endif
+#endif /* ! NOSTREAMIO */
+#endif
+#include <NTL/mat_ZZ.h>
+#include <NTL/mat_lzz_p.h>
+
+#ifdef NTL_CLIENT               // in <NTL/tools.h>: using of name space NTL
+NTL_CLIENT
+#endif
+
+#endif
 
 //#include "polys.h"
@@ -1971 +1994 @@
 #endif
 
+#ifdef HAVE_NTL
+matrix singntl_rref(matrix  m, const ring R)
+{
+  int r=m->rows();
+  int c=m->cols();
+  int i,j;
+  matrix M=mpNew(r,c);
+  if (rField_is_Zp(R))
+  {
+    zz_p::init(rChar(R));
+    mat_zz_p *NTLM=new mat_zz_p;
+    NTLM->SetDims(r,c);
+    for(i=r;i>0;i--)
+    {
+      for(j=c;j>0;j--)
+      {
+        poly h=MATELEM(m,i,j);
+        if ((h!=NULL)
+        && (p_Totaldegree(h,R)==0))
+        {
+          (*NTLM)(i,j)=(long)p_GetCoeff(h,R);
+        }
+        else
+        {
+          WerrorS("smatrix for rref is not constant");
+          return M;
+        }
+      }
+    }
+    gauss(*NTLM);
+    for(i=r;i>0;i--)
+    {
+      for(j=c;j>0;j--)
+      {
+        number n=n_Init(rep((*NTLM)(i,j)),R->cf);
+        if(!n_IsZero(n,R->cf))
+        {
+          poly p=p_NSet(n,R);
+          p_SetComp(p,i,R);
+          M->m[j]=p_Add_q(M->m[j],p,R);
+        }
+      }
+    }
+    delete NTLM;
+  }
+  else
+  {
+    WerrorS("not implemented for these coefficients");
+  }
+  return M;
+}
+#endif
+
+#ifdef HAVE_NTL
+ideal singntl_rref(ideal  m, const ring R) /*assume smatrix m*/
+{
+  int r=m->rank;
+  int c=m->ncols;
+  int i,j;
+  ideal M=idInit(c,r);
+  if (rField_is_Zp(R))
+  {
+    zz_p::init(rChar(R));
+    mat_zz_p *NTLM=new mat_zz_p;
+    NTLM->SetDims(r,c);
+    for(j=c;j>0;j--)
+    {
+      poly h=m->m[j];
+      while(h!=NULL)
+      {
+        i=p_GetComp(h,R);
+        if (p_Totaldegree(h,R)==0)
+          (*NTLM)(i,j)=(long)p_GetCoeff(h,R);
+        else
+        {
+          WerrorS("smatrix for rref is not constant");
+          return M;
+        }
+        pIter(h);
+      }
+    }
+    gauss(*NTLM);
+    for(i=r;i>0;i--)
+    {
+      for(j=c;j>0;j--)
+      {
+        number n=n_Init(rep((*NTLM)(i,j)),R->cf);
+        if(!n_IsZero(n,R->cf))
+        {
+          poly p=p_NSet(n,R);
+          p_SetComp(p,i,R);
+          M->m[j]=p_Add_q(M->m[j],p,R);
+        }
+      }
+    }
+    delete NTLM;
+  }
+  else
+  {
+    WerrorS("not implemented for these coefficients");
+  }
+  return M;
+}
+#endif
+
 #if defined(HAVE_NTL) || defined(HAVE_FLINT)
 ideal singclap_absFactorize ( poly f, ideal & mipos, intvec ** exps, int & numFactors, const ring r)

libpolys/polys/clapsing.h

@@ -44 +44 @@
 intvec* singntl_LLL(intvec* A);
 
+#ifdef HAVE_NTL
+ideal  singntl_rref(ideal  m, const ring R) /*assume smatrix m*/;
+matrix singntl_rref(matrix  m, const ring R);
+#endif
+
 ideal singclap_absFactorize ( poly f, ideal & mipos, intvec ** exps, int & n, const ring r);
 

libpolys/polys/flintconv.cc

@@ -391 +391 @@
       {
         poly h=MATELEM(m,i,j);
-        if (h==NULL)
-          convSingNFlintN(fmpq_mat_entry(FLINTM,i-1,j-1),n,R->cf);
-        else if (p_Totaldegree(h,R)==0)
+        if ((h!=NULL)
+        && (p_Totaldegree(h,R)==0))
           convSingNFlintN(fmpq_mat_entry(FLINTM,i-1,j-1),pGetCoeff(h),R->cf);
-        else
-        {
-          WerrorS("matrix for rref is nor constant");
-          return M;
-        }
+        else
+        {
+          WerrorS("matrix for rref is not constant");
+          return M;
+        }
       }
     }
@@ -447 +446 @@
 }
 
+ideal singflint_rref(ideal  m, const ring R) /*assume smatrix m*/
+{
+  int r=m->rank;
+  int c=m->ncols;
+  int i,j;
+  ideal M=idInit(c,r);
+  if (rField_is_Q(R))
+  {
+    fmpq_mat_t FLINTM;
+    fmpq_mat_init(FLINTM,r,c);
+    for(j=c;j>0;j--)
+    {
+      poly h=m->m[j];
+      while(h!=NULL)
+      {
+        i=p_GetComp(h,R);
+        if (p_Totaldegree(h,R)==0)
+          convSingNFlintN(fmpq_mat_entry(FLINTM,i-1,j-1),p_GetCoeff(h,R),R->cf);
+        else
+        {
+          WerrorS("smatrix for rref is not constant");
+          return M;
+        }
+        pIter(h);
+      }
+    }
+    fmpq_mat_rref(FLINTM,FLINTM);
+    for(i=r;i>0;i--)
+    {
+      for(j=c;j>0;j--)
+      {
+        number n=convFlintNSingN(fmpq_mat_entry(FLINTM,i-1,j-1),R->cf);
+        if(!n_IsZero(n,R->cf))
+        {
+          poly p=p_NSet(n,R);
+          p_SetComp(p,i,R);
+          M->m[j]=p_Add_q(M->m[j],p,R);
+        }
+      }
+    }
+    fmpq_mat_clear(FLINTM);
+  }
+  else if (rField_is_Zp(R))
+  {
+    nmod_mat_t FLINTM;
+    nmod_mat_init(FLINTM,r,c,rChar(R));
+    for(j=c;j>0;j--)
+    {
+      poly h=m->m[j];
+      while(h!=NULL)
+      {
+        i=p_GetComp(h,R);
+        if (p_Totaldegree(h,R)==0)
+          nmod_mat_entry(FLINTM,i-1,j-1)=(long)p_GetCoeff(h,R);
+        else
+        {
+          WerrorS("smatrix for rref is not constant");
+          return M;
+        }
+        pIter(h);
+      }
+    }
+    nmod_mat_rref(FLINTM);
+    for(i=r;i>0;i--)
+    {
+      for(j=c;j>0;j--)
+      {
+        number n=n_Init(nmod_mat_entry(FLINTM,i-1,j-1),R->cf);
+        if(!n_IsZero(n,R->cf))
+        {
+          poly p=p_NSet(n,R);
+          p_SetComp(p,i,R);
+          M->m[j]=p_Add_q(M->m[j],p,R);
+        }
+      }
+    }
+    nmod_mat_clear(FLINTM);
+  }
+  else
+  {
+    #if 0
+    fmpz_t p;
+    convSingIFlintI(p,rChar(currRing));
+    fq_nmod_ctx_init(ctx,p,1,"t");
+    fq_nmod_mat_t FLINTM;
+    // convert matrix
+    convSingMFlintFq_nmod_mat(M,FLINTM,ctx,currRing);
+    // rank
+    long rk= fq_nmod_mat_rref (FLINTM,ctx);
+    res->data=(void*)convFlintFq_nmod_matSingM(FLINTM,ctx,currRing);
+    // clean up
+    fq_nmod_mat_clear (FLINTM,ctx);
+    fq_nmod_ctx_clear(ctx);
+    fmpz_clear(p);
+    #endif
+    WerrorS("not implemented for these coefficients");
+  }
+  return M;
+}
+
 bigintmat* singflint_LLL(bigintmat*  m, bigintmat* T)
 {

libpolys/polys/flintconv.h

@@ -67 +67 @@
 
 matrix singflint_rref(matrix  m, const ring R);
+ideal  singflint_rref(ideal  m, const ring R);
 #if __FLINT_RELEASE >= 20500
 void convSingPFlintnmod_poly_t(nmod_poly_t result, const poly p, const ring r);