source: git/kernel/polys.cc @ a241bd

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