source: git/kernel/polys.cc

spielwiese
Last change on this file was e4e64b, checked in by Hans Schoenemann <hannes@…>, 5 weeks ago
fix: save options with rChangeCurrRing
  • 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 implementations.
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  if (currRing!=NULL)
18    currRing->options = si_opt_1 & TEST_RINGDEP_OPTS;
19  //------------ set global ring vars --------------------------------
20  currRing = r;
21  if( r != NULL )
22  {
23    rTest(r);
24    //------------ global variables related to coefficients ------------
25    assume( r->cf!= NULL );
26    nSetChar(r->cf);
27    //------------ global variables related to polys
28    p_SetGlobals(r); // also setting TEST_RINGDEP_OPTS
29    //------------ global variables related to factory -----------------
30  }
31}
32
33poly p_Divide(poly p, poly q, const ring r)
34{
35  assume(q!=NULL);
36  if (q==NULL)
37  {
38    WerrorS("div. by 0");
39    return NULL;
40  }
41  if (p==NULL)
42  {
43    p_Delete(&q,r);
44    return NULL;
45  }
46  if ((pNext(q)!=NULL)||rIsPluralRing(r))
47  { /* This means that q != 0 consists of at least two terms*/
48    if(p_GetComp(p,r)==0)
49    {
50      if((rFieldType(r)==n_transExt)
51      &&(convSingTrP(p,r))
52      &&(convSingTrP(q,r))
53      &&(!rIsNCRing(r)))
54      {
55        poly res=singclap_pdivide(p, q, r);
56        p_Delete(&p,r);
57        p_Delete(&q,r);
58        return res;
59      }
60      else if ((r->cf->convSingNFactoryN!=ndConvSingNFactoryN)
61      &&(!rField_is_Ring(r))
62      &&(!rIsNCRing(r)))
63      {
64        poly res=singclap_pdivide(p, q, r);
65        p_Delete(&p,r);
66        p_Delete(&q,r);
67        return res;
68      }
69      else
70      {
71        ideal vi=idInit(1,1); vi->m[0]=q;
72        ideal ui=idInit(1,1); ui->m[0]=p;
73        ideal R; matrix U;
74        ring save_ring=currRing;
75        if (r!=currRing) rChangeCurrRing(r);
76        int save_opt;
77        SI_SAVE_OPT1(save_opt);
78        si_opt_1 &= ~(Sy_bit(OPT_PROT));
79        ideal m = idLift(vi,ui,&R, FALSE,TRUE,TRUE,&U);
80        SI_RESTORE_OPT1(save_opt);
81        if (r!=save_ring) rChangeCurrRing(save_ring);
82        p=m->m[0]; m->m[0]=NULL;
83        id_Delete(&m,r);
84        p_SetCompP(p,0,r);
85        id_Delete((ideal *)&U,r);
86        id_Delete(&R,r);
87        //vi->m[0]=NULL; ui->m[0]=NULL;
88        id_Delete(&vi,r);
89        id_Delete(&ui,r);
90        return p;
91      }
92    }
93    else
94    {
95      int comps=p_MaxComp(p,r);
96      ideal I=idInit(comps,1);
97      poly h;
98      int i;
99      // conversion to a list of polys:
100      while (p!=NULL)
101      {
102        i=p_GetComp(p,r)-1;
103        h=pNext(p);
104        pNext(p)=NULL;
105        p_SetComp(p,0,r);
106        I->m[i]=p_Add_q(I->m[i],p,r);
107        p=h;
108      }
109      // division and conversion to vector:
110      h=NULL;
111      p=NULL;
112      for(i=comps-1;i>=0;i--)
113      {
114        if (I->m[i]!=NULL)
115        {
116          if((rFieldType(r)==n_transExt)
117          &&(convSingTrP(I->m[i],r))
118          &&(convSingTrP(q,r))
119          &&(!rIsNCRing(r)))
120          {
121            h=singclap_pdivide(I->m[i],q,r);
122          }
123          else if ((r->cf->convSingNFactoryN!=ndConvSingNFactoryN)
124          &&(!rField_is_Ring(r))
125          &&(!rIsNCRing(r)))
126            h=singclap_pdivide(I->m[i],q,r);
127          else
128          {
129            ideal vi=idInit(1,1); vi->m[0]=q;
130            ideal ui=idInit(1,1); ui->m[0]=I->m[i];
131            ideal R; matrix U;
132            ring save_ring=currRing;
133            if (r!=currRing) rChangeCurrRing(r);
134            int save_opt;
135            SI_SAVE_OPT1(save_opt);
136            si_opt_1 &= ~(Sy_bit(OPT_PROT));
137            ideal m = idLift(vi,ui,&R, FALSE,TRUE,TRUE,&U);
138            SI_RESTORE_OPT1(save_opt);
139            if (r!=save_ring) rChangeCurrRing(save_ring);
140            if (idIs0(R))
141            {
142              matrix T = id_Module2formatedMatrix(m,1,1,r);
143              p=MATELEM(T,1,1); MATELEM(T,1,1)=NULL;
144              id_Delete((ideal *)&T,r);
145            }
146            else p=NULL;
147            id_Delete((ideal*)&U,r);
148            id_Delete(&R,r);
149            vi->m[0]=NULL; ui->m[0]=NULL;
150            id_Delete(&vi,r);
151            id_Delete(&ui,r);
152          }
153          p_SetCompP(h,i+1,r);
154          p=p_Add_q(p,h,r);
155        }
156      }
157      id_Delete(&I,r);
158      p_Delete(&q,r);
159      return p;
160    }
161  }
162  else
163  { /* This means that q != 0 consists of just one term, or LetterPlace */
164#ifdef HAVE_RINGS
165    if (pNext(q)!=NULL)
166    {
167      WerrorS("division over a coefficient domain only implemented for terms");
168      return NULL;
169    }
170#endif
171    return p_DivideM(p,q,r);
172  }
173  return NULL;
174}
175
176poly pp_Divide(poly p, poly q, const ring r)
177{
178  if (q==NULL)
179  {
180    WerrorS("div. by 0");
181    return NULL;
182  }
183  if (p==NULL)
184  {
185    return NULL;
186  }
187  if ((pNext(q)!=NULL)||rIsPluralRing(r))
188  { /* This means that q != 0 consists of at least two terms*/
189    if(p_GetComp(p,r)==0)
190    {
191      if((rFieldType(r)==n_transExt)
192      &&(convSingTrP(p,r))
193      &&(convSingTrP(q,r))
194      &&(!rIsNCRing(r)))
195      {
196        poly res=singclap_pdivide(p, q, r);
197        return res;
198      }
199      else if ((r->cf->convSingNFactoryN!=ndConvSingNFactoryN)
200      &&(!rField_is_Ring(r))
201      &&(!rIsNCRing(r)))
202      {
203        poly res=singclap_pdivide(p, q, r);
204        return res;
205      }
206      else
207      {
208        ideal vi=idInit(1,1); vi->m[0]=p_Copy(q,r);
209        ideal ui=idInit(1,1); ui->m[0]=p_Copy(p,r);
210        ideal R; matrix U;
211        ring save_ring=currRing;
212        if (r!=currRing) rChangeCurrRing(r);
213        int save_opt;
214        SI_SAVE_OPT1(save_opt);
215        si_opt_1 &= ~(Sy_bit(OPT_PROT));
216        ideal m = idLift(vi,ui,&R, FALSE,TRUE,TRUE,&U);
217        SI_RESTORE_OPT1(save_opt);
218        if (r!=save_ring) rChangeCurrRing(save_ring);
219        matrix T = id_Module2formatedMatrix(m,1,1,r);
220        p=MATELEM(T,1,1); MATELEM(T,1,1)=NULL;
221        id_Delete((ideal *)&T,r);
222        id_Delete((ideal *)&U,r);
223        id_Delete(&R,r);
224        //vi->m[0]=NULL; ui->m[0]=NULL;
225        id_Delete(&vi,r);
226        id_Delete(&ui,r);
227        return p;
228      }
229    }
230    else
231    {
232      p=p_Copy(p,r);
233      int comps=p_MaxComp(p,r);
234      ideal I=idInit(comps,1);
235      poly h;
236      int i;
237      // conversion to a list of polys:
238      while (p!=NULL)
239      {
240        i=p_GetComp(p,r)-1;
241        h=pNext(p);
242        pNext(p)=NULL;
243        p_SetComp(p,0,r);
244        I->m[i]=p_Add_q(I->m[i],p,r);
245        p=h;
246      }
247      // division and conversion to vector:
248      h=NULL;
249      p=NULL;
250      q=p_Copy(q,r);
251      for(i=comps-1;i>=0;i--)
252      {
253        if (I->m[i]!=NULL)
254        {
255          if((rFieldType(r)==n_transExt)
256          &&(convSingTrP(I->m[i],r))
257          &&(convSingTrP(q,r))
258          &&(!rIsNCRing(r)))
259          {
260            h=singclap_pdivide(I->m[i],q,r);
261          }
262          else if ((r->cf->convSingNFactoryN!=ndConvSingNFactoryN)
263          &&(!rField_is_Ring(r))
264          &&(!rIsNCRing(r)))
265            h=singclap_pdivide(I->m[i],q,r);
266          else
267          {
268            ideal vi=idInit(1,1); vi->m[0]=q;
269            ideal ui=idInit(1,1); ui->m[0]=I->m[i];
270            ideal R; matrix U;
271            ring save_ring=currRing;
272            if (r!=currRing) rChangeCurrRing(r);
273            int save_opt;
274            SI_SAVE_OPT1(save_opt);
275            si_opt_1 &= ~(Sy_bit(OPT_PROT));
276            ideal m = idLift(vi,ui,&R, FALSE,TRUE,TRUE,&U);
277            SI_RESTORE_OPT1(save_opt);
278            if (r!=save_ring) rChangeCurrRing(save_ring);
279            if (idIs0(R))
280            {
281              matrix T = id_Module2formatedMatrix(m,1,1,r);
282              p=MATELEM(T,1,1); MATELEM(T,1,1)=NULL;
283              id_Delete((ideal *)&T,r);
284            }
285            else p=NULL;
286            id_Delete((ideal*)&U,r);
287            id_Delete(&R,r);
288            vi->m[0]=NULL; ui->m[0]=NULL;
289            id_Delete(&vi,r);
290            id_Delete(&ui,r);
291          }
292          p_SetCompP(h,i+1,r);
293          p=p_Add_q(p,h,r);
294        }
295      }
296      id_Delete(&I,r);
297      p_Delete(&q,r);
298      return p;
299    }
300  }
301  else
302  { /* This means that q != 0 consists of just one term,
303       or that r is over a coefficient ring. */
304#ifdef HAVE_RINGS
305    if (pNext(q)!=NULL)
306    {
307      WerrorS("division over a coefficient domain only implemented for terms");
308      return NULL;
309    }
310#endif
311    return pp_DivideM(p,q,r);
312  }
313  return NULL;
314}
315
316poly p_DivRem(poly p, poly q, poly &rest, const ring r)
317{
318  assume(q!=NULL);
319  rest=NULL;
320  if (q==NULL)
321  {
322    WerrorS("div. by 0");
323    return NULL;
324  }
325  if (p==NULL)
326  {
327    p_Delete(&q,r);
328    return NULL;
329  }
330  if(p_GetComp(p,r)==0)
331  {
332    if((rFieldType(r)==n_transExt)
333    &&(convSingTrP(p,r))
334    &&(convSingTrP(q,r))
335    &&(!rIsNCRing(r)))
336    {
337      poly res=singclap_pdivide(p, q, r);
338      rest=singclap_pmod(p,q,r);
339      p_Delete(&p,r);
340      p_Delete(&q,r);
341      return res;
342    }
343    else if ((r->cf->convSingNFactoryN!=ndConvSingNFactoryN)
344    &&(!rField_is_Ring(r))
345    &&(!rIsNCRing(r)))
346    {
347      poly res=singclap_pdivide(p, q, r);
348      rest=singclap_pmod(p,q,r);
349      p_Delete(&p,r);
350      p_Delete(&q,r);
351      return res;
352    }
353    else
354    {
355      ideal vi=idInit(1,1); vi->m[0]=q;
356      ideal ui=idInit(1,1); ui->m[0]=p;
357      ideal R; matrix U;
358      ring save_ring=currRing;
359      if (r!=currRing) rChangeCurrRing(r);
360      int save_opt;
361      SI_SAVE_OPT1(save_opt);
362      si_opt_1 &= ~(Sy_bit(OPT_PROT));
363      ideal m = idLift(vi,ui,&R, FALSE,TRUE,TRUE,&U);
364      SI_RESTORE_OPT1(save_opt);
365      if (r!=save_ring) rChangeCurrRing(save_ring);
366      p=m->m[0]; m->m[0]=NULL;
367      id_Delete(&m,r);
368      p_SetCompP(p,0,r);
369      rest=R->m[0]; R->m[0]=NULL;
370      id_Delete(&R,r);
371      p_SetCompP(rest,0,r);
372      id_Delete((ideal *)&U,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 if (!rField_is_Ring(r))  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 if (!rField_is_Ring(r))  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    if (nCoeff_is_Ring(r->cf)) 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.