Changeset 4b1d39 in git
 Timestamp:
 Jan 20, 2002, 12:44:48 PM (22 years ago)
 Branches:
 (u'spielwiese', 'd1b01e9d51ade4b46b745d3bada5c5f3696be3a8')
 Children:
 d704faffc1735628a9173c79d65d0b7331bcc012
 Parents:
 b8490e74a539f406ed4238becdea28121774e338
 Location:
 Singular
 Files:

 2 deleted
 2 edited
Legend:
 Unmodified
 Added
 Removed

Singular/extra.cc
rb8490e7 r4b1d39 2 2 * Computer Algebra System SINGULAR * 3 3 *****************************************/ 4 /* $Id: extra.cc,v 1.17 3 20020119 17:20:58 obachmanExp $ */4 /* $Id: extra.cc,v 1.174 20020120 11:44:47 Singular Exp $ */ 5 5 /* 6 6 * ABSTRACT: general interface to internals of Singular ("system" command) … … 88 88 #include "silink.h" 89 89 #include "walk.h" 90 91 #include "fast_maps.h" 90 92 91 93 /* … … 651 653 652 654 #include "mod_raw.h" 653 // #include "fast_maps.cc"654 655 655 656 static BOOLEAN jjEXTENDED_SYSTEM(leftv res, leftv h) … … 783 784 } 784 785 else 785 /*==================== ring debug ==================================*/786 if(strcmp(sys_cmd,"map")==0)787 {788 ring image_r = currRing;789 map theMap = (map)h>Data();790 ideal image_id = (ideal) theMap;791 ring map_r = IDRING(idroot>get(theMap>preimage, myynest));792 ideal map_id = IDIDEAL(map_r>idroot>get(h>Next()>Name(), myynest));793 #if 0794 ring src_r, dest_r;795 maMap_CreateRings(map_id, map_r, image_id, image_r, src_r, dest_r);796 mapoly mp;797 maideal mideal;798 799 maMap_CreatePolyIdeal(map_id, map_r, src_r, dest_r, mp, mideal);800 maPoly_Out(mp, src_r);801 #endif802 return FALSE;803 }804 else805 786 #endif 806 787 /*==================== mtrack ==================================*/ 
Singular/fast_maps.cc
rb8490e7 r4b1d39 7 7 * Author: obachman (Olaf Bachmann) 8 8 * Created: 02/01 9 * Version: $Id: fast_maps.cc,v 1.2 7 20020120 10:01:48 Singular Exp $9 * Version: $Id: fast_maps.cc,v 1.28 20020120 11:44:48 Singular Exp $ 10 10 *******************************************************************/ 11 11 #include "mod2.h" … … 30 30 */ 31 31 // return maximal monomial if max_map_monomials are substituted into pi_m 32 static poly maGetMaxExpP(poly* max_map_monomials, 33 int n_max_map_monomials, ring map_r, 32 static poly maGetMaxExpP(poly* max_map_monomials, 33 int n_max_map_monomials, ring map_r, 34 34 poly pi_m, ring pi_r) 35 35 { … … 38 38 Exponent_t e_i, e_j; 39 39 poly m_i, map_j = p_Init(map_r); 40 40 41 41 for (i=1; i <= n; i++) 42 42 { … … 59 59 60 60 // returns maximal exponent if map_id is applied to pi_id 61 static Exponent_t maGetMaxExp(ideal map_id, ring map_r, ideal pi_id, ring pi_r)61 static Exponent_t maGetMaxExp(ideal pi_id, ring pi_r, ideal map_id, ring map_r) 62 62 { 63 63 Exponent_t max=0; 64 64 poly* max_map_monomials = (poly*) omAlloc(IDELEMS(map_id)*sizeof(poly)); 65 65 poly max_pi_i, max_map_i; 66 66 67 67 int i, j; 68 68 for (i=0; i<IDELEMS(map_id); i++) … … 70 70 max_map_monomials[i] = p_GetMaxExpP(map_id>m[i], map_r); 71 71 } 72 72 73 73 for (i=0; i<IDELEMS(pi_id); i++) 74 74 { 75 75 max_pi_i = p_GetMaxExpP(pi_id>m[i], pi_r); 76 max_map_i = maGetMaxExpP(max_map_monomials, IDELEMS(map_id), map_r, 76 max_map_i = maGetMaxExpP(max_map_monomials, IDELEMS(map_id), map_r, 77 77 max_pi_i, pi_r); 78 78 Exponent_t temp = p_GetMaxExp(max_map_i, map_r); … … 102 102 p_wrp(monomial>dest, dest_r); 103 103 } 104 if (monomial>f1!=NULL) { printf(" f1:%x", (unsigned int)monomial>f1>src); 104 if (monomial>f1!=NULL) { printf(" f1:%x", (unsigned int)monomial>f1>src); 105 105 // p_wrp(monomial>f1>src, src_r); 106 107 if (monomial>f2!=NULL) { printf(" f2:%x",(unsigned int)monomial>f2>src); 106 } 107 if (monomial>f2!=NULL) { printf(" f2:%x",(unsigned int)monomial>f2>src); 108 108 // p_wrp(monomial>f2>src, src_r); 109 109 } 110 110 printf("\n"); 111 111 fflush(stdout); … … 161 161 } 162 162 while (next != NULL); 163 if (mp>dest != NULL) 163 if (mp>dest != NULL) 164 164 { 165 165 assume(dest_r != NULL); … … 182 182 return what; 183 183 } 184 184 185 185 mapoly iter = into; 186 186 mapoly prev = NULL; 187 187 188 188 Top: 189 189 p_LmCmpAction(iter>src, what>src, src_r, goto Equal, goto Greater, goto Smaller); 190 190 191 191 Greater: 192 192 if (iter>next == NULL) … … 198 198 iter = iter>next; 199 199 goto Top; 200 200 201 201 Smaller: 202 202 if (prev == NULL) … … 209 209 what>next = iter; 210 210 return what; 211 211 212 212 Equal: 213 213 iter>ref += what>ref; 214 214 macoeff coeff = what>coeff; 215 if (coeff != NULL) 215 if (coeff != NULL) 216 216 { 217 217 while (coeff>next != NULL) coeff = coeff>next; … … 232 232 { 233 233 poly next; 234 234 235 235 while (what != NULL) 236 236 { … … 254 254 int i; 255 255 mp = NULL; 256 256 257 257 for (i=0; i<mideal>n; i++) 258 258 { … … 260 260 { 261 261 mideal>buckets[i] = sBucketCreate(dest_r); 262 maPoly_InsertPoly(mp, 262 maPoly_InsertPoly(mp, 263 263 #ifdef PDEBUG 264 264 prShallowCopyR(map_id>m[i], map_r, src_r), … … 266 266 prShallowCopyR_NoSort(map_id>m[i], map_r, src_r), 267 267 #endif 268 src_r, 268 src_r, 269 269 mideal>buckets[i]); 270 270 } … … 273 273 274 274 275 void maMap_CreateRings(ideal map_id, ring map_r, 276 ideal image_id, ring image_r, 275 void maMap_CreateRings(ideal map_id, ring map_r, 276 ideal image_id, ring image_r, 277 277 ring &src_r, ring &dest_r, BOOLEAN &simple) 278 278 { … … 294 294 Exponent_t maxExp = maGetMaxExp(map_id, map_r, image_id, image_r); 295 295 if (maxExp <= 1) maxExp = 2; 296 else if (maxExp > (Exponent_t) image_r>bitmask) 296 else if (maxExp > (Exponent_t) image_r>bitmask) 297 297 maxExp = (Exponent_t) image_r>bitmask; 298 298 dest_r = rModifyRing_Simple(image_r, TRUE, TRUE, maxExp, simple); … … 314 314 *F misc . . . . . . . . . . . . . . . . . . . . . . . . . . . . misc stuff 315 315 */ 316 316 317 317 ideal maIdeal_2_Ideal(maideal m_id, ring dest_r) 318 318 { 319 319 ideal res = idInit(m_id>n, 1); 320 320 int l; 321 321 322 322 for (int i= 0; i < m_id>n; i++) 323 323 { … … 338 338 } 339 339 340 340 341 341 /******************************************************************************* 342 342 ** … … 351 351 BOOLEAN no_sort; 352 352 353 // construct rings we work in: 353 // construct rings we work in: 354 354 // src_r: Wp with Weights set to length of poly in image_id 355 355 // dest_r: Simple ring without degree ordering and short exponents 356 356 maMap_CreateRings(map_id, map_r, image_id, image_r, src_r, dest_r, no_sort); 357 357 358 358 // construct dest_id 359 359 if (dest_r != image_r) … … 366 366 maideal mideal; 367 367 maMap_CreatePolyIdeal(map_id, map_r, src_r, dest_r, mp, mideal); 368 368 369 369 if (TEST_OPT_PROT) 370 370 { … … 372 372 Print("map[%d:%d]{%d:", dest_r>bitmask, dest_r>ExpL_Size, length); 373 373 } 374 374 375 375 // do the optimization step 376 376 #if HAVE_MAP_OPTIMIZE > 0 … … 382 382 Print("%d}", length); 383 383 } 384 384 385 385 // do the actual evaluation 386 386 maPoly_Eval(mp, src_r, dest_id, dest_r, length); 387 387 if (TEST_OPT_PROT) Print("."); 388 388 389 389 // collect the results back into an ideal 390 390 ideal res_dest_id = maIdeal_2_Ideal(mideal, dest_r); … … 397 397 if (no_sort) 398 398 res_image_id = idrShallowCopyR_NoSort(res_dest_id, dest_r, image_r); 399 else 399 else 400 400 res_image_id = idrShallowCopyR(res_dest_id, dest_r, image_r); 401 401 id_ShallowDelete(&res_dest_id, dest_r); … … 403 403 else 404 404 res_image_id = res_dest_id; 405 405 406 406 if (TEST_OPT_PROT) Print("."); 407 407 // cleanup the rings … … 415 415 416 416 417 418 419 420 421 422 423 424 425 426 427 428 429 430 431 432 433 434 435 436 437 438 417 /********************************************************************** 418 * Evaluation stuff * 419 **********************************************************************/ 420 421 // substitute p everywhere the monomial occours, 422 // return the number of substitutions 423 static int maPoly_Substitute(macoeff c, poly p, ring dest_r) 424 { 425 // substitute the monomial: go through macoeff 426 int len=pLength(p); 427 int done=0; 428 while (c!=NULL) 429 { 430 done++; 431 poly t=pp_Mult_nn(p,c>n,dest_r); 432 sBucket_Add_p(c>bucket, t, len); 433 c=c>next; 434 } 435 return done; 436 } 437 438 static poly maPoly_EvalMon(poly src, ring src_r, poly* dest_id, ring dest_r) 439 { 440 int i; 441 int e; 442 poly p=NULL; 443 poly pp; 444 for(i=1;i<=src_r>N;i++) 445 { 446 e=p_GetExp(src,i,src_r); 447 if (e>0) 448 { 449 pp=dest_id[i1]; 450 if (pp==NULL) 451 { 452 p_Delete(&p,dest_r); 453 return NULL; 454 } 455 if ((p==NULL) /* && (e>0)*/) 456 { 457 p=p_Copy(pp /*dest_id[i1]*/,dest_r); 458 e; 459 } 460 while (e>0) 461 { 462 p=p_Mult_q(p,p_Copy(pp /*dest_id[i1]*/,dest_r),dest_r); 463 e; 464 } 465 } 466 } 467 if (p==NULL) p=p_ISet(1,dest_r); 468 return p; 469 } 470 471 void maPoly_Eval(mapoly root, ring src_r, ideal dest_id, ring dest_r, int total_cost) 472 { 473 // invert the list rooted at root: 474 if ((root!=NULL) && (root>next!=NULL)) 475 { 476 mapoly q=root>next; 477 mapoly qn; 478 root>next=NULL; 479 do 480 { 481 qn=q>next; 482 q>next=root; 483 root=q; 484 q=qn; 485 } 486 while (qn !=NULL); 487 } 488 489 total_cost /= 10; 490 int next_print_cost = total_cost; 491 492 // the evaluation  493 mapoly p=root; 494 int cost = 0; 495 496 while (p!=NULL) 497 { 498 // look at each mapoly: compute its value in >dest 499 assume (p>dest==NULL); 500 { 501 if ((p>f1!=NULL)&&(p>f2!=NULL)) 502 { 503 #if 0 504 printf("found prod:"); 505 p_wrp(p>src,src_r);printf("="); 506 p_wrp(p>f1>src,src_r);printf(" * "); 507 p_wrp(p>f2>src,src_r);printf("\n"); 508 #endif 509 poly f1=p>f1>dest; 510 poly f2=p>f2>dest; 511 if (p>f1>ref>0) f1=p_Copy(f1,dest_r); 512 else 513 { 514 // we own p>f1>dest now (in f1) 515 p>f1>dest=NULL; 516 } 517 if (p>f2>ref>0) f2=p_Copy(f2,dest_r); 518 else 519 { 520 // we own p>f2>dest now (in f2) 521 p>f2>dest=NULL; 522 } 523 maMonomial_Free(p>f1,src_r, dest_r); 524 maMonomial_Free(p>f2,src_r, dest_r); 525 p>dest=p_Mult_q(f1,f2,dest_r); 526 } /* factors : 2 */ 527 else 528 { 529 //printf("compute "); p_wrp(p>src,src_r);printf("\n"); 530 assume((p>f1==NULL) && (p>f2==NULL)); 531 //if(p>f1!=NULL) { printf(" but f1="); p_wrp(p>f1>src,src_r);printf("\n"); } 532 //if(p>f2!=NULL) { printf(" but f2="); p_wrp(p>f2>src,src_r);printf("\n"); } 533 // no factorization provided, use the classical method: 534 p>dest=maPoly_EvalMon(p>src,src_r,dest_id>m,dest_r); 535 //p_wrp(p>dest, dest_r); printf(" is p>dest\n"); 536 } 537 } /* p>dest==NULL */ 538 // substitute the monomial: go through macoeff 539 p>ref = maPoly_Substitute(p>coeff, p>dest, dest_r); 540 //printf("subst done\n"); 541 if (total_cost) 542 { 543 assume(TEST_OPT_PROT); 544 cost++; 545 if (cost > next_print_cost) 546 { 547 Print(""); 548 next_print_cost += total_cost; 549 } 550 } 551 552 mapoly pp=p; 553 p=p>next; 554 //p_wrp(pp>src, src_r); 555 if (pp>ref<=0) 556 { 557 //printf(" (%x) killed\n",pp); 558 maMonomial_Destroy(pp, src_r, dest_r); 559 } 560 //else 561 // printf(" (%x) not killed, ref=%d\n",pp,pp>ref); 562 } 563 } 564 565 566 /******************************************************************************* 567 ** 568 *F maEggt . . . . . . . . . . . . . . . . . . . . . . . . returns extended ggt 569 */ 570 // return NULL if deg(ggt(m1, m2)) < 2 571 // else return m = ggT(m1, m2) and q1, q2 such that m1 = q1*m m2 = q2*m 572 static poly maEggT(const poly m1, const poly m2, poly &q1, poly &q2,const ring r) 573 { 574 575 int i; 576 int dg = 0; 577 poly ggt = NULL; 578 q1 = p_Init(r); 579 q2 = p_Init(r); 580 ggt=p_Init(r); 581 582 for (i=1;i<=r>N;i++) { 583 Exponent_t e1 = p_GetExp(m1, i, r); 584 Exponent_t e2 = p_GetExp(m2, i, r); 585 if (e1 > 0 && e2 > 0){ 586 Exponent_t em = (e1 > e2 ? e2 : e1); 587 dg += em; 588 p_SetExp(ggt, i, em, r); 589 p_SetExp(q1, i, e1  em, r); 590 p_SetExp(q2, i, e2  em, r); 591 } 592 else { 593 p_SetExp(q1, i, e1, r); 594 p_SetExp(q2, i, e2, r); 595 } 596 } 597 if (dg>1) 598 { 599 p_Setm(ggt, r); 600 p_Setm(q1, r); 601 p_Setm(q2, r); 602 603 604 } 605 else { 606 p_LmFree(ggt, r); 607 p_LmFree(q1, r); 608 p_LmFree(q2, r); 609 ggt = NULL; 610 } 611 return ggt; 612 } 613 614 /******************************************************************** 615 ** * 616 * maFindBestggT * 617 * finds ggT with the highest cost * 618 *******************************************************************/ 619 620 static mapoly maFindBestggT(mapoly mp, mapoly & choice, mapoly & fp, mapoly & fq,const ring r) 621 { 622 int ggt_deg = 0; 623 poly p = mp>src; 624 mapoly iter = choice; 625 poly ggT = NULL; 626 fp = NULL; 627 fq = NULL; 628 poly fp_p=NULL; 629 poly fq_p=NULL; 630 choice=NULL; 631 while ((iter != NULL) && (pDeg(iter>src, r) > ggt_deg)) 632 { 633 // maMonomial_Out(iter, r, NULL); 634 poly q1, q2, q; 635 636 q = maEggT(p, iter>src, q1, q2,r); 637 if (q != NULL) 638 { 639 assume((q1!=NULL)&&(q2!=NULL)); 640 if (pDeg(q,r) > ggt_deg) 641 { 642 choice=iter; 643 if (ggT != NULL) 644 { 645 p_LmFree(ggT, r); 646 p_LmFree(fp_p, r); 647 p_LmFree(fq_p, r); 648 649 } 650 ggt_deg = pDeg(q, r); 651 ggT = q; 652 fp_p = q1; 653 fq_p = q2; 654 } 655 else 656 { 657 p_LmFree(q, r); 658 p_LmFree(q1, r); 659 p_LmFree(q2, r); 660 } 661 } 662 iter=iter>next; 663 } 664 if(ggT!=NULL){ 665 666 int dq =pTotaldegree(fq_p,r); 667 if (dq!=0){ 668 fq=maPoly_InsertMonomial(mp, fq_p, r, NULL); 669 fp=maPoly_InsertMonomial(mp, fp_p, r, NULL); 670 return maPoly_InsertMonomial(mp, ggT, r, NULL); 671 } 672 else { 673 fq=NULL; 674 p_LmFree(fq_p, r); 675 p_LmFree(ggT, r); 676 fp=maPoly_InsertMonomial(mp, fp_p, r, NULL); 677 choice>ref++; 678 return choice; 679 } 680 } 681 else { 682 return NULL; 683 } 684 685 686 } 687 688 /******************************************************************** 689 ** * 690 * maPoly_Optimize * 691 * adds and integrates subexpressions * 692 *******************************************************************/ 693 694 void maPoly_Optimize(mapoly mpoly, ring src_r){ 695 assume(mpoly!=NULL && mpoly>src!=NULL); 696 mapoly iter = mpoly; 697 mapoly choice; 698 mapoly ggT=NULL; 699 mapoly fp=NULL; 700 mapoly fq=NULL; 701 while (iter>next!=NULL){ 702 choice=iter>next; 703 if ((iter>f1==NULL)){ 704 ggT=maFindBestggT(iter, choice, fp, fq,src_r); 705 if (choice!=NULL){ 706 iter>f1=fp; 707 iter>f2=ggT; 708 if (fq!=NULL){ 709 ggT>ref++; 710 choice>f1=fq; 711 choice>f2=ggT; 712 } 713 } 714 715 } 716 iter=iter>next; 717 } 718 }
Note: See TracChangeset
for help on using the changeset viewer.