Changeset cec9624 in git

Timestamp:

Jul 11, 2019, 2:59:38 PM (5 years ago)

Author:

Hans Schoenemann <hannes@…>

Branches:

(u'spielwiese', 'fe61d9c35bf7c61f2b6cbf1b56e25e2f08d536cc')

Children:

dffcf06de72c5d8f3347abee6740a159f6cabd41

Parents:

18296cd9190bb2f8c652446aa5cf3a783ad2cb66

Message:

add: pp_DivideM, pp_Divide

Files:

• kernel/numeric/mpr_base.cc (modified) (2 diffs)
• kernel/polys.cc (modified) (1 diff)
• kernel/polys.h (modified) (1 diff)
• libpolys/polys/monomials/p_polys.cc (modified) (1 diff)
• libpolys/polys/monomials/p_polys.h (modified) (1 diff)

kernel/numeric/mpr_base.cc

@@ -2398 +2398 @@
       // compute row poly
       poly pi= ppMult_qq( (gls->m)[ resVectorList[k].elementOfS ] , resVectorList[k].mon );
-      pi= pDivideM( pCopy( pi ), pCopy( resVectorList[k].dividedBy ) );
+      pi= pp_DivideM( pi, resVectorList[k].dividedBy, currRing  );
 
       // fill in "matrix"
@@ -2430 +2430 @@
 
       pi= (gls->m)[ resVectorList[k].elementOfS ];
-      factor= pDivideM( pCopy( resVectorList[k].mon ), pCopy( resVectorList[k].dividedBy ) );
+      factor= pp_DivideM( resVectorList[k].mon, resVectorList[k].dividedBy, currRing );
 
       j=0;

kernel/polys.cc

@@ -179 +179 @@
 }
 
