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);
81  p=MATELEM(T,1,1); MATELEM(T,1,1)=NULL;
82  id_Delete((ideal *)&T,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  assume(q!=NULL);
177  if (q==NULL)
178  {
179  WerrorS("div. by 0");
180  return NULL;
181  }
182  if (p==NULL)
183  {
184  return NULL;
185  }
186  if ((pNext(q)!=NULL)||rIsPluralRing(r))
187  { /* This means that q != 0 consists of at least two terms*/
188  if(p_GetComp(p,r)==0)
189  {
190  if((rFieldType(r)==n_transExt)
191  &&(convSingTrP(p,r))
192  &&(convSingTrP(q,r))
193  &&(!rIsNCRing(r)))
194  {
195  poly res=singclap_pdivide(p, q, r);
196  return res;
197  }
198  else if ((r->cf->convSingNFactoryN!=ndConvSingNFactoryN)
199  &&(!rField_is_Ring(r))
200  &&(!rIsNCRing(r)))
201  {
202  poly res=singclap_pdivide(p, q, r);
203  return res;
204  }
205  else
206  {
207  ideal vi=idInit(1,1); vi->m[0]=p_Copy(q,r);
208  ideal ui=idInit(1,1); ui->m[0]=p_Copy(p,r);
209  ideal R; matrix U;
210  ring save_ring=currRing;
211  if (r!=currRing) rChangeCurrRing(r);
212  int save_opt;
213  SI_SAVE_OPT1(save_opt);
214  si_opt_1 &= ~(Sy_bit(OPT_PROT));
215  ideal m = idLift(vi,ui,&R, FALSE,TRUE,TRUE,&U);
216  SI_RESTORE_OPT1(save_opt);
217  if (r!=save_ring) rChangeCurrRing(save_ring);
219  p=MATELEM(T,1,1); MATELEM(T,1,1)=NULL;
220  id_Delete((ideal *)&T,r);
221  id_Delete((ideal *)&U,r);
222  id_Delete(&R,r);
223  //vi->m[0]=NULL; ui->m[0]=NULL;
224  id_Delete(&vi,r);
225  id_Delete(&ui,r);
226  return p;
227  }
228  }
229  else
230  {
231  p=p_Copy(p,r);
232  int comps=p_MaxComp(p,r);
233  ideal I=idInit(comps,1);
234  poly h;
235  int i;
236  // conversion to a list of polys:
237  while (p!=NULL)
238  {
239  i=p_GetComp(p,r)-1;
240  h=pNext(p);
241  pNext(p)=NULL;
242  p_SetComp(p,0,r);
243  I->m[i]=p_Add_q(I->m[i],p,r);
244  p=h;
245  }
246  // division and conversion to vector:
247  h=NULL;
248  p=NULL;
249  q=p_Copy(q,r);
250  for(i=comps-1;i>=0;i--)
251  {
252  if (I->m[i]!=NULL)
253  {
254  if((rFieldType(r)==n_transExt)
255  &&(convSingTrP(I->m[i],r))
256  &&(convSingTrP(q,r))
257  &&(!rIsNCRing(r)))
258  {
259  h=singclap_pdivide(I->m[i],q,r);
260  }
261  else if ((r->cf->convSingNFactoryN!=ndConvSingNFactoryN)
262  &&(!rField_is_Ring(r))
263  &&(!rIsNCRing(r)))
264  h=singclap_pdivide(I->m[i],q,r);
265  else
266  {
267  ideal vi=idInit(1,1); vi->m[0]=q;
268  ideal ui=idInit(1,1); ui->m[0]=I->m[i];
269  ideal R; matrix U;
270  ring save_ring=currRing;
271  if (r!=currRing) rChangeCurrRing(r);
272  int save_opt;
273  SI_SAVE_OPT1(save_opt);
274  si_opt_1 &= ~(Sy_bit(OPT_PROT));
275  ideal m = idLift(vi,ui,&R, FALSE,TRUE,TRUE,&U);
276  SI_RESTORE_OPT1(save_opt);
277  if (r!=save_ring) rChangeCurrRing(save_ring);
278  if (idIs0(R))
279  {
281  p=MATELEM(T,1,1); MATELEM(T,1,1)=NULL;
282  id_Delete((ideal *)&T,r);
283  }
284  else p=NULL;
285  id_Delete((ideal*)&U,r);
286  id_Delete(&R,r);
287  vi->m[0]=NULL; ui->m[0]=NULL;
288  id_Delete(&vi,r);
289  id_Delete(&ui,r);
290  }
291  p_SetCompP(h,i+1,r);
292  p=p_Add_q(p,h,r);
293  }
294  }
295  id_Delete(&I,r);
296  p_Delete(&q,r);
297  return p;
298  }
299  }
300  else
301  { /* This means that q != 0 consists of just one term,
302  or that r is over a coefficient ring. */
303 #ifdef HAVE_RINGS
304  if (pNext(q)!=NULL)
305  {
306  WerrorS("division over a coefficient domain only implemented for terms");
307  return NULL;
308  }
309 #endif
310  return pp_DivideM(p,q,r);
311  }
312  return NULL;
313 }
314 
315 poly p_DivRem(poly p, poly q, poly &rest, const ring r)
316 {
317  assume(q!=NULL);
318  rest=NULL;
319  if (q==NULL)
320  {
321  WerrorS("div. by 0");
322  return NULL;
323  }
324  if (p==NULL)
325  {
326  p_Delete(&q,r);
327  return NULL;
328  }
329  if(p_GetComp(p,r)==0)
330  {
331  if((rFieldType(r)==n_transExt)
332  &&(convSingTrP(p,r))
333  &&(convSingTrP(q,r))
334  &&(!rIsNCRing(r)))
335  {
336  poly res=singclap_pdivide(p, q, r);
337  rest=singclap_pmod(p,q,r);
338  p_Delete(&p,r);
339  p_Delete(&q,r);
340  return res;
341  }
342  else if ((r->cf->convSingNFactoryN!=ndConvSingNFactoryN)
343  &&(!rField_is_Ring(r))
344  &&(!rIsNCRing(r)))
345  {
346  poly res=singclap_pdivide(p, q, r);
347  rest=singclap_pmod(p,q,r);
348  p_Delete(&p,r);
349  p_Delete(&q,r);
350  return res;
351  }
352  else
353  {
354  ideal vi=idInit(1,1); vi->m[0]=q;
355  ideal ui=idInit(1,1); ui->m[0]=p;
356  ideal R; matrix U;
357  ring save_ring=currRing;
358  if (r!=currRing) rChangeCurrRing(r);
359  int save_opt;
360  SI_SAVE_OPT1(save_opt);
361  si_opt_1 &= ~(Sy_bit(OPT_PROT));
362  ideal m = idLift(vi,ui,&R, FALSE,TRUE,TRUE,&U);
363  SI_RESTORE_OPT1(save_opt);
364  if (r!=save_ring) rChangeCurrRing(save_ring);
366  p=MATELEM(T,1,1); MATELEM(T,1,1)=NULL;
367  id_Delete((ideal *)&T,r);
368  T = id_Module2formatedMatrix(R,1,1,r);
369  rest=MATELEM(T,1,1); MATELEM(T,1,1)=NULL;
370  id_Delete((ideal *)&T,r);
371  id_Delete((ideal *)&U,r);
372  id_Delete(&R,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 
382 poly 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 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 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  p_Content(res,r);
436  return res;
437  }
438  p_Delete(&f, r);
439  p_Delete(&g, r);
440  return res;
441 }
#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:4080
g
Definition: cfModGcd.cc:4092
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:668
poly singclap_pdivide(poly f, poly g, const ring r)
Definition: clapsing.cc:590
poly singclap_gcd_r(poly f, poly g, const ring r)
Definition: clapsing.cc:45
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:39
static FORCE_INLINE void nSetChar(const coeffs r)
initialisations after each ring change
Definition: coeffs.h:437
#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)
Definition: ideals.cc:1099
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:3565
#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:275
#define NULL
Definition: omList.c:12
VAR unsigned si_opt_1
Definition: options.c:5
#define OPT_PROT
Definition: options.h:74
#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:1565
void p_Content(poly ph, const ring r)
Definition: p_polys.cc:2282
poly pp_DivideM(poly a, poly b, const ring r)
Definition: p_polys.cc:1620
void p_Norm(poly p1, const ring r)
Definition: p_polys.cc:3789
poly p_Cleardenom(poly p, const ring r)
Definition: p_polys.cc:2900
poly p_One(const ring r)
Definition: p_polys.cc:1308
static poly p_Add_q(poly p, poly q, const ring r)
Definition: p_polys.h:896
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:1971
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:861
static poly p_Copy(poly p, const ring r)
returns a copy of p
Definition: p_polys.h:812
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:382
poly p_DivRem(poly p, poly q, poly &rest, const ring r)
Definition: polys.cc:315
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:3368
static BOOLEAN rField_is_Ring(const ring r)
Definition: ring.h:489
static BOOLEAN rField_is_Zp(const ring r)
Definition: ring.h:505
static n_coeffType rFieldType(const ring r)
the type of the coefficient filed of r (n_Zp, n_Q, etc)
Definition: ring.h:561
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:790
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:43