source: git/kernel/polys.cc @ c858487

spielwiese
Last change on this file since c858487 was c858487, checked in by Frédéric Chapoton <chapoton@…>, 18 months ago
fix many typos in kernel/
  • Property mode set to 100644
File size: 11.3 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  //------------ 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        p=m->m[0]; m->m[0]=NULL;
81        id_Delete(&m,r);
82        p_SetCompP(p,0,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  if (q==NULL)
177  {
178    WerrorS("div. by 0");
179    return NULL;
180  }
181  if (p==NULL)
182  {
183    return NULL;
184  }
185  if ((pNext(q)!=NULL)||rIsPluralRing(r))
186  { /* This means that q != 0 consists of at least two terms*/
187    if(p_GetComp(p,r)==0)
188    {
189      if((rFieldType(r)==n_transExt)
190      &&(convSingTrP(p,r))
191      &&(convSingTrP(q,r))
192      &&(!rIsNCRing(r)))
193      {
194        poly res=singclap_pdivide(p, q, r);
195        return res;
196      }
197      else if ((r->cf->convSingNFactoryN!=ndConvSingNFactoryN)
198      &&(!rField_is_Ring(r))
199      &&(!rIsNCRing(r)))
200      {
201        poly res=singclap_pdivide(p, q, r);
202        return res;
203      }
204      else
205      {
206        ideal vi=idInit(1,1); vi->m[0]=p_Copy(q,r);
207        ideal ui=idInit(1,1); ui->m[0]=p_Copy(p,r);
208        ideal R; matrix U;
209        ring save_ring=currRing;
210        if (r!=currRing) rChangeCurrRing(r);
211        int save_opt;
212        SI_SAVE_OPT1(save_opt);
213        si_opt_1 &= ~(Sy_bit(OPT_PROT));
214        ideal m = idLift(vi,ui,&R, FALSE,TRUE,TRUE,&U);
215        SI_RESTORE_OPT1(save_opt);
216        if (r!=save_ring) rChangeCurrRing(save_ring);
217        matrix T = id_Module2formatedMatrix(m,1,1,r);
218        p=MATELEM(T,1,1); MATELEM(T,1,1)=NULL;
219        id_Delete((ideal *)&T,r);
220        id_Delete((ideal *)&U,r);
221        id_Delete(&R,r);
222        //vi->m[0]=NULL; ui->m[0]=NULL;
223        id_Delete(&vi,r);
224        id_Delete(&ui,r);
225        return p;
226      }
227    }
228    else
229    {
230      p=p_Copy(p,r);
231      int comps=p_MaxComp(p,r);
232      ideal I=idInit(comps,1);
233      poly h;
234      int i;
235      // conversion to a list of polys:
236      while (p!=NULL)
237      {
238        i=p_GetComp(p,r)-1;
239        h=pNext(p);
240        pNext(p)=NULL;
241        p_SetComp(p,0,r);
242        I->m[i]=p_Add_q(I->m[i],p,r);
243        p=h;
244      }
245      // division and conversion to vector:
246      h=NULL;
247      p=NULL;
248      q=p_Copy(q,r);
249      for(i=comps-1;i>=0;i--)
250      {
251        if (I->m[i]!=NULL)
252        {
253          if((rFieldType(r)==n_transExt)
254          &&(convSingTrP(I->m[i],r))
255          &&(convSingTrP(q,r))
256          &&(!rIsNCRing(r)))
257          {
258            h=singclap_pdivide(I->m[i],q,r);
259          }
260          else if ((r->cf->convSingNFactoryN!=ndConvSingNFactoryN)
261          &&(!rField_is_Ring(r))
262          &&(!rIsNCRing(r)))
263            h=singclap_pdivide(I->m[i],q,r);
264          else
265          {
266            ideal vi=idInit(1,1); vi->m[0]=q;
267            ideal ui=idInit(1,1); ui->m[0]=I->m[i];
268            ideal R; matrix U;
269            ring save_ring=currRing;
270            if (r!=currRing) rChangeCurrRing(r);
271            int save_opt;
272            SI_SAVE_OPT1(save_opt);
273            si_opt_1 &= ~(Sy_bit(OPT_PROT));
274            ideal m = idLift(vi,ui,&R, FALSE,TRUE,TRUE,&U);
275            SI_RESTORE_OPT1(save_opt);
276            if (r!=save_ring) rChangeCurrRing(save_ring);
277            if (idIs0(R))
278            {
279              matrix T = id_Module2formatedMatrix(m,1,1,r);
280              p=MATELEM(T,1,1); MATELEM(T,1,1)=NULL;
281              id_Delete((ideal *)&T,r);
282            }
283            else p=NULL;
284            id_Delete((ideal*)&U,r);
285            id_Delete(&R,r);
286            vi->m[0]=NULL; ui->m[0]=NULL;
287            id_Delete(&vi,r);
288            id_Delete(&ui,r);
289          }
290          p_SetCompP(h,i+1,r);
291          p=p_Add_q(p,h,r);
292        }
293      }
294      id_Delete(&I,r);
295      p_Delete(&q,r);
296      return p;
297    }
298  }
299  else
300  { /* This means that q != 0 consists of just one term,
301       or that r is over a coefficient ring. */
302#ifdef HAVE_RINGS
303    if (pNext(q)!=NULL)
304    {
305      WerrorS("division over a coefficient domain only implemented for terms");
306      return NULL;
307    }
308#endif
309    return pp_DivideM(p,q,r);
310  }
311  return NULL;
312}
313
314poly p_DivRem(poly p, poly q, poly &rest, const ring r)
315{
316  assume(q!=NULL);
317  rest=NULL;
318  if (q==NULL)
319  {
320    WerrorS("div. by 0");
321    return NULL;
322  }
323  if (p==NULL)
324  {
325    p_Delete(&q,r);
326    return NULL;
327  }
328  if(p_GetComp(p,r)==0)
329  {
330    if((rFieldType(r)==n_transExt)
331    &&(convSingTrP(p,r))
332    &&(convSingTrP(q,r))
333    &&(!rIsNCRing(r)))
334    {
335      poly res=singclap_pdivide(p, q, r);
336      rest=singclap_pmod(p,q,r);
337      p_Delete(&p,r);
338      p_Delete(&q,r);
339      return res;
340    }
341    else if ((r->cf->convSingNFactoryN!=ndConvSingNFactoryN)
342    &&(!rField_is_Ring(r))
343    &&(!rIsNCRing(r)))
344    {
345      poly res=singclap_pdivide(p, q, r);
346      rest=singclap_pmod(p,q,r);
347      p_Delete(&p,r);
348      p_Delete(&q,r);
349      return res;
350    }
351    else
352    {
353      ideal vi=idInit(1,1); vi->m[0]=q;
354      ideal ui=idInit(1,1); ui->m[0]=p;
355      ideal R; matrix U;
356      ring save_ring=currRing;
357      if (r!=currRing) rChangeCurrRing(r);
358      int save_opt;
359      SI_SAVE_OPT1(save_opt);
360      si_opt_1 &= ~(Sy_bit(OPT_PROT));
361      ideal m = idLift(vi,ui,&R, FALSE,TRUE,TRUE,&U);
362      SI_RESTORE_OPT1(save_opt);
363      if (r!=save_ring) rChangeCurrRing(save_ring);
364      p=m->m[0]; m->m[0]=NULL;
365      id_Delete(&m,r);
366      p_SetCompP(p,0,r);
367      rest=R->m[0]; R->m[0]=NULL;
368      id_Delete(&R,r);
369      p_SetCompP(rest,0,r);
370      id_Delete((ideal *)&U,r);
371      //vi->m[0]=NULL; ui->m[0]=NULL;
372      id_Delete(&vi,r);
373      id_Delete(&ui,r);
374      return p;
375    }
376  }
377  return NULL;
378}
379
380poly singclap_gcd ( poly f, poly g, const ring r )
381{
382  poly res=NULL;
383
384  if (f!=NULL)
385  {
386    //if (r->cf->has_simple_Inverse) p_Norm(f,r);
387    if (rField_is_Zp(r))          p_Norm(f,r);
388    else if (!rField_is_Ring(r))  p_Cleardenom(f, r);
389  }
390  if (g!=NULL)
391  {
392    //if (r->cf->has_simple_Inverse) p_Norm(g,r);
393    if (rField_is_Zp(r))          p_Norm(g,r);
394    else if (!rField_is_Ring(r))  p_Cleardenom(g, r);
395  }
396  else         return f; // g==0 => gcd=f (but do a p_Cleardenom/pNorm)
397  if (f==NULL) return g; // f==0 => gcd=g (but do a p_Cleardenom/pNorm)
398  if(!rField_is_Ring(r)
399  && (p_IsConstant(f,r)
400  ||p_IsConstant(g,r)))
401  {
402    res=p_One(r);
403  }
404  else if (r->cf->convSingNFactoryN!=ndConvSingNFactoryN)
405  {
406    res=singclap_gcd_r(f,g,r);
407  }
408  else
409  {
410    ideal I=idInit(2,1);
411    I->m[0]=f;
412    I->m[1]=p_Copy(g,r);
413    intvec *w=NULL;
414    ring save_ring=currRing;
415    if (r!=currRing) rChangeCurrRing(r);
416    int save_opt;
417    SI_SAVE_OPT1(save_opt);
418    si_opt_1 &= ~(Sy_bit(OPT_PROT));
419    ideal S1=idSyzygies(I,testHomog,&w);
420    if (w!=NULL) delete w;
421    // expect S1->m[0]=(-g/gcd,f/gcd)
422    if (IDELEMS(S1)!=1) WarnS("error in syzygy computation for GCD");
423    int lp;
424    p_TakeOutComp(&S1->m[0],1,&res,&lp,r);
425    p_Delete(&S1->m[0],r);
426    // GCD is g divided iby (-g/gcd):
427    res=p_Divide(g,res,r);
428    // restore, r, opt:
429    SI_RESTORE_OPT1(save_opt);
430    if (r!=save_ring) rChangeCurrRing(save_ring);
431    // clean the result
432    res=p_Cleardenom(res,r);
433    if (nCoeff_is_Ring(r->cf)) p_Content(res,r);
434    return res;
435  }
436  p_Delete(&f, r);
437  p_Delete(&g, r);
438  return res;
439}
Note: See TracBrowser for help on using the repository browser.