+poly pp_Divide(poly p, poly q, const ring r)
+{
+  assume(q!=NULL);
+  if (q==NULL)
+  {
+    WerrorS("div. by 0");
+    return NULL;
+  }
+  if (p==NULL)
+  {
+    return NULL;
+  }
+  if (pNext(q)!=NULL)
+  { /* This means that q != 0 consists of at least two terms*/
+    if (rIsLPRing(r))
+    {
+      WerrorS("not implemented for letterplace rings");
+      return NULL;
+    }
+    if(p_GetComp(p,r)==0)
+    {
+      if((rFieldType(r)==n_transExt)
+      &&(convSingTrP(p,r))
+      &&(convSingTrP(q,r)))
+      {
+        poly res=singclap_pdivide(p, q, r);
+        return res;
+      }
+      else if ((r->cf->convSingNFactoryN!=ndConvSingNFactoryN)
+      &&(!rField_is_Ring(r)))
+      {
+        poly res=singclap_pdivide(p, q, r);
+        return res;
+      }
+      else
+      {
+        ideal vi=idInit(1,1); vi->m[0]=p_Copy(q,r);
+        ideal ui=idInit(1,1); ui->m[0]=p_Copy(p,r);
+        ideal R; matrix U;
+        ring save_ring=currRing;
+        if (r!=currRing) rChangeCurrRing(r);
+        int save_opt;
+        SI_SAVE_OPT1(save_opt);
+        si_opt_1 &= ~(Sy_bit(OPT_PROT));
+        ideal m = idLift(vi,ui,&R, FALSE,TRUE,TRUE,&U);
+        SI_RESTORE_OPT1(save_opt);
+        if (r!=save_ring) rChangeCurrRing(save_ring);
+        matrix T = id_Module2formatedMatrix(m,1,1,r);
+        p=MATELEM(T,1,1); MATELEM(T,1,1)=NULL;
+        id_Delete((ideal *)&T,r);
+        id_Delete((ideal *)&U,r);
+        id_Delete(&R,r);
+        //vi->m[0]=NULL; ui->m[0]=NULL;
+        id_Delete(&vi,r);
+        id_Delete(&ui,r);
+        return p;
+      }
+    }
+    else
+    {
+      p=p_Copy(p,r);
+      int comps=p_MaxComp(p,r);
+      ideal I=idInit(comps,1);
+      poly h;
+      int i;
+      // conversion to a list of polys:
+      while (p!=NULL)
+      {
+        i=p_GetComp(p,r)-1;
+        h=pNext(p);
+        pNext(p)=NULL;
+        p_SetComp(p,0,r);
+        I->m[i]=p_Add_q(I->m[i],p,r);
+        p=h;
+      }
+      // division and conversion to vector:
+      h=NULL;
+      p=NULL;
+      q=p_Copy(q,r);
+      for(i=comps-1;i>=0;i--)
+      {
+        if (I->m[i]!=NULL)
+        {
+          if((rFieldType(r)==n_transExt)
+          &&(convSingTrP(I->m[i],r))
+          &&(convSingTrP(q,r)))
+          {
+            h=singclap_pdivide(I->m[i],q,r);
+          }
+          else if ((r->cf->convSingNFactoryN!=ndConvSingNFactoryN)
+          &&(!rField_is_Ring(r)))
+            h=singclap_pdivide(I->m[i],q,r);
+          else
+          {
+            ideal vi=idInit(1,1); vi->m[0]=q;
+            ideal ui=idInit(1,1); ui->m[0]=I->m[i];
+            ideal R; matrix U;
+            ring save_ring=currRing;
+            if (r!=currRing) rChangeCurrRing(r);
+            int save_opt;
+            SI_SAVE_OPT1(save_opt);
+            si_opt_1 &= ~(Sy_bit(OPT_PROT));
+            ideal m = idLift(vi,ui,&R, FALSE,TRUE,TRUE,&U);
+            SI_RESTORE_OPT1(save_opt);
+            if (r!=save_ring) rChangeCurrRing(save_ring);
+            if (idIs0(R))
+            {
+              matrix T = id_Module2formatedMatrix(m,1,1,r);
+              p=MATELEM(T,1,1); MATELEM(T,1,1)=NULL;
+              id_Delete((ideal *)&T,r);
+            }
+            else p=NULL;
+            id_Delete((ideal*)&U,r);
+            id_Delete(&R,r);
+            vi->m[0]=NULL; ui->m[0]=NULL;
+            id_Delete(&vi,r);
+            id_Delete(&ui,r);
+          }
+          p_SetCompP(h,i+1,r);
+          p=p_Add_q(p,h,r);
+        }
+      }
+      id_Delete(&I,r);
+      p_Delete(&q,r);
+      return p;
+    }
+  }
+  else
+  { /* This means that q != 0 consists of just one term,
+       or that r is over a coefficient ring. */
+#ifdef HAVE_RINGS
+    if (!rField_is_Domain(r))
+    {
+      WerrorS("division only defined over coefficient domains");
+      return NULL;
+    }
+    if (pNext(q)!=NULL)
+    {
+      WerrorS("division over a coefficient domain only implemented for terms");
+      return NULL;
+    }
+#endif
+    return pp_DivideM(p,q,r);
+  }
+  return NULL;
+}
+
 poly p_DivRem(poly p, poly q, poly &rest, const ring r)
 {

kernel/polys.h

@@ -164 +164 @@
 /// destroyes a,b
 poly p_Divide(poly a, poly b, const ring r);
+/// polynomial division a/b, ignoring the rest
+/// via singclap_pdivide resp. idLift
+/// does not destroy a,b
+poly pp_Divide(poly a, poly b, const ring r);
 poly p_DivRem(poly a, poly b, poly &rest, const ring r); /*julia*/
 

libpolys/polys/monomials/p_polys.cc

@@ -1610 +1610 @@
   }
   p_Delete(&b, r);
+  return result;
+}
+
+poly pp_DivideM(poly a, poly b, const ring r)
+{
+  if (a==NULL) { return NULL; }
+  poly result=a;
+
+  if(!p_IsConstant(b,r))
+  {
+    if (rIsLPRing(r))
+    {
+      WerrorS("not implemented for letterplace rings");
+      return NULL;
+    }
+    poly prev=NULL;
+    a=p_Copy(a,r);
+    while (a!=NULL)
+    {
+      if (p_DivisibleBy(b,a,r))
+      {
+        p_ExpVectorSub(a,b,r);
+        prev=a;
+        pIter(a);
+      }
+      else
+      {
+        if (prev==NULL)
+        {
+          p_LmDelete(&result,r);
+          a=result;
+        }
+        else
+        {
+          p_LmDelete(&pNext(prev),r);
+          a=pNext(prev);
+        }
+      }
+    }
+  }
+  if (result!=NULL)
+  {
+    number inv=pGetCoeff(b);
+    //if ((!rField_is_Ring(r)) || n_IsUnit(inv,r->cf))
+    if (rField_is_Zp(r))
+    {
+      inv = n_Invers(inv,r->cf);
+      __p_Mult_nn(result,inv,r);
+      n_Delete(&inv, r->cf);
+    }
+    else
+    {
+      result = p_Div_nn(result,inv,r);
+    }
+  }
   return result;
 }

libpolys/polys/monomials/p_polys.h

@@ -1990 +1990 @@
 poly      p_MDivide(poly a, poly b, const ring r);
 poly      p_DivideM(poly a, poly b, const ring r);
+poly      pp_DivideM(poly a, poly b, const ring r);
 poly      p_Div_nn(poly p, const number n, const ring r);