My Project
polys.cc
Go to the documentation of this file.
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.
13 VAR ring currRing = NULL;
14 
15 void 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 
31 poly 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  {
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 
174 poly 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);
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  {
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 
314 poly 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 
380 poly 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 }
#define TRUE
Definition: auxiliary.h:100
#define FALSE
Definition: auxiliary.h:96
int m
Definition: cfEzgcd.cc:128
int i
Definition: cfEzgcd.cc:132
int p
Definition: cfModGcd.cc:4078
g
Definition: cfModGcd.cc:4090
FILE * f
Definition: checklibs.c:9
BOOLEAN convSingTrP(poly p, const ring r)
Definition: clapconv.cc:352
poly singclap_pmod(poly f, poly g, const ring r)
Definition: clapsing.cc:702
poly singclap_pdivide(poly f, poly g, const ring r)
Definition: clapsing.cc:624
poly singclap_gcd_r(poly f, poly g, const ring r)
Definition: clapsing.cc:68
Definition: intvec.h:23
@ n_transExt
used for all transcendental extensions, i.e., the top-most extension in an extension tower is transce...
Definition: coeffs.h:38
static FORCE_INLINE void nSetChar(const coeffs r)
initialisations after each ring change
Definition: coeffs.h:436
static FORCE_INLINE BOOLEAN nCoeff_is_Ring(const coeffs r)
Definition: coeffs.h:730
#define WarnS
Definition: emacs.cc:78
CanonicalForm res
Definition: facAbsFact.cc:60
const CanonicalForm & w
Definition: facAbsFact.cc:51
void WerrorS(const char *s)
Definition: feFopen.cc:24
#define VAR
Definition: globaldefs.h:5
ideal idSyzygies(ideal h1, tHomog h, intvec **w, BOOLEAN setSyzComp, BOOLEAN setRegularity, int *deg, GbVariant alg)
Definition: ideals.cc:830
ideal idLift(ideal mod, ideal submod, ideal *rest, BOOLEAN goodShape, BOOLEAN isSB, BOOLEAN divide, matrix *unit, GbVariant alg)
represents the generators of submod in terms of the generators of mod (Matrix(SM)*U-Matrix(rest)) = M...
Definition: ideals.cc:1105
BOOLEAN idIs0(ideal h)
returns true if h is the zero ideal
STATIC_VAR jList * T
Definition: janet.cc:30
STATIC_VAR Poly * h
Definition: janet.cc:971
void p_TakeOutComp(poly *p, long comp, poly *q, int *lq, const ring r)
Definition: p_polys.cc:3574
#define MATELEM(mat, i, j)
1-based access to matrix
Definition: matpol.h:29
#define assume(x)
Definition: mod2.h:387
#define p_GetComp(p, r)
Definition: monomials.h:64
#define pNext(p)
Definition: monomials.h:36
CanonicalForm ndConvSingNFactoryN(number, BOOLEAN, const coeffs)
Definition: numbers.cc:276
#define NULL
Definition: omList.c:12
VAR unsigned si_opt_1
Definition: options.c:5
#define OPT_PROT
Definition: options.h:75
#define SI_SAVE_OPT1(A)
Definition: options.h:21
#define SI_RESTORE_OPT1(A)
Definition: options.h:24
#define Sy_bit(x)
Definition: options.h:31
poly p_DivideM(poly a, poly b, const ring r)
Definition: p_polys.cc:1574
void p_Content(poly ph, const ring r)
Definition: p_polys.cc:2291
poly pp_DivideM(poly a, poly b, const ring r)
Definition: p_polys.cc:1629
void p_Norm(poly p1, const ring r)
Definition: p_polys.cc:3797
poly p_Cleardenom(poly p, const ring r)
Definition: p_polys.cc:2910
poly p_One(const ring r)
Definition: p_polys.cc:1313
static poly p_Add_q(poly p, poly q, const ring r)
Definition: p_polys.h:936
static void p_SetCompP(poly p, int i, ring r)
Definition: p_polys.h:254
static unsigned long p_SetComp(poly p, unsigned long c, ring r)
Definition: p_polys.h:247
static BOOLEAN p_IsConstant(const poly p, const ring r)
Definition: p_polys.h:2011
static long p_MaxComp(poly p, ring lmRing, ring tailRing)
Definition: p_polys.h:292
static void p_Delete(poly *p, const ring r)
Definition: p_polys.h:901
static poly p_Copy(poly p, const ring r)
returns a copy of p
Definition: p_polys.h:846
void rChangeCurrRing(ring r)
Definition: polys.cc:15
poly p_Divide(poly p, poly q, const ring r)
polynomial division a/b, ignoring the rest via singclap_pdivide resp. idLift destroyes a,...
Definition: polys.cc:31
VAR ring currRing
Widely used global variable which specifies the current polynomial ring for Singular interpreter and ...
Definition: polys.cc:13
poly pp_Divide(poly p, poly q, const ring r)
polynomial division a/b, ignoring the rest via singclap_pdivide resp. idLift does not destroy a,...
Definition: polys.cc:174
poly singclap_gcd(poly f, poly g, const ring r)
polynomial gcd via singclap_gcd_r resp. idSyzygies destroys f and g
Definition: polys.cc:380
poly p_DivRem(poly p, poly q, poly &rest, const ring r)
Definition: polys.cc:314
Compatiblity layer for legacy polynomial operations (over currRing)
void p_SetGlobals(const ring r, BOOLEAN complete)
set all properties of a new ring - also called by rComplete
Definition: ring.cc:3457
static BOOLEAN rField_is_Zp(const ring r)
Definition: ring.h:501
static n_coeffType rFieldType(const ring r)
the type of the coefficient filed of r (n_Zp, n_Q, etc)
Definition: ring.h:557
static BOOLEAN rIsPluralRing(const ring r)
we must always have this test!
Definition: ring.h:400
static BOOLEAN rIsNCRing(const ring r)
Definition: ring.h:421
#define rTest(r)
Definition: ring.h:786
#define rField_is_Ring(R)
Definition: ring.h:486
ideal idInit(int idsize, int rank)
initialise an ideal / module
Definition: simpleideals.cc:35
void id_Delete(ideal *h, ring r)
deletes an ideal/module/matrix
matrix id_Module2formatedMatrix(ideal mod, int rows, int cols, const ring R)
#define IDELEMS(i)
Definition: simpleideals.h:23
#define R
Definition: sirandom.c:27
@ testHomog
Definition: structs.h:38