My Project
flint_mpoly.cc
Go to the documentation of this file.
1 // emacs edit mode for this file is -*- C++ -*-
2 /****************************************
3 * Computer Algebra System SINGULAR *
4 ****************************************/
5 /*
6 * ABSTRACT: flint mpoly
7 */
8 
9 #include "misc/auxiliary.h"
10 #include "flintconv.h"
11 #include "flint_mpoly.h"
12 
13 #ifdef HAVE_FLINT
14 #if __FLINT_RELEASE >= 20503
15 #include "coeffs/coeffs.h"
16 #include "coeffs/longrat.h"
18 
19 #include <vector>
20 
21 /****** ring conversion ******/
22 
23 BOOLEAN convSingRFlintR(fmpq_mpoly_ctx_t ctx, const ring r)
24 {
25  if (rRing_ord_pure_dp(r))
26  {
27  fmpq_mpoly_ctx_init(ctx,r->N,ORD_DEGREVLEX);
28  return FALSE;
29  }
30  else if (rRing_ord_pure_Dp(r))
31  {
32  fmpq_mpoly_ctx_init(ctx,r->N,ORD_DEGLEX);
33  return FALSE;
34  }
35  else if (rRing_ord_pure_lp(r))
36  {
37  fmpq_mpoly_ctx_init(ctx,r->N,ORD_LEX);
38  return FALSE;
39  }
40  return TRUE;
41 }
42 
43 BOOLEAN convSingRFlintR(nmod_mpoly_ctx_t ctx, const ring r)
44 {
45  if (rRing_ord_pure_dp(r))
46  {
47  nmod_mpoly_ctx_init(ctx,r->N,ORD_DEGREVLEX,r->cf->ch);
48  return FALSE;
49  }
50  else if (rRing_ord_pure_Dp(r))
51  {
52  nmod_mpoly_ctx_init(ctx,r->N,ORD_DEGLEX,r->cf->ch);
53  return FALSE;
54  }
55  else if (rRing_ord_pure_lp(r))
56  {
57  nmod_mpoly_ctx_init(ctx,r->N,ORD_LEX,r->cf->ch);
58  return FALSE;
59  }
60  return TRUE;
61 }
62 
63 BOOLEAN convSingRFlintR(fmpz_mpoly_ctx_t ctx, const ring r)
64 {
65  if (rRing_ord_pure_dp(r))
66  {
67  fmpz_mpoly_ctx_init(ctx,r->N,ORD_DEGREVLEX);
68  return FALSE;
69  }
70  else if (rRing_ord_pure_Dp(r))
71  {
72  fmpz_mpoly_ctx_init(ctx,r->N,ORD_DEGLEX);
73  return FALSE;
74  }
75  else if (rRing_ord_pure_lp(r))
76  {
77  fmpz_mpoly_ctx_init(ctx,r->N,ORD_LEX);
78  return FALSE;
79  }
80  return TRUE;
81 }
82 
83 /******** polynomial conversion ***********/
84 
85 // memory allocation is not thread safe; singular polynomials must be constructed in serial
86 
87 /*
88  We agree the that result of a singular -> fmpq_mpoly conversion is
89  readonly. This restricts the usage of the result in flint functions to
90  const arguments. However, the real readonly conversion is currently only
91  implemented in the threaded conversion below since it requires a scan of
92  all coefficients anyways. The _fmpq_mpoly_clear_readonly_sing needs to
93  be provided for a consistent interface in the polynomial operations.
94 */
95 static void _fmpq_mpoly_clear_readonly_sing(fmpq_mpoly_t a, fmpq_mpoly_ctx_t ctx)
96 {
97  fmpq_mpoly_clear(a, ctx);
98 }
99 
100 void convSingPFlintMP(fmpq_mpoly_t res, fmpq_mpoly_ctx_t ctx, poly p, int lp, const ring r)
101 {
102  fmpq_mpoly_init2(res, lp, ctx);
103  ulong* exp=(ulong*)omAlloc((r->N+1)*sizeof(ulong));
104  while(p!=NULL)
105  {
106  number n=pGetCoeff(p);
107  fmpq_t c;
108  convSingNFlintN_QQ(c,n);
109  #if SIZEOF_LONG==8
110  p_GetExpVL(p,(int64*)exp,r);
111  fmpq_mpoly_push_term_fmpq_ui(res, c, exp, ctx);
112  #else
113  p_GetExpV(p,(int*)exp,r);
114  fmpq_mpoly_push_term_fmpq_ui(res, c, &(exp[1]), ctx);
115  #endif
116  fmpq_clear(c);
117  pIter(p);
118  }
119  fmpq_mpoly_reduce(res, ctx); // extra step for QQ ensures res has content canonically factored out
120  omFreeSize(exp,(r->N+1)*sizeof(ulong));
121 }
122 
123 poly convFlintMPSingP(fmpq_mpoly_t f, fmpq_mpoly_ctx_t ctx, const ring r)
124 {
125  int d=fmpq_mpoly_length(f,ctx)-1;
126  poly p=NULL;
127  ulong* exp=(ulong*)omAlloc0((r->N+1)*sizeof(ulong));
128  fmpq_t c;
129  fmpq_init(c);
130  for(int i=d; i>=0; i--)
131  {
132  fmpq_mpoly_get_term_coeff_fmpq(c,f,i,ctx);
133  poly pp=p_Init(r);
134  #if SIZEOF_LONG==8
135  fmpq_mpoly_get_term_exp_ui(exp,f,i,ctx);
136  p_SetExpVL(pp,(int64*)exp,r);
137  #else
138  fmpq_mpoly_get_term_exp_ui(&(exp[1]),f,i,ctx);
139  p_SetExpV(pp,(int*)exp,r);
140  #endif
141  p_Setm(pp,r);
142  number n=convFlintNSingN_QQ(c,r->cf);
143  pSetCoeff0(pp,n);
144  pNext(pp)=p;
145  p=pp;
146  }
147  fmpq_clear(c);
148  omFreeSize(exp,(r->N+1)*sizeof(ulong));
149  p_Test(p,r);
150  return p;
151 }
152 
153 void convSingPFlintMP(fmpz_mpoly_t res, fmpz_mpoly_ctx_t ctx, poly p, int lp, const ring r)
154 {
155  fmpz_mpoly_init2(res, lp, ctx);
156  ulong* exp=(ulong*)omAlloc((r->N+1)*sizeof(ulong));
157  while(p!=NULL)
158  {
159  number n=pGetCoeff(p);
160  fmpz_t c;
161  convSingNFlintN(c,n);
162  #if SIZEOF_LONG==8
163  p_GetExpVL(p,(int64*)exp,r);
164  fmpz_mpoly_push_term_fmpz_ui(res, c, exp, ctx);
165  #else
166  p_GetExpV(p,(int*)exp,r);
167  fmpz_mpoly_push_term_fmpz_ui(res, c, &(exp[1]), ctx);
168  #endif
169  fmpz_clear(c);
170  pIter(p);
171  }
172  omFreeSize(exp,(r->N+1)*sizeof(ulong));
173 }
174 
175 poly convFlintMPSingP(fmpz_mpoly_t f, fmpz_mpoly_ctx_t ctx, const ring r)
176 {
177  int d=fmpz_mpoly_length(f,ctx)-1;
178  poly p=NULL;
179  ulong* exp=(ulong*)omAlloc0((r->N+1)*sizeof(ulong));
180  fmpz_t c;
181  fmpz_init(c);
182  for(int i=d; i>=0; i--)
183  {
184  fmpz_mpoly_get_term_coeff_fmpz(c,f,i,ctx);
185  poly pp=p_Init(r);
186  #if SIZEOF_LONG==8
187  fmpz_mpoly_get_term_exp_ui(exp,f,i,ctx);
188  p_SetExpVL(pp,(int64*)exp,r);
189  #else
190  fmpz_mpoly_get_term_exp_ui(&(exp[1]),f,i,ctx);
191  p_SetExpV(pp,(int*)exp,r);
192  #endif
193  p_Setm(pp,r);
194  number n=convFlintNSingN(c,r->cf);
195  pSetCoeff0(pp,n);
196  pNext(pp)=p;
197  p=pp;
198  }
199  fmpz_clear(c);
200  omFreeSize(exp,(r->N+1)*sizeof(ulong));
201  p_Test(p,r);
202  return p;
203 }
204 
205 poly convFlintMPSingP(nmod_mpoly_t f, nmod_mpoly_ctx_t ctx, const ring r)
206 {
207  int d=nmod_mpoly_length(f,ctx)-1;
208  poly p=NULL;
209  ulong* exp=(ulong*)omAlloc0((r->N+1)*sizeof(ulong));
210  for(int i=d; i>=0; i--)
211  {
212  ulong c=nmod_mpoly_get_term_coeff_ui(f,i,ctx);
213  poly pp=p_Init(r);
214  #if SIZEOF_LONG==8
215  nmod_mpoly_get_term_exp_ui(exp,f,i,ctx);
216  p_SetExpVL(pp,(int64*)exp,r);
217  #else
218  nmod_mpoly_get_term_exp_ui(&(exp[1]),f,i,ctx);
219  p_SetExpV(pp,(int*)exp,r);
220  #endif
221  p_Setm(pp,r);
222  pSetCoeff0(pp,(number)c);
223  pNext(pp)=p;
224  p=pp;
225  }
226  omFreeSize(exp,(r->N+1)*sizeof(ulong));
227  p_Test(p,r);
228  return p;
229 }
230 
231 void convSingPFlintMP(nmod_mpoly_t res, nmod_mpoly_ctx_t ctx, poly p, int lp,const ring r)
232 {
233  nmod_mpoly_init2(res, lp, ctx);
234  ulong* exp=(ulong*)omAlloc((r->N+1)*sizeof(ulong));
235  while(p!=NULL)
236  {
237  number n=pGetCoeff(p);
238  #if SIZEOF_LONG==8
239  p_GetExpVL(p,(int64*)exp,r);
240  nmod_mpoly_push_term_ui_ui(res, (ulong)n, exp, ctx);
241  #else
242  p_GetExpV(p,(int*)exp,r);
243  nmod_mpoly_push_term_ui_ui(res, (ulong)n, &(exp[1]), ctx);
244  #endif
245  pIter(p);
246  }
247  omFreeSize(exp,(r->N+1)*sizeof(ulong));
248 }
249 
250 /****** polynomial operations ***********/
251 
252 poly Flint_Mult_MP(poly p,int lp, poly q, int lq, fmpq_mpoly_ctx_t ctx, const ring r)
253 {
254  fmpq_mpoly_t pp,qq,res;
255  convSingPFlintMP(pp,ctx,p,lp,r); // pp read only
256  convSingPFlintMP(qq,ctx,q,lq,r); // qq read only
257  fmpq_mpoly_init(res,ctx);
258  fmpq_mpoly_mul(res,pp,qq,ctx);
259  poly pres=convFlintMPSingP(res,ctx,r);
260  fmpq_mpoly_clear(res,ctx);
261  _fmpq_mpoly_clear_readonly_sing(pp,ctx);
262  _fmpq_mpoly_clear_readonly_sing(qq,ctx);
263  fmpq_mpoly_ctx_clear(ctx);
264  p_Test(pres,r);
265  return pres;
266 }
267 
268 poly Flint_Mult_MP(poly p,int lp, poly q, int lq, nmod_mpoly_ctx_t ctx, const ring r)
269 {
270  nmod_mpoly_t pp,qq,res;
271  convSingPFlintMP(pp,ctx,p,lp,r);
272  convSingPFlintMP(qq,ctx,q,lq,r);
273  nmod_mpoly_init(res,ctx);
274  nmod_mpoly_mul(res,pp,qq,ctx);
275  poly pres=convFlintMPSingP(res,ctx,r);
276  nmod_mpoly_clear(res,ctx);
277  nmod_mpoly_clear(pp,ctx);
278  nmod_mpoly_clear(qq,ctx);
279  nmod_mpoly_ctx_clear(ctx);
280  p_Test(pres,r);
281  return pres;
282 }
283 
284 poly Flint_Mult_MP(poly p,int lp, poly q, int lq, fmpz_mpoly_ctx_t ctx, const ring r)
285 {
286  fmpz_mpoly_t pp,qq,res;
287  convSingPFlintMP(pp,ctx,p,lp,r); // pp read only
288  convSingPFlintMP(qq,ctx,q,lq,r); // qq read only
289  fmpz_mpoly_init(res,ctx);
290  fmpz_mpoly_mul(res,pp,qq,ctx);
291  poly pres=convFlintMPSingP(res,ctx,r);
292  fmpz_mpoly_clear(res,ctx);
293  fmpz_mpoly_clear(pp,ctx);
294  fmpz_mpoly_clear(qq,ctx);
295  fmpz_mpoly_ctx_clear(ctx);
296  p_Test(pres,r);
297  return pres;
298 }
299 
300 // Zero will be returned if the division is not exact
301 poly Flint_Divide_MP(poly p,int lp, poly q, int lq, fmpq_mpoly_ctx_t ctx, const ring r)
302 {
303  fmpq_mpoly_t pp,qq,res;
304  convSingPFlintMP(pp,ctx,p,lp,r); // pp read only
305  convSingPFlintMP(qq,ctx,q,lq,r); // qq read only
306  fmpq_mpoly_init(res,ctx);
307  fmpq_mpoly_divides(res,pp,qq,ctx);
308  poly pres = convFlintMPSingP(res,ctx,r);
309  fmpq_mpoly_clear(res,ctx);
310  _fmpq_mpoly_clear_readonly_sing(pp,ctx);
311  _fmpq_mpoly_clear_readonly_sing(qq,ctx);
312  fmpq_mpoly_ctx_clear(ctx);
313  p_Test(pres,r);
314  return pres;
315 }
316 
317 poly Flint_Divide_MP(poly p,int lp, poly q, int lq, nmod_mpoly_ctx_t ctx, const ring r)
318 {
319  nmod_mpoly_t pp,qq,res;
320  convSingPFlintMP(pp,ctx,p,lp,r);
321  convSingPFlintMP(qq,ctx,q,lq,r);
322  nmod_mpoly_init(res,ctx);
323  nmod_mpoly_divides(res,pp,qq,ctx);
324  poly pres=convFlintMPSingP(res,ctx,r);
325  nmod_mpoly_clear(res,ctx);
326  nmod_mpoly_clear(pp,ctx);
327  nmod_mpoly_clear(qq,ctx);
328  nmod_mpoly_ctx_clear(ctx);
329  p_Test(pres,r);
330  return pres;
331 }
332 
333 poly Flint_GCD_MP(poly p,int lp,poly q,int lq,nmod_mpoly_ctx_t ctx,const ring r)
334 {
335  nmod_mpoly_t pp,qq,res;
336  convSingPFlintMP(pp,ctx,p,lp,r);
337  convSingPFlintMP(qq,ctx,q,lq,r);
338  nmod_mpoly_init(res,ctx);
339  int ok=nmod_mpoly_gcd(res,pp,qq,ctx);
340  poly pres;
341  if (ok)
342  {
343  pres=convFlintMPSingP(res,ctx,r);
344  p_Test(pres,r);
345  }
346  else
347  {
348  pres=p_One(r);
349  }
350  nmod_mpoly_clear(res,ctx);
351  nmod_mpoly_clear(pp,ctx);
352  nmod_mpoly_clear(qq,ctx);
353  nmod_mpoly_ctx_clear(ctx);
354  return pres;
355 }
356 
357 poly Flint_GCD_MP(poly p,int lp,poly q,int lq,fmpq_mpoly_ctx_t ctx,const ring r)
358 {
359  fmpq_mpoly_t pp,qq,res;
360  convSingPFlintMP(pp,ctx,p,lp,r); // pp read only
361  convSingPFlintMP(qq,ctx,q,lq,r); // qq read only
362  fmpq_mpoly_init(res,ctx);
363  int ok=fmpq_mpoly_gcd(res,pp,qq,ctx);
364  poly pres;
365  if (ok)
366  {
367  // Flint normalizes the gcd to be monic.
368  // Singular wants a gcd defined over ZZ that is primitive and has a positive leading coeff.
369  if (!fmpq_mpoly_is_zero(res, ctx))
370  {
371  fmpq_t content;
372  fmpq_init(content);
373  fmpq_mpoly_content(content, res, ctx);
374  fmpq_mpoly_scalar_div_fmpq(res, res, content, ctx);
375  fmpq_clear(content);
376  }
377  pres=convFlintMPSingP(res,ctx,r);
378  p_Test(pres,r);
379  }
380  else
381  {
382  pres=p_One(r);
383  }
384  fmpq_mpoly_clear(res,ctx);
385  _fmpq_mpoly_clear_readonly_sing(pp,ctx);
386  _fmpq_mpoly_clear_readonly_sing(qq,ctx);
387  fmpq_mpoly_ctx_clear(ctx);
388  return pres;
389 }
390 
391 #endif
392 #endif
All the auxiliary stuff.
long int64
Definition: auxiliary.h:68
int BOOLEAN
Definition: auxiliary.h:87
#define TRUE
Definition: auxiliary.h:100
#define FALSE
Definition: auxiliary.h:96
CanonicalForm pp(const CanonicalForm &)
CanonicalForm pp ( const CanonicalForm & f )
Definition: cf_gcd.cc:676
CanonicalForm content(const CanonicalForm &)
CanonicalForm content ( const CanonicalForm & f )
Definition: cf_gcd.cc:603
int i
Definition: cfEzgcd.cc:132
int p
Definition: cfModGcd.cc:4080
FILE * f
Definition: checklibs.c:9
Coefficient rings, fields and other domains suitable for Singular polynomials.
CanonicalForm res
Definition: facAbsFact.cc:60
This file is work in progress and currently not part of the official Singular.
void convSingNFlintN(fmpz_t f, mpz_t z)
void convSingNFlintN_QQ(fmpq_t f, number n)
void convFlintNSingN(mpz_t z, fmpz_t f)
number convFlintNSingN_QQ(fmpq_t f, const coeffs cf)
#define pIter(p)
Definition: monomials.h:37
#define pNext(p)
Definition: monomials.h:36
static number & pGetCoeff(poly p)
return an alias to the leading coefficient of p assumes that p != NULL NOTE: not copy
Definition: monomials.h:44
#define pSetCoeff0(p, n)
Definition: monomials.h:59
gmp_float exp(const gmp_float &a)
Definition: mpr_complex.cc:357
Definition: lq.h:40
#define omFreeSize(addr, size)
Definition: omAllocDecl.h:260
#define omAlloc(size)
Definition: omAllocDecl.h:210
#define omAlloc0(size)
Definition: omAllocDecl.h:211
#define NULL
Definition: omList.c:12
poly p_One(const ring r)
Definition: p_polys.cc:1308
static void p_SetExpVL(poly p, int64 *ev, const ring r)
Definition: p_polys.h:1513
static void p_SetExpV(poly p, int *ev, const ring r)
Definition: p_polys.h:1504
static void p_Setm(poly p, const ring r)
Definition: p_polys.h:233
static void p_GetExpVL(poly p, int64 *ev, const ring r)
Definition: p_polys.h:1489
static void p_GetExpV(poly p, int *ev, const ring r)
Definition: p_polys.h:1480
static poly p_Init(const ring r, omBin bin)
Definition: p_polys.h:1280
#define p_Test(p, r)
Definition: p_polys.h:162
BOOLEAN rRing_ord_pure_Dp(const ring r)
Definition: ring.cc:5174
BOOLEAN rRing_ord_pure_dp(const ring r)
Definition: ring.cc:5164
BOOLEAN rRing_ord_pure_lp(const ring r)
Definition: ring.cc:5184