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