Changeset 690c98 in git
 Timestamp:
 Jul 24, 2019, 4:31:32 PM (4 years ago)
 Branches:
 (u'jengelhdatetime', 'ceac47cbc86fe4a15902392bdbb9bd2ae0ea02c6')(u'spielwiese', '0604212ebb110535022efecad887940825b97c3f')
 Children:
 3f94cbdb24129f0002a999ab5fb63a34f0fc933d
 Parents:
 fa1cd304b94fb2782b47874564b731a57670c34c
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

libpolys/polys/flint_mpoly.cc
rfa1cd3 r690c98 65 65 #if 1 66 66 // 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 */ 76 static void _fmpq_mpoly_clear_readonly_sing(fmpq_mpoly_t a, fmpq_mpoly_ctx_t ctx) 77 { 78 fmpq_mpoly_clear(a, ctx); 79 } 67 80 68 81 void convSingPFlintMP(fmpq_mpoly_t res, fmpq_mpoly_ctx_t ctx, poly p, int lp, const ring r) … … 186 199 } 187 200 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 212 static 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 188 230 /* singular > fmpq_mpoly conversion */ 189 231 … … 197 239 std::vector<poly> markers; 198 240 ring r; 241 fmpq_t content; 199 242 200 243 convert_sing_to_fmpq_mpoly_base(slong num_threads_, fmpq_mpoly_struct * res_, … … 205 248 r(r_) 206 249 { 250 fmpq_t c; 251 fmpq_init(c); 252 fmpq_init(content); 253 fmpq_zero(content); 254 207 255 length = 0; 208 256 while (1) 209 257 { 210 258 if ((length % 4096) == 0) 259 { 260 my_convSingNFlintN_QQ(c, number(pGetCoeff(p))); 261 fmpq_gcd(content, content, c); 211 262 markers.push_back(p); 263 } 212 264 if (p == NULL) 213 265 return; … … 215 267 pIter(p); 216 268 } 269 270 fmpq_clear(c); 271 } 272 273 ~convert_sing_to_fmpq_mpoly_base() 274 { 275 fmpq_clear(content); 217 276 } 218 277 }; … … 248 307 249 308 flint_bitcnt_t required_bits = MPOLY_MIN_BITS; 250 fmpq_ zero(arg>content);309 fmpq_set(arg>content, base>content); 251 310 252 311 while (idx < arg>end_idx) 253 312 { 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 } 256 324 257 325 #if SIZEOF_LONG==8 … … 293 361 } 294 362 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 295 368 while (idx < arg>end_idx) 296 369 { 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 singularsmall 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 singularlarge integer, n may be flintsmall 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 304 408 #if SIZEOF_LONG==8 305 409 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); 307 411 #else 308 412 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); 310 414 #endif 311 415 … … 340 444 } 341 445 446 /* the constructor works out the length of p and sets some markers */ 342 447 convert_sing_to_fmpq_mpoly_base base(num_handles + 1, res, ctx, r, p); 343 448 … … 348 453 { 349 454 slong next_idx = i + 1 < base.num_threads ? (i + 1)*base.length/base.num_threads : base.length; 455 FLINT_ASSERT(cur_idx <= base.length); 350 456 next_idx = FLINT_MAX(next_idx, cur_idx); 351 457 next_idx = FLINT_MIN(next_idx, base.length); … … 367 473 required_bits = FLINT_MAX(required_bits, args[i].required_bits); 368 474 369 /* initialize res with optimal bits */370 fmpq_mpoly_init3(res, base.length, mpoly_fix_bits(required_bits, ctx>zctx>minfo), ctx);371 372 475 /* 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); 374 479 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); 376 481 if (p != NULL) 377 482 { … … 380 485 my_convSingNFlintN_QQ(c, number(pGetCoeff(p))); 381 486 if (fmpq_sgn(c) < 0) 382 fmpq_neg( base.res>content, base.res>content);487 fmpq_neg(content, content); 383 488 fmpq_clear(c); 384 489 } 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); 385 523 386 524 /* fill in res>zpoly */ … … 492 630 slong next_idx = i + 1 < base.num_threads ? (i + 1)*base.f>zpoly>length/base.num_threads 493 631 : base.f>zpoly>length; 632 FLINT_ASSERT(cur_idx <= base.f>zpoly>length); 494 633 next_idx = FLINT_MAX(next_idx, cur_idx); 495 634 next_idx = FLINT_MIN(next_idx, base.f>zpoly>length); … … 624 763 } 625 764 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 626 770 while (idx < arg>end_idx) 627 771 { 628 slong N = mpoly_words_per_exp(base>res>bits, base>ctx>minfo);629 772 #if SIZEOF_LONG==8 630 773 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); 632 775 #else 633 776 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); 635 778 #endif 636 779 637 base>res>coeffs[idx] = (ulong)(number(pGetCoeff(p)));780 res_coeffs[idx] = (ulong)(number(pGetCoeff(p))); 638 781 639 782 pIter(p); … … 674 817 slong next_idx = i + 1 < base.num_threads ? (i + 1)*base.length/base.num_threads 675 818 : base.length; 819 FLINT_ASSERT(cur_idx <= base.length); 676 820 next_idx = FLINT_MAX(next_idx, cur_idx); 677 821 next_idx = FLINT_MIN(next_idx, base.length); … … 801 945 slong next_idx = i + 1 < base.num_threads ? (i + 1)*base.f>length/base.num_threads 802 946 : base.f>length; 947 FLINT_ASSERT(cur_idx <= base.f>length); 803 948 next_idx = FLINT_MAX(next_idx, cur_idx); 804 949 next_idx = FLINT_MIN(next_idx, base.f>length); … … 845 990 { 846 991 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 849 994 fmpq_mpoly_init(res,ctx); 850 995 fmpq_mpoly_mul(res,pp,qq,ctx); 851 996 poly pres=convFlintMPSingP(res,ctx,r); 852 997 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); 855 1000 fmpq_mpoly_ctx_clear(ctx); 856 1001 p_Test(pres,r); … … 878 1023 { 879 1024 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 882 1027 fmpq_mpoly_init(res,ctx); 883 1028 fmpq_mpoly_divides(res,pp,qq,ctx); 884 1029 poly pres = convFlintMPSingP(res,ctx,r); 885 1030 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); 888 1033 fmpq_mpoly_ctx_clear(ctx); 889 1034 p_Test(pres,r); … … 934 1079 { 935 1080 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 938 1083 fmpq_mpoly_init(res,ctx); 939 1084 int ok=fmpq_mpoly_gcd(res,pp,qq,ctx); … … 959 1104 } 960 1105 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); 963 1108 fmpq_mpoly_ctx_clear(ctx); 964 1109 return pres;
Note: See TracChangeset
for help on using the changeset viewer.