source: git/kernel/polys.cc @ a1ef3a2

fieker-DuValspielwiese
Last change on this file since a1ef3a2 was 83aaa6, checked in by Hans Schoenemann <hannes@…>, 3 years ago
format
  • Property mode set to 100644
File size: 11.4 KB
Line 
1#include "kernel/mod2.h"
2
3#include "misc/options.h"
4
5#include "polys.h"
6#include "kernel/ideals.h"
7#include "kernel/ideals.h"
8#include "polys/clapsing.h"
9#include "polys/clapconv.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.
13VAR ring  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)||rIsPluralRing(r))
45  { /* This means that q != 0 consists of at least two terms*/
46    if(p_GetComp(p,r)==0)
47    {
48      if((rFieldType(r)==n_transExt)
49      &&(convSingTrP(p,r))
50      &&(convSingTrP(q,r))
51      &&(!rIsNCRing(r)))
52      {
53        poly res=singclap_pdivide(p, q, r);
54        p_Delete(&p,r);
55        p_Delete(&q,r);
56        return res;
57      }
58      else if ((r->cf->convSingNFactoryN!=ndConvSingNFactoryN)
59      &&(!rField_is_Ring(r))
60      &&(!rIsNCRing(r)))
61      {
62        poly res=singclap_pdivide(p, q, r);
63        p_Delete(&p,r);
64        p_Delete(&q,r);
65        return res;
66      }
67      else
68      {
69        ideal vi=idInit(1,1); vi->m[0]=q;
70        ideal ui=idInit(1,1); ui->m[0]=p;
71        ideal R; matrix U;
72        ring save_ring=currRing;
73        if (r!=currRing) rChangeCurrRing(r);
74        int save_opt;
75        SI_SAVE_OPT1(save_opt);
76        si_opt_1 &= ~(Sy_bit(OPT_PROT));
77        ideal m = idLift(vi,ui,&R, FALSE,TRUE,TRUE,&U);
78        SI_RESTORE_OPT1(save_opt);
79        if (r!=save_ring) rChangeCurrRing(save_ring);
80        matrix T = id_Module2formatedMatrix(m,1,1,r);
81        p=MATELEM(T,1,1); MATELEM(T,1,1)=NULL;
82        id_Delete((ideal *)&T,r);
83        id_Delete((ideal *)&U,r);
84        id_Delete(&R,r);
85        //vi->m[0]=NULL; ui->m[0]=NULL;
86        id_Delete(&vi,r);
87        id_Delete(&ui,r);
88        return p;
89      }
90    }
91    else
92    {
93      int comps=p_MaxComp(p,r);
94      ideal I=idInit(comps,1);
95      poly h;
96      int i;
97      // conversion to a list of polys:
98      while (p!=NULL)
99      {
100        i=p_GetComp(p,r)-1;
101        h=pNext(p);
102        pNext(p)=NULL;
103        p_SetComp(p,0,r);
104        I->m[i]=p_Add_q(I->m[i],p,r);
105        p=h;
106      }
107      // division and conversion to vector:
108      h=NULL;
109      p=NULL;
110      for(i=comps-1;i>=0;i--)
111      {
112        if (I->m[i]!=NULL)
113        {
114          if((rFieldType(r)==n_transExt)
115          &&(convSingTrP(I->m[i],r))
116          &&(convSingTrP(q,r))
117          &&(!rIsNCRing(r)))
118          {
119            h=singclap_pdivide(I->m[i],q,r);
120          }
121          else if ((r->cf->convSingNFactoryN!=ndConvSingNFactoryN)
122          &&(!rField_is_Ring(r))
123          &&(!rIsNCRing(r)))
124            h=singclap_pdivide(I->m[i],q,r);
125          else
126          {
127            ideal vi=idInit(1,1); vi->m[0]=q;
128            ideal ui=idInit(1,1); ui->m[0]=I->m[i];
129            ideal R; matrix U;
130            ring save_ring=currRing;
131            if (r!=currRing) rChangeCurrRing(r);
132            int save_opt;
133            SI_SAVE_OPT1(save_opt);
134            si_opt_1 &= ~(Sy_bit(OPT_PROT));
135            ideal m = idLift(vi,ui,&R, FALSE,TRUE,TRUE,&U);
136            SI_RESTORE_OPT1(save_opt);
137            if (r!=save_ring) rChangeCurrRing(save_ring);
138            if (idIs0(R))
139            {
140              matrix T = id_Module2formatedMatrix(m,1,1,r);
141              p=MATELEM(T,1,1); MATELEM(T,1,1)=NULL;
142              id_Delete((ideal *)&T,r);
143            }
144            else p=NULL;
145            id_Delete((ideal*)&U,r);
146            id_Delete(&R,r);
147            vi->m[0]=NULL; ui->m[0]=NULL;
148            id_Delete(&vi,r);
149            id_Delete(&ui,r);
150          }
151          p_SetCompP(h,i+1,r);
152          p=p_Add_q(p,h,r);
153        }
154      }
155      id_Delete(&I,r);
156      p_Delete(&q,r);
157      return p;
158    }
159  }
160  else
161  { /* This means that q != 0 consists of just one term, or LetterPlace */
162#ifdef HAVE_RINGS
163    if (pNext(q)!=NULL)
164    {
165      WerrorS("division over a coefficient domain only implemented for terms");
166      return NULL;
167    }
168#endif
169    return p_DivideM(p,q,r);
170  }
171  return NULL;
172}
173
174poly pp_Divide(poly p, poly q, const ring r)
175{
176  assume(q!=NULL);
177  if (q==NULL)
178  {
179    WerrorS("div. by 0");
180    return NULL;
181  }
182  if (p==NULL)
183  {
184    return NULL;
185  }
186  if ((pNext(q)!=NULL)||rIsPluralRing(r))
187  { /* This means that q != 0 consists of at least two terms*/
188    if(p_GetComp(p,r)==0)
189    {
190      if((rFieldType(r)==n_transExt)
191      &&(convSingTrP(p,r))
192      &&(convSingTrP(q,r))
193      &&(!rIsNCRing(r)))
194      {
195        poly res=singclap_pdivide(p, q, r);
196        return res;
197      }
198      else if ((r->cf->convSingNFactoryN!=ndConvSingNFactoryN)
199      &&(!rField_is_Ring(r))
200      &&(!rIsNCRing(r)))
201      {
202        poly res=singclap_pdivide(p, q, r);
203        return res;
204      }
205      else
206      {
207        ideal vi=idInit(1,1); vi->m[0]=p_Copy(q,r);
208        ideal ui=idInit(1,1); ui->m[0]=p_Copy(p,r);
209        ideal R; matrix U;
210        ring save_ring=currRing;
211        if (r!=currRing) rChangeCurrRing(r);
212        int save_opt;
213        SI_SAVE_OPT1(save_opt);
214        si_opt_1 &= ~(Sy_bit(OPT_PROT));
215        ideal m = idLift(vi,ui,&R, FALSE,TRUE,TRUE,&U);
216        SI_RESTORE_OPT1(save_opt);
217        if (r!=save_ring) rChangeCurrRing(save_ring);
218        matrix T = id_Module2formatedMatrix(m,1,1,r);
219        p=MATELEM(T,1,1); MATELEM(T,1,1)=NULL;
220        id_Delete((ideal *)&T,r);
221        id_Delete((ideal *)&U,r);
222        id_Delete(&R,r);
223        //vi->m[0]=NULL; ui->m[0]=NULL;
224        id_Delete(&vi,r);
225        id_Delete(&ui,r);
226        return p;
227      }
228    }
229    else
230    {
231      p=p_Copy(p,r);
232      int comps=p_MaxComp(p,r);
233      ideal I=idInit(comps,1);
234      poly h;
235      int i;
236      // conversion to a list of polys:
237      while (p!=NULL)
238      {
239        i=p_GetComp(p,r)-1;
240        h=pNext(p);
241        pNext(p)=NULL;
242        p_SetComp(p,0,r);
243        I->m[i]=p_Add_q(I->m[i],p,r);
244        p=h;
245      }
246      // division and conversion to vector:
247      h=NULL;
248      p=NULL;
249      q=p_Copy(q,r);
250      for(i=comps-1;i>=0;i--)
251      {
252        if (I->m[i]!=NULL)
253        {
254          if((rFieldType(r)==n_transExt)
255          &&(convSingTrP(I->m[i],r))
256          &&(convSingTrP(q,r))
257          &&(!rIsNCRing(r)))
258          {
259            h=singclap_pdivide(I->m[i],q,r);
260          }
261          else if ((r->cf->convSingNFactoryN!=ndConvSingNFactoryN)
262          &&(!rField_is_Ring(r))
263          &&(!rIsNCRing(r)))
264            h=singclap_pdivide(I->m[i],q,r);
265          else
266          {
267            ideal vi=idInit(1,1); vi->m[0]=q;
268            ideal ui=idInit(1,1); ui->m[0]=I->m[i];
269            ideal R; matrix U;
270            ring save_ring=currRing;
271            if (r!=currRing) rChangeCurrRing(r);
272            int save_opt;
273            SI_SAVE_OPT1(save_opt);
274            si_opt_1 &= ~(Sy_bit(OPT_PROT));
275            ideal m = idLift(vi,ui,&R, FALSE,TRUE,TRUE,&U);
276            SI_RESTORE_OPT1(save_opt);
277            if (r!=save_ring) rChangeCurrRing(save_ring);
278            if (idIs0(R))
279            {
280              matrix T = id_Module2formatedMatrix(m,1,1,r);
281              p=MATELEM(T,1,1); MATELEM(T,1,1)=NULL;
282              id_Delete((ideal *)&T,r);
283            }
284            else p=NULL;
285            id_Delete((ideal*)&U,r);
286            id_Delete(&R,r);
287            vi->m[0]=NULL; ui->m[0]=NULL;
288            id_Delete(&vi,r);
289            id_Delete(&ui,r);
290          }
291          p_SetCompP(h,i+1,r);
292          p=p_Add_q(p,h,r);
293        }
294      }
295      id_Delete(&I,r);
296      p_Delete(&q,r);
297      return p;
298    }
299  }
300  else
301  { /* This means that q != 0 consists of just one term,
302       or that r is over a coefficient ring. */
303#ifdef HAVE_RINGS
304    if (pNext(q)!=NULL)
305    {
306      WerrorS("division over a coefficient domain only implemented for terms");
307      return NULL;
308    }
309#endif
310    return pp_DivideM(p,q,r);
311  }
312  return NULL;
313}
314
315poly p_DivRem(poly p, poly q, poly &rest, const ring r)
316{
317  assume(q!=NULL);
318  rest=NULL;
319  if (q==NULL)
320  {
321    WerrorS("div. by 0");
322    return NULL;
323  }
324  if (p==NULL)
325  {
326    p_Delete(&q,r);
327    return NULL;
328  }
329  if(p_GetComp(p,r)==0)
330  {
331    if((rFieldType(r)==n_transExt)
332    &&(convSingTrP(p,r))
333    &&(convSingTrP(q,r))
334    &&(!rIsNCRing(r)))
335    {
336      poly res=singclap_pdivide(p, q, r);
337      rest=singclap_pmod(p,q,r);
338      p_Delete(&p,r);
339      p_Delete(&q,r);
340      return res;
341    }
342    else if ((r->cf->convSingNFactoryN!=ndConvSingNFactoryN)
343    &&(!rField_is_Ring(r))
344    &&(!rIsNCRing(r)))
345    {
346      poly res=singclap_pdivide(p, q, r);
347      rest=singclap_pmod(p,q,r);
348      p_Delete(&p,r);
349      p_Delete(&q,r);
350      return res;
351    }
352    else
353    {
354      ideal vi=idInit(1,1); vi->m[0]=q;
355      ideal ui=idInit(1,1); ui->m[0]=p;
356      ideal R; matrix U;
357      ring save_ring=currRing;
358      if (r!=currRing) rChangeCurrRing(r);
359      int save_opt;
360      SI_SAVE_OPT1(save_opt);
361      si_opt_1 &= ~(Sy_bit(OPT_PROT));
362      ideal m = idLift(vi,ui,&R, FALSE,TRUE,TRUE,&U);
363      SI_RESTORE_OPT1(save_opt);
364      if (r!=save_ring) rChangeCurrRing(save_ring);
365      matrix T = id_Module2formatedMatrix(m,1,1,r);
366      p=MATELEM(T,1,1); MATELEM(T,1,1)=NULL;
367      id_Delete((ideal *)&T,r);
368      T = id_Module2formatedMatrix(R,1,1,r);
369      rest=MATELEM(T,1,1); MATELEM(T,1,1)=NULL;
370      id_Delete((ideal *)&T,r);
371      id_Delete((ideal *)&U,r);
372      id_Delete(&R,r);
373      //vi->m[0]=NULL; ui->m[0]=NULL;
374      id_Delete(&vi,r);
375      id_Delete(&ui,r);
376      return p;
377    }
378  }
379  return NULL;
380}
381
382poly singclap_gcd ( poly f, poly g, const ring r )
383{
384  poly res=NULL;
385
386  if (f!=NULL)
387  {
388    //if (r->cf->has_simple_Inverse) p_Norm(f,r);
389    if (rField_is_Zp(r)) p_Norm(f,r);
390    else                 p_Cleardenom(f, r);
391  }
392  if (g!=NULL)
393  {
394    //if (r->cf->has_simple_Inverse) p_Norm(g,r);
395    if (rField_is_Zp(r)) p_Norm(g,r);
396    else                 p_Cleardenom(g, r);
397  }
398  else         return f; // g==0 => gcd=f (but do a p_Cleardenom/pNorm)
399  if (f==NULL) return g; // f==0 => gcd=g (but do a p_Cleardenom/pNorm)
400  if(!rField_is_Ring(r)
401  && (p_IsConstant(f,r)
402  ||p_IsConstant(g,r)))
403  {
404    res=p_One(r);
405  }
406  else if (r->cf->convSingNFactoryN!=ndConvSingNFactoryN)
407  {
408    res=singclap_gcd_r(f,g,r);
409  }
410  else
411  {
412    ideal I=idInit(2,1);
413    I->m[0]=f;
414    I->m[1]=p_Copy(g,r);
415    intvec *w=NULL;
416    ring save_ring=currRing;
417    if (r!=currRing) rChangeCurrRing(r);
418    int save_opt;
419    SI_SAVE_OPT1(save_opt);
420    si_opt_1 &= ~(Sy_bit(OPT_PROT));
421    ideal S1=idSyzygies(I,testHomog,&w);
422    if (w!=NULL) delete w;
423    // expect S1->m[0]=(-g/gcd,f/gcd)
424    if (IDELEMS(S1)!=1) WarnS("error in syzygy computation for GCD");
425    int lp;
426    p_TakeOutComp(&S1->m[0],1,&res,&lp,r);
427    p_Delete(&S1->m[0],r);
428    // GCD is g divided iby (-g/gcd):
429    res=p_Divide(g,res,r);
430    // restore, r, opt:
431    SI_RESTORE_OPT1(save_opt);
432    if (r!=save_ring) rChangeCurrRing(save_ring);
433    // clean the result
434    res=p_Cleardenom(res,r);
435    p_Content(res,r);
436    return res;
437  }
438  p_Delete(&f, r);
439  p_Delete(&g, r);
440  return res;
441}
Note: See TracBrowser for help on using the repository browser.