Changeset 690c98 in git


Ignore:
Timestamp:
Jul 24, 2019, 4:31:32 PM (5 years ago)
Author:
tthsqe12 <tthsqe12@…>
Branches:
(u'fieker-DuVal', '117eb8c30fc9e991c4decca4832b1d19036c4c65')(u'spielwiese', 'b4f17ed1d25f93d46dbe29e4b499baecc2fd51bb')
Children:
3f94cbdb24129f0002a999ab5fb63a34f0fc933d
Parents:
fa1cd304b94fb2782b47874564b731a57670c34c
Message:
opt: let fmpq_mpoly borrow in some cases
File:
1 edited

Legend:

Unmodified
Added
Removed

  • libpolys/polys/flint_mpoly.cc

    rfa1cd3 r690c98  
    6565#if 1
    6666// memory allocation is not thread safe; singular polynomials must be constructed in serial
     67
     68/*
     69    We agree the that result of a singular -> fmpq_mpoly conversion is
     70    readonly. This restricts the usage of the result in flint functions to
     71    const arguments. However, the real readonly conversion is currently only
     72    implemented in the threaded conversion below since it requires a scan of
     73    all coefficients anyways. The _fmpq_mpoly_clear_readonly_sing needs to
     74    be provided for a consistent interface in the polynomial operations.
     75*/
     76static void _fmpq_mpoly_clear_readonly_sing(fmpq_mpoly_t a, fmpq_mpoly_ctx_t ctx)
     77{
     78    fmpq_mpoly_clear(a, ctx);
     79}
    6780
    6881void convSingPFlintMP(fmpq_mpoly_t res, fmpq_mpoly_ctx_t ctx, poly p, int lp, const ring r)
    … …  
    186199}
    187200
     201
     202/*
     203    In order that flint may sometimes borrow the large integer coeffs of
     204    polynomials over QQ (borrow means: simply point to the same GMP structs
     205    that singular has already allocated), we define the result of a
     206    singular -> fmpq_mpoly conversion to be readonly. This means we agree
     207    that
     208        - it can only be used as an const argument to a flint function
     209        - singular must not mutate the original coeffs while the readonly object is in use
     210*/
     211
     212static void _fmpq_mpoly_clear_readonly_sing(fmpq_mpoly_t a, fmpq_mpoly_ctx_t ctx)
     213{
     214    if (fmpq_is_one(a->content))
     215    {
     216        if (a->zpoly->alloc > 0)
     217        {
     218            flint_free(a->zpoly->coeffs);
     219            flint_free(a->zpoly->exps);
     220        }
     221
     222        fmpq_clear(a->content);
     223    }
     224    else
     225    {
     226        fmpq_mpoly_clear(a, ctx);
     227    }
     228}
     229
    188230/* singular -> fmpq_mpoly conversion */
    189231
    … …  
    197239    std::vector<poly> markers;
    198240    ring r;
     241    fmpq_t content;
    199242
    200243    convert_sing_to_fmpq_mpoly_base(slong num_threads_, fmpq_mpoly_struct * res_,
    … …  
    205248        r(r_)
    206249    {
     250        fmpq_t c;
     251        fmpq_init(c);
     252        fmpq_init(content);
     253        fmpq_zero(content);
     254
    207255        length = 0;
    208256        while (1)
    209257        {
    210258            if ((length % 4096) == 0)
     259            {
     260                my_convSingNFlintN_QQ(c, number(pGetCoeff(p)));
     261                fmpq_gcd(content, content, c);
    211262                markers.push_back(p);
     263            }
    212264            if (p == NULL)
    213265                return;
    … …  
    215267            pIter(p);
    216268        }
     269
     270        fmpq_clear(c);
     271    }
     272
     273    ~convert_sing_to_fmpq_mpoly_base()
     274    {
     275        fmpq_clear(content);
    217276    }
    218277};
    … …  
    248307
    249308    flint_bitcnt_t required_bits = MPOLY_MIN_BITS;
    250     fmpq_zero(arg->content);
     309    fmpq_set(arg->content, base->content);
    251310
    252311    while (idx < arg->end_idx)
    253312    {
    254         my_convSingNFlintN_QQ(c, number(pGetCoeff(p)));
    255         fmpq_gcd(arg->content, arg->content, c);
     313        number n = number(pGetCoeff(p));
     314
     315        if (fmpq_is_one(arg->content) && (SR_HDL(n)&SR_INT || n->s >= 3))
     316        {
     317            /* content is 1 and n is an integer, nothing to do */
     318        }
     319        else
     320        {
     321            my_convSingNFlintN_QQ(c, n);
     322            fmpq_gcd(arg->content, arg->content, c);
     323        }
    256324
    257325        #if SIZEOF_LONG==8
    … …  
    293361    }
    294362
     363    slong N = mpoly_words_per_exp(base->res->zpoly->bits, base->ctx->zctx->minfo);
     364    fmpz * res_coeffs = base->res->zpoly->coeffs;
     365    ulong * res_exps = base->res->zpoly->exps;
     366    flint_bitcnt_t res_bits = base->res->zpoly->bits;
     367
    295368    while (idx < arg->end_idx)
    296369    {
    297         my_convSingNFlintN_QQ(c, number(pGetCoeff(p)));
    298         FLINT_ASSERT(!fmpq_is_zero(base->res->content));
    299         fmpq_div(t, c, base->res->content);
    300         FLINT_ASSERT(fmpz_is_one(fmpq_denref(t)));
    301 
    302         slong N = mpoly_words_per_exp(base->res->zpoly->bits, base->ctx->zctx->minfo);
    303         fmpz_swap(base->res->zpoly->coeffs + idx, fmpq_numref(t));
     370        if (fmpq_is_one(base->res->content))
     371        {
     372            // borrowing singular integers
     373            // the entry res_coeffs[idx] is junk, we should just overwrite it
     374
     375            number n = number(pGetCoeff(p));
     376
     377            if (SR_HDL(n)&SR_INT)
     378            {
     379                // n is a singular-small integer
     380                res_coeffs[idx] = SR_TO_INT(n);
     381            }
     382            else if (n->s<3)
     383            {
     384                // n is an element of QQ \ ZZ, should not happen
     385                assume(false);
     386            }
     387            else
     388            {
     389                // n is a singular-large integer, n may be flint-small
     390                res_coeffs[idx] = PTR_TO_COEFF(n->z);
     391                if (fmpz_fits_si(res_coeffs + idx))
     392                {
     393                    slong val = fmpz_get_si(res_coeffs + idx);
     394                    if (val >= COEFF_MIN && val <= COEFF_MAX)
     395                        res_coeffs[idx] = val;
     396                }
     397            }
     398        }
     399        else
     400        {
     401            my_convSingNFlintN_QQ(c, number(pGetCoeff(p)));
     402            FLINT_ASSERT(!fmpq_is_zero(base->res->content));
     403            fmpq_div(t, c, base->res->content);
     404            FLINT_ASSERT(fmpz_is_one(fmpq_denref(t)));
     405            fmpz_swap(res_coeffs + idx, fmpq_numref(t));
     406        }
     407
    304408        #if SIZEOF_LONG==8
    305409        p_GetExpVL(p, (int64*)exp, base->r);
    306         mpoly_set_monomial_ui(base->res->zpoly->exps + N*idx, exp, base->res->zpoly->bits, base->ctx->zctx->minfo);
     410        mpoly_set_monomial_ui(res_exps + N*idx, exp, res_bits, base->ctx->zctx->minfo);
    307411        #else
    308412        p_GetExpV(p, (int*)exp, base->r);
    309         mpoly_set_monomial_ui(base->res->zpoly->exps + N*idx, &(exp[1]), base->res->zpoly->bits, base->ctx->minfo);
     413        mpoly_set_monomial_ui(res_exps + N*idx, &(exp[1]), res_bits, base->ctx->minfo);
    310414        #endif
    311415
    … …  
    340444    }
    341445
     446    /* the constructor works out the length of p and sets some markers */
    342447    convert_sing_to_fmpq_mpoly_base base(num_handles + 1, res, ctx, r, p);
    343448
    … …  
    348453    {
    349454        slong next_idx = i + 1 < base.num_threads ? (i + 1)*base.length/base.num_threads : base.length;
     455        FLINT_ASSERT(cur_idx <= base.length);
    350456        next_idx = FLINT_MAX(next_idx, cur_idx);
    351457        next_idx = FLINT_MIN(next_idx, base.length);
    … …  
    367473        required_bits = FLINT_MAX(required_bits, args[i].required_bits);
    368474
    369     /* initialize res with optimal bits */
    370     fmpq_mpoly_init3(res, base.length, mpoly_fix_bits(required_bits, ctx->zctx->minfo), ctx);
    371 
    372475    /* sign of content should match sign of first coeff */
    373     fmpq_zero(base.res->content);
     476    fmpq_t content;
     477    fmpq_init(content);
     478    fmpq_zero(content);
    374479    for (slong i = 0; i < base.num_threads; i++)
    375         fmpq_gcd(base.res->content, base.res->content, args[i].content);
     480        fmpq_gcd(content, content, args[i].content);
    376481    if (p != NULL)
    377482    {
    … …  
    380485        my_convSingNFlintN_QQ(c, number(pGetCoeff(p)));
    381486        if (fmpq_sgn(c) < 0)
    382             fmpq_neg(base.res->content, base.res->content);
     487            fmpq_neg(content, content);
    383488        fmpq_clear(c);
    384489    }
     490
     491    /* initialize res with optimal bits */
     492    required_bits = mpoly_fix_bits(required_bits, ctx->zctx->minfo);
     493    if (fmpq_is_one(content))
     494    {
     495        /* initialize borrowed coeffs */
     496        slong N = mpoly_words_per_exp(required_bits, ctx->zctx->minfo);
     497        slong alloc = base.length;
     498        if (alloc != 0)
     499        {
     500            res->zpoly->coeffs = (fmpz *) flint_malloc(alloc*sizeof(fmpz));
     501            res->zpoly->exps   = (ulong *) flint_malloc(alloc*N*sizeof(ulong));
     502        }
     503        else
     504        {
     505            res->zpoly->coeffs = NULL;
     506            res->zpoly->exps = NULL;
     507        }
     508        res->zpoly->alloc = alloc;
     509        res->zpoly->length = 0;
     510        res->zpoly->bits = required_bits;
     511
     512        fmpq_init(res->content);
     513        fmpq_one(res->content);
     514    }
     515    else
     516    {
     517        /* initialize coeffs that will be created and destroyed */
     518        fmpq_mpoly_init3(res, base.length, required_bits, ctx);
     519        fmpq_swap(res->content, content);
     520    }
     521
     522    fmpq_clear(content);
    385523
    386524    /* fill in res->zpoly */
    … …  
    492630        slong next_idx = i + 1 < base.num_threads ? (i + 1)*base.f->zpoly->length/base.num_threads
    493631                                                  : base.f->zpoly->length;
     632        FLINT_ASSERT(cur_idx <= base.f->zpoly->length);
    494633        next_idx = FLINT_MAX(next_idx, cur_idx);
    495634        next_idx = FLINT_MIN(next_idx, base.f->zpoly->length);
    … …  
    624763    }
    625764
     765    slong N = mpoly_words_per_exp(base->res->bits, base->ctx->minfo);
     766    ulong * res_coeffs = base->res->coeffs;
     767    ulong * res_exps = base->res->exps;
     768    flint_bitcnt_t res_bits = base->res->bits;
     769
    626770    while (idx < arg->end_idx)
    627771    {
    628         slong N = mpoly_words_per_exp(base->res->bits, base->ctx->minfo);
    629772        #if SIZEOF_LONG==8
    630773        p_GetExpVL(p, (int64*)exp, base->r);
    631         mpoly_set_monomial_ui(base->res->exps + N*idx, exp, base->res->bits, base->ctx->minfo);
     774        mpoly_set_monomial_ui(res_exps + N*idx, exp, res_bits, base->ctx->minfo);
    632775        #else
    633776        p_GetExpV(p, (int*)exp, base->r);
    634         mpoly_set_monomial_ui(base->res->exps + N*idx, &(exp[1]), base->res->bits, base->ctx->minfo);
     777        mpoly_set_monomial_ui(res_exps + N*idx, &(exp[1]), res_bits, base->ctx->minfo);
    635778        #endif
    636779
    637         base->res->coeffs[idx] = (ulong)(number(pGetCoeff(p)));
     780        res_coeffs[idx] = (ulong)(number(pGetCoeff(p)));
    638781
    639782        pIter(p);
    … …  
    674817        slong next_idx = i + 1 < base.num_threads ? (i + 1)*base.length/base.num_threads
    675818                                                  : base.length;
     819        FLINT_ASSERT(cur_idx <= base.length);
    676820        next_idx = FLINT_MAX(next_idx, cur_idx);
    677821        next_idx = FLINT_MIN(next_idx, base.length);
    … …  
    801945        slong next_idx = i + 1 < base.num_threads ? (i + 1)*base.f->length/base.num_threads
    802946                                                  : base.f->length;
     947        FLINT_ASSERT(cur_idx <= base.f->length);
    803948        next_idx = FLINT_MAX(next_idx, cur_idx);
    804949        next_idx = FLINT_MIN(next_idx, base.f->length);
    … …  
    845990{
    846991  fmpq_mpoly_t pp,qq,res;
    847   convSingPFlintMP(pp,ctx,p,lp,r);
    848   convSingPFlintMP(qq,ctx,q,lq,r);
     992  convSingPFlintMP(pp,ctx,p,lp,r); // pp read only
     993  convSingPFlintMP(qq,ctx,q,lq,r); // qq read only
    849994  fmpq_mpoly_init(res,ctx);
    850995  fmpq_mpoly_mul(res,pp,qq,ctx);
    851996  poly pres=convFlintMPSingP(res,ctx,r);
    852997  fmpq_mpoly_clear(res,ctx);
    853   fmpq_mpoly_clear(pp,ctx);
    854   fmpq_mpoly_clear(qq,ctx);
     998  _fmpq_mpoly_clear_readonly_sing(pp,ctx);
     999  _fmpq_mpoly_clear_readonly_sing(qq,ctx);
    8551000  fmpq_mpoly_ctx_clear(ctx);
    8561001  p_Test(pres,r);
    … …  
    8781023{
    8791024  fmpq_mpoly_t pp,qq,res;
    880   convSingPFlintMP(pp,ctx,p,lp,r);
    881   convSingPFlintMP(qq,ctx,q,lq,r);
     1025  convSingPFlintMP(pp,ctx,p,lp,r); // pp read only
     1026  convSingPFlintMP(qq,ctx,q,lq,r); // qq read only
    8821027  fmpq_mpoly_init(res,ctx);
    8831028  fmpq_mpoly_divides(res,pp,qq,ctx);
    8841029  poly pres = convFlintMPSingP(res,ctx,r);
    8851030  fmpq_mpoly_clear(res,ctx);
    886   fmpq_mpoly_clear(pp,ctx);
    887   fmpq_mpoly_clear(qq,ctx);
     1031  _fmpq_mpoly_clear_readonly_sing(pp,ctx);
     1032  _fmpq_mpoly_clear_readonly_sing(qq,ctx);
    8881033  fmpq_mpoly_ctx_clear(ctx);
    8891034  p_Test(pres,r);
    … …  
    9341079{
    9351080  fmpq_mpoly_t pp,qq,res;
    936   convSingPFlintMP(pp,ctx,p,lp,r);
    937   convSingPFlintMP(qq,ctx,q,lq,r);
     1081  convSingPFlintMP(pp,ctx,p,lp,r); // pp read only
     1082  convSingPFlintMP(qq,ctx,q,lq,r); // qq read only
    9381083  fmpq_mpoly_init(res,ctx);
    9391084  int ok=fmpq_mpoly_gcd(res,pp,qq,ctx);
    … …  
    9591104  }
    9601105  fmpq_mpoly_clear(res,ctx);
    961   fmpq_mpoly_clear(pp,ctx);
    962   fmpq_mpoly_clear(qq,ctx);
     1106  _fmpq_mpoly_clear_readonly_sing(pp,ctx);
     1107  _fmpq_mpoly_clear_readonly_sing(qq,ctx);
    9631108  fmpq_mpoly_ctx_clear(ctx);
    9641109  return pres;
Note: See TracChangeset for help on using the changeset viewer.