Changeset 1e351b3 in git
- Timestamp:
- Jul 8, 2019, 2:07:08 PM (5 years ago)
- 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
- Files:
-
- 7 edited
Legend:
- Unmodified
- Added
- Removed
-
Singular/LIB/algebra.lib
r84f0d8 r1e351b3 278 278 //with the ring k[y(0),...,y(n)] is computed, the result is ker 279 279 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;");281 280 if (mp!="0") 282 281 { execute ("minpoly=number("+mp+");"); } -
Singular/LIB/gitfan.lib
r84f0d8 r1e351b3 266 266 newVars = newVars+string(var(k)); 267 267 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"); 270 270 ideal I = satstd(imap(origin,I)); 271 271 I = simplify(I,2+4+32); -
Singular/feOpt.cc
rbc4eb7 r1e351b3 23 23 #include "fehelp.h" 24 24 25 #ifdef HAVE_FLINT 26 #include <flint/flint.h> 27 #endif 25 28 26 29 const char SHORT_OPTS_STRING[] = "bdhpqstvxec:r:u:"; … … 309 312 } 310 313 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 311 331 default: 312 332 return NULL; -
Singular/feOptTab.h
rbc4eb7 r1e351b3 3 3 4 4 #define LONG_OPTION_RETURN 13 5 6 #ifdef HAVE_FLINT 7 #include <flint/flint.h> 8 #endif 5 9 6 10 // Define here which cmd-line options are recognized … … 144 148 "#threads", "maximal number of CPUs to use for threads", feOptInt, (void*)2, 0}, 145 149 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 146 157 147 158 {"MPport", required_argument, LONG_OPTION_RETURN, -
libpolys/polys/flint_mpoly.cc
rbc4eb7 r1e351b3 17 17 #include "polys/monomials/p_polys.h" 18 18 19 #include <vector> 20 21 /****** ring conversion ******/ 22 19 23 BOOLEAN convSingRFlintR(fmpq_mpoly_ctx_t ctx, const ring r) 20 24 { … … 37 41 } 38 42 43 BOOLEAN 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 39 68 void convSingPFlintMP(fmpq_mpoly_t res, fmpq_mpoly_ctx_t ctx, poly p, int lp, const ring r) 40 69 { 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); 44 71 ulong* exp=(ulong*)omAlloc((r->N+1)*sizeof(ulong)); 45 int i=0;46 72 while(p!=NULL) 47 73 { … … 51 77 #if SIZEOF_LONG==8 52 78 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); 54 80 #else 55 81 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); 57 83 #endif 58 fmpq_mpoly_set_term_coeff_fmpq(res,i,c,ctx);59 84 fmpq_clear(c); 60 i++;61 85 pIter(p); 62 86 } 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 64 88 omFreeSize(exp,(r->N+1)*sizeof(ulong)); 65 89 } … … 70 94 poly p=NULL; 71 95 ulong* exp=(ulong*)omAlloc0((r->N+1)*sizeof(ulong)); 96 fmpq_t c; 97 fmpq_init(c); 72 98 for(int i=d; i>=0; i--) 73 99 { 74 fmpq_t c;75 100 fmpq_mpoly_get_term_coeff_fmpq(c,f,i,ctx); 76 101 poly pp=p_Init(r); … … 84 109 p_Setm(pp,r); 85 110 number n=convFlintNSingN_QQ(c,r->cf); 86 //fmpq_clear(c); // LEAK?87 111 pSetCoeff0(pp,n); 88 112 pNext(pp)=p; 89 113 p=pp; 90 114 } 115 fmpq_clear(c); 116 omFreeSize(exp,(r->N+1)*sizeof(ulong)); 91 117 p_Test(p,r); 92 118 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;129 119 } 130 120 … … 150 140 p=pp; 151 141 } 142 omFreeSize(exp,(r->N+1)*sizeof(ulong)); 152 143 p_Test(p,r); 153 144 return p; … … 156 147 void convSingPFlintMP(nmod_mpoly_t res, nmod_mpoly_ctx_t ctx, poly p, int lp,const ring r) 157 148 { 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); 161 150 ulong* exp=(ulong*)omAlloc((r->N+1)*sizeof(ulong)); 162 int i=0;163 151 while(p!=NULL) 164 152 { … … 166 154 #if SIZEOF_LONG==8 167 155 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); 169 157 #else 170 158 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); 172 160 #endif 173 nmod_mpoly_set_term_coeff_ui(res,i,(ulong)n,ctx);174 i++;175 161 pIter(p); 176 162 } 177 //nmod_mpoly_print_pretty(res,r->names,ctx);PrintLn();178 163 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 170 static 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 190 class convert_sing_to_fmpq_mpoly_base 191 { 192 public: 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 228 class convert_sing_to_fmpq_mpoly_arg 229 { 230 public: 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 239 static 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 267 static 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!!!! */ 312 void 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 390 class convert_fmpq_mpoly_to_sing_base 391 { 392 public: 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 408 class convert_fmpq_mpoly_to_sing_arg 409 { 410 public: 411 poly poly_start, poly_end; 412 slong start_idx, end_idx; 413 convert_fmpq_mpoly_to_sing_base* base; 414 }; 415 416 static 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 451 poly 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 519 class convert_sing_to_nmod_mpoly_base 520 { 521 public: 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 557 class convert_sing_to_nmod_mpoly_arg 558 { 559 public: 560 slong start_idx, end_idx; 561 convert_sing_to_nmod_mpoly_base* base; 562 }; 563 564 static 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!!!! */ 600 void 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 657 class convert_nmod_mpoly_to_sing_base 658 { 659 public: 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 675 class convert_nmod_mpoly_to_sing_arg 676 { 677 public: 678 poly poly_start, poly_end; 679 slong start_idx, end_idx; 680 convert_nmod_mpoly_to_sing_base* base; 681 }; 682 683 static 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 714 poly 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 783 poly 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; 179 797 } 180 798 … … 184 802 convSingPFlintMP(pp,ctx,p,lp,r); 185 803 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); 188 805 nmod_mpoly_mul(res,pp,qq,ctx); 189 806 poly pres=convFlintMPSingP(res,ctx,r); … … 195 812 return pres; 196 813 } 814 197 815 poly Flint_GCD_MP(poly p,int lp,poly q,int lq,nmod_mpoly_ctx_t ctx,const ring r) 198 816 { … … 200 818 convSingPFlintMP(pp,ctx,p,lp,r); 201 819 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); 204 821 int ok=nmod_mpoly_gcd(res,pp,qq,ctx); 205 822 poly pres; … … 225 842 convSingPFlintMP(pp,ctx,p,lp,r); 226 843 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); 229 845 int ok=fmpq_mpoly_gcd(res,pp,qq,ctx); 230 846 poly pres; 231 847 if (ok) 232 848 { 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 } 233 859 pres=convFlintMPSingP(res,ctx,r); 234 860 p_Test(pres,r); … … 244 870 return pres; 245 871 } 872 246 873 #endif 247 874 #endif -
libpolys/polys/flintconv.cc
rbc4eb7 r1e351b3 226 226 int d=fmpq_poly_length(f); 227 227 poly p=NULL; 228 fmpq_t c; 229 fmpq_init(c); 228 230 for(int i=0; i<=d; i++) 229 231 { 230 fmpq_t c;231 fmpq_init(c);232 232 fmpq_poly_get_coeff_fmpq(c,f,i); 233 233 number n=convFlintNSingN(c,r->cf); … … 238 238 p=p_Add_q(p,pp,r); 239 239 } 240 fmpq_clear(c); 240 241 p_Test(p,r); 241 242 return p; -
libpolys/polys/monomials/ring.cc
r84f0d8 r1e351b3 5816 5816 R->names=(char**)omReallocSize(R->names,r->N*sizeof(char_ptr),R->N*sizeof(char_ptr)); 5817 5817 } 5818 elsei--;5818 i--; 5819 5819 } 5820 5820 R->block1[p]=R->N; 5821 rComplete(R );5821 rComplete(R,1); 5822 5822 return R; 5823 5823 }
Note: See TracChangeset
for help on using the changeset viewer.