source: git/factory/FLINTconvert.cc @ cc5a86

spielwiese
Last change on this file since cc5a86 was cc5a86, checked in by Hans Schoenemann <hannes@…>, 3 years ago
fix: order of args.
  • Property mode set to 100644
File size: 22.1 KB
Line 
1/*****************************************************************************\
2 * Computer Algebra System SINGULAR
3\*****************************************************************************/
4/** @file FLINTconvert.cc
5 *
6 * This file implements functions for conversion to FLINT (www.flintlib.org)
7 * and back.
8 *
9 * @author Martin Lee
10 *
11 **/
12/*****************************************************************************/
13
14
15
16#include <config.h>
17
18
19#include "canonicalform.h"
20#include "fac_util.h"
21#include "cf_iter.h"
22#include "cf_factory.h"
23#include "gmpext.h"
24#include "singext.h"
25#include "cf_algorithm.h"
26
27#ifdef HAVE_OMALLOC
28#define Alloc(L) omAlloc(L)
29#define Free(A,L) omFreeSize(A,L)
30#else
31#define Alloc(L) malloc(L)
32#define Free(A,L) free(A)
33#endif
34
35#ifdef HAVE_FLINT
36#ifdef HAVE_CSTDIO
37#include <cstdio>
38#else
39#include <stdio.h>
40#endif
41#ifdef __cplusplus
42extern "C"
43{
44#endif
45#ifndef __GMP_BITS_PER_MP_LIMB
46#define __GMP_BITS_PER_MP_LIMB GMP_LIMB_BITS
47#endif
48#include <flint/fmpz.h>
49#include <flint/fmpq.h>
50#include <flint/fmpz_poly.h>
51#include <flint/fmpz_mod_poly.h>
52#include <flint/nmod_poly.h>
53#include <flint/fmpq_poly.h>
54#include <flint/nmod_mat.h>
55#include <flint/fmpz_mat.h>
56#if ( __FLINT_RELEASE >= 20400)
57#include <flint/fq.h>
58#include <flint/fq_poly.h>
59#include <flint/fq_nmod.h>
60#include <flint/fq_nmod_poly.h>
61#include <flint/fq_nmod_mat.h>
62#endif
63#if ( __FLINT_RELEASE >= 20503)
64#include <flint/fmpq_mpoly.h>
65
66#if (__FLINT_RELEASE < 20601)
67// helper for fq_nmod_t -> nmod_poly_t
68static void fq_nmod_get_nmod_poly(nmod_poly_t a, const fq_nmod_t b, const fq_nmod_ctx_t ctx)
69{
70    FLINT_ASSERT(b->mod.n == ctx->modulus->mod.n);
71    a->mod = ctx->modulus->mod;
72    nmod_poly_set(a, b);
73}
74
75// helper for nmod_poly_t -> fq_nmod_t
76static void fq_nmod_set_nmod_poly(fq_nmod_t a, const nmod_poly_t b, const fq_nmod_ctx_t ctx)
77{
78    FLINT_ASSERT(a->mod.n == b->mod.n);
79    FLINT_ASSERT(a->mod.n == ctx->modulus->mod.n);
80    nmod_poly_set(a, b);
81    fq_nmod_reduce(a, ctx);
82}
83#endif
84
85
86#endif
87#ifdef __cplusplus
88}
89#endif
90
91#include "FLINTconvert.h"
92
93void convertCF2Fmpz (fmpz_t result, const CanonicalForm& f)
94{
95  if (f.isImm())
96    fmpz_set_si (result, f.intval());
97  else
98  {
99    mpz_t gmp_val;
100    f.mpzval(gmp_val);
101    fmpz_set_mpz (result, gmp_val);
102    mpz_clear (gmp_val);
103  }
104}
105
106void convertFacCF2Fmpz_poly_t (fmpz_poly_t result, const CanonicalForm& f)
107{
108  fmpz_poly_init2 (result, degree (f)+1);
109  _fmpz_poly_set_length(result, degree(f)+1);
110  for (CFIterator i= f; i.hasTerms(); i++)
111    convertCF2Fmpz (fmpz_poly_get_coeff_ptr(result, i.exp()), i.coeff());
112}
113
114CanonicalForm convertFmpz2CF (const fmpz_t coefficient)
115{
116  if (fmpz_cmp_si (coefficient, MINIMMEDIATE) >= 0 &&
117      fmpz_cmp_si (coefficient, MAXIMMEDIATE) <= 0)
118  {
119    long coeff= fmpz_get_si (coefficient);
120    return CanonicalForm (coeff);
121  }
122  else
123  {
124    mpz_t gmp_val;
125    mpz_init (gmp_val);
126    fmpz_get_mpz (gmp_val, coefficient);
127    CanonicalForm result= CanonicalForm (CFFactory::basic (gmp_val));
128    return result;
129  }
130}
131
132CanonicalForm
133convertFmpz_poly_t2FacCF (const fmpz_poly_t poly, const Variable& x)
134{
135  CanonicalForm result= 0;
136  fmpz* coeff;
137  for (int i= 0; i < fmpz_poly_length (poly); i++)
138  {
139    coeff= fmpz_poly_get_coeff_ptr (poly, i);
140    if (!fmpz_is_zero (coeff))
141      result += convertFmpz2CF (coeff)*power (x,i);
142  }
143  return result;
144}
145
146void
147convertFacCF2nmod_poly_t (nmod_poly_t result, const CanonicalForm& f)
148{
149  bool save_sym_ff= isOn (SW_SYMMETRIC_FF);
150  if (save_sym_ff) Off (SW_SYMMETRIC_FF);
151  nmod_poly_init2 (result, getCharacteristic(), degree (f)+1);
152  for (CFIterator i= f; i.hasTerms(); i++)
153  {
154    CanonicalForm c= i.coeff();
155    if (!c.isImm()) c=c.mapinto(); //c%= getCharacteristic();
156    if (!c.isImm())
157    {  //This case will never happen if the characteristic is in fact a prime
158       // number, since all coefficients are represented as immediates
159       printf("convertCF2nmod_poly_t: coefficient not immediate!, char=%d\n",
160              getCharacteristic());
161    }
162    else
163      nmod_poly_set_coeff_ui (result, i.exp(), c.intval());
164  }
165  if (save_sym_ff) On (SW_SYMMETRIC_FF);
166}
167
168CanonicalForm
169convertnmod_poly_t2FacCF (const nmod_poly_t poly, const Variable& x)
170{
171  CanonicalForm result= 0;
172  for (int i= 0; i < nmod_poly_length (poly); i++)
173  {
174    ulong coeff= nmod_poly_get_coeff_ui (poly, i);
175    if (coeff != 0)
176      result += CanonicalForm ((long)coeff)*power (x,i);
177  }
178  return result;
179}
180
181void convertCF2Fmpq (fmpq_t result, const CanonicalForm& f)
182{
183  //ASSERT (isOn (SW_RATIONAL), "expected rational");
184  if (f.isImm ())
185  {
186    fmpq_set_si (result, f.intval(), 1);
187  }
188  else if(f.inQ())
189  {
190    mpz_t gmp_val;
191    gmp_numerator (f, gmp_val);
192    fmpz_set_mpz (fmpq_numref (result), gmp_val);
193    mpz_clear (gmp_val);
194    gmp_denominator (f, gmp_val);
195    fmpz_set_mpz (fmpq_denref (result), gmp_val);
196    mpz_clear (gmp_val);
197  }
198  else if(f.inZ())
199  {
200    mpz_t gmp_val;
201    f.mpzval(gmp_val);
202    fmpz_set_mpz (fmpq_numref (result), gmp_val);
203    mpz_clear (gmp_val);
204    fmpz_one(fmpq_denref(result));
205  }
206  else
207  {
208    printf("wrong type\n");
209  }
210}
211
212CanonicalForm convertFmpq_t2CF (const fmpq_t q)
213{
214  bool isRat= isOn (SW_RATIONAL);
215  if (!isRat)
216    On (SW_RATIONAL);
217
218  CanonicalForm num, den;
219  mpz_t nnum, nden;
220  mpz_init (nnum);
221  mpz_init (nden);
222  fmpz_get_mpz (nnum, fmpq_numref (q));
223  fmpz_get_mpz (nden, fmpq_denref (q));
224
225  CanonicalForm result;
226  if (mpz_is_imm (nden))
227  {
228    if (mpz_is_imm(nnum))
229    {
230      num= CanonicalForm (mpz_get_si(nnum));
231      den= CanonicalForm (mpz_get_si(nden));
232      mpz_clear (nnum);
233      mpz_clear (nden);
234      result= num/den;
235    }
236    else if (mpz_cmp_si(nden,1)==0)
237    {
238      result= CanonicalForm( CFFactory::basic(nnum));
239      mpz_clear (nden);
240    }
241    else
242      result= CanonicalForm( CFFactory::rational( nnum, nden, false));
243  }
244  else
245  {
246    result= CanonicalForm( CFFactory::rational( nnum, nden, false));
247  }
248  if (!isRat)
249    Off (SW_RATIONAL);
250  return result;
251}
252
253CanonicalForm
254convertFmpq_poly_t2FacCF (const fmpq_poly_t p, const Variable& x)
255{
256  CanonicalForm result= 0;
257  fmpq_t coeff;
258  long n= p->length;
259  for (long i= 0; i < n; i++)
260  {
261    fmpq_init (coeff);
262    fmpq_poly_get_coeff_fmpq (coeff, p, i);
263    if (fmpq_is_zero (coeff))
264    {
265      fmpq_clear (coeff);
266      continue;
267    }
268    result += convertFmpq_t2CF (coeff)*power (x, i);
269    fmpq_clear (coeff);
270  }
271  return result;
272}
273
274void convertFacCF2Fmpz_array (fmpz* result, const CanonicalForm& f)
275{
276  for (CFIterator i= f; i.hasTerms(); i++)
277    convertCF2Fmpz (&result[i.exp()], i.coeff());
278}
279
280void convertFacCF2Fmpq_poly_t (fmpq_poly_t result, const CanonicalForm& f)
281{
282  bool isRat= isOn (SW_RATIONAL);
283  if (!isRat)
284    On (SW_RATIONAL);
285
286  fmpq_poly_init2 (result, degree (f)+1);
287  _fmpq_poly_set_length (result, degree (f) + 1);
288  CanonicalForm den= bCommonDen (f);
289  convertFacCF2Fmpz_array (fmpq_poly_numref (result), f*den);
290  convertCF2Fmpz (fmpq_poly_denref (result), den);
291
292  if (!isRat)
293    Off (SW_RATIONAL);
294}
295
296CFFList
297convertFLINTnmod_poly_factor2FacCFFList (const nmod_poly_factor_t fac,
298                                          const mp_limb_t leadingCoeff,
299                                          const Variable& x
300                                         )
301{
302  CFFList result;
303  if (leadingCoeff != 1)
304    result.insert (CFFactor (CanonicalForm ((long) leadingCoeff), 1));
305
306  long i;
307
308  for (i = 0; i < fac->num; i++)
309    result.append (CFFactor (convertnmod_poly_t2FacCF (
310                             (nmod_poly_t &)fac->p[i],x),
311                             fac->exp[i]));
312  return result;
313}
314
315#if __FLINT_RELEASE >= 20503
316CFFList
317convertFLINTfmpz_poly_factor2FacCFFList (
318                   const fmpz_poly_factor_t fac, ///< [in] a nmod_poly_factor_t
319                   const Variable& x       ///< [in] variable the result should
320                                           ///< have
321                                        )
322
323{
324  CFFList result;
325  long i;
326
327  result.append (CFFactor(convertFmpz2CF(&fac->c),1));
328
329  for (i = 0; i < fac->num; i++)
330    result.append (CFFactor (convertFmpz_poly_t2FacCF (
331                             (fmpz_poly_t &)fac->p[i],x),
332                             fac->exp[i]));
333  return result;
334}
335#endif
336
337#if __FLINT_RELEASE >= 20400
338CFFList
339convertFLINTFq_nmod_poly_factor2FacCFFList (const fq_nmod_poly_factor_t fac,
340                                       const Variable& x, const Variable& alpha,
341                                       const fq_nmod_ctx_t fq_con
342                                         )
343{
344  CFFList result;
345
346  long i;
347
348  for (i = 0; i < fac->num; i++)
349    result.append (CFFactor (convertFq_nmod_poly_t2FacCF (
350                             (fq_nmod_poly_t &)fac->poly[i], x, alpha, fq_con),
351                             fac->exp[i]));
352  return result;
353}
354#endif
355
356void
357convertFacCF2Fmpz_mod_poly_t (fmpz_mod_poly_t result, const CanonicalForm& f,
358                              const fmpz_t p)
359{
360  fmpz_mod_poly_init2 (result, p, degree (f) + 1);
361  fmpz_poly_t buf;
362  convertFacCF2Fmpz_poly_t (buf, f);
363  fmpz_mod_poly_set_fmpz_poly (result, buf);
364  fmpz_poly_clear (buf);
365}
366
367CanonicalForm
368convertFmpz_mod_poly_t2FacCF (const fmpz_mod_poly_t poly, const Variable& x,
369                              const modpk& b)
370{
371  fmpz_poly_t buf;
372  fmpz_poly_init (buf);
373  fmpz_mod_poly_get_fmpz_poly (buf, poly);
374  CanonicalForm result= convertFmpz_poly_t2FacCF (buf, x);
375  fmpz_poly_clear (buf);
376  return b (result);
377}
378
379#if __FLINT_RELEASE >= 20400
380void
381convertFacCF2Fq_nmod_t (fq_nmod_t result, const CanonicalForm& f,
382                        const fq_nmod_ctx_t ctx)
383{
384  bool save_sym_ff= isOn (SW_SYMMETRIC_FF);
385  if (save_sym_ff) Off (SW_SYMMETRIC_FF);
386  nmod_poly_t res;
387  nmod_poly_init(res,getCharacteristic());
388  for (CFIterator i= f; i.hasTerms(); i++)
389  {
390    CanonicalForm c= i.coeff();
391    if (!c.isImm()) c=c.mapinto(); //c%= getCharacteristic();
392    if (!c.isImm())
393    {  //This case will never happen if the characteristic is in fact a prime
394       // number, since all coefficients are represented as immediates
395       printf("convertFacCF2Fq_nmod_t: coefficient not immediate!, char=%d\n",
396              getCharacteristic());
397    }
398    else
399    {
400      STICKYASSERT (i.exp() <= fq_nmod_ctx_degree(ctx), "convertFacCF2Fq_nmod_t: element is not reduced");
401      nmod_poly_set_coeff_ui (res, i.exp(), c.intval());
402    }
403  }
404  fq_nmod_init(result,ctx);
405  fq_nmod_set_nmod_poly(result,res,ctx);
406  if (save_sym_ff) On (SW_SYMMETRIC_FF);
407}
408
409CanonicalForm
410convertFq_nmod_t2FacCF (const fq_nmod_t poly, const Variable& alpha, const fq_nmod_ctx_t ctx)
411{
412  return convertnmod_poly_t2FacCF (poly, alpha);
413}
414
415void
416convertFacCF2Fq_t (fq_t result, const CanonicalForm& f, const fq_ctx_t ctx)
417{
418  fmpz_poly_init2 (result, fq_ctx_degree(ctx));
419  ASSERT (degree (f) < fq_ctx_degree (ctx), "input is not reduced");
420  _fmpz_poly_set_length(result, degree(f)+1);
421  for (CFIterator i= f; i.hasTerms(); i++)
422    convertCF2Fmpz (fmpz_poly_get_coeff_ptr(result, i.exp()), i.coeff());
423  _fmpz_vec_scalar_mod_fmpz (result->coeffs, result->coeffs, degree (f) + 1,
424                             &ctx->p);
425  _fmpz_poly_normalise (result);
426}
427
428CanonicalForm
429convertFq_t2FacCF (const fq_t poly, const Variable& alpha)
430{
431  return convertFmpz_poly_t2FacCF (poly, alpha);
432}
433
434void
435convertFacCF2Fq_poly_t (fq_poly_t result, const CanonicalForm& f,
436                        const fq_ctx_t ctx)
437{
438  fq_poly_init2 (result, degree (f)+1, ctx);
439  _fq_poly_set_length (result, degree (f) + 1, ctx);
440  fmpz_poly_t buf;
441  for (CFIterator i= f; i.hasTerms(); i++)
442  {
443    convertFacCF2Fmpz_poly_t (buf, i.coeff());
444    _fmpz_vec_scalar_mod_fmpz (buf->coeffs, buf->coeffs, degree (i.coeff()) + 1,
445                               &ctx->p);
446    _fmpz_poly_normalise (buf);
447    fq_poly_set_coeff (result, i.exp(), buf, ctx);
448    fmpz_poly_clear (buf);
449  }
450}
451
452void
453convertFacCF2Fq_nmod_poly_t (fq_nmod_poly_t result, const CanonicalForm& f,
454                             const fq_nmod_ctx_t ctx)
455{
456  fq_nmod_poly_init2 (result, degree (f)+1, ctx);
457  _fq_nmod_poly_set_length (result, degree (f) + 1, ctx);
458  fq_nmod_t buf;
459  fq_nmod_init2 (buf, ctx);
460  for (CFIterator i= f; i.hasTerms(); i++)
461  {
462    convertFacCF2Fq_nmod_t (buf, i.coeff(), ctx);
463    fq_nmod_poly_set_coeff (result, i.exp(), buf, ctx);
464    fq_nmod_zero (buf, ctx);
465  }
466  fq_nmod_clear (buf, ctx);
467}
468
469CanonicalForm
470convertFq_poly_t2FacCF (const fq_poly_t p, const Variable& x,
471                        const Variable& alpha, const fq_ctx_t ctx)
472{
473  CanonicalForm result= 0;
474  fq_t coeff;
475  long n= fq_poly_length (p, ctx);
476  fq_init2 (coeff, ctx);
477  for (long i= 0; i < n; i++)
478  {
479    fq_poly_get_coeff (coeff, p, i, ctx);
480    if (fq_is_zero (coeff, ctx))
481      continue;
482    result += convertFq_t2FacCF (coeff, alpha)*power (x, i);
483    fq_zero (coeff, ctx);
484  }
485  fq_clear (coeff, ctx);
486
487  return result;
488}
489
490CanonicalForm
491convertFq_nmod_poly_t2FacCF (const fq_nmod_poly_t p, const Variable& x,
492                             const Variable& alpha, const fq_nmod_ctx_t ctx)
493{
494  CanonicalForm result= 0;
495  fq_nmod_t coeff;
496  long n= fq_nmod_poly_length (p, ctx);
497  fq_nmod_init2 (coeff, ctx);
498  for (long i= 0; i < n; i++)
499  {
500    fq_nmod_poly_get_coeff (coeff, p, i, ctx);
501    if (fq_nmod_is_zero (coeff, ctx))
502      continue;
503    result += convertFq_nmod_t2FacCF (coeff, alpha, ctx)*power (x, i);
504    fq_nmod_zero (coeff, ctx);
505  }
506  fq_nmod_clear (coeff, ctx);
507
508  return result;
509}
510#endif
511
512void convertFacCFMatrix2Fmpz_mat_t (fmpz_mat_t M, const CFMatrix &m)
513{
514  fmpz_mat_init (M, (long) m.rows(), (long) m.columns());
515
516  int i,j;
517  for(i=m.rows();i>0;i--)
518  {
519    for(j=m.columns();j>0;j--)
520    {
521      convertCF2Fmpz (fmpz_mat_entry (M,i-1,j-1), m(i,j));
522    }
523  }
524}
525CFMatrix* convertFmpz_mat_t2FacCFMatrix(const fmpz_mat_t m)
526{
527  CFMatrix *res=new CFMatrix(fmpz_mat_nrows (m),fmpz_mat_ncols (m));
528  int i,j;
529  for(i=res->rows();i>0;i--)
530  {
531    for(j=res->columns();j>0;j--)
532    {
533      (*res)(i,j)=convertFmpz2CF(fmpz_mat_entry (m,i-1,j-1));
534    }
535  }
536  return res;
537}
538
539void convertFacCFMatrix2nmod_mat_t (nmod_mat_t M, const CFMatrix &m)
540{
541  nmod_mat_init (M, (long) m.rows(), (long) m.columns(), getCharacteristic());
542
543  bool save_sym_ff= isOn (SW_SYMMETRIC_FF);
544  if (save_sym_ff) Off (SW_SYMMETRIC_FF);
545  int i,j;
546  for(i=m.rows();i>0;i--)
547  {
548    for(j=m.columns();j>0;j--)
549    {
550      if(!(m(i,j)).isImm()) printf("convertFacCFMatrix2FLINTmat_zz_p: not imm.\n");
551      nmod_mat_entry (M,i-1,j-1)= (m(i,j)).intval();
552    }
553  }
554  if (save_sym_ff) On (SW_SYMMETRIC_FF);
555}
556
557CFMatrix* convertNmod_mat_t2FacCFMatrix(const nmod_mat_t m)
558{
559  CFMatrix *res=new CFMatrix(nmod_mat_nrows (m), nmod_mat_ncols (m));
560  int i,j;
561  for(i=res->rows();i>0;i--)
562  {
563    for(j=res->columns();j>0;j--)
564    {
565      (*res)(i,j)=CanonicalForm((long) nmod_mat_entry (m, i-1, j-1));
566    }
567  }
568  return res;
569}
570
571#if __FLINT_RELEASE >= 20400
572void
573convertFacCFMatrix2Fq_nmod_mat_t (fq_nmod_mat_t M,
574                                  const fq_nmod_ctx_t fq_con, const CFMatrix &m)
575{
576  fq_nmod_mat_init (M, (long) m.rows(), (long) m.columns(), fq_con);
577  int i,j;
578  for(i=m.rows();i>0;i--)
579  {
580    for(j=m.columns();j>0;j--)
581    {
582      convertFacCF2nmod_poly_t (M->rows[i-1]+j-1, m (i,j));
583    }
584  }
585}
586
587CFMatrix*
588convertFq_nmod_mat_t2FacCFMatrix(const fq_nmod_mat_t m,
589                                 const fq_nmod_ctx_t& fq_con,
590                                 const Variable& alpha)
591{
592  CFMatrix *res=new CFMatrix(fq_nmod_mat_nrows (m, fq_con),
593                             fq_nmod_mat_ncols (m, fq_con));
594  int i,j;
595  for(i=res->rows();i>0;i--)
596  {
597    for(j=res->columns();j>0;j--)
598    {
599      (*res)(i,j)=convertFq_nmod_t2FacCF (fq_nmod_mat_entry (m, i-1, j-1),
600                                          alpha, fq_con);
601    }
602  }
603  return res;
604}
605#endif
606#if __FLINT_RELEASE >= 20503
607static void convFlint_RecPP ( const CanonicalForm & f, ulong * exp, nmod_mpoly_t result, nmod_mpoly_ctx_t ctx, int N )
608{
609  // assume f!=0
610  if ( ! f.inCoeffDomain() )
611  {
612    int l = f.level();
613    for ( CFIterator i = f; i.hasTerms(); i++ )
614    {
615      exp[N-l] = i.exp();
616      convFlint_RecPP( i.coeff(), exp, result, ctx, N );
617    }
618    exp[N-l] = 0;
619  }
620  else
621  {
622    int c=f.intval(); // with Off(SW_SYMMETRIC_FF): 0<=c<p
623    nmod_mpoly_push_term_ui_ui(result,c,exp,ctx);
624  }
625}
626
627static void convFlint_RecPP ( const CanonicalForm & f, ulong * exp, fmpq_mpoly_t result, fmpq_mpoly_ctx_t ctx, int N )
628{
629  // assume f!=0
630  if ( ! f.inBaseDomain() )
631  {
632    int l = f.level();
633    for ( CFIterator i = f; i.hasTerms(); i++ )
634    {
635      exp[N-l] = i.exp();
636      convFlint_RecPP( i.coeff(), exp, result, ctx, N );
637    }
638    exp[N-l] = 0;
639  }
640  else
641  {
642    fmpq_t c;
643    fmpq_init(c);
644    convertCF2Fmpq(c,f);
645    fmpq_mpoly_push_term_fmpq_ui(result,c,exp,ctx);
646    fmpq_clear(c);
647  }
648}
649
650void convFactoryPFlintMP ( const CanonicalForm & f, nmod_mpoly_t res, nmod_mpoly_ctx_t ctx, int N )
651{
652  if (f.isZero()) return;
653  ulong * exp = (ulong*)Alloc(N*sizeof(ulong));
654  memset(exp,0,N*sizeof(ulong));
655  bool save_sym_ff= isOn (SW_SYMMETRIC_FF);
656  if (save_sym_ff) Off (SW_SYMMETRIC_FF);
657  convFlint_RecPP( f, exp, res, ctx, N );
658  if (save_sym_ff) On(SW_SYMMETRIC_FF);
659  Free(exp,N*sizeof(ulong));
660}
661
662void convFactoryPFlintMP ( const CanonicalForm & f, fmpq_mpoly_t res, fmpq_mpoly_ctx_t ctx, int N )
663{
664  if (f.isZero()) return;
665  ulong * exp = (ulong*)Alloc(N*sizeof(ulong));
666  memset(exp,0,N*sizeof(ulong));
667  convFlint_RecPP( f, exp, res, ctx, N );
668  fmpq_mpoly_reduce(res,ctx);
669  Free(exp,N*sizeof(ulong));
670}
671
672CanonicalForm convFlintMPFactoryP(nmod_mpoly_t f, nmod_mpoly_ctx_t ctx, int N)
673{
674  CanonicalForm result;
675  int d=nmod_mpoly_length(f,ctx)-1;
676  ulong* exp=(ulong*)Alloc(N*sizeof(ulong));
677  for(int i=d; i>=0; i--)
678  {
679    ulong c=nmod_mpoly_get_term_coeff_ui(f,i,ctx);
680    nmod_mpoly_get_term_exp_ui(exp,f,i,ctx);
681    CanonicalForm term=(int)c;
682    for ( int i = 0; i <N; i++ )
683    {
684      if (exp[i]!=0) term*=CanonicalForm( Variable( N-i ), exp[i] );
685    }
686    result+=term;
687  }
688  Free(exp,N*sizeof(ulong));
689  return result;
690}
691
692CanonicalForm convFlintMPFactoryP(fmpq_mpoly_t f, fmpq_mpoly_ctx_t ctx, int N)
693{
694  CanonicalForm result;
695  int d=fmpq_mpoly_length(f,ctx)-1;
696  ulong* exp=(ulong*)Alloc(N*sizeof(ulong));
697  fmpq_t c;
698  fmpq_init(c);
699  for(int i=d; i>=0; i--)
700  {
701    fmpq_mpoly_get_term_coeff_fmpq(c,f,i,ctx);
702    fmpq_mpoly_get_term_exp_ui(exp,f,i,ctx);
703    CanonicalForm term=convertFmpq_t2CF(c);
704    for ( int i = 0; i <N; i++ )
705    {
706      if (exp[i]!=0) term*=CanonicalForm( Variable( N-i ), exp[i] );
707    }
708    result+=term;
709  }
710  fmpq_clear(c);
711  Free(exp,N*sizeof(ulong));
712  return result;
713}
714
715// stolen from:
716// https://graphics.stanford.edu/~seander/bithacks.html#IntegerLog
717static inline int SI_LOG2(int v)
718{
719  const unsigned int b[] = {0x2, 0xC, 0xF0, 0xFF00, 0xFFFF0000};
720  const unsigned int S[] = {1, 2, 4, 8, 16};
721
722  unsigned int r = 0; // result of log2(v) will go here
723  if (v & b[4]) { v >>= S[4]; r |= S[4]; }
724  if (v & b[3]) { v >>= S[3]; r |= S[3]; }
725  if (v & b[2]) { v >>= S[2]; r |= S[2]; }
726  if (v & b[1]) { v >>= S[1]; r |= S[1]; }
727  if (v & b[0]) { v >>= S[0]; r |= S[0]; }
728  return (int)r;
729}
730
731CanonicalForm mulFlintMP_Zp(const CanonicalForm& F,int lF, const CanonicalForm& G, int lG,int m)
732{
733  int bits=SI_LOG2(m)+1;
734  int N=F.level();
735  nmod_mpoly_ctx_t ctx;
736  nmod_mpoly_ctx_init(ctx,N,ORD_LEX,getCharacteristic());
737  nmod_mpoly_t f,g,res;
738  nmod_mpoly_init3(f,lF,bits,ctx);
739  nmod_mpoly_init3(g,lG,bits,ctx);
740  convFactoryPFlintMP(F,f,ctx,N);
741  convFactoryPFlintMP(G,g,ctx,N);
742  nmod_mpoly_init(res,ctx);
743  nmod_mpoly_mul(res,f,g,ctx);
744  nmod_mpoly_clear(g,ctx);
745  nmod_mpoly_clear(f,ctx);
746  CanonicalForm RES=convFlintMPFactoryP(res,ctx,N);
747  nmod_mpoly_clear(res,ctx);
748  nmod_mpoly_ctx_clear(ctx);
749  return RES;
750}
751
752CanonicalForm mulFlintMP_QQ(const CanonicalForm& F,int lF, const CanonicalForm& G, int lG, int m)
753{
754  int bits=SI_LOG2(m)+1;
755  int N=F.level();
756  fmpq_mpoly_ctx_t ctx;
757  fmpq_mpoly_ctx_init(ctx,N,ORD_LEX);
758  fmpq_mpoly_t f,g,res;
759  fmpq_mpoly_init3(f,lF,bits,ctx);
760  fmpq_mpoly_init3(g,lG,bits,ctx);
761  convFactoryPFlintMP(F,f,ctx,N);
762  convFactoryPFlintMP(G,g,ctx,N);
763  fmpq_mpoly_init(res,ctx);
764  fmpq_mpoly_mul(res,f,g,ctx);
765  fmpq_mpoly_clear(g,ctx);
766  fmpq_mpoly_clear(f,ctx);
767  CanonicalForm RES=convFlintMPFactoryP(res,ctx,N);
768  fmpq_mpoly_clear(res,ctx);
769  fmpq_mpoly_ctx_clear(ctx);
770  return RES;
771}
772
773CanonicalForm gcdFlintMP_Zp(const CanonicalForm& F, const CanonicalForm& G)
774{
775  int N=F.level();
776  int lf,lg,m=1<<MPOLY_MIN_BITS;
777  lf=size_maxexp(F,m);
778  lg=size_maxexp(G,m);
779  int bits=SI_LOG2(m)+1;
780  nmod_mpoly_ctx_t ctx;
781  nmod_mpoly_ctx_init(ctx,N,ORD_LEX,getCharacteristic());
782  nmod_mpoly_t f,g,res;
783  nmod_mpoly_init3(f,lf,bits,ctx);
784  nmod_mpoly_init3(g,lg,bits,ctx);
785  convFactoryPFlintMP(F,f,ctx,N);
786  convFactoryPFlintMP(G,g,ctx,N);
787  nmod_mpoly_init(res,ctx);
788  int ok=nmod_mpoly_gcd(res,f,g,ctx);
789  nmod_mpoly_clear(g,ctx);
790  nmod_mpoly_clear(f,ctx);
791  CanonicalForm RES=1;
792  if (ok)
793  {
794    RES=convFlintMPFactoryP(res,ctx,N);
795  }
796  nmod_mpoly_clear(res,ctx);
797  nmod_mpoly_ctx_clear(ctx);
798  return RES;
799}
800
801static CanonicalForm b_content ( const CanonicalForm & f )
802{
803    if ( f.inCoeffDomain() )
804        return f;
805    else
806    {
807        CanonicalForm result = 0;
808        CFIterator i;
809        for ( i = f; i.hasTerms() && (!result.isOne()); i++ )
810            result=bgcd( b_content(i.coeff()) , result );
811        return result;
812    }
813}
814
815
816CanonicalForm gcdFlintMP_QQ(const CanonicalForm& F, const CanonicalForm& G)
817{
818  int N=F.level();
819  fmpq_mpoly_ctx_t ctx;
820  fmpq_mpoly_ctx_init(ctx,N,ORD_LEX);
821  fmpq_mpoly_t f,g,res;
822  fmpq_mpoly_init(f,ctx);
823  fmpq_mpoly_init(g,ctx);
824  convFactoryPFlintMP(F,f,ctx,N);
825  convFactoryPFlintMP(G,g,ctx,N);
826  fmpq_mpoly_init(res,ctx);
827  int ok=fmpq_mpoly_gcd(res,f,g,ctx);
828  fmpq_mpoly_clear(g,ctx);
829  fmpq_mpoly_clear(f,ctx);
830  CanonicalForm RES=1;
831  if (ok)
832  {
833    // Flint normalizes the gcd to be monic.
834    // Singular wants a gcd defined over ZZ that is primitive and has a positive leading coeff.
835    if (!fmpq_mpoly_is_zero(res, ctx))
836    {
837      fmpq_t content;
838      fmpq_init(content);
839      fmpq_mpoly_content(content, res, ctx);
840      fmpq_mpoly_scalar_div_fmpq(res, res, content, ctx);
841      fmpq_clear(content);
842    }
843    RES=convFlintMPFactoryP(res,ctx,N);
844    // gcd(2x,4x) should be 2x, so RES should also have the gcd(lc(F),lc(G))
845    RES*=bgcd(b_content(F),b_content(G));
846  }
847  fmpq_mpoly_clear(res,ctx);
848  fmpq_mpoly_ctx_clear(ctx);
849  return RES;
850}
851
852#endif
853
854#endif
855
856
Note: See TracBrowser for help on using the repository browser.