source: git/kernel/polys.cc @ 49c522

spielwiese
Last change on this file since 49c522 was 49c522, checked in by Hans Schoenemann <hannes@…>, 6 years ago
fix: Singular.jl/issues/23: polynomial division by a monomial
  • Property mode set to 100644
File size: 4.2 KB
Line 
1#include "kernel/mod2.h"
2
3#include "omalloc/omalloc.h"
4#include "misc/options.h"
5
6#include "polys.h"
7#include "kernel/ideals.h"
8#include "kernel/ideals.h"
9#include "polys/clapsing.h"
10
11/// Widely used global variable which specifies the current polynomial ring for Singular interpreter and legacy implementatins.
12/// @Note: one should avoid using it in newer designs, for example due to possible problems in parallelization with threads.
13ring  currRing = NULL;
14
15void rChangeCurrRing(ring r)
16{
17  //------------ set global ring vars --------------------------------
18  currRing = r;
19  if( r != NULL )
20  {
21    rTest(r);
22    //------------ global variables related to coefficients ------------
23    assume( r->cf!= NULL );
24    nSetChar(r->cf);
25    //------------ global variables related to polys
26    p_SetGlobals(r); // also setting TEST_RINGDEP_OPTS
27    //------------ global variables related to factory -----------------
28  }
29}
30
31poly p_Divide(poly p, poly q, const ring r)
32{
33  assume(q!=NULL);
34  if (q==NULL)
35  {
36    WerrorS("div. by 0");
37    return NULL;
38  }
39  if (p==NULL)
40  {
41    p_Delete(&q,r);
42    return NULL;
43  }
44  if (pNext(q)!=NULL)
45  { /* This means that q != 0 consists of at least two terms*/
46    if(p_GetComp(p,r)==0)
47    {
48      if ((r->cf->convSingNFactoryN!=ndConvSingNFactoryN)
49      &&(!rField_is_Ring(r)))
50      {
51        poly res=singclap_pdivide(p, q, r);
52        p_Delete(&p,r);
53        p_Delete(&q,r);
54        return res;
55      }
56      else
57      {
58        ideal vi=idInit(1,1); vi->m[0]=q;
59        ideal ui=idInit(1,1); ui->m[0]=p;
60        ideal R; matrix U;
61        ring save_ring=currRing;
62        if (r!=currRing) rChangeCurrRing(r);
63        int save_opt;
64        SI_SAVE_OPT1(save_opt);
65        si_opt_1 &= ~(Sy_bit(OPT_PROT));
66        ideal m = idLift(vi,ui,&R, FALSE,TRUE,TRUE,&U);
67        SI_RESTORE_OPT1(save_opt);
68        if (r!=save_ring) rChangeCurrRing(save_ring);
69        matrix T = id_Module2formatedMatrix(m,1,1,r);
70        p=MATELEM(T,1,1); MATELEM(T,1,1)=NULL;
71        id_Delete((ideal *)&T,r);
72        id_Delete((ideal *)&U,r);
73        id_Delete(&R,r);
74        //vi->m[0]=NULL; ui->m[0]=NULL;
75        id_Delete(&vi,r);
76        id_Delete(&ui,r);
77        return p;
78      }
79    }
80    else
81    {
82      int comps=p_MaxComp(p,r);
83      ideal I=idInit(comps,1);
84      poly h;
85      int i;
86      // conversion to a list of polys:
87      while (p!=NULL)
88      {
89        i=p_GetComp(p,r)-1;
90        h=pNext(p);
91        pNext(p)=NULL;
92        p_SetComp(p,0,r);
93        I->m[i]=p_Add_q(I->m[i],p,r);
94        p=h;
95      }
96      // division and conversion to vector:
97      h=NULL;
98      p=NULL;
99      for(i=comps-1;i>=0;i--)
100      {
101        if (I->m[i]!=NULL)
102        {
103          if ((r->cf->convSingNFactoryN!=ndConvSingNFactoryN)
104          &&(!rField_is_Ring(r)))
105            h=singclap_pdivide(I->m[i],q,r);
106          else
107          {
108            ideal vi=idInit(1,1); vi->m[0]=q;
109            ideal ui=idInit(1,1); ui->m[0]=I->m[i];
110            ideal R; matrix U;
111            ring save_ring=currRing;
112            if (r!=currRing) rChangeCurrRing(r);
113            int save_opt;
114            SI_SAVE_OPT1(save_opt);
115            si_opt_1 &= ~(Sy_bit(OPT_PROT));
116            ideal m = idLift(vi,ui,&R, FALSE,TRUE,TRUE,&U);
117            SI_RESTORE_OPT1(save_opt);
118            if (r!=save_ring) rChangeCurrRing(save_ring);
119            matrix T = id_Module2formatedMatrix(m,1,1,r);
120            h=MATELEM(T,1,1); MATELEM(T,1,1)=NULL;
121            id_Delete((ideal*)&T,r);
122            id_Delete((ideal*)&U,r);
123            id_Delete(&R,r);
124            vi->m[0]=NULL; ui->m[0]=NULL;
125            id_Delete(&vi,r);
126            id_Delete(&ui,r);
127          }
128          p_SetCompP(h,i+1,r);
129          p=p_Add_q(p,h,r);
130        }
131      }
132      id_Delete(&I,r);
133      p_Delete(&q,r);
134      return p;
135    }
136  }
137  else
138  { /* This means that q != 0 consists of just one term,
139       or that r is over a coefficient ring. */
140#ifdef HAVE_RINGS
141    if (!rField_is_Domain(r))
142    {
143      WerrorS("division only defined over coefficient domains");
144      return NULL;
145    }
146    if (pNext(q)!=NULL)
147    {
148      WerrorS("division over a coefficient domain only implemented for terms");
149      return NULL;
150    }
151#endif
152    return p_DivideM(p,q,r);
153  }
154  return FALSE;
155}
Note: See TracBrowser for help on using the repository browser.