source: git/kernel/polys.cc @ cec9624

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