Changeset 1e351b3 in git


Ignore:
Timestamp:
Jul 8, 2019, 2:07:08 PM (5 years ago)
Author:
Hans Schoenemann <hannes@…>
Branches:
(u'spielwiese', '5b153614cbc72bfa198d75b1e9e33dab2645d9fe')
Children:
11cc5e9e0459f60e3740564f9982a3180ad57c07
Parents:
bc4eb794179fe80f037ef650a09e4a41d316208e84f0d81502771d37e63ee8041e613ac62d48668d
git-author:
Hans Schoenemann <hannes@mathematik.uni-kl.de>2019-07-08 14:07:08+02:00
git-committer:
GitHub <noreply@github.com>2019-07-08 14:07:08+02:00
Message:
Merge pull request #938 from tthsqe12/parallel_flint_conversion

Parallel flint conversion
Files:
7 edited

Legend:

Unmodified
Added
Removed
  • Singular/LIB/algebra.lib

    r84f0d8 r1e351b3  
    278278  //with the ring k[y(0),...,y(n)] is computed, the result is ker
    279279   execute ("ring r1= ("+charstr(basering)+"),(x(1..n),y(0..z)),lp;");
    280  //  execute ("ring r1= ("+charstr(basering)+"),(y(0..z),x(1..n)),dp;");
    281280  if (mp!="0")
    282281  { execute ("minpoly=number("+mp+");"); }
  • Singular/LIB/gitfan.lib

    r84f0d8 r1e351b3  
    266266  newVars = newVars+string(var(k));
    267267  newWeightVector[n]=weightVector[k];
    268   execute("ring ringForSaturation = ("+charstr(origin)+"),("+newVars+"),(wp(newWeightVector));");
    269 
     268  ring ringForSaturation = create_ring(ringlist(origin)[1], "("+newVars+")",
     269      list(list("wp", newWeightVector)),"no_minpoly");
    270270  ideal I = satstd(imap(origin,I));
    271271  I = simplify(I,2+4+32);
  • Singular/feOpt.cc

    rbc4eb7 r1e351b3  
    2323#include "fehelp.h"
    2424
     25#ifdef HAVE_FLINT
     26#include <flint/flint.h>
     27#endif
    2528
    2629const char SHORT_OPTS_STRING[] = "bdhpqstvxec:r:u:";
     
    309312      }
    310313
     314      #ifdef HAVE_FLINT
     315      #if __FLINT_RELEASE >= 20503
     316      case FE_OPT_FLINT_THREADS:
     317      {
     318        slong nthreads = (slong)feOptSpec[FE_OPT_FLINT_THREADS].value;
     319        nthreads = FLINT_MAX(nthreads, WORD(1));
     320        flint_set_num_threads(nthreads);
     321        int * cpu_affinities = new int[nthreads];
     322        for (slong i = 0; i < nthreads; i++)
     323          cpu_affinities[i] = (int)i;
     324        flint_set_thread_affinity(cpu_affinities, nthreads);
     325        delete[] cpu_affinities;
     326        return NULL;
     327      }
     328      #endif
     329      #endif
     330
    311331      default:
    312332        return NULL;
  • Singular/feOptTab.h

    rbc4eb7 r1e351b3  
    33
    44#define LONG_OPTION_RETURN 13
     5
     6#ifdef HAVE_FLINT
     7#include <flint/flint.h>
     8#endif
    59
    610// Define here which cmd-line options are recognized
     
    144148   "#threads", "maximal number of CPUs to use for threads",            feOptInt,    (void*)2,      0},
    145149
     150#ifdef HAVE_FLINT
     151#if __FLINT_RELEASE >= 20503
     152  {"flint-threads",   required_argument,    LONG_OPTION_RETURN,
     153   "#flintthreads", "maximal number of threads to use in flint library", feOptInt,    (void*)1,      0},
     154#endif
     155#endif
     156
    146157
    147158  {"MPport",           required_argument,   LONG_OPTION_RETURN,
  • libpolys/polys/flint_mpoly.cc

    rbc4eb7 r1e351b3  
    1717#include "polys/monomials/p_polys.h"
    1818
     19#include <vector>
     20
     21/****** ring conversion ******/
     22
    1923BOOLEAN convSingRFlintR(fmpq_mpoly_ctx_t ctx, const ring r)
    2024{
     
    3741}
    3842
     43BOOLEAN 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/******** polynomial conversion ***********/
     64
     65#if 1
     66// malloc is not thread safe; singular polynomials must be constructed in serial
     67
    3968void convSingPFlintMP(fmpq_mpoly_t res, fmpq_mpoly_ctx_t ctx, poly p, int lp, const ring r)
    4069{
    41   int bits=SI_LOG2(r->bitmask);
    42   fmpq_mpoly_init3(res,lp,bits,ctx);
    43   fmpq_mpoly_resize(res,lp,ctx);
     70  fmpq_mpoly_init2(res, lp, ctx);
    4471  ulong* exp=(ulong*)omAlloc((r->N+1)*sizeof(ulong));
    45   int i=0;
    4672  while(p!=NULL)
    4773  {
     
    5177    #if SIZEOF_LONG==8
    5278    p_GetExpVL(p,(int64*)exp,r);
    53     fmpq_mpoly_set_term_exp_ui(res,i,exp,ctx);
     79    fmpq_mpoly_push_term_fmpq_ui(res, c, exp, ctx);
    5480    #else
    5581    p_GetExpV(p,(int*)exp,r);
    56     fmpq_mpoly_set_term_exp_ui(res,i,&(exp[1]),ctx);
     82    fmpq_mpoly_push_term_fmpq_ui(res, c, &(exp[1]), ctx);
    5783    #endif
    58     fmpq_mpoly_set_term_coeff_fmpq(res,i,c,ctx);
    5984    fmpq_clear(c);
    60     i++;
    6185    pIter(p);
    6286  }
    63   //fmpq_mpoly_print_pretty(res,r->names,ctx);PrintLn();
     87  fmpq_mpoly_reduce(res, ctx); // extra step for QQ ensures res has content canonically factored out
    6488  omFreeSize(exp,(r->N+1)*sizeof(ulong));
    6589}
     
    7094  poly p=NULL;
    7195  ulong* exp=(ulong*)omAlloc0((r->N+1)*sizeof(ulong));
     96  fmpq_t c;
     97  fmpq_init(c);
    7298  for(int i=d; i>=0; i--)
    7399  {
    74     fmpq_t c;
    75100    fmpq_mpoly_get_term_coeff_fmpq(c,f,i,ctx);
    76101    poly pp=p_Init(r);
     
    84109    p_Setm(pp,r);
    85110    number n=convFlintNSingN_QQ(c,r->cf);
    86     //fmpq_clear(c); // LEAK?
    87111    pSetCoeff0(pp,n);
    88112    pNext(pp)=p;
    89113    p=pp;
    90114  }
     115  fmpq_clear(c);
     116  omFreeSize(exp,(r->N+1)*sizeof(ulong));
    91117  p_Test(p,r);
    92118  return p;
    93 }
    94 
    95 poly Flint_Mult_MP(poly p,int lp, poly q, int lq, fmpq_mpoly_ctx_t ctx, const ring r)
    96 {
    97   fmpq_mpoly_t pp,qq,res;
    98   convSingPFlintMP(pp,ctx,p,lp,r);
    99   convSingPFlintMP(qq,ctx,q,lq,r);
    100   int bits=SI_LOG2(r->bitmask);
    101   fmpq_mpoly_init3(res,lp*lq,bits,ctx);
    102   fmpq_mpoly_mul(res,pp,qq,ctx);
    103   poly pres=convFlintMPSingP(res,ctx,r);
    104   fmpq_mpoly_clear(res,ctx);
    105   fmpq_mpoly_clear(pp,ctx);
    106   fmpq_mpoly_clear(qq,ctx);
    107   fmpq_mpoly_ctx_clear(ctx);
    108   p_Test(pres,r);
    109   return pres;
    110 }
    111 BOOLEAN convSingRFlintR(nmod_mpoly_ctx_t ctx, const ring r)
    112 {
    113   if (rRing_ord_pure_dp(r))
    114   {
    115     nmod_mpoly_ctx_init(ctx,r->N,ORD_DEGREVLEX,r->cf->ch);
    116     return FALSE;
    117   }
    118   else if (rRing_ord_pure_Dp(r))
    119   {
    120     nmod_mpoly_ctx_init(ctx,r->N,ORD_DEGLEX,r->cf->ch);
    121     return FALSE;
    122   }
    123   else if (rRing_ord_pure_lp(r))
    124   {
    125     nmod_mpoly_ctx_init(ctx,r->N,ORD_LEX,r->cf->ch);
    126     return FALSE;
    127   }
    128   return TRUE;
    129119}
    130120
     
    150140    p=pp;
    151141  }
     142  omFreeSize(exp,(r->N+1)*sizeof(ulong));
    152143  p_Test(p,r);
    153144  return p;
     
    156147void convSingPFlintMP(nmod_mpoly_t res, nmod_mpoly_ctx_t ctx, poly p, int lp,const ring r)
    157148{
    158   int bits=SI_LOG2(r->bitmask);
    159   nmod_mpoly_init3(res,lp,bits,ctx);
    160   nmod_mpoly_resize(res,lp,ctx);
     149  nmod_mpoly_init2(res, lp, ctx);
    161150  ulong* exp=(ulong*)omAlloc((r->N+1)*sizeof(ulong));
    162   int i=0;
    163151  while(p!=NULL)
    164152  {
     
    166154    #if SIZEOF_LONG==8
    167155    p_GetExpVL(p,(int64*)exp,r);
    168     nmod_mpoly_set_term_exp_ui(res,i,exp,ctx);
     156    nmod_mpoly_push_term_ui_ui(res, (ulong)n, exp, ctx);
    169157    #else
    170158    p_GetExpV(p,(int*)exp,r);
    171     nmod_mpoly_set_term_exp_ui(res,i,&(exp[1]),ctx);
     159    nmod_mpoly_push_term_ui_ui(res, (ulong)n, &(exp[1]), ctx);
    172160    #endif
    173     nmod_mpoly_set_term_coeff_ui(res,i,(ulong)n,ctx);
    174     i++;
    175161    pIter(p);
    176162  }
    177   //nmod_mpoly_print_pretty(res,r->names,ctx);PrintLn();
    178163  omFreeSize(exp,(r->N+1)*sizeof(ulong));
     164}
     165
     166#else
     167// malloc is thread safe; singular polynomials may be constructed in serial
     168
     169// like convSingNFlintN_QQ but it takes an initialized fmpq_t f
     170static void my_convSingNFlintN_QQ(fmpq_t f, number n)
     171{
     172    if (SR_HDL(n)&SR_INT)
     173    {
     174        fmpq_set_si(f,SR_TO_INT(n),1);
     175    }
     176    else if (n->s<3)
     177    {
     178        fmpz_set_mpz(fmpq_numref(f), n->z);
     179        fmpz_set_mpz(fmpq_denref(f), n->n);
     180    }
     181    else
     182    {
     183        fmpz_set_mpz(fmpq_numref(f), n->z);
     184        fmpz_one(fmpq_denref(f));
     185    }
     186}
     187
     188/* singular -> fmpq_mpoly conversion */
     189
     190class convert_sing_to_fmpq_mpoly_base
     191{
     192public:
     193    slong num_threads;
     194    slong length;
     195    fmpq_mpoly_struct * res;
     196    const fmpq_mpoly_ctx_struct * ctx;
     197    std::vector<poly> markers;
     198    ring r;
     199
     200    convert_sing_to_fmpq_mpoly_base(slong num_threads_, fmpq_mpoly_struct * res_,
     201                             const fmpq_mpoly_ctx_struct * ctx_, const ring r_)
     202      : num_threads(num_threads_),
     203        res(res_),
     204        ctx(ctx_),
     205        r(r_)
     206    {
     207    }
     208
     209    void init_poly(poly p)
     210    {
     211        length = 0;
     212        while (p != NULL)
     213        {
     214            markers.push_back(p);
     215            for (int i = 4096; i > 0; i--)
     216            {
     217                pIter(p);
     218                length++;
     219                if (p == NULL)
     220                    break;
     221            }
     222        }
     223
     224        fmpq_mpoly_init3(res, length, SI_LOG2(r->bitmask), ctx);
     225    }
     226};
     227
     228class convert_sing_to_fmpq_mpoly_arg
     229{
     230public:
     231    fmpq_t content;
     232    slong start_idx, end_idx;
     233    convert_sing_to_fmpq_mpoly_base* base;
     234
     235    convert_sing_to_fmpq_mpoly_arg() {fmpq_init(content);}
     236    ~convert_sing_to_fmpq_mpoly_arg() {fmpq_clear(content);}
     237};
     238
     239static void convert_sing_to_fmpq_mpoly_content_worker(void * varg)
     240{
     241    convert_sing_to_fmpq_mpoly_arg* arg = (convert_sing_to_fmpq_mpoly_arg*) varg;
     242    fmpq_t c;
     243    fmpq_init(c);
     244
     245    slong idx = arg->start_idx/4096;
     246    poly p = arg->base->markers[idx];
     247    idx *= 4096;
     248    while (idx < arg->start_idx)
     249    {
     250        pIter(p);
     251        idx++;
     252    }
     253
     254    /* first find gcd of coeffs */
     255    fmpq_zero(arg->content);
     256    while (idx < arg->end_idx)
     257    {
     258        my_convSingNFlintN_QQ(c, number(pGetCoeff(p)));
     259        fmpq_gcd(arg->content, arg->content, c);
     260        pIter(p);
     261        idx++;
     262    }
     263
     264    fmpq_clear(c);
     265}
     266
     267static void convert_sing_to_fmpq_mpoly_zpoly_worker(void * varg)
     268{
     269    convert_sing_to_fmpq_mpoly_arg * arg = (convert_sing_to_fmpq_mpoly_arg *) varg;
     270    convert_sing_to_fmpq_mpoly_base * base = arg->base;
     271    ulong * exp = (ulong*) flint_malloc((base->r->N + 1)*sizeof(ulong));
     272    fmpq_t c, t;
     273    fmpq_init(c);
     274    fmpq_init(t);
     275
     276    slong idx = arg->start_idx/4096;
     277    poly p = base->markers[idx];
     278    idx *= 4096;
     279    while (idx < arg->start_idx)
     280    {
     281        pIter(p);
     282        idx++;
     283    }
     284
     285    while (idx < arg->end_idx)
     286    {
     287        my_convSingNFlintN_QQ(c, number(pGetCoeff(p)));
     288        FLINT_ASSERT(!fmpq_is_zero(base->res->content));
     289        fmpq_div(t, c, base->res->content);
     290        FLINT_ASSERT(fmpz_is_one(fmpq_denref(t)));
     291
     292        slong N = mpoly_words_per_exp(base->res->zpoly->bits, base->ctx->zctx->minfo);
     293        fmpz_swap(base->res->zpoly->coeffs + idx, fmpq_numref(t));
     294        #if SIZEOF_LONG==8
     295        p_GetExpVL(p, (int64*)exp, base->r);
     296        mpoly_set_monomial_ui(base->res->zpoly->exps + N*idx, exp, base->res->zpoly->bits, base->ctx->zctx->minfo);
     297        #else
     298        p_GetExpV(p, (int*)exp, base->r);
     299        mpoly_set_monomial_ui(base->res->zpoly->exps + N*idx, &(exp[1]), base->res->zpoly->bits, base->ctx->minfo);
     300        #endif
     301
     302        pIter(p);
     303        idx++;
     304    }
     305
     306    flint_free(exp);
     307    fmpq_clear(c);
     308    fmpq_clear(t);
     309}
     310
     311/* res comes in unitialized!!!! */
     312void convSingPFlintMP(fmpq_mpoly_t res, fmpq_mpoly_ctx_t ctx, poly p, int lp, const ring r)
     313{
     314    thread_pool_handle * handles;
     315    slong num_handles;
     316    slong thread_limit = 1000;
     317
     318    /* get workers */
     319    handles = NULL;
     320    num_handles = 0;
     321    if (thread_limit > 1 && global_thread_pool_initialized)
     322    {
     323        slong max_num_handles = thread_pool_get_size(global_thread_pool);
     324        max_num_handles = FLINT_MIN(thread_limit - 1, max_num_handles);
     325        if (max_num_handles > 0)
     326        {
     327            handles = (thread_pool_handle *) flint_malloc(max_num_handles*sizeof(thread_pool_handle));
     328            num_handles = thread_pool_request(global_thread_pool, handles, max_num_handles);
     329        }
     330    }
     331
     332    convert_sing_to_fmpq_mpoly_base base(num_handles + 1, res, ctx, r);
     333    base.init_poly(p);
     334
     335    convert_sing_to_fmpq_mpoly_arg * args = new convert_sing_to_fmpq_mpoly_arg[base.num_threads];
     336    slong cur_idx = 0;
     337    for (slong i = 0; i < base.num_threads; i++)
     338    {
     339        slong next_idx = i + 1 < base.num_threads ? (i + 1)*base.length/base.num_threads : base.length;
     340        next_idx = FLINT_MAX(next_idx, cur_idx);
     341        next_idx = FLINT_MIN(next_idx, base.length);
     342        args[i].base = &base;
     343        args[i].start_idx = cur_idx;
     344        args[i].end_idx   = next_idx;
     345        cur_idx = next_idx;
     346    }
     347
     348    /* get content */
     349    for (slong i = 0; i < num_handles; i++)
     350        thread_pool_wake(global_thread_pool, handles[i], convert_sing_to_fmpq_mpoly_content_worker, args + i);
     351    convert_sing_to_fmpq_mpoly_content_worker(args + num_handles);
     352    for (slong i = 0; i < num_handles; i++)
     353    {
     354        thread_pool_wait(global_thread_pool, handles[i]);
     355    }
     356
     357    /* sign of content should match sign of first coeff */
     358    fmpq_zero(base.res->content);
     359    for (slong i = 0; i < base.num_threads; i++)
     360        fmpq_gcd(base.res->content, base.res->content, args[i].content);
     361    if (p != NULL)
     362    {
     363        fmpq_t c;
     364        fmpq_init(c);
     365        my_convSingNFlintN_QQ(c, number(pGetCoeff(p)));
     366        if (fmpq_sgn(c) < 0)
     367            fmpq_neg(base.res->content, base.res->content);
     368        fmpq_clear(c);
     369    }
     370
     371    /* fill in res->zpoly */
     372    for (slong i = 0; i < num_handles; i++)
     373        thread_pool_wake(global_thread_pool, handles[i], convert_sing_to_fmpq_mpoly_zpoly_worker, args + i);
     374    convert_sing_to_fmpq_mpoly_zpoly_worker(args + num_handles);
     375    for (slong i = 0; i < num_handles; i++)
     376        thread_pool_wait(global_thread_pool, handles[i]);
     377
     378    base.res->zpoly->length = base.length;
     379
     380    delete[] args;
     381
     382    for (slong i = 0; i < num_handles; i++)
     383        thread_pool_give_back(global_thread_pool, handles[i]);
     384    if (handles)
     385        flint_free(handles);
     386}
     387
     388/* fmpq_mpoly -> singular conversion */
     389
     390class convert_fmpq_mpoly_to_sing_base
     391{
     392public:
     393    slong num_threads;
     394    fmpq_mpoly_struct * f;
     395    const fmpq_mpoly_ctx_struct * ctx;
     396    ring r;
     397
     398    convert_fmpq_mpoly_to_sing_base(slong num_threads_, fmpq_mpoly_struct * f_,
     399                            const fmpq_mpoly_ctx_struct * ctx_, const ring r_)
     400      : num_threads(num_threads_),
     401        f(f_),
     402        ctx(ctx_),
     403        r(r_)
     404    {
     405    }
     406};
     407
     408class convert_fmpq_mpoly_to_sing_arg
     409{
     410public:
     411    poly poly_start, poly_end;
     412    slong start_idx, end_idx;
     413    convert_fmpq_mpoly_to_sing_base* base;
     414};
     415
     416static void convert_fmpq_mpoly_to_sing_worker(void * varg)
     417{
     418    convert_fmpq_mpoly_to_sing_arg * arg = (convert_fmpq_mpoly_to_sing_arg *) varg;
     419    convert_fmpq_mpoly_to_sing_base * base = arg->base;
     420    ulong* exp = (ulong*) flint_calloc(base->r->N + 1, sizeof(ulong));
     421    fmpq_t c;
     422    fmpq_init(c);
     423
     424    for (slong idx = arg->end_idx - 1; idx >= arg->start_idx; idx--)
     425    {
     426        poly pp = p_Init(base->r);
     427
     428        #if SIZEOF_LONG==8
     429        fmpq_mpoly_get_term_exp_ui(exp, base->f, idx, base->ctx);
     430        p_SetExpVL(pp, (int64*)exp, base->r);
     431        #else
     432        fmpq_mpoly_get_term_exp_ui(&(exp[1]), base->f, idx, base->ctx);
     433        p_SetExpV(pp, (int*)exp, base->r);
     434        #endif
     435        p_Setm(pp, base->r);
     436
     437        fmpq_mpoly_get_term_coeff_fmpq(c, base->f, idx, base->ctx);
     438        pSetCoeff0(pp, number(convFlintNSingN_QQ(c, base->r->cf)));
     439
     440        if (idx == arg->end_idx - 1)
     441            arg->poly_end = pp;
     442
     443        pNext(pp) = arg->poly_start;
     444        arg->poly_start = pp;
     445    }
     446
     447    flint_free(exp);
     448    fmpq_clear(c);
     449}
     450
     451poly convFlintMPSingP(fmpq_mpoly_t f, fmpq_mpoly_ctx_t ctx, const ring r)
     452{
     453    thread_pool_handle * handles;
     454    slong num_handles;
     455    slong thread_limit = 1000;
     456
     457    /* get workers */
     458    handles = NULL;
     459    num_handles = 0;
     460    if (thread_limit > 1 && global_thread_pool_initialized)
     461    {
     462        slong max_num_handles = thread_pool_get_size(global_thread_pool);
     463        max_num_handles = FLINT_MIN(thread_limit - 1, max_num_handles);
     464        if (max_num_handles > 0)
     465        {
     466            handles = (thread_pool_handle *) flint_malloc(max_num_handles*sizeof(thread_pool_handle));
     467            num_handles = thread_pool_request(global_thread_pool, handles, max_num_handles);
     468        }
     469    }
     470
     471    convert_fmpq_mpoly_to_sing_base base(num_handles + 1, f, ctx, r);
     472
     473    convert_fmpq_mpoly_to_sing_arg * args = new convert_fmpq_mpoly_to_sing_arg[base.num_threads];
     474    slong cur_idx = 0;
     475    for (slong i = 0; i < base.num_threads; i++)
     476    {
     477        slong next_idx = i + 1 < base.num_threads ? (i + 1)*base.f->zpoly->length/base.num_threads
     478                                                  : base.f->zpoly->length;
     479        next_idx = FLINT_MAX(next_idx, cur_idx);
     480        next_idx = FLINT_MIN(next_idx, base.f->zpoly->length);
     481        args[i].base = &base;
     482        args[i].start_idx = cur_idx;
     483        args[i].end_idx   = next_idx;
     484        cur_idx = next_idx;
     485        args[i].poly_end = NULL;
     486        args[i].poly_start = NULL;
     487    }
     488
     489    /* construct pieces */
     490    for (slong i = 0; i < num_handles; i++)
     491        thread_pool_wake(global_thread_pool, handles[i], convert_fmpq_mpoly_to_sing_worker, args + i);
     492    convert_fmpq_mpoly_to_sing_worker(args + num_handles);
     493    for (slong i = 0; i < num_handles; i++)
     494        thread_pool_wait(global_thread_pool, handles[i]);
     495
     496    /* join pieces */
     497    poly p = NULL;
     498    for (slong i = num_handles; i >= 0; i--)
     499    {
     500        if (args[i].poly_end != NULL)
     501        {
     502            pNext(args[i].poly_end) = p;
     503            p = args[i].poly_start;
     504        }
     505    }
     506
     507    for (slong i = 0; i < num_handles; i++)
     508        thread_pool_give_back(global_thread_pool, handles[i]);
     509    if (handles)
     510        flint_free(handles);
     511
     512    p_Test(p,r);
     513    return p;
     514}
     515
     516
     517/* singular -> nmod_mpoly conversion */
     518
     519class convert_sing_to_nmod_mpoly_base
     520{
     521public:
     522    slong num_threads;
     523    slong length;
     524    nmod_mpoly_struct * res;
     525    const nmod_mpoly_ctx_struct * ctx;
     526    std::vector<poly> markers;
     527    ring r;
     528
     529    convert_sing_to_nmod_mpoly_base(slong num_threads_, nmod_mpoly_struct * res_,
     530                            const nmod_mpoly_ctx_struct * ctx_, const ring r_)
     531      : num_threads(num_threads_),
     532        res(res_),
     533        ctx(ctx_),
     534        r(r_)
     535    {
     536    }
     537
     538    void init_poly(poly p)
     539    {
     540        length = 0;
     541        while (p != NULL)
     542        {
     543            markers.push_back(p);
     544            for (int i = 4096; i > 0; i--)
     545            {
     546                pIter(p);
     547                length++;
     548                if (p == NULL)
     549                    break;
     550            }
     551        }
     552
     553        nmod_mpoly_init3(res, length, SI_LOG2(r->bitmask), ctx);
     554    }
     555};
     556
     557class convert_sing_to_nmod_mpoly_arg
     558{
     559public:
     560    slong start_idx, end_idx;
     561    convert_sing_to_nmod_mpoly_base* base;
     562};
     563
     564static void convert_sing_to_nmod_mpoly_worker(void * varg)
     565{
     566    convert_sing_to_nmod_mpoly_arg * arg = (convert_sing_to_nmod_mpoly_arg *) varg;
     567    convert_sing_to_nmod_mpoly_base * base = arg->base;
     568    ulong * exp = (ulong*) flint_malloc((base->r->N + 1)*sizeof(ulong));
     569
     570    slong idx = arg->start_idx/4096;
     571    poly p = base->markers[idx];
     572    idx *= 4096;
     573    while (idx < arg->start_idx)
     574    {
     575        pIter(p);
     576        idx++;
     577    }
     578
     579    while (idx < arg->end_idx)
     580    {
     581        slong N = mpoly_words_per_exp(base->res->bits, base->ctx->minfo);
     582        #if SIZEOF_LONG==8
     583        p_GetExpVL(p, (int64*)exp, base->r);
     584        mpoly_set_monomial_ui(base->res->exps + N*idx, exp, base->res->bits, base->ctx->minfo);
     585        #else
     586        p_GetExpV(p, (int*)exp, base->r);
     587        mpoly_set_monomial_ui(base->res->exps + N*idx, &(exp[1]), base->res->bits, base->ctx->minfo);
     588        #endif
     589
     590        base->res->coeffs[idx] = (ulong)(number(pGetCoeff(p)));
     591
     592        pIter(p);
     593        idx++;
     594    }
     595
     596    flint_free(exp);
     597}
     598
     599/* res comes in unitialized!!!! */
     600void convSingPFlintMP(nmod_mpoly_t res, nmod_mpoly_ctx_t ctx, poly p, int lp, const ring r)
     601{
     602    thread_pool_handle * handles;
     603    slong num_handles;
     604    slong thread_limit = 1000;
     605
     606    /* get workers */
     607    handles = NULL;
     608    num_handles = 0;
     609    if (thread_limit > 1 && global_thread_pool_initialized)
     610    {
     611        slong max_num_handles = thread_pool_get_size(global_thread_pool);
     612        max_num_handles = FLINT_MIN(thread_limit - 1, max_num_handles);
     613        if (max_num_handles > 0)
     614        {
     615            handles = (thread_pool_handle *) flint_malloc(max_num_handles*sizeof(thread_pool_handle));
     616            num_handles = thread_pool_request(global_thread_pool, handles, max_num_handles);
     617        }
     618    }
     619
     620    convert_sing_to_nmod_mpoly_base base(num_handles + 1, res, ctx, r);
     621    base.init_poly(p);
     622
     623    convert_sing_to_nmod_mpoly_arg * args = new convert_sing_to_nmod_mpoly_arg[base.num_threads];
     624    slong cur_idx = 0;
     625    for (slong i = 0; i < base.num_threads; i++)
     626    {
     627        slong next_idx = i + 1 < base.num_threads ? (i + 1)*base.length/base.num_threads
     628                                                  : base.length;
     629        next_idx = FLINT_MAX(next_idx, cur_idx);
     630        next_idx = FLINT_MIN(next_idx, base.length);
     631        args[i].base = &base;
     632        args[i].start_idx = cur_idx;
     633        args[i].end_idx   = next_idx;
     634        cur_idx = next_idx;
     635    }
     636
     637    /* fill in res */
     638    for (slong i = 0; i < num_handles; i++)
     639        thread_pool_wake(global_thread_pool, handles[i], convert_sing_to_nmod_mpoly_worker, args + i);
     640    convert_sing_to_nmod_mpoly_worker(args + num_handles);
     641    for (slong i = 0; i < num_handles; i++)
     642        thread_pool_wait(global_thread_pool, handles[i]);
     643
     644    base.res->length = base.length;
     645
     646    delete[] args;
     647
     648    for (slong i = 0; i < num_handles; i++)
     649        thread_pool_give_back(global_thread_pool, handles[i]);
     650    if (handles)
     651        flint_free(handles);
     652}
     653
     654
     655/* nmod_mpoly -> singular conversion */
     656
     657class convert_nmod_mpoly_to_sing_base
     658{
     659public:
     660    slong num_threads;
     661    nmod_mpoly_struct * f;
     662    const nmod_mpoly_ctx_struct * ctx;
     663    ring r;
     664
     665    convert_nmod_mpoly_to_sing_base(slong num_threads_, nmod_mpoly_struct * f_,
     666                             const nmod_mpoly_ctx_struct * ctx_, const ring r_)
     667      : num_threads(num_threads_),
     668        f(f_),
     669        ctx(ctx_),
     670        r(r_)
     671    {
     672    }
     673};
     674
     675class convert_nmod_mpoly_to_sing_arg
     676{
     677public:
     678    poly poly_start, poly_end;
     679    slong start_idx, end_idx;
     680    convert_nmod_mpoly_to_sing_base* base;
     681};
     682
     683static void convert_nmod_mpoly_to_sing_worker(void * varg)
     684{
     685    convert_nmod_mpoly_to_sing_arg * arg = (convert_nmod_mpoly_to_sing_arg *) varg;
     686    convert_nmod_mpoly_to_sing_base * base = arg->base;
     687    ulong* exp = (ulong*) flint_calloc(base->r->N + 1, sizeof(ulong));
     688
     689    for (slong idx = arg->end_idx - 1; idx >= arg->start_idx; idx--)
     690    {
     691        poly pp = p_Init(base->r);
     692
     693        #if SIZEOF_LONG==8
     694        nmod_mpoly_get_term_exp_ui(exp, base->f, idx, base->ctx);
     695        p_SetExpVL(pp, (int64*)exp, base->r);
     696        #else
     697        nmod_mpoly_get_term_exp_ui(&(exp[1]), base->f, idx, base->ctx);
     698        p_SetExpV(pp, (int*)exp, base->r);
     699        #endif
     700        p_Setm(pp, base->r);
     701
     702        pSetCoeff0(pp, number(nmod_mpoly_get_term_coeff_ui(base->f, idx, base->ctx)));
     703
     704        if (idx == arg->end_idx - 1)
     705            arg->poly_end = pp;
     706
     707        pNext(pp) = arg->poly_start;
     708        arg->poly_start = pp;
     709    }
     710
     711    flint_free(exp);
     712}
     713
     714poly convFlintMPSingP(nmod_mpoly_t f, nmod_mpoly_ctx_t ctx, const ring r)
     715{
     716    thread_pool_handle * handles;
     717    slong num_handles;
     718    slong thread_limit = 1000;
     719
     720    /* get workers */
     721    handles = NULL;
     722    num_handles = 0;
     723    if (thread_limit > 1 && global_thread_pool_initialized)
     724    {
     725        slong max_num_handles = thread_pool_get_size(global_thread_pool);
     726        max_num_handles = FLINT_MIN(thread_limit - 1, max_num_handles);
     727        if (max_num_handles > 0)
     728        {
     729            handles = (thread_pool_handle *) flint_malloc(max_num_handles*sizeof(thread_pool_handle));
     730            num_handles = thread_pool_request(global_thread_pool, handles, max_num_handles);
     731        }
     732    }
     733
     734    convert_nmod_mpoly_to_sing_base base(num_handles + 1, f, ctx, r);
     735
     736    convert_nmod_mpoly_to_sing_arg * args = new convert_nmod_mpoly_to_sing_arg[base.num_threads];
     737    slong cur_idx = 0;
     738    for (slong i = 0; i < base.num_threads; i++)
     739    {
     740        slong next_idx = i + 1 < base.num_threads ? (i + 1)*base.f->length/base.num_threads
     741                                                  : base.f->length;
     742        next_idx = FLINT_MAX(next_idx, cur_idx);
     743        next_idx = FLINT_MIN(next_idx, base.f->length);
     744        args[i].base = &base;
     745        args[i].start_idx = cur_idx;
     746        args[i].end_idx   = next_idx;
     747        cur_idx = next_idx;
     748        args[i].poly_end = NULL;
     749        args[i].poly_start = NULL;
     750    }
     751
     752    /* construct pieces */
     753    for (slong i = 0; i < num_handles; i++)
     754        thread_pool_wake(global_thread_pool, handles[i], convert_nmod_mpoly_to_sing_worker, args + i);
     755    convert_nmod_mpoly_to_sing_worker(args + num_handles);
     756    for (slong i = 0; i < num_handles; i++)
     757        thread_pool_wait(global_thread_pool, handles[i]);
     758
     759    /* join pieces */
     760    poly p = NULL;
     761    for (slong i = num_handles; i >= 0; i--)
     762    {
     763        if (args[i].poly_end != NULL)
     764        {
     765            pNext(args[i].poly_end) = p;
     766            p = args[i].poly_start;
     767        }
     768    }
     769
     770    for (slong i = 0; i < num_handles; i++)
     771        thread_pool_give_back(global_thread_pool, handles[i]);
     772    if (handles)
     773        flint_free(handles);
     774
     775    p_Test(p,r);
     776    return p;
     777}
     778
     779#endif
     780
     781/****** polynomial operations ***********/
     782
     783poly Flint_Mult_MP(poly p,int lp, poly q, int lq, fmpq_mpoly_ctx_t ctx, const ring r)
     784{
     785  fmpq_mpoly_t pp,qq,res;
     786  convSingPFlintMP(pp,ctx,p,lp,r);
     787  convSingPFlintMP(qq,ctx,q,lq,r);
     788  fmpq_mpoly_init(res,ctx);
     789  fmpq_mpoly_mul(res,pp,qq,ctx);
     790  poly pres=convFlintMPSingP(res,ctx,r);
     791  fmpq_mpoly_clear(res,ctx);
     792  fmpq_mpoly_clear(pp,ctx);
     793  fmpq_mpoly_clear(qq,ctx);
     794  fmpq_mpoly_ctx_clear(ctx);
     795  p_Test(pres,r);
     796  return pres;
    179797}
    180798
     
    184802  convSingPFlintMP(pp,ctx,p,lp,r);
    185803  convSingPFlintMP(qq,ctx,q,lq,r);
    186   int bits=SI_LOG2(r->bitmask);
    187   nmod_mpoly_init3(res,lp*lq,bits,ctx);
     804  nmod_mpoly_init(res,ctx);
    188805  nmod_mpoly_mul(res,pp,qq,ctx);
    189806  poly pres=convFlintMPSingP(res,ctx,r);
     
    195812  return pres;
    196813}
     814
    197815poly Flint_GCD_MP(poly p,int lp,poly q,int lq,nmod_mpoly_ctx_t ctx,const ring r)
    198816{
     
    200818  convSingPFlintMP(pp,ctx,p,lp,r);
    201819  convSingPFlintMP(qq,ctx,q,lq,r);
    202   int bits=SI_LOG2(r->bitmask);
    203   nmod_mpoly_init3(res,si_max(lp,lq),bits,ctx);
     820  nmod_mpoly_init(res,ctx);
    204821  int ok=nmod_mpoly_gcd(res,pp,qq,ctx);
    205822  poly pres;
     
    225842  convSingPFlintMP(pp,ctx,p,lp,r);
    226843  convSingPFlintMP(qq,ctx,q,lq,r);
    227   int bits=SI_LOG2(r->bitmask);
    228   fmpq_mpoly_init3(res,si_max(lp,lq),bits,ctx);
     844  fmpq_mpoly_init(res,ctx);
    229845  int ok=fmpq_mpoly_gcd(res,pp,qq,ctx);
    230846  poly pres;
    231847  if (ok)
    232848  {
     849    // Flint normalizes the gcd to be monic.
     850    // Singular wants a gcd defined over ZZ that is primitive and has a positive leading coeff.
     851    if (!fmpq_mpoly_is_zero(res, ctx))
     852    {
     853      fmpq_t content;
     854      fmpq_init(content);
     855      fmpq_mpoly_content(content, res, ctx);
     856      fmpq_mpoly_scalar_div_fmpq(res, res, content, ctx);
     857      fmpq_clear(content);
     858    }
    233859    pres=convFlintMPSingP(res,ctx,r);
    234860    p_Test(pres,r);
     
    244870  return pres;
    245871}
     872
    246873#endif
    247874#endif
  • libpolys/polys/flintconv.cc

    rbc4eb7 r1e351b3  
    226226  int d=fmpq_poly_length(f);
    227227  poly p=NULL;
     228  fmpq_t c;
     229  fmpq_init(c);
    228230  for(int i=0; i<=d; i++)
    229231  {
    230     fmpq_t c;
    231     fmpq_init(c);
    232232    fmpq_poly_get_coeff_fmpq(c,f,i);
    233233    number n=convFlintNSingN(c,r->cf);
     
    238238    p=p_Add_q(p,pp,r);
    239239  }
     240  fmpq_clear(c);
    240241  p_Test(p,r);
    241242  return p;
  • libpolys/polys/monomials/ring.cc

    r84f0d8 r1e351b3  
    58165816      R->names=(char**)omReallocSize(R->names,r->N*sizeof(char_ptr),R->N*sizeof(char_ptr));
    58175817    }
    5818     else i--;
     5818    i--;
    58195819  }
    58205820  R->block1[p]=R->N;
    5821   rComplete(R);
     5821  rComplete(R,1);
    58225822  return R;
    58235823}
Note: See TracChangeset for help on using the changeset viewer.