Changeset e92c6a in git
- Timestamp:
- Apr 19, 2013, 11:52:27 AM (11 years ago)
- Branches:
- (u'spielwiese', '5b153614cbc72bfa198d75b1e9e33dab2645d9fe')
- Children:
- 4d998e87e374f7d2eeb8255db038e35cecbd06f6
- Parents:
- 8659a91de7a2bca7e3533e0b2270e4ecd33c1809
- git-author:
- Stephan Oberfranz <oberfran@pipo.mathematik.uni-kl.de>2013-04-19 11:52:27+02:00
- git-committer:
- Hans Schoenemann <hannes@mathematik.uni-kl.de>2013-04-19 11:53:54+02:00
- Location:
- Singular
- Files:
-
- 2 added
- 3 edited
Legend:
- Unmodified
- Added
- Removed
-
Singular/extra.cc
r8659a9 re92c6a 2021 2021 } 2022 2022 else 2023 if (strcmp(sys_cmd, "Mrwalk") == 0) 2024 { // Random Walk 2025 if (h == NULL || h->Typ() != IDEAL_CMD || 2026 h->next == NULL || h->next->Typ() != INTVEC_CMD || 2027 h->next->next == NULL || h->next->next->Typ() != INTVEC_CMD || 2028 h->next->next->next == NULL || h->next->next->next->Typ() != INT_CMD || 2029 h->next->next->next->next == NULL || h->next->next->next->next->Typ() != INT_CMD) 2030 { 2031 Werror("system(\"Mrwalk\", ideal, intvec, intvec, int, int) expected"); 2032 return TRUE; 2033 } 2034 2035 if (((intvec*) h->next->Data())->length() != currRing->N && 2036 ((intvec*) h->next->next->Data())->length() != currRing->N ) 2037 { 2038 Werror("system(\"Mrwalk\" ...) intvecs not of length %d\n", 2039 currRing->N); 2040 return TRUE; 2041 } 2042 ideal arg1 = (ideal) h->Data(); 2043 intvec* arg2 = (intvec*) h->next->Data(); 2044 intvec* arg3 = (intvec*) h->next->next->Data(); 2045 int arg4 = (int)(long) h->next->next->next->Data(); 2046 int arg5 = (int)(long) h->next->next->next->next->Data(); 2047 2048 2049 ideal result = (ideal) Mrwalk(arg1, arg2, arg3, arg4, arg5); 2050 2051 res->rtyp = IDEAL_CMD; 2052 res->data = result; 2053 2054 return FALSE; 2055 } 2056 else 2023 2057 if (strcmp(sys_cmd, "MAltwalk1") == 0) 2024 2058 { … … 2119 2153 } 2120 2154 else 2155 if (strcmp(sys_cmd, "Mfrwalk") == 0) 2156 { 2157 if (h == NULL || h->Typ() != IDEAL_CMD || 2158 h->next == NULL || h->next->Typ() != INTVEC_CMD || 2159 h->next->next == NULL || h->next->next->Typ() != INTVEC_CMD || 2160 h->next->next->next == NULL || h->next->next->next->Typ() != INT_CMD) 2161 { 2162 Werror("system(\"Mfrwalk\", ideal, intvec, intvec, int) expected"); 2163 return TRUE; 2164 } 2165 2166 if (((intvec*) h->next->Data())->length() != currRing->N && 2167 ((intvec*) h->next->next->Data())->length() != currRing->N ) 2168 { 2169 Werror("system(\"Mfrwalk\" ...) intvecs not of length %d\n", 2170 currRing->N); 2171 return TRUE; 2172 } 2173 ideal arg1 = (ideal) h->Data(); 2174 intvec* arg2 = (intvec*) h->next->Data(); 2175 intvec* arg3 = (intvec*) h->next->next->Data(); 2176 int arg4 = (int)(long) h->next->next->next->Data(); 2177 2178 ideal result = (ideal) Mfrwalk(arg1, arg2, arg3, arg4); 2179 2180 res->rtyp = IDEAL_CMD; 2181 res->data = result; 2182 2183 return FALSE; 2184 } 2185 else 2121 2186 2122 2187 #ifdef TRAN_Orig … … 2213 2278 } 2214 2279 else 2280 if (strcmp(sys_cmd, "TranMrImprovwalk") == 0) 2281 { 2282 if (h == NULL || h->Typ() != IDEAL_CMD || 2283 h->next == NULL || h->next->Typ() != INTVEC_CMD || 2284 h->next->next == NULL || h->next->next->Typ() != INTVEC_CMD || 2285 h->next->next->next == NULL || h->next->next->next->Typ() != INT_CMD || 2286 h->next->next->next == NULL || h->next->next->next->next->Typ() != INT_CMD || 2287 h->next->next->next == NULL || h->next->next->next->next->next->Typ() != INT_CMD) 2288 { 2289 Werror("system(\"TranMrImprovwalk\", ideal, intvec, intvec) expected"); 2290 return TRUE; 2291 } 2292 2293 if (((intvec*) h->next->Data())->length() != currRing->N && 2294 ((intvec*) h->next->next->Data())->length() != currRing->N ) 2295 { 2296 Werror("system(\"TranMrImprovwalk\" ...) intvecs not of length %d\n", currRing->N); 2297 return TRUE; 2298 } 2299 ideal arg1 = (ideal) h->Data(); 2300 intvec* arg2 = (intvec*) h->next->Data(); 2301 intvec* arg3 = (intvec*) h->next->next->Data(); 2302 int arg4 = (int)(long) h->next->next->next->Data(); 2303 int arg5 = (int)(long) h->next->next->next->next->Data(); 2304 int arg6 = (int)(long) h->next->next->next->next->next->Data(); 2305 2306 ideal result = (ideal) TranMrImprovwalk(arg1, arg2, arg3, arg4, arg5, arg6); 2307 2308 res->rtyp = IDEAL_CMD; 2309 res->data = result; 2310 2311 return FALSE; 2312 } 2313 else 2314 2215 2315 #endif 2216 2316 /*================= Extended system call ========================*/ -
Singular/walk.cc
r8659a9 re92c6a 2 2 * Computer Algebra System SINGULAR * 3 3 *****************************************/ 4 /* $Id$ */ 4 5 /* 5 6 * ABSTRACT: Implementation of the Groebner walk … … 30 31 31 32 //#define TIME_TEST // print the used time of each subroutine 32 //#define ENDWALKS //print the size of the last omega-homogenoues Gr ï¿œbner basis33 //#define ENDWALKS //print the size of the last omega-homogenoues Groebner basis 33 34 34 35 /* includes */ … … 38 39 #include <time.h> 39 40 #include <sys/time.h> 41 #include <math.h> 40 42 #include <sys/stat.h> 41 43 #include <unistd.h> 42 #include <stdio.h>43 44 #include <float.h> 44 45 #include <misc/mylimits.h> … … 106 107 clock_t xftostd, xtextra, xftinput, to; 107 108 108 /* 2109 *utilities for TSet, LSet 110 */109 /**************************** 110 * utilities for TSet, LSet * 111 ****************************/ 111 112 inline static intset initec (int maxnr) 112 113 { … … 123 124 } 124 125 125 #if 0 /*unused*/ 126 /*2 127 *construct the set s from F u {P} 128 */ 126 /************************************ 127 * construct the set s from F u {P} * 128 ************************************/ 129 // unused 130 #if 0 129 131 static void initSSpecialCC (ideal F, ideal Q, ideal P,kStrategy strat) 130 132 { … … 141 143 strat->S=strat->Shdl->m; 142 144 143 / *- put polys into S -*/145 // - put polys into S - 144 146 if (Q!=NULL) 145 147 { … … 183 185 } 184 186 } 185 / *- put polys into S -*/187 //- put polys into S - 186 188 for (i=0; i<IDELEMS(F); i++) 187 189 { … … 276 278 #endif 277 279 278 /* 2279 *interreduces F 280 */280 /***************** 281 *interreduce F * 282 *****************/ 281 283 static ideal kInterRedCC(ideal F, ideal Q) 282 284 { … … 295 297 initBuchMoraCrit(strat); 296 298 strat->NotUsedAxis = (BOOLEAN *)omAlloc((currRing->N+1)*sizeof(BOOLEAN)); 297 for (j=currRing->N; j>0; j--) strat->NotUsedAxis[j] = TRUE; 299 for(j=currRing->N; j>0; j--) 300 { 301 strat->NotUsedAxis[j] = TRUE; 302 } 298 303 strat->enterS = enterSBba; 299 304 strat->posInT = posInT0; … … 305 310 strat->R = initR(); 306 311 strat->sevT = initsevT(); 307 if (rHasLocalOrMixedOrdering_currRing()) strat->honey = TRUE; 308 312 if(rHasLocalOrMixedOrdering_currRing()) 313 { 314 strat->honey = TRUE; 315 } 309 316 310 317 //initSCC(F,Q,strat); … … 316 323 tininitS=tininitS+clock()-timetmp;//22.01.02 317 324 */ 318 if (TEST_OPT_REDSB) 325 if(TEST_OPT_REDSB) 326 { 319 327 strat->noTailReduction=FALSE; 320 328 } 321 329 updateS(TRUE,strat); 322 330 323 if (TEST_OPT_REDSB && TEST_OPT_INTSTRATEGY) 331 if(TEST_OPT_REDSB && TEST_OPT_INTSTRATEGY) 332 { 324 333 completeReduce(strat); 334 } 325 335 pDelete(&strat->kHEdge); 326 336 omFreeSize((ADDRESS)strat->T,strat->tmax*sizeof(TObject)); … … 332 342 omfree(strat->R); 333 343 334 if (strat->fromQ) 335 { 336 for (j=0;j<IDELEMS(strat->Shdl);j++) 337 { 338 if(strat->fromQ[j]) pDelete(&strat->Shdl->m[j]); 344 if(strat->fromQ) 345 { 346 for(j=0; j<IDELEMS(strat->Shdl); j++) 347 { 348 if(strat->fromQ[j]) 349 { 350 pDelete(&strat->Shdl->m[j]); 351 } 339 352 } 340 353 omFreeSize((ADDRESS)strat->fromQ,IDELEMS(strat->Shdl)*sizeof(int)); 341 strat->fromQ =NULL;354 strat->fromQ = NULL; 342 355 } 343 356 // if (TEST_OPT_PROT) … … 353 366 } 354 367 355 #if 0 /*unused*/ 368 //unused 369 #if 0 356 370 static void TimeString(clock_t tinput, clock_t tostd, clock_t tif,clock_t tstd, 357 371 clock_t tlf,clock_t tred, clock_t tnw, int step) 358 372 { 359 373 double totm = ((double) (clock() - tinput))/1000000; 360 double ostd,mostd, mif, mstd, m extra, mlf, mred, mnw, mxif,mxstd,mxlf,mxred,mxnw,tot;361 374 double ostd,mostd, mif, mstd, mlf, mred, mnw, mxif,mxstd,mxlf,mxred,mxnw,tot; 375 // double mextra 362 376 Print("\n// total time = %.2f sec", totm); 363 377 Print("\n// tostd = %.2f sec = %.2f", ostd=((double) tostd)/1000000, … … 393 407 #endif 394 408 395 #if 0 /*unused*/ 409 //unused 410 #if 0 396 411 static void TimeStringFractal(clock_t tinput, clock_t tostd, clock_t tif,clock_t tstd, 397 412 clock_t textra, clock_t tlf,clock_t tred, clock_t tnw) … … 428 443 Print("\n// ideal %s = ", st); 429 444 for(i=0; i<nL-1; i++) 445 { 430 446 Print(" %s, ", pString(L->m[i])); 431 447 } 432 448 Print(" %s;", pString(L->m[nL-1])); 433 449 } 434 450 435 #if 0 /*unused*/ 451 //unused 452 //#if 0 436 453 static void headidString(ideal L, char* st) 437 454 { … … 440 457 Print("\n// ideal %s = ", st); 441 458 for(i=0; i<nL-1; i++) 459 { 442 460 Print(" %s, ", pString(pHead(L->m[i]))); 443 461 } 444 462 Print(" %s;", pString(pHead(L->m[nL-1]))); 445 463 } 446 #endif 447 448 #if 0 /*unused*/ 464 //#endif 465 466 //unused 467 //#if 0 449 468 static void idElements(ideal L, char* st) 450 469 { … … 454 473 Print("\n// #monoms of %s = ", st); 455 474 for(i=0; i<nL; i++) 475 { 456 476 K[i] = pLength(L->m[i]); 457 458 int j, nsame, nk=0; 477 } 478 int j, nsame; 479 // int nk=0; 459 480 for(i=0; i<nL; i++) 460 481 { … … 462 483 { 463 484 nsame = 1; 464 for(j=i+1; j<nL; j++){ 465 if(K[j]==K[i]){ 485 for(j=i+1; j<nL; j++) 486 { 487 if(K[j]==K[i]) 488 { 466 489 nsame ++; 467 490 K[j]=0; … … 469 492 } 470 493 if(nsame == 1) 494 { 471 495 Print("%d, ",K[i]); 496 } 472 497 else 498 { 473 499 Print("%d[%d], ", K[i], nsame); 500 } 474 501 } 475 502 } 476 503 omFree(K); 477 504 } 478 #endif505 //#endif 479 506 480 507 … … 482 509 { 483 510 int nV = iv->length()-1; 484 //Print("\n// vector %s = (", ch);485 511 Print("\n// intvec %s = ", ch); 486 512 487 513 for(int i=0; i<nV; i++) 514 { 488 515 Print("%d, ", (*iv)[i]); 516 } 489 517 Print("%d;", (*iv)[nV]); 490 518 } 491 519 492 #if 0 /*unused*/ 520 //unused 521 //#if 0 493 522 static void MivString(intvec* iva, intvec* ivb, intvec* ivc) 494 523 { … … 497 526 PrintS("\n// ("); 498 527 for(i=0; i<nV; i++) 528 { 499 529 Print("%d, ", (*iva)[i]); 530 } 500 531 Print("%d) ==> (", (*iva)[nV]); 501 502 532 for(i=0; i<nV; i++) 533 { 503 534 Print("%d, ", (*ivb)[i]); 535 } 504 536 Print("%d) := (", (*ivb)[nV]); 505 537 506 538 for(i=0; i<nV; i++) 539 { 507 540 Print("%d, ", (*ivc)[i]); 541 } 508 542 Print("%d)", (*ivc)[nV]); 509 543 } 510 #endif 511 512 // returns gcd of integers a and b 544 //#endif 545 546 /******************************************************************** 547 * returns gcd of integers a and b * 548 ********************************************************************/ 513 549 static inline long gcd(const long a, const long b) 514 550 { … … 516 552 //assume(p0 >= 0 && p1 >= 0); 517 553 if(p0 < 0) 554 { 518 555 p0 = -p0; 519 556 } 520 557 if(p1 < 0) 558 { 521 559 p1 = -p1; 522 560 } 523 561 while(p1 != 0) 524 562 { … … 530 568 } 531 569 532 // cancel gcd of integers zaehler and nenner 570 /********************************************* 571 * cancel gcd of integers zaehler and nenner * 572 *********************************************/ 533 573 static void cancel(mpz_t zaehler, mpz_t nenner) 534 574 { … … 544 584 } 545 585 546 #if 0 /*unused*/ 547 /* 23.07.03 */ 586 //unused 587 #if 0 548 588 static int isVectorNeg(intvec* omega) 549 589 { … … 551 591 552 592 for(i=omega->length(); i>=0; i--) 593 { 553 594 if((*omega)[i]<0) 595 { 554 596 return 1; 555 597 } 598 } 556 599 return 0; 557 600 } … … 587 630 if(mpz_cmp(zsum, sing_int)>0) 588 631 { 589 if(Overflow_Error == FALSE) { 632 if(Overflow_Error == FALSE) 633 { 590 634 PrintLn(); 591 635 PrintS("\n// ** OVERFLOW in \"MwalkInitialForm\": "); … … 596 640 } 597 641 642 mpz_clear(zmul); 643 mpz_clear(zvec); 644 mpz_clear(zsum); 645 mpz_clear(sing_int); 646 598 647 return wgrad; 599 648 } … … 613 662 614 663 if (maxtemp > max) 664 { 615 665 max = maxtemp; 666 } 616 667 } 617 668 return max; … … 620 671 621 672 /******************************************************************** 622 * compute a weight degree of a monomial p w.r.t. a weight_vector *673 * compute a weight degree of a monomial p w.r.t. a weight_vector * 623 674 ********************************************************************/ 624 675 static void MLmWeightedDegree_gmp(mpz_t result, const poly p, intvec* weight) … … 628 679 mpz_init_set_ui(sing_int, 2147483647); 629 680 630 int i , wgrad;681 int i; 631 682 632 683 mpz_t zmul; … … 657 708 { 658 709 if(g == NULL) 710 { 659 711 return NULL; 660 712 } 661 713 mpz_t max; mpz_init(max); 662 714 mpz_t maxtmp; mpz_init(maxtmp); … … 676 728 in_w_g = pHead(hg); 677 729 } 678 else { 730 else 731 { 679 732 if(mpz_cmp(maxtmp, max)==0) 733 { 680 734 in_w_g = pAdd(in_w_g, pHead(hg)); 735 } 681 736 } 682 737 } … … 696 751 697 752 for(i=nG-1; i>=0; i--) 753 { 698 754 Gomega->m[i] = MpolyInitialForm(G->m[i], ivw); 699 755 } 700 756 if(Overflow_Error == FALSE) 757 { 701 758 Overflow_Error = nError; 702 759 } 703 760 return Gomega; 704 761 } … … 706 763 /************************************************************************ 707 764 * test whether the weight vector iv is in the cone of the ideal G * 708 * i.e. are in(in_w(g)) =? in(g),for all g in G *765 * i.e. test whether in(in_w(g)) = in(g) for all g in G * 709 766 ************************************************************************/ 710 767 … … 731 788 { 732 789 pDelete(&mi); 733 734 790 if(Overflow_Error == FALSE) 791 { 735 792 Overflow_Error = nError; 736 793 } 737 794 return 0; 738 795 } … … 740 797 { 741 798 pDelete(&mi); 742 743 799 if(Overflow_Error == FALSE) 800 { 744 801 Overflow_Error = nError; 745 802 } 746 803 return 0; 747 804 } 748 749 805 pDelete(&mi); 750 806 } 751 807 752 808 if(Overflow_Error == FALSE) 809 { 753 810 Overflow_Error = nError; 754 811 } 755 812 return 1; 756 813 } 757 814 758 759 //compute a least common multiple of two integers 815 /*************************************************** 816 * compute a least common multiple of two integers * 817 ***************************************************/ 760 818 static inline long Mlcm(long &i1, long &i2) 761 819 { … … 775 833 776 834 for(i=n-1; i>=0; i--) 835 { 777 836 result += (*a)[i] * (*b)[i]; 778 837 } 779 838 return result; 780 839 } 781 840 782 841 /***************************************************** 842 * Substract two given intvecs componentwise * 843 *****************************************************/ 783 844 static intvec* MivSub(intvec* a, intvec* b) 784 845 { … … 788 849 789 850 for(i=n-1; i>=0; i--) 851 { 790 852 (*result)[i] = (*a)[i] - (*b)[i]; 791 853 } 792 854 return result; 793 855 } 794 856 795 /** 21.10.00*******************************************857 /***************************************************** 796 858 * return the "intvec" lead exponent of a polynomial * 797 859 *****************************************************/ … … 802 864 803 865 for(i=nR-1; i>=0; i--) 866 { 804 867 (*result)[i] = pGetExp(f,i+1); 805 868 } 806 869 return result; 807 870 } 808 871 809 /* return 1, if two given intvecs are the same, otherwise 0*/ 872 /***************************************************** 873 * Compare two given intvecs and return 1, if they * 874 * are the same, otherwise 0 * 875 *****************************************************/ 810 876 int MivSame(intvec* u , intvec* v) 811 877 { … … 815 881 816 882 for (i=0; i<niv; i++) 883 { 817 884 if ((*u)[i] != (*v)[i]) 885 { 818 886 return 0; 819 887 } 888 } 820 889 return 1; 821 890 } 822 891 892 /****************************************************** 893 * Compare 3 given intvecs and return 0, if the first * 894 * and the second are the same. Return 1, if the * 895 * the second and the third are the same, otherwise 2 * 896 ******************************************************/ 823 897 int M3ivSame(intvec* temp, intvec* u , intvec* v) 824 898 { … … 826 900 827 901 if((MivSame(temp, u)) == 1) 902 { 828 903 return 0; 829 904 } 830 905 if((MivSame(temp, v)) == 1) 906 { 831 907 return 1; 832 908 } 833 909 return 2; 834 910 } 835 911 836 837 /* compute a Groebner basis of an ideal */ 912 /***************************************************** 913 * compute a Groebner basis of an ideal * 914 *****************************************************/ 838 915 static ideal MstdCC(ideal G) 839 916 { … … 848 925 } 849 926 850 851 /* compute a Groebner basis of a homogenoues ideal */ 927 /***************************************************** 928 * compute a Groebner basis of an homogeneous ideal * 929 *****************************************************/ 852 930 static ideal MstdhomCC(ideal G) 853 931 { … … 868 946 intvec* MivMatrixOrder(intvec* iv) 869 947 { 870 int i,j, nR = iv->length(); 948 int i, nR = iv->length(); 949 871 950 intvec* ivm = new intvec(nR*nR); 872 951 873 952 for(i=0; i<nR; i++) 953 { 874 954 (*ivm)[i] = (*iv)[i]; 875 955 } 876 956 for(i=1; i<nR; i++) 957 { 877 958 (*ivm)[i*nR+i-1] = 1; 878 959 } 879 960 return ivm; 880 961 } 881 962 882 /* return intvec = (1, ..., 1) */ 963 /******************************* 964 * return intvec = (1, ..., 1) * 965 *******************************/ 883 966 intvec* Mivdp(int nR) 884 967 { … … 887 970 888 971 for(i=nR-1; i>=0; i--) 972 { 889 973 (*ivm)[i] = 1; 890 974 } 891 975 return ivm; 892 976 } 893 977 894 /* return intvvec = (1,0, ..., 0) */ 978 /********************************** 979 * return intvvec = (1,0, ..., 0) * 980 **********************************/ 895 981 intvec* Mivlp(int nR) 896 982 { 897 int i;898 983 intvec* ivm = new intvec(nR); 899 984 (*ivm)[0] = 1; … … 902 987 } 903 988 904 /**** 28.10.02 print the max total degree and the max coefficient of G***/ 905 #if 0 /*unused*/ 989 //unused 990 /***************************************************************************** 991 * print the max total degree and the max coefficient of G * 992 *****************************************************************************/ 993 #if 0 906 994 static void checkComplexity(ideal G, char* cG) 907 995 { 908 996 int nV = currRing->N; 909 997 int nG = IDELEMS(G); 910 intvec* ivUnit = Mivdp(nV); //19.02911 int i, j,tmpdeg, maxdeg=0;998 intvec* ivUnit = Mivdp(nV); 999 int i, tmpdeg, maxdeg=0; 912 1000 number tmpcoeff , maxcoeff=currRing->cf->nNULL; 913 1001 poly p; … … 915 1003 { 916 1004 tmpdeg = MwalkWeightDegree(G->m[i], ivUnit); 917 if (tmpdeg > maxdeg ) 1005 if(tmpdeg > maxdeg ) 1006 { 918 1007 maxdeg = tmpdeg; 1008 } 919 1009 } 920 1010 … … 927 1017 tmpcoeff = pGetCoeff(p); 928 1018 if(nGreater(tmpcoeff,maxcoeff)) 1019 { 929 1020 maxcoeff = nCopy(tmpcoeff); 1021 } 930 1022 pIter(p); 931 1023 } … … 934 1026 p = pNSet(maxcoeff); 935 1027 char* pStr = pString(p); 1028 delete ivUnit; 936 1029 Print("// max total degree of %s = %d\n",cG, maxdeg); 937 1030 Print("// max coefficient of %s = %s", cG, pStr);//ing(p)); … … 940 1033 } 941 1034 #endif 942 943 1035 944 1036 /***************************************************************************** … … 953 1045 * degree which smaller than the numbers of variables * 954 1046 ******************************************************************************/ 955 /* ivtarget is a matrix order of a degree reverse lex. order */956 1047 intvec* MPertVectors(ideal G, intvec* ivtarget, int pdeg) 957 1048 { 1049 // ivtarget is a matrix order of a degree reverse lex. order 958 1050 int nV = currRing->N; 959 1051 //assume(pdeg <= nV && pdeg >= 0); … … 963 1055 964 1056 965 // Checkingthat the perturbed degree is valid1057 // Check that the perturbed degree is valid 966 1058 if(pdeg > nV || pdeg <= 0) 967 1059 { … … 972 1064 973 1065 if(pdeg == 1) 1066 { 974 1067 return ivtarget; 975 976 mpz_t *pert_vector=(mpz_t*)omAlloc(nV*sizeof(mpz_t)); 1068 } 1069 mpz_t *pert_vector = (mpz_t*)omAlloc(nV*sizeof(mpz_t)); 1070 //mpz_t *pert_vector1 = (mpz_t*)omAlloc(nV*sizeof(mpz_t)); 977 1071 978 1072 for(i=0; i<nV; i++) 1073 { 979 1074 mpz_init_set_si(pert_vector[i], (*ivtarget)[i]); 980 981 1075 // mpz_init_set_si(pert_vector1[i], (*ivtarget)[i]); 1076 } 982 1077 // Calculate max1 = Max(A2)+Max(A3)+...+Max(Apdeg), 983 1078 // where the Ai are the i-te rows of the matrix target_ord. 984 985 1079 int ntemp, maxAi, maxA=0; 986 1080 for(i=1; i<pdeg; i++) 987 1081 { 988 1082 maxAi = (*ivtarget)[i*nV]; 989 if(maxAi<0) maxAi = -maxAi; 990 1083 if(maxAi<0) 1084 { 1085 maxAi = -maxAi; 1086 } 991 1087 for(j=i*nV+1; j<(i+1)*nV; j++) 992 1088 { 993 1089 ntemp = (*ivtarget)[j]; 994 if(ntemp < 0) ntemp = -ntemp; 995 1090 if(ntemp < 0) 1091 { 1092 ntemp = -ntemp; 1093 } 996 1094 if(ntemp > maxAi) 1095 { 997 1096 maxAi = ntemp; 1097 } 998 1098 } 999 1099 maxA += maxAi; … … 1011 1111 for(i=nG-1; i>=0; i--) 1012 1112 { 1013 mpz_set_ui(maxdeg, MwalkWeightDegree(G->m[i], ivUnit)); 1014 if (mpz_cmp(maxdeg, tot_deg) > 0 ) 1015 mpz_set(tot_deg, maxdeg); 1113 mpz_set_ui(maxdeg, MwalkWeightDegree(G->m[i], ivUnit)); 1114 if (mpz_cmp(maxdeg, tot_deg) > 0 ) 1115 { 1116 mpz_set(tot_deg, maxdeg); 1117 } 1016 1118 } 1017 1119 … … 1021 1123 1022 1124 1023 // xx1.06.02takes "small" inveps1125 // takes "small" inveps 1024 1126 #ifdef INVEPS_SMALL_IN_MPERTVECTOR 1025 1127 if(mpz_cmp_ui(inveps, pdeg)>0 && pdeg > 3) 1026 1128 { 1027 /* 1028 Print("\n// choose the\"small\" inverse epsilon := %d / %d = ", 1029 mpz_get_si(inveps), pdeg); 1030 */ 1129 // Print("\n// choose the\"small\" inverse epsilon := %d / %d = ", mpz_get_si(inveps), pdeg); 1031 1130 mpz_fdiv_q_ui(inveps, inveps, pdeg); 1032 // mpz_out_str(stdout, 10, inveps);1131 // mpz_out_str(stdout, 10, inveps); 1033 1132 } 1034 1133 #else 1035 // PrintS("\n// the \"big\" inverse epsilon: ");1134 // PrintS("\n// the \"big\" inverse epsilon: "); 1036 1135 mpz_out_str(stdout, 10, inveps); 1037 1136 #endif … … 1039 1138 // pert(A1) = inveps^(pdeg-1)*A1 + inveps^(pdeg-2)*A2+...+A_pdeg, 1040 1139 // pert_vector := A1 1041 for ( i=1; i < pdeg; i++ ) 1140 for( i=1; i < pdeg; i++ ) 1141 { 1042 1142 for(j=0; j<nV; j++) 1043 1143 { 1044 1144 mpz_mul(pert_vector[j], pert_vector[j], inveps); 1045 1145 if((*ivtarget)[i*nV+j]<0) 1146 { 1046 1147 mpz_sub_ui(pert_vector[j], pert_vector[j],-(*ivtarget)[i*nV+j]); 1148 } 1047 1149 else 1150 { 1048 1151 mpz_add_ui(pert_vector[j], pert_vector[j],(*ivtarget)[i*nV+j]); 1049 } 1050 1152 } 1153 } 1154 } 1051 1155 mpz_t ztemp; 1052 1156 mpz_init(ztemp); … … 1056 1160 mpz_gcd(ztemp, ztemp, pert_vector[i]); 1057 1161 if(mpz_cmp_si(ztemp, 1) == 0) 1162 { 1058 1163 break; 1164 } 1059 1165 } 1060 1166 if(mpz_cmp_si(ztemp, 1) != 0) 1167 { 1061 1168 for(i=0; i<nV; i++) 1169 { 1062 1170 mpz_divexact(pert_vector[i], pert_vector[i], ztemp); 1063 1171 } 1172 } 1173 1174 intvec *pert_vector1= new intvec(nV); 1175 j = 0; 1176 for(i=0; i<nV; i++) 1177 { 1178 (* pert_vector1)[i] = mpz_get_si(pert_vector[i]); 1179 (* pert_vector1)[i] = 0.1*(* pert_vector1)[i]; 1180 (* pert_vector1)[i] = floor((* pert_vector1)[i] + 0.5); 1181 if((* pert_vector1)[i] == 0) 1182 { 1183 j++; 1184 } 1185 } 1186 if(j > nV - 1) 1187 { 1188 // Print("\n// MPertVectors: geaenderter vector gleich Null! \n"); 1189 delete pert_vector1; 1190 goto CHECK_OVERFLOW; 1191 } 1192 1193 // check that the perturbed weight vector lies in the Groebner cone 1194 if(test_w_in_ConeCC(G,pert_vector1) != 0) 1195 { 1196 // Print("\n// MPertVectors: geaenderter vector liegt in Groebnerkegel! \n"); 1197 for(i=0; i<nV; i++) 1198 { 1199 mpz_set_si(pert_vector[i], (*pert_vector1)[i]); 1200 } 1201 } 1202 else 1203 { 1204 //Print("\n// MpertVectors: geaenderter vector liegt nicht in Groebnerkegel! \n"); 1205 } 1206 delete pert_vector1; 1207 1208 CHECK_OVERFLOW: 1064 1209 intvec* result = new intvec(nV); 1210 1065 1211 /* 2147483647 is max. integer representation in SINGULAR */ 1066 1212 mpz_t sing_int; … … 1071 1217 { 1072 1218 (*result)[i] = mpz_get_si(pert_vector[i]); 1073 1074 1219 if(mpz_cmp(pert_vector[i], sing_int)>=0) 1075 1220 { … … 1095 1240 mpz_clear(sing_int); 1096 1241 omFree(pert_vector); 1242 //omFree(pert_vector1); 1243 mpz_clear(tot_deg); 1244 mpz_clear(maxdeg); 1245 mpz_clear(inveps); 1246 1247 rComplete(currRing); 1248 for(j=0; j<IDELEMS(G); j++) 1249 { 1250 poly p=G->m[j]; 1251 while(p!=NULL) 1252 { 1253 p_Setm(p,currRing); pIter(p); 1254 } 1255 } 1097 1256 return result; 1098 1257 } 1099 1258 1100 1101 /* ivtarget is a matrix order of the lex. order */ 1259 /***************************************************************************** 1260 * The following procedure returns * 1261 * Pert(A1) = 1/eps^(pdeg-1)*A_1 + 1/eps^(pdeg-2)*A_2+...+A_pdeg, * 1262 * where the A_i are the i-th rows of the matrix target_ord and * 1263 * 1/eps > deg(p)*(max(A_2) + max(A_3)+...+max(A_pdeg)) * 1264 *****************************************************************************/ 1102 1265 intvec* MPertVectorslp(ideal G, intvec* ivtarget, int pdeg) 1103 1266 { 1267 // ivtarget is a matrix order of the lex. order 1104 1268 int nV = currRing->N; 1105 1269 //assume(pdeg <= nV && pdeg >= 0); … … 1115 1279 } 1116 1280 for(i=0; i<nV; i++) 1281 { 1117 1282 (*pert_vector)[i]=(*ivtarget)[i]; 1118 1283 } 1119 1284 if(pdeg == 1) 1285 { 1120 1286 return pert_vector; 1121 1287 } 1122 1288 // Calculate max1 = Max(A2)+Max(A3)+...+Max(Apdeg), 1123 1289 // where the Ai are the i-te rows of the matrix target_ord. … … 1130 1296 ntemp = (*ivtarget)[j]; 1131 1297 if(ntemp > maxAi) 1298 { 1132 1299 maxAi = ntemp; 1300 } 1133 1301 } 1134 1302 maxA += maxAi; … … 1141 1309 for(i=nG-1; i>=0; i--) 1142 1310 { 1143 // maxdeg = pTotaldegree(G->m[i], currRing); //it's wrong for ex1,2,rose1311 // maxdeg = pTotaldegree(G->m[i], currRing); //it's wrong for ex1,2,rose 1144 1312 maxdeg = MwalkWeightDegree(G->m[i], ivUnit); 1145 1313 if (maxdeg > tot_deg ) 1314 { 1146 1315 tot_deg = maxdeg; 1316 } 1147 1317 } 1148 1318 delete ivUnit; … … 1150 1320 inveps = (tot_deg * maxA) + 1; 1151 1321 1152 //9.10.011153 1322 #ifdef INVEPS_SMALL_IN_FRACTAL 1154 /* 1155 Print("\n// choose the\"small\" inverse epsilon := %d / %d = ", 1156 inveps, pdeg); 1157 */ 1323 // Print("\n// choose the\"small\" inverse epsilon := %d / %d = ", inveps, pdeg); 1158 1324 if(inveps > pdeg && pdeg > 3) 1325 { 1159 1326 inveps = inveps / pdeg; 1160 1161 // Print(" %d", inveps);1327 } 1328 // Print(" %d", inveps); 1162 1329 #else 1163 1330 PrintS("\n// the \"big\" inverse epsilon %d", inveps); 1164 1331 #endif 1165 1332 1166 // Pert(A1) = inveps^(pdeg-1)*A1 + inveps^(pdeg-2)*A2+...+A_pdeg ,1333 // Pert(A1) = inveps^(pdeg-1)*A1 + inveps^(pdeg-2)*A2+...+A_pdeg 1167 1334 for ( i=1; i < pdeg; i++ ) 1335 { 1168 1336 for(j=0; j<nV; j++) 1337 { 1169 1338 (*pert_vector)[j] = inveps*((*pert_vector)[j]) + (*ivtarget)[i*nV+j]; 1339 } 1340 } 1170 1341 1171 1342 int temp = (*pert_vector)[0]; … … 1174 1345 temp = gcd(temp, (*pert_vector)[i]); 1175 1346 if(temp == 1) 1347 { 1176 1348 break; 1349 } 1177 1350 } 1178 1351 if(temp != 1) 1352 { 1179 1353 for(i=0; i<nV; i++) 1354 { 1180 1355 (*pert_vector)[i] = (*pert_vector)[i] / temp; 1356 } 1357 } 1181 1358 1182 1359 intvec* result = pert_vector; 1183 pert_vector = NULL;1360 delete pert_vector; 1184 1361 return result; 1185 1362 } 1186 1363 1187 1188 /* define a lexicographic order matrix as intvec */ 1364 /***************************************************************************** 1365 * define a lexicographic order matrix as intvec * 1366 *****************************************************************************/ 1189 1367 intvec* MivMatrixOrderlp(int nV) 1190 1368 { … … 1193 1371 1194 1372 for(i=0; i<nV; i++) 1373 { 1195 1374 (*ivM)[i*nV + i] = 1; 1196 1375 } 1197 1376 return(ivM); 1198 1377 } 1199 1378 1200 /* define a rlex order (dp) matrix as intvec */ 1379 1380 /***************************************************************************** 1381 * define a reverse lexicographic order (dp) matrix as intvec * 1382 *****************************************************************************/ 1201 1383 intvec* MivMatrixOrderdp(int nV) 1202 1384 { … … 1205 1387 1206 1388 for(i=0; i<nV; i++) 1389 { 1207 1390 (*ivM)[i] = 1; 1208 1391 } 1209 1392 for(i=1; i<nV; i++) 1393 { 1210 1394 (*ivM)[(i+1)*nV - i] = -1; 1211 1395 } 1212 1396 return(ivM); 1213 1397 } 1214 1398 1215 //creates an intvec of the monomial order Wp(ivstart) 1399 /***************************************************************************** 1400 * creates an intvec of the monomial order Wp(ivstart) * 1401 *****************************************************************************/ 1216 1402 intvec* MivWeightOrderlp(intvec* ivstart) 1217 1403 { … … 1221 1407 1222 1408 for(i=0; i<nV; i++) 1409 { 1223 1410 (*ivM)[i] = (*ivstart)[i]; 1224 1411 } 1225 1412 for(i=1; i<nV; i++) 1413 { 1226 1414 (*ivM)[i*nV + i-1] = 1; 1227 1415 } 1228 1416 return(ivM); 1229 1417 } 1230 1418 1419 /***************************************************************************** 1420 * creates an intvec of the monomial order dp(ivstart) * 1421 *****************************************************************************/ 1231 1422 intvec* MivWeightOrderdp(intvec* ivstart) 1232 1423 { … … 1236 1427 1237 1428 for(i=0; i<nV; i++) 1429 { 1238 1430 (*ivM)[i] = (*ivstart)[i]; 1239 1431 } 1240 1432 for(i=0; i<nV; i++) 1433 { 1241 1434 (*ivM)[nV+i] = 1; 1242 1435 } 1243 1436 for(i=2; i<nV; i++) 1437 { 1244 1438 (*ivM)[(i+1)*nV - i] = -1; 1245 1439 } 1246 1440 return(ivM); 1247 1441 } 1248 1442 1249 #if 0 /*unused*/ 1443 //unused 1444 #if 0 1250 1445 static intvec* MatrixOrderdp(int nV) 1251 1446 { … … 1254 1449 1255 1450 for(i=0; i<nV; i++) 1451 { 1256 1452 (*ivM)[i] = 1; 1257 1453 } 1258 1454 for(i=1; i<nV; i++) 1455 { 1259 1456 (*ivM)[(i+1)*nV - i] = -1; 1260 1457 } 1261 1458 return(ivM); 1262 1459 } … … 1267 1464 int i; 1268 1465 intvec* ivM = new intvec(nV); 1269 1270 1466 for(i=nV-1; i>=0; i--) 1467 { 1271 1468 (*ivM)[i] = 1; 1272 1469 } 1273 1470 return(ivM); 1274 1471 } … … 1279 1476 *************************************************************************/ 1280 1477 int Xnlev; 1281 1282 1478 intvec* Mfpertvector(ideal G, intvec* ivtarget) 1283 1479 { … … 1286 1482 int niv = nV*nV; 1287 1483 1288 // Calculate max1 = Max(A2) + Max(A3) + ... + Max(AnV), 1484 1485 // Calculate maxA = Max(A2) + Max(A3) + ... + Max(AnV), 1289 1486 // where the Ai are the i-te rows of the matrix 'targer_ord'. 1290 1487 int ntemp, maxAi, maxA=0; … … 1292 1489 { 1293 1490 maxAi = (*ivtarget)[i*nV]; 1294 if(maxAi<0) maxAi = -maxAi; 1295 1491 if(maxAi<0) 1492 { 1493 maxAi = -maxAi; 1494 } 1296 1495 for(j=i*nV+1; j<(i+1)*nV; j++) 1297 1496 { 1298 1497 ntemp = (*ivtarget)[j]; 1299 if(ntemp < 0) ntemp = -ntemp; 1300 1498 if(ntemp < 0) 1499 { 1500 ntemp = -ntemp; 1501 } 1301 1502 if(ntemp > maxAi) 1503 { 1302 1504 maxAi = ntemp; 1505 } 1303 1506 } 1304 1507 maxA = maxA + maxAi; … … 1306 1509 intvec* ivUnit = Mivdp(nV); 1307 1510 1308 // Calculate inveps = 1/eps, where 1/eps > deg(p)*max 1for all p in G.1511 // Calculate inveps = 1/eps, where 1/eps > deg(p)*maxA for all p in G. 1309 1512 mpz_t tot_deg; mpz_init(tot_deg); 1310 1513 mpz_t maxdeg; mpz_init(maxdeg); … … 1316 1519 mpz_set_ui(maxdeg, MwalkWeightDegree(G->m[i], ivUnit)); 1317 1520 if (mpz_cmp(maxdeg, tot_deg) > 0 ) 1521 { 1318 1522 mpz_set(tot_deg, maxdeg); 1523 } 1319 1524 } 1320 1525 … … 1324 1529 mpz_add_ui(inveps, inveps, 1); 1325 1530 1326 // xx1.06.02takes "small" inveps1531 // takes "small" inveps 1327 1532 #ifdef INVEPS_SMALL_IN_FRACTAL 1328 1533 if(mpz_cmp_ui(inveps, nV)>0 && nV > 3) 1534 { 1329 1535 mpz_cdiv_q_ui(inveps, inveps, nV); 1330 1536 } 1331 1537 //PrintS("\n// choose the \"small\" inverse epsilon!"); 1332 1538 #endif … … 1338 1544 mpz_t *pert_vector=(mpz_t *)omAlloc(niv*sizeof(mpz_t)); 1339 1545 1340 for(i=0; i <nV; i++)1546 for(i=0; i < nV; i++) 1341 1547 { 1342 1548 mpz_init_set_si(ivtemp[i], (*ivtarget)[i]); … … 1345 1551 1346 1552 mpz_t ztmp; mpz_init(ztmp); 1347 BOOLEAN isneg = FALSE;1553 // BOOLEAN isneg = FALSE; 1348 1554 1349 1555 for(i=1; i<nV; i++) … … 1352 1558 { 1353 1559 mpz_mul(ztmp, inveps, ivtemp[j]); 1354 1355 1560 if((*ivtarget)[i*nV+j]<0) 1561 { 1356 1562 mpz_sub_ui(ivtemp[j], ztmp, -(*ivtarget)[i*nV+j]); 1563 } 1357 1564 else 1565 { 1358 1566 mpz_add_ui(ivtemp[j], ztmp,(*ivtarget)[i*nV+j]); 1567 } 1359 1568 } 1360 1569 1361 1570 for(j=0; j<nV; j++) 1571 { 1362 1572 mpz_init_set(pert_vector[i*nV+j],ivtemp[j]); 1573 } 1363 1574 } 1364 1575 … … 1368 1579 1369 1580 intvec* result = new intvec(niv); 1581 intvec* result1 = new intvec(niv); 1370 1582 BOOLEAN nflow = FALSE; 1371 1583 … … 1376 1588 mpz_gcd(ztmp, ztmp, pert_vector[i]); 1377 1589 if(mpz_cmp_si(ztmp, 1)==0) 1590 { 1378 1591 break; 1592 } 1379 1593 } 1380 1594 … … 1383 1597 mpz_divexact(pert_vector[i], pert_vector[i], ztmp); 1384 1598 (* result)[i] = mpz_get_si(pert_vector[i]); 1385 1599 } 1600 1601 j = 0; 1602 for(i=0; i<nV; i++) 1603 { 1604 (* result1)[i] = mpz_get_si(pert_vector[i]); 1605 (* result1)[i] = 0.1*(* result1)[i]; 1606 (* result1)[i] = floor((* result1)[i] + 0.5); 1607 if((* result1)[i] == 0) 1608 { 1609 j++; 1610 } 1611 } 1612 if(j > nV - 1) 1613 { 1614 // Print("\n// MfPertwalk: geaenderter vector gleich Null! \n"); 1615 delete result1; 1616 goto CHECK_OVERFLOW; 1617 } 1618 1619 // check that the perturbed weight vector lies in the Groebner cone 1620 if(test_w_in_ConeCC(G,result1) != 0) 1621 { 1622 // Print("\n// MfPertwalk: geaenderter vector liegt in Groebnerkegel! \n"); 1623 delete result; 1624 result = result1; 1625 for(i=0; i<nV; i++) 1626 { 1627 mpz_set_si(pert_vector[i], (*result1)[i]); 1628 } 1629 } 1630 else 1631 { 1632 delete result1; 1633 // Print("\n// Mfpertwalk: geaenderter vector liegt nicht in Groebnerkegel! \n"); 1634 } 1635 1636 CHECK_OVERFLOW: 1637 1638 for(i=0; i<niv; i++) 1639 { 1386 1640 if(mpz_cmp(pert_vector[i], sing_int)>0) 1641 { 1387 1642 if(nflow == FALSE) 1388 1643 { … … 1390 1645 nflow = TRUE; 1391 1646 Overflow_Error = TRUE; 1392 1393 1647 Print("\n// Xlev = %d and the %d-th element is", Xnlev, i+1); 1394 1648 PrintS("\n// ** OVERFLOW in \"Mfpertvector\": "); … … 1397 1651 Print("\n// So vector[%d] := %d is wrong!!", i+1, (*result)[i]); 1398 1652 } 1399 }1400 1653 } 1654 } 1401 1655 if(Overflow_Error == TRUE) 1656 { 1402 1657 ivString(result, "new_vector"); 1403 1658 } 1404 1659 omFree(pert_vector); 1405 1660 omFree(ivtemp); 1406 1661 mpz_clear(ztmp); 1407 1662 mpz_clear(tot_deg); 1663 mpz_clear(maxdeg); 1664 mpz_clear(inveps); 1665 mpz_clear(sing_int); 1666 1667 rComplete(currRing); 1668 for(j=0; j<IDELEMS(G); j++) 1669 { 1670 poly p=G->m[j]; 1671 while(p!=NULL) 1672 { 1673 p_Setm(p,currRing); pIter(p); 1674 } 1675 } 1408 1676 return result; 1409 1677 } 1410 1678 1411 1679 /**************************************************************** 1412 * Multipli kation of two ideals by elementwise*1680 * Multiplication of two ideals element by element * 1413 1681 * i.e. Let be A := (a_i) and B := (b_i), return C := (a_i*b_i) * 1414 1682 * destroy A, keeps B * … … 1419 1687 1420 1688 if(A==NULL || B==NULL) 1689 { 1421 1690 return NULL; 1422 1691 } 1423 1692 if(mB < mA) 1693 { 1424 1694 mA = mB; 1425 1695 } 1426 1696 ideal result = idInit(mA, 1); 1427 1697 … … 1431 1701 result->m[k] = pMult(A->m[i], pCopy(B->m[i])); 1432 1702 A->m[i]=NULL; 1433 if (result->m[k]!=NULL) k++; 1703 if (result->m[k]!=NULL) 1704 { 1705 k++; 1706 } 1434 1707 } 1435 1708 … … 1451 1724 ideal Mtmp = idLift(Gw, M, NULL, FALSE, TRUE, TRUE, NULL); 1452 1725 1453 // 3.12.02 Note: if Gw is a GB, then isSB = TRUE, otherwise FALSE1454 // So, it is better, if one tests whether Gw is a GB1455 // in ideals.cc:1456 // idLift (ideal mod, ideal submod,ideal * rest, BOOLEAN goodShape,1726 // If Gw is a GB, then isSB = TRUE, otherwise FALSE 1727 // So, it is better, if one tests whether Gw is a GB 1728 // in ideals.cc: 1729 // idLift (ideal mod, ideal submod,ideal * rest, BOOLEAN goodShape, 1457 1730 // BOOLEAN isSB,BOOLEAN divide,matrix * unit) 1458 1731 1459 / *Let be Mtmp = {m1,...,ms}, where mi=sum hij.in_gj, for all i=1,...,s1460 We compute F = {f1,...,fs}, where fi=sum hij.gj */1732 // Let be Mtmp = {m1,...,ms}, where mi=sum hij.in_gj, for all i=1,...,s 1733 // We compute F = {f1,...,fs}, where fi=sum hij.gj 1461 1734 int i, j, nM = IDELEMS(Mtmp); 1462 1735 ideal idpol, idLG; … … 1480 1753 } 1481 1754 1482 1483 #if 0 /*unused*/1755 //unused 1756 #if 0 1484 1757 static void checkidealCC(ideal G, char* Ch) 1485 1758 { … … 1495 1768 1496 1769 if(i != n) 1770 { 1497 1771 Print("%d, ", ntmp); 1772 } 1498 1773 else 1774 { 1499 1775 Print(" bzw. %d ", ntmp); 1776 } 1500 1777 } 1501 1778 PrintS(" Monomen.\n"); … … 1505 1782 #endif 1506 1783 1507 #if 0 /*unused*/ 1784 //unused 1785 #if 0 1508 1786 static void HeadidString(ideal L, char* st) 1509 1787 { … … 1512 1790 Print("// The head terms of the ideal %s = ", st); 1513 1791 for(i=0; i<nL; i++) 1792 { 1514 1793 Print(" %s, ", pString(pHead(L->m[i]))); 1515 1794 } 1516 1795 Print(" %s;\n", pString(pHead(L->m[nL]))); 1517 1796 } … … 1522 1801 assume(iva->length() == ivb->length()); 1523 1802 int i; 1524 1525 1803 for(i=iva->length()-1; i>=0; i--) 1804 { 1526 1805 if((*iva)[i] - (*ivb)[i] != 0) 1806 { 1527 1807 return 0; 1528 1808 } 1809 } 1529 1810 return 1; 1530 1811 } 1531 1812 1532 1533 /* 1534 compute a next weight vector between curr_weight and target_weight 1535 with respect to an ideal G. 1536 */ 1813 /********************************************** 1814 * Look for the smallest absolut value in vec * 1815 **********************************************/ 1816 static int MivAbsMax(intvec* vec) 1817 { 1818 int i,k; 1819 if((*vec)[0] < 0) 1820 { 1821 k = -(*vec)[0]; 1822 } 1823 else 1824 { 1825 k = (*vec)[0]; 1826 } 1827 for(i=1; i < (vec->length()); i++) 1828 { 1829 if((*vec)[i] < 0) 1830 { 1831 if(-(*vec)[i] > k) 1832 { 1833 k = -(*vec)[i]; 1834 } 1835 } 1836 else 1837 { 1838 if((*vec)[i] > k) 1839 { 1840 k = (*vec)[i]; 1841 } 1842 } 1843 } 1844 return k; 1845 } 1846 1847 /********************************************************************** 1848 * Compute a next weight vector between curr_weight and target_weight * 1849 * with respect to an ideal <G>. * 1850 **********************************************************************/ 1537 1851 static intvec* MwalkNextWeightCC(intvec* curr_weight, intvec* target_weight, 1538 1852 ideal G) … … 1545 1859 1546 1860 int nRing = currRing->N; 1547 int j, nG = IDELEMS(G);1861 int checkRed, j, kkk, nG = IDELEMS(G); 1548 1862 intvec* ivtemp; 1549 1863 … … 1558 1872 mpz_init(MwWd); 1559 1873 1560 1874 mpz_t sing_int; 1875 mpz_init(sing_int); 1876 mpz_set_si(sing_int, 2147483647); 1877 1878 mpz_t sing_int_half; 1879 mpz_init(sing_int_half); 1880 mpz_set_si(sing_int_half, 3*(1073741824/2)); 1881 1561 1882 mpz_t deg_w0_p1, deg_d0_p1; 1562 1883 mpz_init(deg_w0_p1); … … 1566 1887 mpz_init(sztn); 1567 1888 mpz_init(sntz); 1889 1568 1890 mpz_t t_null; 1569 1891 mpz_init(t_null); … … 1572 1894 mpz_init(ggt); 1573 1895 1574 int tn0, tn1, tz1, ncmp, gcd_tmp, ntmp; 1896 mpz_t dcw; 1897 mpz_init(dcw); 1898 1899 //int tn0, tn1, tz1, ncmp, gcd_tmp, ntmp; 1900 int gcd_tmp; 1575 1901 intvec* diff_weight = MivSub(target_weight, curr_weight); 1576 1902 1577 poly g, gw; 1903 intvec* diff_weight1 = MivSub(target_weight, curr_weight); 1904 poly g; 1905 //poly g, gw; 1578 1906 for (j=0; j<nG; j++) 1579 1907 { … … 1611 1939 } 1612 1940 1613 //compute a simpl yfraction of s1941 //compute a simple fraction of s 1614 1942 cancel(s_zaehler, s_nenner); 1615 1943 … … 1637 1965 } 1638 1966 } 1639 1967 //Print("\n// Alloc Size = %d \n", nRing*sizeof(mpz_t)); 1640 1968 mpz_t *vec=(mpz_t*)omAlloc(nRing*sizeof(mpz_t)); 1641 1642 /* there is no 0<t<1 and define the next weight vector that is equal to 1643 the current weight vector */1969 1970 1971 // there is no 0<t<1 and define the next weight vector that is equal to the current weight vector 1644 1972 if(mpz_cmp(t_nenner, t_null) == 0) 1645 1973 { 1974 Print("\n//MwalkNextWeightCC: t_nenner ist Null!"); 1646 1975 delete diff_weight; 1647 1976 diff_weight = ivCopy(curr_weight);//take memory … … 1649 1978 } 1650 1979 1651 / * define the target vector as the next weight vector, if t = 1 */1980 // define the target vector as the next weight vector, if t = 1 1652 1981 if(mpz_cmp_si(t_nenner, 1)==0 && mpz_cmp_si(t_zaehler,1)==0) 1653 1982 { … … 1657 1986 } 1658 1987 1659 1660 //14.08.03 simplify the both vectors curr_weight and diff_weight (C-int) 1988 checkRed = 0; 1989 1990 SIMPLIFY_GCD: 1991 1992 // simplify the vectors curr_weight and diff_weight (C-int) 1661 1993 gcd_tmp = (*curr_weight)[0]; 1662 1994 … … 1665 1997 gcd_tmp = gcd(gcd_tmp, (*curr_weight)[j]); 1666 1998 if(gcd_tmp == 1) 1999 { 1667 2000 break; 1668 }1669 2001 } 2002 } 1670 2003 if(gcd_tmp != 1) 2004 { 1671 2005 for (j=0; j<nRing; j++) 1672 2006 { 1673 2007 gcd_tmp = gcd(gcd_tmp, (*diff_weight)[j]); 1674 2008 if(gcd_tmp == 1) 2009 { 1675 2010 break; 1676 } 1677 2011 } 2012 } 2013 } 1678 2014 if(gcd_tmp != 1) 2015 { 1679 2016 for (j=0; j<nRing; j++) 1680 2017 { … … 1682 2019 (*diff_weight)[j] = (*diff_weight)[j]/gcd_tmp; 1683 2020 } 2021 } 2022 if(checkRed > 0) 2023 { 2024 for (j=0; j<nRing; j++) 2025 { 2026 mpz_set_si(vec[j], (*diff_weight)[j]); 2027 } 2028 goto TEST_OVERFLOW; 2029 } 2030 1684 2031 #ifdef NEXT_VECTORS_CC 1685 2032 Print("\n// gcd of the weight vectors (current and target) = %d", gcd_tmp); … … 1691 2038 #endif 1692 2039 1693 mpz_t ddf; mpz_init(ddf); 1694 mpz_t dcw; mpz_init(dcw); 1695 BOOLEAN isdwpos; 2040 // BOOLEAN isdwpos; 1696 2041 1697 2042 // construct a new weight vector … … 1702 2047 1703 2048 if( (*diff_weight)[j]>0) 2049 { 1704 2050 mpz_mul_ui(s_zaehler, t_zaehler, (*diff_weight)[j]); 2051 } 1705 2052 else 1706 2053 { … … 1708 2055 mpz_neg(s_zaehler, s_zaehler); 1709 2056 } 1710 1711 2057 mpz_add(sntz, s_nenner, s_zaehler); 1712 1713 2058 mpz_init_set(vec[j], sntz); 1714 2059 1715 2060 #ifdef NEXT_VECTORS_CC 1716 2061 Print("\n// j = %d ==> ", j); … … 1728 2073 1729 2074 if(j==0) 2075 { 1730 2076 mpz_set(ggt, sntz); 2077 } 1731 2078 else 2079 { 1732 2080 if(mpz_cmp_si(ggt,1) != 0) 2081 { 1733 2082 mpz_gcd(ggt, ggt, sntz); 1734 2083 } 2084 } 1735 2085 } 1736 2086 … … 1740 2090 #endif 1741 2091 1742 mpz_t omega; 1743 mpz_t sing_int; 1744 mpz_init_set_ui(sing_int, 2147483647); 1745 1746 /* construct a new weight vector and check whether vec[j] is overflow!! 1747 i.e. vec[j] > 2^31. 1748 If vec[j] doesn't overflow, define a weight vector 1749 otherwise, report that overflow appears. 1750 In the second case test whether the defined new vector correct is 1751 plays an important rolle */ 1752 2092 /********************************************************************** 2093 * construct a new weight vector and check whether vec[j] is overflow, * 2094 * i.e. vec[j] > 2^31. * 2095 * If vec[j] doesn't overflow, define a weight vector. Otherwise, * 2096 * report that overflow appears. In the second case, test whether the * 2097 * the correctness of the new vector plays an important role * 2098 **********************************************************************/ 2099 kkk=0; 2100 for(j=0; j<nRing; j++) 2101 { 2102 if(mpz_cmp(vec[j], sing_int_half) >= 0) 2103 { 2104 goto REDUCTION; 2105 } 2106 } 2107 checkRed = 1; 1753 2108 for (j=0; j<nRing; j++) 1754 { 1755 if(mpz_cmp_si(ggt,1)==0) 2109 { 1756 2110 (*diff_weight)[j] = mpz_get_si(vec[j]); 2111 } 2112 goto SIMPLIFY_GCD; 2113 2114 REDUCTION: 2115 for (j=0; j<nRing; j++) 2116 { 2117 (*diff_weight)[j] = mpz_get_si(vec[j]); 2118 } 2119 while(MivAbsMax(diff_weight) >= 5) 2120 { 2121 for (j=0; j<nRing; j++) 2122 { 2123 if(mpz_cmp_si(ggt,1)==0) 2124 { 2125 (*diff_weight1)[j] = floor(0.1*(*diff_weight)[j] + 0.5); 2126 // Print("\n// vector[%d] = %d \n",j+1, (*diff_weight1)[j]); 2127 } 2128 else 2129 { 2130 mpz_divexact(vec[j], vec[j], ggt); 2131 (*diff_weight1)[j] = floor(0.1*(*diff_weight)[j] + 0.5); 2132 // Print("\n// vector[%d] = %d \n",j+1, (*diff_weight1)[j]); 2133 } 2134 /* 2135 if((*diff_weight1)[j] == 0) 2136 { 2137 kkk = kkk + 1; 2138 } 2139 */ 2140 } 2141 2142 2143 /* 2144 if(kkk > nRing - 1) 2145 { 2146 // diff_weight was reduced to zero 2147 // Print("\n // MwalkNextWeightCC: geaenderter Vector gleich Null! \n"); 2148 goto TEST_OVERFLOW; 2149 } 2150 */ 2151 2152 if(test_w_in_ConeCC(G,diff_weight1) != 0) 2153 { 2154 Print("\n// MwalkNextWeightCC: geaenderter vector liegt in Groebnerkegel! \n"); 2155 for (j=0; j<nRing; j++) 2156 { 2157 (*diff_weight)[j] = (*diff_weight1)[j]; 2158 } 2159 if(MivAbsMax(diff_weight) < 5) 2160 { 2161 checkRed = 1; 2162 goto SIMPLIFY_GCD; 2163 } 2164 } 1757 2165 else 1758 2166 { 1759 mpz_divexact(vec[j], vec[j], ggt); 1760 (*diff_weight)[j] = mpz_get_si(vec[j]); 1761 } 1762 2167 // Print("\n// MwalkNextWeightCC: geaenderter vector liegt nicht in Groebnerkegel! \n"); 2168 break; 2169 } 2170 } 2171 2172 TEST_OVERFLOW: 2173 2174 for (j=0; j<nRing; j++) 2175 { 1763 2176 if(mpz_cmp(vec[j], sing_int)>=0) 2177 { 1764 2178 if(Overflow_Error == FALSE) 1765 2179 { 1766 2180 Overflow_Error = TRUE; 1767 1768 PrintS("\n// ** OVERFLOW in \"NextVector\": "); 2181 PrintS("\n// ** OVERFLOW in \"MwalkNextWeightCC\": "); 1769 2182 mpz_out_str( stdout, 10, vec[j]); 1770 2183 PrintS(" is greater than 2147483647 (max. integer representation)"); 1771 Print("\n// So vector[%d] := %d is wrong!!\n",j+1, (*diff_weight)[j]); 1772 } 2184 Print("\n// So vector[%d] := %d is wrong!!\n",j+1, vec[j]); 2185 } 2186 } 1773 2187 } 1774 2188 1775 2189 FINISH: 1776 1777 mpz_clear(ggt); 2190 delete diff_weight1; 1778 2191 mpz_clear(t_zaehler); 1779 2192 mpz_clear(t_nenner); 2193 mpz_clear(s_zaehler); 2194 mpz_clear(s_nenner); 1780 2195 mpz_clear(sntz); 1781 2196 mpz_clear(sztn); … … 1784 2199 mpz_clear(deg_w0_p1); 1785 2200 mpz_clear(deg_d0_p1); 2201 mpz_clear(ggt); 1786 2202 omFree(vec); 1787 2203 mpz_clear(sing_int_half); 2204 mpz_clear(sing_int); 2205 mpz_clear(dcw); 2206 mpz_clear(t_null); 2207 2208 2209 1788 2210 if(Overflow_Error == FALSE) 2211 { 1789 2212 Overflow_Error = nError; 1790 1791 return diff_weight; 1792 } 1793 1794 /* 1795 compute an intermediate weight vector from iva to ivb w.r.t. 1796 the reduced Groebner basis G. 1797 Return NULL, if it is equal to iva or iva = avb. 1798 */ 2213 } 2214 rComplete(currRing); 2215 for(kkk=0; kkk<IDELEMS(G);kkk++) 2216 { 2217 poly p=G->m[kkk]; 2218 while(p!=NULL) 2219 { 2220 p_Setm(p,currRing); 2221 pIter(p); 2222 } 2223 } 2224 return diff_weight; 2225 } 2226 2227 /********************************************************************** 2228 * Compute an intermediate weight vector from iva to ivb w.r.t. * 2229 * the reduced Groebner basis G. * 2230 * Return NULL, if it is equal to iva or iva = avb. * 2231 **********************************************************************/ 1799 2232 intvec* MkInterRedNextWeight(intvec* iva, intvec* ivb, ideal G) 1800 2233 { … … 1803 2236 1804 2237 if(G == NULL) 2238 { 1805 2239 return tmp; 1806 2240 } 1807 2241 if(MivComp(iva, ivb) == 1) 2242 { 1808 2243 return tmp; 1809 2244 } 1810 2245 result = MwalkNextWeightCC(iva, ivb, G); 1811 2246 … … 1820 2255 } 1821 2256 1822 /* 01.11.01 */ 1823 /* define and execute a new ring which order is (a(va),lp,C) */ 1824 static void VMrDefault(intvec* va) 2257 /************************************************************** 2258 * define and execute a new ring which order is (a(vb),a(va),lp,C) * 2259 * ************************************************************/ 2260 static void VMrHomogeneous(intvec* va, intvec* vb) 1825 2261 { 1826 2262 1827 2263 if ((currRing->ppNoether)!=NULL) 2264 { 1828 2265 pDelete(&(currRing->ppNoether)); 1829 2266 } 1830 2267 if (((sLastPrinted.rtyp>BEGIN_RING) && (sLastPrinted.rtyp<END_RING)) || 1831 2268 ((sLastPrinted.rtyp==LIST_CMD)&&(lRingDependend((lists)sLastPrinted.data)))) 1832 1833 2269 { 1834 2270 sLastPrinted.CleanUp(); … … 1842 2278 int nb = 4; 1843 2279 1844 /*names*/ 1845 char* Q; //30.10.01 to avoid the corrupted memory, NOT change!! 2280 2281 //names 2282 char* Q; // In order to avoid the corrupted memory, do not change. 2283 r->names = (char **) omAlloc0(nv * sizeof(char_ptr)); 2284 for(i=0; i<nv; i++) 2285 { 2286 Q = currRing->names[i]; 2287 r->names[i] = omStrDup(Q); 2288 } 2289 2290 //weights: entries for 3 blocks: NULL Made:??? 2291 r->wvhdl = (int **)omAlloc0(nb * sizeof(int_ptr)); 2292 r->wvhdl[0] = (int*) omAlloc(nv*sizeof(int)); 2293 r->wvhdl[1] = (int*) omAlloc((nv-1)*sizeof(int)); 2294 2295 for(i=0; i<nv-1; i++) 2296 { 2297 r->wvhdl[1][i] = (*vb)[i]; 2298 r->wvhdl[0][i] = (*va)[i]; 2299 } 2300 r->wvhdl[0][nv] = (*va)[nv]; 2301 2302 // order: (1..1),a,lp,C 2303 r->order = (int *) omAlloc(nb * sizeof(int *)); 2304 r->block0 = (int *)omAlloc0(nb * sizeof(int *)); 2305 r->block1 = (int *)omAlloc0(nb * sizeof(int *)); 2306 2307 // ringorder a for the first block: var 1..nv 2308 r->order[0] = ringorder_a; 2309 r->block0[0] = 1; 2310 r->block1[0] = nv; 2311 2312 // ringorder a for the second block: var 2..nv 2313 r->order[1] = ringorder_a; 2314 r->block0[1] = 2; 2315 r->block1[1] = nv; 2316 2317 // ringorder lp for the third block: var 2..nv 2318 r->order[2] = ringorder_lp; 2319 r->block0[2] = 2; 2320 r->block1[2] = nv; 2321 2322 // ringorder C for the 4th block 2323 // it is very important within "idLift", 2324 // especially, by ring syz_ring=rCurrRingAssure_SyzComp(); 2325 // therefore, nb must be (nBlocks(currRing) + 1) 2326 r->order[3] = ringorder_C; 2327 2328 // polynomial ring 2329 r->OrdSgn = 1; 2330 2331 // complete ring intializations 2332 2333 rComplete(r); 2334 2335 rChangeCurrRing(r); 2336 } 2337 2338 2339 /************************************************************** 2340 * define and execute a new ring which order is (a(va),lp,C) * 2341 * ************************************************************/ 2342 static void VMrDefault(intvec* va) 2343 { 2344 2345 if ((currRing->ppNoether)!=NULL) 2346 { 2347 pDelete(&(currRing->ppNoether)); 2348 } 2349 if (((sLastPrinted.rtyp>BEGIN_RING) && (sLastPrinted.rtyp<END_RING)) || 2350 ((sLastPrinted.rtyp==LIST_CMD)&&(lRingDependend((lists)sLastPrinted.data)))) 2351 { 2352 sLastPrinted.CleanUp(); 2353 } 2354 2355 ring r = (ring) omAlloc0Bin(sip_sring_bin); 2356 int i, nv = currRing->N; 2357 2358 r->cf = currRing->cf; 2359 r->N = currRing->N; 2360 2361 int nb = 4; 2362 2363 //names 2364 char* Q; // In order to avoid the corrupted memory, do not change. 1846 2365 r->names = (char **) omAlloc0(nv * sizeof(char_ptr)); 1847 2366 for(i=0; i<nv; i++) … … 1862 2381 r->block1 = (int *)omAlloc0(nb * sizeof(int *)); 1863 2382 1864 / * ringorder a for the first block: var 1..nv */2383 // ringorder a for the first block: var 1..nv 1865 2384 r->order[0] = ringorder_a; 1866 2385 r->block0[0] = 1; 1867 2386 r->block1[0] = nv; 1868 2387 1869 / * ringorder lp for the second block: var 1..nv */2388 // ringorder lp for the second block: var 1..nv 1870 2389 r->order[1] = ringorder_lp; 1871 2390 r->block0[1] = 1; 1872 2391 r->block1[1] = nv; 1873 2392 1874 / * ringorder C for the third block */2393 // ringorder C for the third block 1875 2394 // it is very important within "idLift", 1876 2395 // especially, by ring syz_ring=rCurrRingAssure_SyzComp(); … … 1878 2397 r->order[2] = ringorder_C; 1879 2398 1880 / * the last block: everything is 0 */2399 // the last block: everything is 0 1881 2400 r->order[3] = 0; 1882 2401 1883 / *polynomial ring*/2402 // polynomial ring 1884 2403 r->OrdSgn = 1; 1885 2404 1886 /* complete ring intializations */ 2405 // complete ring intializations 2406 1887 2407 rComplete(r); 1888 2408 … … 1890 2410 } 1891 2411 1892 /* 03.11.01 */ 1893 /* define and execute a new ring which order is a lexicographic order */ 2412 /********************************************************************** 2413 * define and execute a new ring which order is a lexicographic order * 2414 ***********************************************************************/ 1894 2415 static void VMrDefaultlp(void) 1895 2416 { 1896 2417 1897 2418 if ((currRing->ppNoether)!=NULL) 2419 { 1898 2420 pDelete(&(currRing->ppNoether)); 1899 2421 } 1900 2422 if (((sLastPrinted.rtyp>BEGIN_RING) && (sLastPrinted.rtyp<END_RING)) || 1901 2423 ((sLastPrinted.rtyp==LIST_CMD)&&(lRingDependend((lists)sLastPrinted.data)))) … … 1910 2432 r->cf = currRing->cf; 1911 2433 r->N = currRing->N; 1912 int nb = rBlocks(currRing) + 1; //31.10.01 (+1)1913 1914 / *names*/1915 char* Q; // 30.10.01 to avoid the corrupted memory, NOTchange!!2434 int nb = rBlocks(currRing) + 1; 2435 2436 // names 2437 char* Q; // to avoid the corrupted memory, do not change!! 1916 2438 r->names = (char **) omAlloc0(nv * sizeof(char_ptr)); 1917 2439 for(i=0; i<nv; i++) … … 1945 2467 1946 2468 /* complete ring intializations */ 2469 1947 2470 rComplete(r); 1948 2471 … … 1967 2490 res->cf = currRing->cf; currRing->cf->ref++; 1968 2491 1969 // if (currRing->cf->extRing!=NULL)1970 // currRing->cf->extRing->ref++;1971 //1972 // if (rParameter (currRing)!=NULL)1973 // {1974 // res->cf->extRing->qideal->m[0]=p_Copy(currRing->cf->extRing->qideal->m[0],currRing->cf->extRing);1975 // int l=rPar(currRing);1976 //1977 // res->cf->extRing->names=(char **)omAlloc(l*sizeof(char_ptr));1978 //1979 // for(i=l-1;i>=0;i--)1980 // rParameter (res)[i]=omStrDup(rParameter (currRing)[i]);1981 // }1982 2492 1983 2493 /*weights: entries for 3 blocks: NULL Made:???*/ … … 1988 2498 1989 2499 /* order: a,lp,C,0 */ 2500 1990 2501 res->order = (int *) omAlloc(nb * sizeof(int *)); 1991 2502 res->block0 = (int *)omAlloc0(nb * sizeof(int *)); 1992 2503 res->block1 = (int *)omAlloc0(nb * sizeof(int *)); 1993 2504 1994 / * ringorder a for the first block: var 1..nv */2505 // ringorder a for the first block: var 1..nv 1995 2506 res->order[0] = ringorder_a; 1996 2507 res->block0[0] = 1; 1997 2508 res->block1[0] = nv; 1998 2509 1999 / * ringorder lp for the second block: var 1..nv */2510 // ringorder lp for the second block: var 1..nv 2000 2511 res->order[1] = ringorder_lp; 2001 2512 res->block0[1] = 1; 2002 2513 res->block1[1] = nv; 2003 2514 2004 / * ringorder C for the third block */2515 // ringorder C for the third block 2005 2516 // it is very important within "idLift", 2006 2517 // especially, by ring syz_ring=rCurrRingAssure_SyzComp(); … … 2008 2519 res->order[2] = ringorder_C; 2009 2520 2010 / * the last block: everything is 0 */2521 // the last block: everything is 0 2011 2522 res->order[3] = 0; 2012 2523 2013 / *polynomial ring*/2524 // polynomial ring 2014 2525 res->OrdSgn = 1; 2015 2526 … … 2017 2528 res->names = (char **)omAlloc0(nv * sizeof(char_ptr)); 2018 2529 for (i=nv-1; i>=0; i--) 2530 { 2019 2531 res->names[i] = omStrDup(currRing->names[i]); 2020 2021 /* complete ring intializations */ 2022 rComplete(res); 2023 2532 } 2533 // complete ring intializations 2534 rComplete(res); 2024 2535 2025 2536 // clean up history … … 2030 2541 2031 2542 2032 / * execute the created ring */2543 // execute the created ring 2033 2544 rChangeCurrRing(res); 2034 2545 } … … 2048 2559 r->cf = currRing->cf; currRing->cf->ref++; 2049 2560 2050 // if (currRing->cf->extRing!=NULL)2051 // currRing->cf->extRing->ref++;2052 //2053 // if (rParameter (currRing)!=NULL)2054 // {2055 // r->cf->extRing->qideal->m[0]=p_Copy(currRing->cf->extRing->qideal->m[0], currRing->cf->extRing);2056 // int l=rPar(currRing);2057 // r->cf->extRing->names=(char **)omAlloc(l*sizeof(char_ptr));2058 //2059 // for(i=l-1;i>=0;i--)2060 // rParameter(r)[i]=omStrDup(rParameter (currRing)[i]);2061 // }2062 2063 2561 2064 2562 r->cf = currRing->cf; 2065 2563 r->N = currRing->N; 2066 int nb = rBlocks(currRing) + 1; //31.10.01 (+1)2067 2068 / *names*/2564 int nb = rBlocks(currRing) + 1; 2565 2566 // names 2069 2567 char* Q; 2070 2568 r->names = (char **) omAlloc0(nv * sizeof(char_ptr)); … … 2106 2604 // 2107 2605 // for(i=l-1;i>=0;i--) 2606 // { 2108 2607 // rParameter(r)[i]=omStrDup(rParameter(currRing)[i]); 2608 // } 2109 2609 // } 2110 2610 2111 /* complete ring intializations */ 2611 // complete ring intializations 2612 2112 2613 rComplete(r); 2113 2614 … … 2118 2619 } 2119 2620 2120 / * execute the created ring */2621 // execute the created ring 2121 2622 rChangeCurrRing(r); 2122 2623 } 2123 2624 2124 #if 0 /*unused*/ 2125 /* check wheather one or more components of a vector are zero */ 2625 //unused 2626 /************************************************************** 2627 * check wheather one or more components of a vector are zero * 2628 **************************************************************/ 2629 //#if 0 2126 2630 static int isNolVector(intvec* hilb) 2127 2631 { 2128 2632 int i; 2129 2633 for(i=hilb->length()-1; i>=0; i--) 2634 { 2130 2635 if((* hilb)[i]==0) 2636 { 2131 2637 return 1; 2132 2638 } 2639 } 2133 2640 return 0; 2134 2641 } 2135 #endif 2136 2642 //#endif 2137 2643 2138 2644 /****************************** Februar 2002 **************************** … … 2141 2647 * its perturbation degree is tp_deg * 2142 2648 * We call the following subfunction LastGB, if * 2143 * the computed intermediate weight vector or * 2144 * the perturbed target weight vector * 2145 * does NOT in the correct cone. * 2649 * the computed intermediate weight vector or * 2650 * if the perturbed target weight vector does NOT lie n the correct cone * 2146 2651 **************************************************************************/ 2147 2652 … … 2162 2667 intvec* ivNull = new intvec(nV); 2163 2668 intvec* extra_curr_weight = new intvec(nV); 2669 intvec* next_weight; 2670 2671 #ifndef BUCHBERGER_ALG 2164 2672 intvec* hilb_func; 2165 intvec* next_weight; 2166 2167 / * to avoid (1,0,...,0) as the target vector */2673 #endif 2674 2675 // to avoid (1,0,...,0) as the target vector 2168 2676 intvec* last_omega = new intvec(nV); 2169 2677 for(i=nV-1; i>0; i--) 2678 { 2170 2679 (*last_omega)[i] = 1; 2680 } 2171 2681 (*last_omega)[0] = 10000; 2172 2682 2173 2683 ring EXXRing = currRing; 2174 2684 2175 / * compute a pertubed weight vector of the target weight vector */2685 // compute a pertubed weight vector of the target weight vector 2176 2686 if(tp_deg > 1 && tp_deg <= nV) 2177 2687 { 2178 2688 //..25.03.03 VMrDefaultlp();// VMrDefault(target_weight); 2179 2689 if (rParameter (currRing) != NULL) 2690 { 2180 2691 DefRingParlp(); 2692 } 2181 2693 else 2694 { 2182 2695 VMrDefaultlp(); 2183 2696 } 2184 2697 TargetRing = currRing; 2185 2698 ssG = idrMoveR(G,EXXRing,currRing); … … 2194 2707 } 2195 2708 else 2709 { 2196 2710 target_weight = Mivlp(nV); 2197 2711 } 2198 2712 //Print("\n// ring r%d_%d = %s;\n", tp_deg, nwalk, rString(currRing)); 2199 2713 … … 2203 2717 nstep++; 2204 2718 to=clock(); 2205 /* compute a next weight vector */2719 // compute a next weight vector 2206 2720 next_weight = MkInterRedNextWeight(curr_weight,target_weight, G); 2207 2721 xtnw=xtnw+clock()-to; 2722 2208 2723 #ifdef PRINT_VECTORS 2209 2724 MivString(curr_weight, target_weight, next_weight); 2210 2725 #endif 2211 2726 2212 if(Overflow_Error == TRUE){ 2727 if(Overflow_Error == TRUE) 2728 { 2213 2729 newRing = currRing; 2214 2730 nnwinC = 0; 2215 2731 if(tp_deg == 1) 2732 { 2216 2733 nlast = 1; 2734 } 2217 2735 delete next_weight; 2218 2736 … … 2223 2741 } 2224 2742 2225 if(MivComp(next_weight, ivNull) == 1){ 2743 if(MivComp(next_weight, ivNull) == 1) 2744 { 2226 2745 //Print("\n// ring r%d_%d = %s;\n", tp_deg, nwalk, rString(currRing)); 2227 2746 newRing = currRing; … … 2234 2753 2235 2754 for(i=nV-1; i>=0; i--) 2236 (*extra_curr_weight)[i] = (*curr_weight)[i]; 2237 2755 { 2756 (*extra_curr_weight)[i] = (*curr_weight)[i]; 2757 } 2238 2758 /* 06.11.01 NOT Changed */ 2239 2759 for(i=nV-1; i>=0; i--) 2760 { 2240 2761 (*curr_weight)[i] = (*next_weight)[i]; 2241 2762 } 2242 2763 oldRing = currRing; 2243 2764 to=clock(); 2244 / * compute an initial form ideal of <G> w.r.t. "curr_vector" */2765 // compute an initial form ideal of <G> w.r.t. "curr_vector" 2245 2766 Gomega = MwalkInitialForm(G, curr_weight); 2246 2767 xtif=xtif+clock()-to; 2247 2768 2248 2769 #ifdef ENDWALKS 2249 if(endwalks == 1){ 2770 if(endwalks == 1) 2771 { 2250 2772 Print("\n// ring r%d_%d = %s;\n", tp_deg, nwalk, rString(currRing)); 2251 2773 idElements(Gomega, "Gw"); … … 2256 2778 #ifndef BUCHBERGER_ALG 2257 2779 if(isNolVector(curr_weight) == 0) 2780 { 2258 2781 hilb_func = hFirstSeries(Gomega,NULL,NULL,curr_weight,currRing); 2782 } 2259 2783 else 2784 { 2260 2785 hilb_func = hFirstSeries(Gomega,NULL,NULL,last_omega,currRing); 2786 } 2261 2787 #endif // BUCHBERGER_ALG 2262 2788 … … 2264 2790 //..25.03.03 VMrDefault(curr_weight); 2265 2791 if (rParameter (currRing) != NULL) 2792 { 2266 2793 DefRingPar(curr_weight); 2794 } 2267 2795 else 2796 { 2268 2797 VMrDefault(curr_weight); 2269 2798 } 2270 2799 newRing = currRing; 2271 2800 Gomega1 = idrMoveR(Gomega, oldRing,currRing); … … 2303 2832 idDelete(&F1); 2304 2833 2305 if(endwalks == 1){ 2834 if(endwalks == 1) 2835 { 2306 2836 //Print("\n// ring r%d_%d = %s;\n", tp_deg, nwalk, rString(currRing)); 2307 2837 break; … … 2317 2847 //..25.03.03 VMrDefaultlp();//define and execute the ring "lp" 2318 2848 if (rParameter (currRing) != NULL) 2849 { 2319 2850 DefRingParlp(); 2851 } 2320 2852 else 2853 { 2321 2854 VMrDefaultlp(); 2322 2855 } 2323 2856 F1 = idrMoveR(G, newRing,currRing); 2324 2857 … … 2358 2891 //..25.03.03 VMrDefaultlp();//define and execute the ring "lp" 2359 2892 if (rParameter (currRing) != NULL) 2893 { 2360 2894 DefRingParlp(); 2895 } 2361 2896 else 2897 { 2362 2898 VMrDefaultlp(); 2363 2899 } 2364 2900 2365 2901 F1 = idrMoveR(G, newRing,currRing); … … 2379 2915 2380 2916 if(Overflow_Error == FALSE) 2917 { 2381 2918 Overflow_Error = nError; 2382 2919 } 2383 2920 return(result); 2384 2921 } 2385 2922 2386 2387 /* check whether a polynomial of G has least 3 monomials */ 2923 /********************************************************** 2924 * check whether a polynomial of G has least 3 monomials * 2925 **********************************************************/ 2388 2926 static int lengthpoly(ideal G) 2389 2927 { 2390 2928 int i; 2391 2929 for(i=IDELEMS(G)-1; i>=0; i--) 2930 { 2392 2931 #if 0 2393 2932 if(pLength(G->m[i])>2) 2933 { 2394 2934 return 1; 2935 } 2395 2936 #else 2396 2937 if((G->m[i]!=NULL) /* len >=0 */ … … 2399 2940 && (G->m[i]->next->next->next!=NULL) /* len >=3 */ 2400 2941 //&& (G->m[i]->next->next->next->next!=NULL) /* len >=4 */ 2401 ) return 1; 2402 #endif 2942 ) 2943 { 2944 return 1; 2945 } 2946 #endif 2947 } 2403 2948 return 0; 2404 2949 } 2405 2950 2406 /* check whether a polynomial of G has least 2 monomials */ 2951 /********************************************************* 2952 * check whether a polynomial of G has least 2 monomials * 2953 **********************************************************/ 2407 2954 static int islengthpoly2(ideal G) 2408 2955 { 2409 2956 int i; 2410 2957 for(i=IDELEMS(G)-1; i>=0; i--) 2958 { 2411 2959 if((G->m[i]!=NULL) /* len >=0 */ 2412 2960 && (G->m[i]->next!=NULL) /* len >=1 */ 2413 2961 && (G->m[i]->next->next!=NULL)) /* len >=2 */ 2962 { 2414 2963 return 1; 2415 2964 } 2965 } 2416 2966 return 0; 2417 2967 } … … 2435 2985 */ 2436 2986 2437 /* 22438 * return the initial term of an ideal 2439 */2987 /*************************************** 2988 * return the initial term of an ideal * 2989 ***************************************/ 2440 2990 static ideal idHeadCC(ideal h) 2441 2991 { … … 2447 2997 { 2448 2998 if (h->m[i]!=NULL) 2999 { 2449 3000 m->m[i]=pHead(h->m[i]); 3001 } 2450 3002 } 2451 3003 return m; 2452 3004 } 2453 3005 2454 /* check whether two head-ideals are the same */ 3006 /********************************************** 3007 * check whether two head-ideals are the same * 3008 **********************************************/ 2455 3009 static inline int test_G_GB_walk(ideal H0, ideal H1) 2456 3010 { … … 2458 3012 2459 3013 if(nG != IDELEMS(H1)) 3014 { 2460 3015 return 0; 2461 3016 } 2462 3017 for(i=nG-1; i>=0; i--) 2463 3018 { … … 2472 3027 #else 2473 3028 if(!pEqualPolys(H0->m[i],H1->m[i])) 3029 { 2474 3030 return 0; 2475 #endif 2476 } 2477 3031 } 3032 #endif 3033 } 2478 3034 return 1; 2479 3035 } 2480 3036 2481 #if 0 /*unused*/ 2482 /* 19.11.01 */ 2483 /* find the maximal total degree of polynomials in G */ 3037 //unused 3038 /***************************************************** 3039 * find the maximal total degree of polynomials in G * 3040 *****************************************************/ 3041 #if 0 2484 3042 static int Trandegreebound(ideal G) 2485 3043 { 2486 3044 int i, nG = IDELEMS(G); 2487 int np=1, nV = currRing->N; 3045 // int np=1; 3046 int nV = currRing->N; 2488 3047 int degtmp, result = 0; 2489 3048 intvec* ivUnit = Mivdp(nV); … … 2491 3050 for(i=nG-1; i>=0; i--) 2492 3051 { 2493 / * find the maximal total degree of the polynomial G[i] */3052 // find the maximal total degree of the polynomial G[i] 2494 3053 degtmp = MwalkWeightDegree(G->m[i], ivUnit); 2495 3054 if(degtmp > result) 3055 { 2496 3056 result = degtmp; 3057 } 2497 3058 } 2498 3059 delete ivUnit; … … 2501 3062 #endif 2502 3063 2503 / * perturb the weight vector iva w.r.t. the ideal G.2504 the monomial order of the current ring is the w_1 weight lex. order. 2505 define w := d^(n-1)w_1+ d^(n-2)w_2, ...+ dw_(n-1)+ w_n2506 where d := 1 + max{totdeg(g):g in G}*m, or2507 d := (2*maxdeg*maxdeg + (nV+1)*maxdeg)*m;2508 */ 2509 2510 #if 0 /*unused*/2511 //GMP 3064 //unused 3065 /************************************************************************ 3066 * perturb the weight vector iva w.r.t. the ideal G. * 3067 * the monomial order of the current ring is the w_1 weight lex. order * 3068 * define w := d^(n-1)w_1+ d^(n-2)w_2, ...+ dw_(n-1)+ w_n * 3069 * where d := 1 + max{totdeg(g):g in G}*m, or * 3070 * d := (2*maxdeg*maxdeg + (nV+1)*maxdeg)*m; * 3071 ************************************************************************/ 3072 #if 0 2512 3073 static intvec* TranPertVector(ideal G, intvec* iva) 2513 3074 { … … 2515 3076 Overflow_Error = FALSE; 2516 3077 2517 int i, j, nG = IDELEMS(G); 3078 int i, j; 3079 // int nG = IDELEMS(G); 2518 3080 int nV = currRing->N; 2519 3081 … … 2527 3089 { 2528 3090 mtmp = (*ivMat)[i]; 2529 2530 if(mtmp <0) mtmp = -mtmp; 2531 3091 if(mtmp <0) 3092 { 3093 mtmp = -mtmp; 3094 } 2532 3095 if(mtmp > m) 3096 { 2533 3097 m = mtmp; 2534 } 2535 2536 /* define the maximal total degree of polynomials of G */ 3098 } 3099 } 3100 3101 // define the maximal total degree of polynomials of G 2537 3102 mpz_t ndeg; 2538 3103 mpz_init(ndeg); … … 2555 3120 mpz_mul_ui(ndeg, ndeg, m); 2556 3121 3122 mpz_clear(ztmp); 3123 2557 3124 //PrintS("\n// with the new upper degree bound (2d^2+(n+1)d)*m "); 2558 3125 //Print("\n// where d = %d, n = %d and bound = %d", maxdeg, nV, ndeg); 2559 3126 #endif //UPPER_BOUND 2560 3127 2561 /* 29.08.03*/2562 3128 #ifdef INVEPS_SMALL_IN_TRAN 2563 if(mpz_cmp_ui(ndeg, nV)>0 && nV > 3) 3129 if(mpz_cmp_ui(ndeg, nV)>0 && nV > 3) 3130 { 2564 3131 mpz_cdiv_q_ui(ndeg, ndeg, nV); 2565 3132 } 2566 3133 //PrintS("\n// choose the \"small\" inverse epsilon:"); 2567 3134 //mpz_out_str(stdout, 10, ndeg); … … 2581 3148 mpz_t *ivtmp=(mpz_t *)omAlloc(nV*sizeof(mpz_t)); 2582 3149 for(i=0; i<nV; i++) 3150 { 2583 3151 mpz_init(ivtmp[i]); 2584 3152 } 2585 3153 mpz_t sing_int; 2586 3154 mpz_init_set_ui(sing_int, 2147483647); … … 2588 3156 intvec* repr_vector = new intvec(nV); 2589 3157 2590 / * define ivtmp := ndeg^(n-1).w_1 + ndeg^(n-2).w_2 + ... + w_n */3158 // define ivtmp := ndeg^(n-1).w_1 + ndeg^(n-2).w_2 + ... + w_n 2591 3159 for(i=0; i<nV; i++) 3160 { 2592 3161 for(j=0; j<nV; j++) 2593 3162 { 2594 3163 if( (*ivMat)[i*nV+j] >= 0 ) 3164 { 2595 3165 mpz_mul_ui(ivres[i], ivres[i], (*ivMat)[i*nV+j]); 3166 } 2596 3167 else 2597 3168 { … … 2601 3172 mpz_add(ivtmp[j], ivtmp[j], ivres[i]); 2602 3173 } 2603 3174 } 2604 3175 delete ivMat; 2605 3176 … … 2633 3204 omFree(ivres); 2634 3205 omFree(ivtmp); 3206 3207 mpz_clear(sing_int); 3208 mpz_clear(deg_tmp); 3209 mpz_clear(ndeg); 3210 2635 3211 return repr_vector; 2636 3212 } 2637 3213 #endif 2638 3214 2639 2640 2641 #if 0 /*unused*/ 3215 //unused 3216 #if 0 2642 3217 static intvec* TranPertVector_lp(ideal G) 2643 3218 { 2644 3219 BOOLEAN nError = Overflow_Error; 2645 3220 Overflow_Error = FALSE; 2646 2647 int i , j, nG = IDELEMS(G);3221 // int j, nG = IDELEMS(G); 3222 int i; 2648 3223 int nV = currRing->N; 2649 3224 2650 / * define the maximal total degree of polynomials of G */3225 // define the maximal total degree of polynomials of G 2651 3226 mpz_t ndeg; 2652 3227 mpz_init(ndeg); … … 2667 3242 mpz_mul_ui(maxdeg, maxdeg, nV+1); 2668 3243 mpz_add(ndeg, ztmp, maxdeg); 2669 /* 2670 PrintS("\n// with the new upper degree bound (2d^2+(n+1)d)*m "); 2671 Print("\n// where d = %d, n = %d and bound = %d", 2672 mpz_get_si(maxdeg), nV, mpz_get_si(ndeg)); 2673 */ 2674 #endif //UPPER_BOUND 3244 // PrintS("\n// with the new upper degree bound (2d^2+(n+1)d)*m "); 3245 // Print("\n// where d = %d, n = %d and bound = %d", 3246 // mpz_get_si(maxdeg), nV, mpz_get_si(ndeg)); 3247 3248 mpz_clear(ztmp); 3249 3250 #endif 2675 3251 2676 3252 #ifdef INVEPS_SMALL_IN_TRAN … … 2725 3301 2726 3302 omFree(ivres); 3303 3304 mpz_clear(ndeg); 3305 mpz_clear(sing_int); 3306 2727 3307 return repr_vector; 2728 3308 } 2729 3309 #endif 2730 3310 2731 2732 #if 0 /*unused*/ 2733 //GMP 3311 //unused 3312 #if 0 2734 3313 static intvec* RepresentationMatrix_Dp(ideal G, intvec* M) 2735 3314 { … … 2745 3324 for(i=IDELEMS(G)-1; i>=0; i--) 2746 3325 { 2747 / * find the maximal total degree of the polynomial G[i] */3326 // find the maximal total degree of the polynomial G[i] 2748 3327 degtmp = MwalkWeightDegree(G->m[i], ivUnit); 2749 3328 if(degtmp > maxdeg) … … 2766 3345 mpz_init(ivtmp[i]); 2767 3346 2768 / * define ivtmp := ndeg^(n-1).w_1 + ndeg^(n-2).w_2 + ... + w_n */3347 // define ivtmp := ndeg^(n-1).w_1 + ndeg^(n-2).w_2 + ... + w_n 2769 3348 for(i=0; i<nV; i++) 2770 3349 for(j=0; j<nV; j++) … … 2780 3359 mpz_add(ivtmp[j], ivtmp[j], ztmp); 2781 3360 } 2782 3361 delete ivres; 2783 3362 mpz_t sing_int; 2784 3363 mpz_init_set_ui(sing_int, 2147483647); … … 2811 3390 Overflow_Error = nError; 2812 3391 3392 mpz_clear(sing_int); 2813 3393 mpz_clear(ztmp); 2814 3394 omFree(ivtmp); … … 2818 3398 #endif 2819 3399 2820 2821 2822 2823 2824 /* The following subroutine is the implementation of our first improved 2825 Groebner walk algorithm, i.e. the first altervative algorithm. 2826 First we use the Grobner walk algorithm and then we call the changed 2827 perturbation walk algorithm with decreased degree, if an intermediate 2828 weight vector is equal to the current target weight vector. 2829 This call will be only repeated until we get the wanted reduced Groebner 2830 basis or n times, where n is the numbers of variables. 2831 */ 2832 2833 // 19 Juni 2003 2834 #if 0 /* unused*/ 3400 /***************************************************************************** 3401 * The following subroutine is the implementation of our first improved * 3402 * Groebner walk algorithm, i.e. the first altervative algorithm. * 3403 * First we use the Grobner walk algorithm and then we call the changed * 3404 * perturbation walk algorithm with decreased degree, if an intermediate * 3405 * weight vector is equal to the current target weight vector. * 3406 * This call will be only repeated until we get the wanted reduced Groebner * 3407 * basis or n times, where n is the numbers of variables. * 3408 *****************************************************************************/ 3409 3410 //unused 3411 #if 0 2835 3412 static int testnegintvec(intvec* v) 2836 3413 { 2837 3414 int n = v->length(); 2838 3415 int i; 2839 for(i=0; i<n; i++){ 3416 for(i=0; i<n; i++) 3417 { 2840 3418 if((*v)[i]<0) 3419 { 2841 3420 return(1); 3421 } 2842 3422 } 2843 3423 return(0); … … 2845 3425 #endif 2846 3426 2847 2848 /* 7 Februar 2002 */ 2849 /* npwinc = 0, if curr_weight doesn't stay in the correct Groebner cone */ 3427 // npwinc = 0, if curr_weight doesn't stay in the correct Groebner cone 2850 3428 static ideal Rec_LastGB(ideal G, intvec* curr_weight, 2851 3429 intvec* orig_target_weight, int tp_deg, int npwinc) … … 2853 3431 BOOLEAN nError = Overflow_Error; 2854 3432 Overflow_Error = FALSE; 2855 BOOLEAN nOverflow_Error = FALSE;3433 // BOOLEAN nOverflow_Error = FALSE; 2856 3434 2857 3435 clock_t tproc=0; … … 2867 3445 intvec* ivNull = new intvec(nV); //define (0,...,0) 2868 3446 ring EXXRing = currRing; 2869 int NEG=0; //19 juni 03 2870 intvec* extra_curr_weight = new intvec(nV); 3447 //int NEG=0; //19 juni 03 2871 3448 intvec* next_weight; 2872 3449 #ifndef BUCHBERGER_ALG 2873 3450 //08 Juli 03 2874 3451 intvec* hilb_func; 2875 /* to avoid (1,0,...,0) as the target vector */ 3452 #endif 3453 // to avoid (1,0,...,0) as the target vector 2876 3454 intvec* last_omega = new intvec(nV); 2877 3455 for(i=nV-1; i>0; i--) … … 2881 3459 BOOLEAN isGB = FALSE; 2882 3460 2883 /* compute a pertubed weight vector of the target weight vector */ 2884 if(tp_deg > 1 && tp_deg <= nV) { 3461 // compute a pertubed weight vector of the target weight vector 3462 if(tp_deg > 1 && tp_deg <= nV) 3463 { 2885 3464 ideal H0 = idHeadCC(G); 2886 3465 2887 3466 if (rParameter (currRing) != NULL) 3467 { 2888 3468 DefRingParlp(); 3469 } 2889 3470 else 3471 { 2890 3472 VMrDefaultlp(); 2891 3473 } 2892 3474 TargetRing = currRing; 2893 3475 ssG = idrMoveR(G,EXXRing,currRing); … … 2896 3478 ideal H1 = idHeadCC(ssG); 2897 3479 2898 / * Apply Lemma 2.2 in Collart et. al (1997) to check whether2899 cone(k-1) is equal to cone(k) */2900 if(test_G_GB_walk(H0_tmp,H1)==1){3480 // Apply Lemma 2.2 in Collart et. al (1997) to check whether cone(k-1) is equal to cone(k) 3481 if(test_G_GB_walk(H0_tmp,H1)==1) 3482 { 2901 3483 idDelete(&H0_tmp); 2902 3484 idDelete(&H1); … … 2907 3489 2908 3490 if(npwinc != 0) 3491 { 2909 3492 goto LastGB_Finish; 2910 else { 3493 } 3494 else 3495 { 2911 3496 isGB = TRUE; 2912 3497 goto KSTD_Finish; … … 2923 3508 G = idrMoveR(ssG, TargetRing,currRing); 2924 3509 2925 if(Overflow_Error == TRUE) { 2926 nOverflow_Error = Overflow_Error; 2927 NEG = 1; 3510 if(Overflow_Error == TRUE) 3511 { 3512 //nOverflow_Error = Overflow_Error; 3513 //NEG = 1; 2928 3514 newRing = currRing; 2929 3515 goto JUNI_STD; … … 2937 3523 2938 3524 if(nwalk==1) 3525 { 2939 3526 goto FIRST_STEP; 2940 3527 } 2941 3528 to=clock(); 2942 / * compute an initial form ideal of <G> w.r.t. "curr_vector" */3529 // compute an initial form ideal of <G> w.r.t. "curr_vector" 2943 3530 Gomega = MwalkInitialForm(G, curr_weight); 2944 3531 xtif=xtif+clock()-to; … … 2946 3533 #ifndef BUCHBERGER_ALG 2947 3534 if(isNolVector(curr_weight) == 0) 3535 { 2948 3536 hilb_func = hFirstSeries(Gomega,NULL,NULL,curr_weight,currRing); 3537 } 2949 3538 else 3539 { 2950 3540 hilb_func = hFirstSeries(Gomega,NULL,NULL,last_omega,currRing); 3541 } 2951 3542 #endif // BUCHBERGER_ALG 2952 3543 2953 3544 oldRing = currRing; 2954 3545 2955 / * defiNe a new ring that its ordering is "(a(curr_weight),lp) */3546 // defiNe a new ring that its ordering is "(a(curr_weight),lp) 2956 3547 if (rParameter(currRing) != NULL) 3548 { 2957 3549 DefRingPar(curr_weight); 3550 } 2958 3551 else 3552 { 2959 3553 VMrDefault(curr_weight); 2960 3554 } 2961 3555 newRing = currRing; 2962 3556 Gomega1 = idrMoveR(Gomega, oldRing,currRing); 2963 3557 to=clock(); 2964 / * compute a reduced Groebner basis of <Gomega> w.r.t. "newRing" */3558 // compute a reduced Groebner basis of <Gomega> w.r.t. "newRing" 2965 3559 #ifdef BUCHBERGER_ALG 2966 3560 M = MstdhomCC(Gomega1); … … 2970 3564 #endif // BUCHBERGER_ALG 2971 3565 xtstd=xtstd+clock()-to; 2972 / * change the ring to oldRing */3566 // change the ring to oldRing 2973 3567 rChangeCurrRing(oldRing); 2974 3568 M1 = idrMoveR(M, newRing,currRing); … … 2976 3570 2977 3571 to=clock(); 2978 /* compute a reduced Groebner basis of <G> w.r.t. "newRing" by the 2979 lifting process*/ 3572 // compute a reduced Groebner basis of <G> w.r.t. "newRing" by the lifting process 2980 3573 F = MLifttwoIdeal(Gomega2, M1, G); 2981 3574 xtlift=xtlift+clock()-to; … … 2984 3577 idDelete(&G); 2985 3578 2986 / * change the ring to newRing */3579 // change the ring to newRing 2987 3580 rChangeCurrRing(newRing); 2988 3581 F1 = idrMoveR(F, oldRing,currRing); 2989 3582 2990 3583 to=clock(); 2991 / * reduce the Groebner basis <G> w.r.t. new ring */3584 // reduce the Groebner basis <G> w.r.t. new ring 2992 3585 G = kInterRedCC(F1, NULL); 2993 3586 xtred=xtred+clock()-to; … … 2995 3588 2996 3589 if(endwalks == 1) 3590 { 2997 3591 break; 2998 3592 } 2999 3593 FIRST_STEP: 3000 3594 to=clock(); 3001 3595 Overflow_Error = FALSE; 3002 / * compute a next weight vector */3596 // compute a next weight vector 3003 3597 next_weight = MkInterRedNextWeight(curr_weight,target_weight, G); 3004 3598 xtnw=xtnw+clock()-to; … … 3007 3601 #endif 3008 3602 3009 if(Overflow_Error == TRUE) { 3603 if(Overflow_Error == TRUE) 3604 { 3010 3605 //PrintS("\n// ** The next vector does NOT stay in Cone!!\n"); 3011 3606 #ifdef TEST_OVERFLOW … … 3015 3610 nnwinC = 0; 3016 3611 if(tp_deg == nV) 3612 { 3017 3613 nlast = 1; 3018 3614 } 3019 3615 delete next_weight; 3020 3616 break; 3021 3617 } 3022 3618 3023 if(MivComp(next_weight, ivNull) == 1) { 3619 if(MivComp(next_weight, ivNull) == 1) 3620 { 3024 3621 //newRing = currRing; 3025 3622 delete next_weight; … … 3027 3624 } 3028 3625 3029 if(MivComp(next_weight, target_weight) == 1) { 3626 if(MivComp(next_weight, target_weight) == 1) 3627 { 3030 3628 if(tp_deg == nV) 3629 { 3031 3630 endwalks = 1; 3032 else { 3033 REC_LAST_GB_ALT2: 3034 nOverflow_Error = Overflow_Error; 3631 } 3632 else 3633 { 3634 // REC_LAST_GB_ALT2: 3635 //nOverflow_Error = Overflow_Error; 3035 3636 tproc=tproc+clock()-tinput; 3036 3637 /* … … 3045 3646 } 3046 3647 3047 for(i=nV-1; i>=0; i--) {3048 //(*extra_curr_weight)[i] = (*curr_weight)[i];3648 for(i=nV-1; i>=0; i--) 3649 { 3049 3650 (*curr_weight)[i] = (*next_weight)[i]; 3050 3651 } … … 3054 3655 delete ivNull; 3055 3656 3056 if(tp_deg != nV) { 3657 if(tp_deg != nV) 3658 { 3057 3659 newRing = currRing; 3058 3660 3059 3661 if (rParameter(currRing) != NULL) 3662 { 3060 3663 DefRingParlp(); 3664 } 3061 3665 else 3666 { 3062 3667 VMrDefaultlp(); 3063 3668 } 3064 3669 F1 = idrMoveR(G, newRing,currRing); 3065 3670 3066 if(nnwinC == 0 || test_w_in_ConeCC(F1, target_weight) != 1 ) { 3067 nOverflow_Error = Overflow_Error; 3671 if(nnwinC == 0 || test_w_in_ConeCC(F1, target_weight) != 1 ) 3672 { 3673 // nOverflow_Error = Overflow_Error; 3068 3674 //Print("\n// takes %d steps and calls \"Rec_LastGB (%d):", tp_deg+1); 3069 3675 tproc=tproc+clock()-tinput; … … 3076 3682 result = idrMoveR(F1, TargetRing,currRing); 3077 3683 } 3078 else { 3079 if(nlast == 1) { 3684 else 3685 { 3686 if(nlast == 1) 3687 { 3080 3688 JUNI_STD: 3081 3689 3082 3690 newRing = currRing; 3083 3691 if (rParameter(currRing) != NULL) 3692 { 3084 3693 DefRingParlp(); 3694 } 3085 3695 else 3696 { 3086 3697 VMrDefaultlp(); 3087 3698 } 3088 3699 KSTD_Finish: 3089 3700 if(isGB == FALSE) 3701 { 3090 3702 F1 = idrMoveR(G, newRing,currRing); 3703 } 3091 3704 else 3705 { 3092 3706 F1 = G; 3707 } 3093 3708 to=clock(); 3094 3709 // Print("\n// apply the Buchberger's alg in ring = %s",rString(currRing)); … … 3108 3723 3109 3724 if(Overflow_Error == FALSE) 3725 { 3110 3726 Overflow_Error=nError; 3111 /* 3112 Print("\n// \"Rec_LastGB\" (%d) took %d steps and %.2f sec.Overflow_Error (%d)", tp_deg, 3113 nwalk, ((double) tproc)/1000000, nOverflow_Error); 3114 */ 3727 } 3728 // Print("\n// \"Rec_LastGB\" (%d) took %d steps and %.2f sec.Overflow_Error (%d)", tp_deg, nwalk, ((double) tproc)/1000000, nOverflow_Error); 3115 3729 return(result); 3116 3730 } … … 3124 3738 basis or n times, where n is the numbers of variables. 3125 3739 */ 3126 /* walk + recursive LastGB */ 3740 3741 /****************************** 3742 * walk + recursive LastGB * 3743 ******************************/ 3127 3744 ideal MAltwalk2(ideal Go, intvec* curr_weight, intvec* target_weight) 3128 3745 { 3129 3746 Set_Error(FALSE); 3130 3747 Overflow_Error = FALSE; 3131 BOOLEAN nOverflow_Error = FALSE;3748 //BOOLEAN nOverflow_Error = FALSE; 3132 3749 //Print("// pSetm_Error = (%d)", ErrorCheck()); 3133 3750 … … 3138 3755 nstep = 0; 3139 3756 int i, nV = currRing->N; 3140 int nwalk=0, endwalks=0, nhilb=1; 3141 3142 ideal Gomega, M, F, Gomega1, Gomega2, M1, F1, G, G1; 3143 ring endRing, newRing, oldRing; 3757 int nwalk=0, endwalks=0; 3758 // int nhilb = 1; 3759 ideal Gomega, M, F, Gomega1, Gomega2, M1, F1, G; 3760 //ideal G1; 3761 //ring endRing; 3762 ring newRing, oldRing; 3144 3763 intvec* ivNull = new intvec(nV); 3145 3764 intvec* next_weight; 3765 #if 0 3146 3766 intvec* extra_curr_weight = new intvec(nV); 3147 intvec* hilb_func; 3767 #endif 3768 //intvec* hilb_func; 3148 3769 intvec* exivlp = Mivlp(nV); 3149 3770 … … 3162 3783 */ 3163 3784 if(currRing->order[0] == ringorder_a) 3785 { 3164 3786 goto NEXT_VECTOR; 3165 3787 } 3166 3788 while(1) 3167 3789 { … … 3185 3807 /* define a new ring that its ordering is "(a(curr_weight),lp) */ 3186 3808 if (rParameter(currRing) != NULL) 3809 { 3187 3810 DefRingPar(curr_weight); 3811 } 3188 3812 else 3813 { 3189 3814 VMrDefault(curr_weight); 3190 3815 } 3191 3816 newRing = currRing; 3192 3817 Gomega1 = idrMoveR(Gomega, oldRing,currRing); … … 3241 3866 #endif 3242 3867 3243 3244 3868 newRing = currRing; 3245 3869 if (rParameter(currRing) != NULL) 3870 { 3246 3871 DefRingPar(target_weight); 3872 } 3247 3873 else 3874 { 3248 3875 VMrDefault(target_weight); 3249 3876 } 3250 3877 F1 = idrMoveR(G, newRing,currRing); 3251 3878 G = MstdCC(F1); … … 3266 3893 if(MivSame(target_weight, exivlp)==1) 3267 3894 { 3268 LAST_GB_ALT2:3269 nOverflow_Error = Overflow_Error;3895 // LAST_GB_ALT2: 3896 //nOverflow_Error = Overflow_Error; 3270 3897 tproc = clock()-xftinput; 3271 3898 //Print("\n// takes %d steps and calls the recursion of level 2:", nwalk); … … 3286 3913 delete next_weight; 3287 3914 } 3915 #ifdef TEST_OVERFLOW 3288 3916 TEST_OVERFLOW_OI: 3917 #endif 3289 3918 rChangeCurrRing(XXRing); 3290 3919 G = idrMoveR(G, newRing,currRing); … … 3293 3922 3294 3923 #ifdef TIME_TEST 3295 Print("\n// \"Main procedure\" took %d steps dnd %.2f sec. Overflow_Error (%d)", nwalk, 3296 ((double) tproc)/1000000, nOverflow_Error); 3924 // Print("\n// \"Main procedure\" took %d steps dnd %.2f sec. Overflow_Error (%d)", nwalk, ((double) tproc)/1000000, nOverflow_Error); 3297 3925 3298 3926 TimeStringFractal(xftinput, tostd, xtif, xtstd, xtextra,xtlift, xtred,xtnw); … … 3307 3935 3308 3936 3309 /* 5.5.02 */ 3310 /* The implementation of the fractal walk algorithmus */ 3311 3312 /* The main procedur Mfwalk calls the recursive Subroutine 3313 rec_fractal_call to compute the wanted Grï¿œbner basis. 3314 At the main procedur we compute the reduced Grï¿œbner basis w.r.t. a "fast" 3315 order, e.g. "dp" and a sequence of weight vectors which are row vectors 3316 of a matrix. This matrix defines the given monomial order, e.g. "lp" 3317 */ 3318 3319 3320 3321 3322 /* perturb the matrix order of "lex" */ 3937 /************************************** 3938 * perturb the matrix order of "lex" * 3939 **************************************/ 3323 3940 static intvec* NewVectorlp(ideal I) 3324 3941 { … … 3337 3954 intvec* Xivlp; 3338 3955 3339 3340 3341 /*********************************************************************** 3342 The procedur REC_GB_Mwalk computes a GB for <G> w.r.t. the weight order 3343 otw, where G is a reduced GB w.r.t. the weight order cw. 3344 The new procedur Mwalk calls REC_GB. 3345 ************************************************************************/ 3956 /******************************** 3957 * compute a next weight vector * 3958 ********************************/ 3959 static intvec* MWalkRandomNextWeight(ideal G, intvec* curr_weight, intvec* target_weight, int weight_rad, int pert_deg) 3960 { 3961 int i, weight_norm; 3962 int nV = currRing->N; 3963 intvec* next_weight2; 3964 intvec* next_weight22 = new intvec(nV); 3965 intvec* next_weight = MwalkNextWeightCC(curr_weight,target_weight, G); 3966 if(MivComp(next_weight, target_weight) == 1) 3967 { 3968 return(next_weight); 3969 } 3970 else 3971 { 3972 //compute a perturbed next weight vector "next_weight1" 3973 intvec* next_weight1 = MkInterRedNextWeight(MPertVectors(G, MivMatrixOrder(curr_weight), pert_deg), target_weight, G); 3974 //Print("\n // size of next_weight1 = %d", sizeof((*next_weight1))); 3975 3976 //compute a random next weight vector "next_weight2" 3977 while(1) 3978 { 3979 weight_norm = 0; 3980 while(weight_norm == 0) 3981 { 3982 for(i=0; i<nV; i++) 3983 { 3984 //Print("\n// next_weight[%d] = %d", i, (*next_weight)[i]); 3985 (*next_weight22)[i] = rand() % 60000 - 30000; 3986 weight_norm = weight_norm + (*next_weight22)[i]*(*next_weight22)[i]; 3987 } 3988 weight_norm = 1 + floor(sqrt(weight_norm)); 3989 } 3990 3991 for(i=nV-1; i>=0; i--) 3992 { 3993 if((*next_weight22)[i] < 0) 3994 { 3995 (*next_weight22)[i] = 1 + (*curr_weight)[i] + floor(weight_rad*(*next_weight22)[i]/weight_norm); 3996 } 3997 else 3998 { 3999 (*next_weight22)[i] = (*curr_weight)[i] + floor(weight_rad*(*next_weight22)[i]/weight_norm); 4000 } 4001 //Print("\n// next_weight22[%d] = %d", i, (*next_weight22)[i]); 4002 } 4003 4004 if(test_w_in_ConeCC(G, next_weight22) == 1) 4005 { 4006 //Print("\n//MWalkRandomNextWeight: next_weight2 im Kegel\n"); 4007 next_weight2 = MkInterRedNextWeight(next_weight22, target_weight, G); 4008 delete next_weight22; 4009 break; 4010 } 4011 } 4012 intvec* result = new intvec(nV); 4013 ideal G_test = MwalkInitialForm(G, next_weight); 4014 ideal G_test1 = MwalkInitialForm(G, next_weight1); 4015 ideal G_test2 = MwalkInitialForm(G, next_weight2); 4016 4017 // compare next_weights 4018 if(IDELEMS(G_test1) < IDELEMS(G_test)) 4019 { 4020 if(IDELEMS(G_test2) <= IDELEMS(G_test1)) // |G_test2| <= |G_test1| < |G_test| 4021 { 4022 for(i=0; i<nV; i++) 4023 { 4024 (*result)[i] = (*next_weight2)[i]; 4025 } 4026 } 4027 else // |G_test1| < |G_test|, |G_test1| < |G_test2| 4028 { 4029 for(i=0; i<nV; i++) 4030 { 4031 (*result)[i] = (*next_weight1)[i]; 4032 } 4033 } 4034 } 4035 else 4036 { 4037 if(IDELEMS(G_test2) <= IDELEMS(G_test)) // |G_test2| <= |G_test| <= |G_test1| 4038 { 4039 for(i=0; i<nV; i++) 4040 { 4041 (*result)[i] = (*next_weight2)[i]; 4042 } 4043 } 4044 else // |G_test| <= |G_test1|, |G_test| < |G_test2| 4045 { 4046 for(i=0; i<nV; i++) 4047 { 4048 (*result)[i] = (*next_weight)[i]; 4049 } 4050 } 4051 } 4052 delete next_weight; 4053 delete next_weight1; 4054 idDelete(&G_test); 4055 idDelete(&G_test1); 4056 idDelete(&G_test2); 4057 if(test_w_in_ConeCC(G, result) == 1) 4058 { 4059 delete next_weight2; 4060 return result; 4061 } 4062 else 4063 { 4064 delete result; 4065 return next_weight2; 4066 } 4067 } 4068 } 4069 4070 4071 /*************************************************************************** 4072 * The procedur REC_GB_Mwalk computes a GB for <G> w.r.t. the weight order * 4073 * otw, where G is a reduced GB w.r.t. the weight order cw. * 4074 * The new procedur Mwalk calls REC_GB. * 4075 ***************************************************************************/ 3346 4076 static ideal REC_GB_Mwalk(ideal G, intvec* curr_weight, intvec* orig_target_weight, 3347 4077 int tp_deg, int npwinc) … … 3350 4080 Overflow_Error = FALSE; 3351 4081 3352 int i, nV = currRing->N , ntwC, npwinC;4082 int i, nV = currRing->N; 3353 4083 int nwalk=0, endwalks=0, nnwinC=1, nlast = 0; 3354 4084 ideal Gomega, M, F, Gomega1, Gomega2, M1,F1,result,ssG; 3355 4085 ring newRing, oldRing, TargetRing; 3356 intvec* iv_M_lp;3357 4086 intvec* target_weight; 3358 4087 intvec* ivNull = new intvec(nV); 4088 #ifndef BUCHBERGER_ALG 3359 4089 intvec* hilb_func; 3360 BOOLEAN isGB = FALSE; 3361 3362 /* to avoid (1,0,...,0) as the target vector */ 4090 // to avoid (1,0,...,0) as the target vector 3363 4091 intvec* last_omega = new intvec(nV); 3364 4092 for(i=nV-1; i>0; i--) 4093 { 3365 4094 (*last_omega)[i] = 1; 4095 } 3366 4096 (*last_omega)[0] = 10000; 4097 #endif 4098 BOOLEAN isGB = FALSE; 3367 4099 3368 4100 ring EXXRing = currRing; 3369 4101 3370 / * compute a pertubed weight vector of the target weight vector */4102 // compute a pertubed weight vector of the target weight vector 3371 4103 if(tp_deg > 1 && tp_deg <= nV) 3372 4104 { 3373 4105 ideal H0 = idHeadCC(G); 3374 4106 if (rParameter(currRing) != NULL) 4107 { 3375 4108 DefRingPar(orig_target_weight); 4109 } 3376 4110 else 4111 { 3377 4112 VMrDefault(orig_target_weight); 3378 4113 } 3379 4114 TargetRing = currRing; 3380 4115 ssG = idrMoveR(G,EXXRing,currRing); … … 3386 4121 if(test_G_GB_walk(H0_tmp,H1)==1) 3387 4122 { 3388 //Print("//input in %d-th recursive is a GB",tp_deg); 3389 idDelete(&H0_tmp);idDelete(&H1); 4123 Print("\n//REC_GB_Mwalk: input in %d-th recursive is a GB!\n",tp_deg); 4124 idDelete(&H0_tmp); 4125 idDelete(&H1); 3390 4126 G = ssG; 3391 4127 ssG = NULL; … … 3398 4134 } 3399 4135 else 4136 { 3400 4137 goto LastGB_Finish; 3401 } 3402 idDelete(&H0_tmp); idDelete(&H1); 3403 3404 intvec* ivlp = Mivlp(nV); 3405 if( MivSame(orig_target_weight, ivlp)==1 ) 3406 iv_M_lp = MivMatrixOrderlp(nV); 3407 else 3408 iv_M_lp = MivMatrixOrder(orig_target_weight); 3409 3410 //target_weight = MPertVectorslp(ssG, iv_M_lp, tp_deg); 3411 target_weight = MPertVectors(ssG, iv_M_lp, tp_deg); 3412 3413 delete ivlp; 3414 delete iv_M_lp; 4138 } 4139 } 4140 idDelete(&H0_tmp); 4141 idDelete(&H1); 4142 4143 target_weight = MPertVectors(ssG, MivMatrixOrder(orig_target_weight), tp_deg); 3415 4144 3416 4145 rChangeCurrRing(EXXRing); … … 3423 4152 nstep++; 3424 4153 if(nwalk == 1) 4154 { 3425 4155 goto NEXT_STEP; 3426 3427 //Print("\n// Entering the %d-th step in the %d-th recursive:",nwalk,tp_deg);4156 } 4157 Print("\n//REC_GB_Mwalk: Entering the %d-th step in the %d-th recursive:\n",nwalk,tp_deg); 3428 4158 to = clock(); 3429 / * compute an initial form ideal of <G> w.r.t. "curr_vector" */4159 // compute an initial form ideal of <G> w.r.t. "curr_vector" 3430 4160 Gomega = MwalkInitialForm(G, curr_weight); 3431 4161 xtif = xtif + clock()-to; … … 3433 4163 #ifndef BUCHBERGER_ALG 3434 4164 if(isNolVector(curr_weight) == 0) 4165 { 3435 4166 hilb_func = hFirstSeries(Gomega,NULL,NULL,curr_weight,currRing); 4167 } 3436 4168 else 4169 { 3437 4170 hilb_func = hFirstSeries(Gomega,NULL,NULL,last_omega,currRing); 3438 #endif // BUCHBERGER_ALG 4171 } 4172 #endif 3439 4173 3440 4174 oldRing = currRing; 3441 4175 3442 / * define a new ring that its ordering is "(a(curr_weight),lp) */4176 // define a new ring with ordering "(a(curr_weight),lp) 3443 4177 if (rParameter(currRing) != NULL) 4178 { 3444 4179 DefRingPar(curr_weight); 4180 } 3445 4181 else 4182 { 3446 4183 VMrDefault(curr_weight); 3447 4184 } 3448 4185 newRing = currRing; 3449 4186 Gomega1 = idrMoveR(Gomega, oldRing,currRing); 3450 4187 3451 4188 to = clock(); 3452 / * compute a reduced Groebner basis of <Gomega> w.r.t. "newRing" */4189 // compute a reduced Groebner basis of <Gomega> w.r.t. "newRing" 3453 4190 #ifdef BUCHBERGER_ALG 3454 4191 M = MstdhomCC(Gomega1); … … 3456 4193 M=kStd(Gomega1,NULL,isHomog,NULL,hilb_func,0,NULL,curr_weight); 3457 4194 delete hilb_func; 3458 #endif // BUCHBERGER_ALG4195 #endif 3459 4196 xtstd = xtstd + clock() - to; 3460 4197 3461 / * change the ring to oldRing */4198 // change the ring to oldRing 3462 4199 rChangeCurrRing(oldRing); 3463 4200 … … 3474 4211 3475 4212 3476 / * change the ring to newRing */4213 // change the ring to newRing 3477 4214 rChangeCurrRing(newRing); 3478 4215 F1 = idrMoveR(F, oldRing,currRing); 3479 4216 3480 4217 to = clock(); 3481 / * reduce the Groebner basis <G> w.r.t. new ring */4218 // reduce the Groebner basis <G> w.r.t. new ring 3482 4219 G = kInterRedCC(F1, NULL); 3483 4220 xtred = xtred + clock() -to; … … 3486 4223 3487 4224 if(endwalks == 1) 4225 { 3488 4226 break; 3489 4227 } 3490 4228 NEXT_STEP: 3491 4229 to = clock(); 3492 / * compute a next weight vector */4230 // compute a next weight vector 3493 4231 intvec* next_weight = MkInterRedNextWeight(curr_weight,target_weight, G); 4232 4233 3494 4234 xtnw = xtnw + clock() - to; 3495 4235 … … 3498 4238 #endif 3499 4239 3500 /*check whether the computed vector does in the correct cone */3501 //ntwC = test_w_in_ConeCC(G, next_weight);3502 //if(ntwC != 1)3503 4240 if(Overflow_Error == TRUE) 3504 4241 { 3505 PrintS("\n// ** The computed vector does NOT stay in Cone!!\n");4242 PrintS("\n//REC_GB_Mwalk: The computed vector does NOT stay in the correct cone!!\n"); 3506 4243 nnwinC = 0; 3507 4244 if(tp_deg == nV) 4245 { 3508 4246 nlast = 1; 4247 } 3509 4248 delete next_weight; 3510 4249 break; … … 3520 4259 { 3521 4260 if(tp_deg == nV) 4261 { 3522 4262 endwalks = 1; 3523 else { 4263 } 4264 else 4265 { 3524 4266 G = REC_GB_Mwalk(G,curr_weight, orig_target_weight, tp_deg+1,nnwinC); 3525 4267 newRing = currRing; … … 3529 4271 } 3530 4272 3531 /* 06.11.01 NOT Changed */3532 4273 for(i=nV-1; i>=0; i--) 4274 { 3533 4275 (*curr_weight)[i] = (*next_weight)[i]; 3534 4276 } 3535 4277 delete next_weight; 3536 } //while4278 } 3537 4279 3538 4280 delete ivNull; … … 3540 4282 if(tp_deg != nV) 3541 4283 { 3542 //28.07.033543 4284 newRing = currRing; 3544 4285 3545 4286 if (rParameter(currRing) != NULL) 3546 // DefRingParlp(); //4287 { 3547 4288 DefRingPar(orig_target_weight); 4289 } 3548 4290 else 4291 { 3549 4292 VMrDefault(orig_target_weight); 3550 3551 4293 } 3552 4294 F1 = idrMoveR(G, newRing,currRing); 3553 4295 3554 4296 if(nnwinC == 0) 4297 { 3555 4298 F1 = REC_GB_Mwalk(F1,curr_weight, orig_target_weight, tp_deg+1,nnwinC); 4299 } 3556 4300 else 4301 { 3557 4302 if(test_w_in_ConeCC(F1, target_weight) != 1) 4303 { 3558 4304 F1 = REC_GB_Mwalk(F1,curr_weight, orig_target_weight,tp_deg+1,nnwinC); 3559 4305 } 4306 } 3560 4307 delete target_weight; 3561 4308 … … 3569 4316 { 3570 4317 if (rParameter(currRing) != NULL) 4318 { 3571 4319 DefRingPar(orig_target_weight); 4320 } 3572 4321 else 4322 { 3573 4323 VMrDefault(orig_target_weight); 3574 3575 4324 } 3576 4325 KSTD_Finish: 3577 4326 if(isGB == FALSE) 4327 { 3578 4328 F1 = idrMoveR(G, newRing,currRing); 4329 } 3579 4330 else 4331 { 3580 4332 F1 = G; 4333 } 3581 4334 to=clock(); 3582 / * apply Buchberger alg to compute a red. GB of F1 */4335 // apply Buchberger alg to compute a red. GB of F1 3583 4336 G = MstdCC(F1); 3584 4337 xtextra=clock()-to; … … 3593 4346 3594 4347 if(Overflow_Error == FALSE) 4348 { 3595 4349 Overflow_Error = nError; 3596 4350 } 4351 #ifndef BUCHBERGER_ALG 4352 delete last_omega; 4353 #endif 3597 4354 return(result); 3598 4355 } 3599 4356 3600 4357 3601 /* 08.09.02 */ 3602 /******** THE NEW GRï¿œBNER WALK ALGORITHM **********/ 3603 /* Grï¿œbnerwalk with a recursive "second" alternative GW, REC_GB_Mwalk 3604 that only computes the last reduced GB */ 4358 // THE NEW GROEBNER WALK ALGORITHM 4359 // Groebnerwalk with a recursive "second" alternative GW, called REC_GB_Mwalk that only computes the last reduced GB 3605 4360 ideal Mwalk(ideal Go, intvec* curr_weight, intvec* target_weight) 3606 4361 { … … 3614 4369 clock_t tim; 3615 4370 nstep=0; 3616 int i, nV = currRing->N; 3617 int nwalk=0, endwalks=0; 3618 3619 ideal Gomega, M, F, Gomega1, Gomega2, M1, F1, G, G1; 3620 ring endRing, newRing, oldRing; 4371 int i; 4372 int nV = currRing->N; 4373 int nwalk=0; 4374 int endwalks=0; 4375 4376 ideal Gomega, M, F, Gomega1, Gomega2, M1, F1, G; 4377 //ideal G1; 4378 //ring endRing; 4379 ring newRing, oldRing; 3621 4380 intvec* ivNull = new intvec(nV); 3622 4381 intvec* exivlp = Mivlp(nV); 4382 #ifndef BUCHBERGER_ALG 3623 4383 intvec* hilb_func; 3624 4384 #endif 3625 4385 intvec* tmp_weight = new intvec(nV); 3626 4386 for(i=nV-1; i>=0; i--) 3627 4387 (*tmp_weight)[i] = (*curr_weight)[i]; 3628 4388 3629 / * to avoid (1,0,...,0) as the target vector */4389 // to avoid (1,0,...,0) as the target vector 3630 4390 intvec* last_omega = new intvec(nV); 3631 4391 for(i=nV-1; i>0; i--) … … 3636 4396 3637 4397 to = clock(); 3638 / * the monomial ordering of this current ring would be "dp" */4398 // the monomial ordering of this current ring would be "dp" 3639 4399 G = MstdCC(Go); 3640 4400 tostd = clock()-to; … … 3648 4408 nstep ++; 3649 4409 to = clock(); 3650 / * compute an initial form ideal of <G> w.r.t. "curr_vector" */4410 // compute an initial form ideal of <G> w.r.t. "curr_vector" 3651 4411 Gomega = MwalkInitialForm(G, curr_weight); 3652 4412 tif = tif + clock()-to; … … 3691 4451 /* create a new ring newRing */ 3692 4452 if (rParameter(currRing) != NULL) 4453 { 3693 4454 DefRingPar(curr_weight); 4455 } 3694 4456 else 4457 { 3695 4458 VMrDefault(curr_weight); 3696 4459 } 3697 4460 newRing = currRing; 3698 4461 F1 = idrMoveR(F, oldRing,currRing); … … 3703 4466 #ifndef BUCHBERGER_ALG 3704 4467 if(isNolVector(curr_weight) == 0) 4468 { 3705 4469 hilb_func = hFirstSeries(Gomega,NULL,NULL,curr_weight,currRing); 4470 } 3706 4471 else 4472 { 3707 4473 hilb_func = hFirstSeries(Gomega,NULL,NULL,last_omega,currRing); 4474 } 3708 4475 #endif // BUCHBERGER_ALG 3709 4476 3710 / * define a new ring that its ordering is "(a(curr_weight),lp) */4477 // define a new ring that its ordering is "(a(curr_weight),lp) 3711 4478 if (rParameter(currRing) != NULL) 4479 { 3712 4480 DefRingPar(curr_weight); 4481 } 3713 4482 else 4483 { 3714 4484 VMrDefault(curr_weight); 3715 4485 } 3716 4486 newRing = currRing; 3717 4487 Gomega1 = idrMoveR(Gomega, oldRing,currRing); 3718 4488 3719 4489 to = clock(); 3720 / * compute a reduced Groebner basis of <Gomega> w.r.t. "newRing" */4490 // compute a reduced Groebner basis of <Gomega> w.r.t. "newRing" 3721 4491 #ifdef BUCHBERGER_ALG 3722 4492 M = MstdhomCC(Gomega1); … … 3727 4497 tstd = tstd + clock() - to; 3728 4498 3729 / * change the ring to oldRing */4499 // change the ring to oldRing 3730 4500 rChangeCurrRing(oldRing); 3731 4501 M1 = idrMoveR(M, newRing,currRing); … … 3733 4503 3734 4504 to = clock(); 3735 /* compute a representation of the generators of submod (M) 3736 with respect to those of mod (Gomega). 3737 Gomega is a reduced Groebner basis w.r.t. the current ring */ 4505 // compute a representation of the generators of submod (M) with respect to those of mod (Gomega). Gomega is a reduced Groebner basis w.r.t. the current ring. 3738 4506 F = MLifttwoIdeal(Gomega2, M1, G); 3739 4507 tlift = tlift + clock() - to; … … 3743 4511 idDelete(&G); 3744 4512 3745 / * change the ring to newRing */4513 // change the ring to newRing 3746 4514 rChangeCurrRing(newRing); 3747 4515 F1 = idrMoveR(F, oldRing,currRing); … … 3749 4517 3750 4518 to = clock(); 3751 / * reduce the Groebner basis <G> w.r.t. new ring */4519 // reduce the Groebner basis <G> w.r.t. new ring 3752 4520 G = kInterRedCC(F1, NULL); 3753 4521 if(endwalks != 1) 4522 { 3754 4523 tred = tred + clock() - to; 4524 } 3755 4525 else 4526 { 3756 4527 xtred = xtred + clock() - to; 3757 4528 } 3758 4529 idDelete(&F1); 3759 4530 if(endwalks == 1) 4531 { 3760 4532 break; 3761 4533 } 3762 4534 NEXT_VECTOR: 3763 4535 to = clock(); 3764 / * compute a next weight vector */4536 // compute a next weight vector 3765 4537 intvec* next_weight = MkInterRedNextWeight(curr_weight,target_weight,G); 3766 4538 tnw = tnw + clock() - to; … … 3776 4548 3777 4549 if (rParameter(currRing) != NULL) 4550 { 3778 4551 DefRingPar(target_weight); 4552 } 3779 4553 else 4554 { 3780 4555 VMrDefault(target_weight); 3781 4556 } 3782 4557 F1 = idrMoveR(G, newRing,currRing); 3783 4558 G = MstdCC(F1); … … 3795 4570 } 3796 4571 if(MivComp(next_weight, target_weight) == 1) 4572 { 3797 4573 endwalks = 1; 3798 3799 for(i=nV-1; i>=0; i--) { 4574 } 4575 for(i=nV-1; i>=0; i--) 4576 { 3800 4577 (*tmp_weight)[i] = (*curr_weight)[i]; 3801 4578 (*curr_weight)[i] = (*next_weight)[i]; … … 3819 4596 } 3820 4597 3821 /* 2.12.02*/ 4598 // 07.11.2012 4599 // THE RANDOM WALK ALGORITHM 4600 ideal Mrwalk(ideal Go, intvec* curr_weight, intvec* target_weight, int weight_rad, int pert_deg) 4601 { 4602 Set_Error(FALSE); 4603 Overflow_Error = FALSE; 4604 //Print("// pSetm_Error = (%d)", ErrorCheck()); 4605 4606 clock_t tinput, tostd, tif=0, tstd=0, tlift=0, tred=0, tnw=0; 4607 xtif=0; xtstd=0; xtlift=0; xtred=0; xtnw=0; 4608 tinput = clock(); 4609 clock_t tim; 4610 nstep=0; 4611 int i,k,nwalk,endwalks = 0; 4612 int nV = currRing->N; 4613 4614 ideal Gomega, M, F, Gomega1, Gomega2, M1, F1, G; 4615 //ideal G1; 4616 //ring endRing; 4617 ring newRing, oldRing; 4618 intvec* ivNull = new intvec(nV); 4619 intvec* exivlp = Mivlp(nV); 4620 intvec* tmp_weight = new intvec(nV); 4621 for(i=nV-1; i>=0; i--) 4622 { 4623 (*tmp_weight)[i] = (*curr_weight)[i]; 4624 } 4625 #ifndef BUCHBERGER_ALG 4626 intvec* hilb_func; 4627 // to avoid (1,0,...,0) as the target vector 4628 intvec* last_omega = new intvec(nV); 4629 for(i=nV-1; i>0; i--) 4630 { 4631 (*last_omega)[i] = 1; 4632 } 4633 (*last_omega)[0] = 10000; 4634 #endif 4635 ring XXRing = currRing; 4636 4637 to = clock(); 4638 G = MstdCC(Go); 4639 tostd = clock()-to; 4640 4641 //if(currRing->order[0] == ringorder_a) 4642 // goto NEXT_VECTOR; 4643 nwalk = 0; 4644 while(1) 4645 { 4646 nwalk ++; 4647 nstep ++; 4648 to = clock(); 4649 #ifdef CHECK_IDEAL_MWALK 4650 idString(G,"G"); 4651 #endif 4652 Gomega = MwalkInitialForm(G, curr_weight);// compute an initial form ideal of <G> w.r.t. "curr_vector" 4653 tif = tif + clock()-to; 4654 oldRing = currRing; 4655 #ifdef CHECK_IDEAL_MWALK 4656 idString(Gomega,"G_omega"); 4657 #endif 4658 if(endwalks == 1) 4659 { 4660 // compute a reduced Groebner basis of Gomega w.r.t. >>_cw by the recursive changed perturbation walk alg. 4661 tim = clock(); 4662 Print("\n//Mrwalk: Groebnerwalk took %d steps.\n", nwalk); 4663 //PrintS("\n//Mrwalk: call the rec. Pert. Walk to compute a red GB of:"); 4664 //idString(Gomega, "G_omega"); 4665 M = REC_GB_Mwalk(idCopy(Gomega), tmp_weight, curr_weight,pert_deg,1); 4666 Print("\n//Mrwalk: time for the last std(Gw) = %.2f sec\n",((double) (clock()-tim)/1000000)); 4667 4668 #ifdef CHECK_IDEAL_MWALK 4669 idString(Gomega, "G_omega"); 4670 idString(M, "M"); 4671 #endif 4672 to = clock(); 4673 F = MLifttwoIdeal(Gomega, M, G); 4674 xtlift = xtlift + clock() - to; 4675 4676 idDelete(&Gomega); 4677 idDelete(&M); 4678 idDelete(&G); 4679 4680 oldRing = currRing; 4681 VMrDefault(curr_weight); //define a new ring with ordering "(a(curr_weight),lp) 4682 4683 /*if (rParameter(currRing) != NULL) 4684 { 4685 DefRingPar(curr_weight); 4686 } 4687 else 4688 { 4689 VMrDefault(curr_weight); 4690 }*/ 4691 newRing = currRing; 4692 F1 = idrMoveR(F, oldRing,currRing); 4693 } 4694 else 4695 { 4696 NORMAL_GW: 4697 #ifndef BUCHBERGER_ALG 4698 if(isNolVector(curr_weight) == 0) 4699 { 4700 hilb_func = hFirstSeries(Gomega,NULL,NULL,curr_weight,currRing); 4701 } 4702 else 4703 { 4704 hilb_func = hFirstSeries(Gomega,NULL,NULL,last_omega,currRing); 4705 } 4706 #endif 4707 /* 4708 if (rParameter(currRing) != NULL) 4709 { 4710 DefRingPar(curr_weight); 4711 } 4712 else 4713 { 4714 VMrDefault(curr_weight); 4715 }*/ 4716 4717 VMrDefault(curr_weight); //define a new ring with ordering "(a(curr_weight),lp) 4718 newRing = currRing; 4719 Gomega1 = idrMoveR(Gomega, oldRing,currRing); 4720 #ifdef CHECK_IDEAL_MWALK 4721 idString(Gomega1, "G_omega1"); 4722 #endif 4723 // compute a reduced Groebner basis of <Gomega> w.r.t. "newRing" 4724 to = clock(); 4725 #ifndef BUCHBERGER_ALG 4726 M=kStd(Gomega1,NULL,isHomog,NULL,hilb_func,0,NULL,curr_weight); 4727 delete hilb_func; 4728 #else 4729 //M = MstdhomCC(Gomega1); 4730 //M = MstdCC(Gomega1); 4731 //M = kStd(Gomega1, NULL, testHomog,NULL, NULL,0,0,curr_weight); 4732 M = kStd(Gomega1,NULL,testHomog,NULL,NULL,0,0,NULL); 4733 #endif 4734 tstd = tstd + clock() - to; 4735 #ifdef CHECK_IDEAL_MWALK 4736 idString(M, "M"); 4737 #endif 4738 //change the ring to oldRing 4739 rChangeCurrRing(oldRing); 4740 M1 = idrMoveR(M, newRing,currRing); 4741 #ifdef CHECK_IDEAL_MWALK 4742 idString(M1, "M1"); 4743 #endif 4744 Gomega2 = idrMoveR(Gomega1, newRing,currRing); 4745 4746 to = clock(); 4747 // compute a representation of the generators of submod (M) with respect to those of mod (Gomega), where Gomega is a reduced Groebner basis w.r.t. the current ring 4748 F = MLifttwoIdeal(Gomega2, M1, G); 4749 #ifdef CHECK_IDEAL_MWALK 4750 idString(F, "F"); 4751 #endif 4752 tlift = tlift + clock() - to; 4753 4754 idDelete(&M1); 4755 idDelete(&Gomega2); 4756 idDelete(&G); 4757 4758 rChangeCurrRing(newRing); // change the ring to newRing 4759 F1 = idrMoveR(F,oldRing,currRing); 4760 } 4761 4762 to = clock(); 4763 G = kInterRedCC(F1, NULL); //reduce the Groebner basis <G> w.r.t. new ring 4764 #ifdef CHECK_IDEAL_MWALK 4765 idString(G, "G"); 4766 #endif 4767 idDelete(&F1); 4768 if(endwalks == 1) 4769 { 4770 xtred = xtred + clock() - to; 4771 break; 4772 } 4773 else 4774 { 4775 tred = tred + clock() - to; 4776 } 4777 4778 NEXT_VECTOR: 4779 to = clock(); 4780 //intvec* next_weight = MkInterRedNextWeight(curr_weight,target_weight,G); 4781 intvec* next_weight = MWalkRandomNextWeight(G, curr_weight, target_weight, weight_rad, pert_deg); 4782 4783 tnw = tnw + clock() - to; 4784 //#ifdef PRINT_VECTORS 4785 MivString(curr_weight, target_weight, next_weight); 4786 //#endif 4787 if(Overflow_Error == TRUE) 4788 { 4789 newRing = currRing; 4790 PrintS("\n//Mrwalk: The computed vector does NOT stay in Cone!!\n"); 4791 4792 if (rParameter(currRing) != NULL) 4793 { 4794 DefRingPar(target_weight); 4795 } 4796 else 4797 { 4798 VMrDefault(target_weight); //define a new ring with ordering "(a(curr_weight),lp) 4799 } 4800 F1 = idrMoveR(G, newRing,currRing); 4801 G = MstdCC(F1); 4802 idDelete(&F1); 4803 4804 newRing = currRing; 4805 break; 4806 } 4807 4808 if(MivComp(next_weight, ivNull) == 1) 4809 { 4810 newRing = currRing; 4811 delete next_weight; 4812 break; 4813 } 4814 if(MivComp(next_weight, target_weight) == 1) 4815 { 4816 endwalks = 1; 4817 } 4818 for(i=nV-1; i>=0; i--) 4819 { 4820 (*tmp_weight)[i] = (*curr_weight)[i]; 4821 (*curr_weight)[i] = (*next_weight)[i]; 4822 } 4823 delete next_weight; 4824 } 4825 rChangeCurrRing(XXRing); 4826 G = idrMoveR(G, newRing,currRing); 4827 4828 delete tmp_weight; 4829 delete ivNull; 4830 delete exivlp; 4831 #ifndef BUCHBERGER_ALG 4832 delete last_omega; 4833 #endif 4834 #ifdef TIME_TEST 4835 TimeString(tinput, tostd, tif, tstd, tlift, tred, tnw, nstep); 4836 4837 Print("\n// pSetm_Error = (%d)", ErrorCheck()); 4838 Print("\n// Overflow_Error? (%d)\n", Overflow_Error); 4839 #endif 4840 return(G); 4841 } 4842 4843 //unused 4844 #if 0 3822 4845 ideal Mwalk_tst(ideal Go, intvec* curr_weight, intvec* target_weight) 3823 4846 { 3824 clock_t tinput=clock();4847 //clock_t tinput=clock(); 3825 4848 //idString(Go,"Ginp"); 3826 4849 int i, nV = currRing->N; 3827 4850 int nwalk=0, endwalks=0; 3828 4851 3829 ideal Gomega, M, F, Gomega1, Gomega2, M1, F1, G, G1; 3830 ring endRing, newRing, oldRing; 4852 ideal Gomega, M, F, Gomega1, Gomega2, M1, F1, G; 4853 // ideal G1; ring endRing; 4854 ring newRing, oldRing; 3831 4855 intvec* ivNull = new intvec(nV); 3832 4856 ring XXRing = currRing; … … 3834 4858 intvec* tmp_weight = new intvec(nV); 3835 4859 for(i=nV-1; i>=0; i--) 4860 { 3836 4861 (*tmp_weight)[i] = (*curr_weight)[i]; 3837 4862 } 3838 4863 /* the monomial ordering of this current ring would be "dp" */ 3839 4864 G = MstdCC(Go); 3840 4865 #ifndef BUCHBERGER_ALG 3841 4866 intvec* hilb_func; 3842 4867 #endif 3843 4868 /* to avoid (1,0,...,0) as the target vector */ 3844 4869 intvec* last_omega = new intvec(nV); … … 3940 4965 return(G); 3941 4966 } 3942 3943 4967 #endif 3944 4968 3945 4969 /**************************************************************/ … … 3955 4979 2) the changed perturbation walk algorithm with a decreased degree. 3956 4980 */ 3957 / * use kStd, if nP = 0, else call LastGB */4981 // use kStd, if nP = 0, else call LastGB 3958 4982 ideal Mpwalk(ideal Go, int op_deg, int tp_deg,intvec* curr_weight, 3959 4983 intvec* target_weight, int nP) … … 3971 4995 3972 4996 nstep = 0; 3973 int i, ntwC=1, ntestw=1, nV = currRing->N , op_tmp = op_deg;3974 int endwalks=0 , nhilb=0, ntestomega=0;4997 int i, ntwC=1, ntestw=1, nV = currRing->N; 4998 int endwalks=0; 3975 4999 3976 5000 ideal Gomega, M, F, G, Gomega1, Gomega2, M1,F1,Eresult,ssG; … … 3983 5007 intvec* ivNull = new intvec(nV); 3984 5008 intvec* iv_dp = MivUnit(nV);// define (1,1,...,1) 3985 intvec* cw_tmp = curr_weight; 5009 #ifndef BUCHBERGER_ALG 3986 5010 intvec* hilb_func; 5011 #endif 3987 5012 intvec* next_weight; 3988 intvec* extra_curr_weight = new intvec(nV); 3989 3990 /* to avoid (1,0,...,0) as the target vector */ 5013 5014 // to avoid (1,0,...,0) as the target vector 3991 5015 intvec* last_omega = new intvec(nV); 3992 5016 for(i=nV-1; i>0; i--) … … 4011 5035 else 4012 5036 { 4013 // ring order := (a(curr_weight),lp);5037 //define ring order := (a(curr_weight),lp); 4014 5038 if (rParameter(currRing) != NULL) 4015 5039 DefRingPar(curr_weight); … … 4054 5078 } 4055 5079 delete iv_M_lp; 4056 pert_target_vector = target_weight; //vor 19. mai 2003//test 19 Junu 035080 pert_target_vector = target_weight; 4057 5081 rChangeCurrRing(HelpRing); 4058 5082 G = idrMoveR(ssG, TargetRing,currRing); … … 4093 5117 oldRing = currRing; 4094 5118 4095 / * define a new ring that its ordering is "(a(curr_weight),lp) */5119 // define a new ring with ordering "(a(curr_weight),lp) 4096 5120 if (rParameter(currRing) != NULL) 4097 5121 DefRingPar(curr_weight); … … 4187 5211 { 4188 5212 ntwC = 0; 4189 ntestomega = 1;5213 //ntestomega = 1; 4190 5214 //Print("\n// ring r%d = %s;\n", nstep, rString(currRing)); 4191 5215 //idElements(G, "G"); … … 4195 5219 if(MivComp(next_weight, ivNull) == 1){ 4196 5220 newRing = currRing; 4197 delete next_weight; //16.03.025221 delete next_weight; 4198 5222 //Print("\n// ring r%d = %s;\n", nstep, rString(currRing)); 4199 5223 break; … … 4210 5234 if(tp_deg != 1) 4211 5235 { 4212 FINISH_160302: //16.03.025236 FINISH_160302: 4213 5237 if(MivSame(orig_target, exivlp) == 1) 4214 5238 if (rParameter(currRing) != NULL) … … 4244 5268 } 4245 5269 */ 4246 / * LastGB is "better" than the kStd subroutine */5270 // LastGB is "better" than the kStd subroutine 4247 5271 to=clock(); 4248 5272 ideal eF1; … … 4296 5320 intvec* XivNull; 4297 5321 4298 / * define a matrix (1 ... 1) */5322 // define a matrix (1 ... 1) 4299 5323 intvec* MMatrixone(int nV) 4300 5324 { … … 4313 5337 int Xngleich; 4314 5338 4315 /* 27.07.03 */ 4316 /* Perturb the start weight vector at the top level, i.e. nlev = 1 */ 5339 /*********************************************************************** 5340 * Perturb the start weight vector at the top level, i.e. nlev = 1 * 5341 ***********************************************************************/ 4317 5342 static ideal rec_fractal_call(ideal G, int nlev, intvec* omtmp) 4318 5343 { … … 4321 5346 4322 5347 int i, nV = currRing->N; 4323 ring new_ring, testring, extoRing; 5348 ring new_ring, testring; 5349 //ring extoRing; 4324 5350 ideal Gomega, Gomega1, Gomega2, F, F1, Gresult, Gresult1, G1, Gt; 4325 5351 int nwalks = 0; 4326 5352 intvec* Mwlp; 5353 #ifndef BUCHBERGER_ALG 4327 5354 intvec* hilb_func; 4328 intvec* extXtau; 5355 #endif 5356 // intvec* extXtau; 4329 5357 intvec* next_vect; 4330 5358 intvec* omega2 = new intvec(nV); 4331 5359 intvec* altomega = new intvec(nV); 4332 5360 4333 BOOLEAN isnewtarget = FALSE;4334 4335 / * to avoid (1,0,...,0) as the target vector (Hans) */5361 //BOOLEAN isnewtarget = FALSE; 5362 5363 // to avoid (1,0,...,0) as the target vector (Hans) 4336 5364 intvec* last_omega = new intvec(nV); 4337 5365 for(i=nV-1; i>0; i--) … … 4438 5466 { 4439 5467 delete next_vect; 4440 4441 OVERFLOW_IN_NEXT_VECTOR:4442 5468 if (rParameter(currRing) != NULL) 5469 { 4443 5470 DefRingPar(omtmp); 5471 } 4444 5472 else 5473 { 4445 5474 VMrDefault(omtmp); 4446 5475 } 4447 5476 #ifdef TEST_OVERFLOW 4448 5477 Gt = idrMoveR(G, oRing,currRing); … … 4654 5683 } 4655 5684 5685 /************************************************************************ 5686 * Perturb the start weight vector at the top level with random element * 5687 ************************************************************************/ 5688 static ideal rec_r_fractal_call(ideal G, int nlev, intvec* omtmp, int weight_rad) 5689 { 5690 Overflow_Error = FALSE; 5691 //Print("\n\n// Entering the %d-th recursion:", nlev); 5692 5693 int i,k,weight_norm; 5694 int nV = currRing->N; 5695 ring new_ring, testring; 5696 //ring extoRing; 5697 ideal Gomega, Gomega1, Gomega2, F, F1, Gresult, Gresult1, G1, Gt; 5698 int nwalks = 0; 5699 intvec* Mwlp; 5700 #ifndef BUCHBERGER_ALG 5701 intvec* hilb_func; 5702 #endif 5703 //intvec* extXtau; 5704 intvec* next_vect; 5705 intvec* omega2 = new intvec(nV); 5706 intvec* altomega = new intvec(nV); 5707 5708 //BOOLEAN isnewtarget = FALSE; 5709 5710 // to avoid (1,0,...,0) as the target vector (Hans) 5711 intvec* last_omega = new intvec(nV); 5712 for(i=nV-1; i>0; i--) 5713 (*last_omega)[i] = 1; 5714 (*last_omega)[0] = 10000; 5715 5716 intvec* omega = new intvec(nV); 5717 for(i=0; i<nV; i++) { 5718 if(Xsigma->length() == nV) 5719 (*omega)[i] = (*Xsigma)[i]; 5720 else 5721 (*omega)[i] = (*Xsigma)[(nV*(nlev-1))+i]; 5722 5723 (*omega2)[i] = (*Xtau)[(nlev-1)*nV+i]; 5724 } 5725 5726 if(nlev == 1) Xcall = 1; 5727 else Xcall = 0; 5728 5729 ring oRing = currRing; 5730 5731 while(1) 5732 { 5733 #ifdef FIRST_STEP_FRACTAL 5734 // perturb the current weight vector only on the top level or 5735 // after perturbation of the both vectors, nlev = 2 as the top level 5736 if((nlev == 1 && Xcall == 0) || (nlev == 2 && Xngleich == 1)) 5737 if(islengthpoly2(G) == 1) 5738 { 5739 Mwlp = MivWeightOrderlp(omega); 5740 Xsigma = Mfpertvector(G, Mwlp); 5741 delete Mwlp; 5742 Overflow_Error = FALSE; 5743 } 5744 #endif 5745 nwalks ++; 5746 NEXT_VECTOR_FRACTAL: 5747 to=clock(); 5748 /* determine the next border */ 5749 next_vect = MkInterRedNextWeight(omega,omega2,G); 5750 xtnw=xtnw+clock()-to; 5751 #ifdef PRINT_VECTORS 5752 MivString(omega, omega2, next_vect); 5753 #endif 5754 oRing = currRing; 5755 5756 /* We only perturb the current target vector at the recursion level 1 */ 5757 if(Xngleich == 0 && nlev == 1) //(ngleich == 0) important, e.g. ex2, ex3 5758 { 5759 if (MivComp(next_vect, omega2) == 1) 5760 { 5761 /* to dispense with taking initial (and lifting/interreducing 5762 after the call of recursion */ 5763 //Print("\n\n// ** Perturb the both vectors with degree %d with",nlev); 5764 //idElements(G, "G"); 5765 5766 Xngleich = 1; 5767 nlev +=1; 5768 5769 if (rParameter(currRing) != NULL) 5770 DefRingPar(omtmp); 5771 else 5772 VMrDefault(omtmp); 5773 5774 testring = currRing; 5775 Gt = idrMoveR(G, oRing,currRing); 5776 5777 /* perturb the original target vector w.r.t. the current GB */ 5778 delete Xtau; 5779 Xtau = NewVectorlp(Gt); 5780 5781 rChangeCurrRing(oRing); 5782 G = idrMoveR(Gt, testring,currRing); 5783 5784 /* perturb the current vector w.r.t. the current GB */ 5785 Mwlp = MivWeightOrderlp(omega); 5786 Xsigma = Mfpertvector(G, Mwlp); 5787 delete Mwlp; 5788 5789 for(i=nV-1; i>=0; i--) { 5790 (*omega2)[i] = (*Xtau)[nV+i]; 5791 (*omega)[i] = (*Xsigma)[nV+i]; 5792 } 5793 5794 delete next_vect; 5795 to=clock(); 5796 5797 /* to avoid the value of Overflow_Error that occur in Mfpertvector*/ 5798 Overflow_Error = FALSE; 5799 5800 next_vect = MkInterRedNextWeight(omega,omega2,G); 5801 xtnw=xtnw+clock()-to; 5802 5803 #ifdef PRINT_VECTORS 5804 MivString(omega, omega2, next_vect); 5805 #endif 5806 } 5807 else 5808 { 5809 // compute a perturbed next weight vector "next_vect1" 5810 intvec* next_vect11 = MPertVectors(G, MivMatrixOrder(omega), nlev); 5811 intvec* next_vect1 = MkInterRedNextWeight(next_vect11, omega2, G); 5812 // Print("\n// size of next_vect = %d", sizeof((*next_vect))); 5813 // Print("\n// size of next_vect1 = %d", sizeof((*next_vect1))); 5814 // Print("\n// size of next_vect11 = %d", sizeof((*next_vect11))); 5815 delete next_vect11; 5816 5817 // compare next_vect and next_vect1 5818 ideal G_test = MwalkInitialForm(G, next_vect); 5819 ideal G_test1 = MwalkInitialForm(G, next_vect1); 5820 // Print("\n// G_test, G_test 1 erzeugt"); 5821 if(IDELEMS(G_test1) <= IDELEMS(G_test)) 5822 { 5823 next_vect = ivCopy(next_vect1); 5824 } 5825 delete next_vect1; 5826 // compute a random next weight vector "next_vect2" 5827 intvec* next_vect22 = ivCopy(omega2); 5828 // Print("\n// size of next_vect22 = %d", sizeof((*next_vect22))); 5829 k = 0; 5830 while(test_w_in_ConeCC(G, next_vect22) == 0) 5831 { 5832 k = k + 1; 5833 if(k>10) 5834 { 5835 break; 5836 } 5837 weight_norm = 0; 5838 while(weight_norm == 0) 5839 { 5840 for(i=nV-1; i>=0; i--) 5841 { 5842 (*next_vect22)[i] = rand() % 60000 - 30000; 5843 weight_norm = weight_norm + (*next_vect22)[i]*(*next_vect22)[i]; 5844 } 5845 weight_norm = 1 + floor(sqrt(weight_norm)); 5846 } 5847 for(i=nV-1; i>=0; i--) 5848 { 5849 if((*next_vect22)[i] < 0) 5850 { 5851 (*next_vect22)[i] = 1 + (*omega)[i] + floor(weight_rad*(*next_vect22)[i]/weight_norm); 5852 } 5853 else 5854 { 5855 (*next_vect22)[i] = (*omega)[i] + floor(weight_rad*(*next_vect22)[i]/weight_norm); 5856 } 5857 } 5858 } 5859 if(test_w_in_ConeCC(G, next_vect22) == 1) 5860 { 5861 //compare next_weight and next_vect2 5862 intvec* next_vect2 = MkInterRedNextWeight(next_vect22, omega2, G); 5863 // Print("\n// size of next_vect2 = %d", sizeof((*next_vect2))); 5864 ideal G_test2 = MwalkInitialForm(G, next_vect2); 5865 if(IDELEMS(G_test2) <= IDELEMS(G_test)) 5866 { 5867 if(IDELEMS(G_test2) <= IDELEMS(G_test1)) 5868 { 5869 next_vect = ivCopy(next_vect2); 5870 } 5871 } 5872 idDelete(&G_test2); 5873 delete next_vect2; 5874 } 5875 delete next_vect22; 5876 idDelete(&G_test); 5877 idDelete(&G_test1); 5878 #ifdef PRINT_VECTORS 5879 MivString(omega, omega2, next_vect); 5880 #endif 5881 } 5882 } 5883 5884 5885 /* check whether the the computed vector is in the correct cone */ 5886 /* If no, the reduced GB of an omega-homogeneous ideal will be 5887 computed by Buchberger algorithm and stop this recursion step*/ 5888 //if(test_w_in_ConeCC(G, next_vect) != 1) //e.g. Example s7, cyc6 5889 if(Overflow_Error == TRUE) 5890 { 5891 delete next_vect; 5892 5893 //OVERFLOW_IN_NEXT_VECTOR: 5894 if (rParameter(currRing) != NULL) 5895 DefRingPar(omtmp); 5896 else 5897 VMrDefault(omtmp); 5898 5899 #ifdef TEST_OVERFLOW 5900 Gt = idrMoveR(G, oRing,currRing); 5901 Gt = NULL; return(Gt); 5902 #endif 5903 5904 //Print("\n\n// apply BB's alg. in ring r = %s;", rString(currRing)); 5905 to=clock(); 5906 Gt = idrMoveR(G, oRing,currRing); 5907 G1 = MstdCC(Gt); 5908 xtextra=xtextra+clock()-to; 5909 Gt = NULL; 5910 5911 delete omega2; 5912 delete altomega; 5913 5914 //Print("\n// Leaving the %d-th recursion with %d steps", nlev, nwalks); 5915 //Print(" ** Overflow_Error? (%d)", Overflow_Error); 5916 nnflow ++; 5917 5918 Overflow_Error = FALSE; 5919 return (G1); 5920 } 5921 5922 5923 /* If the perturbed target vector stays in the correct cone, 5924 return the current GB, 5925 otherwise, return the computed GB by the Buchberger-algorithm. 5926 Then we update the perturbed target vectors w.r.t. this GB. */ 5927 5928 /* the computed vector is equal to the origin vector, since 5929 t is not defined */ 5930 if (MivComp(next_vect, XivNull) == 1) 5931 { 5932 if (rParameter(currRing) != NULL) 5933 DefRingPar(omtmp); 5934 else 5935 VMrDefault(omtmp); 5936 5937 testring = currRing; 5938 Gt = idrMoveR(G, oRing,currRing); 5939 5940 if(test_w_in_ConeCC(Gt, omega2) == 1) { 5941 delete omega2; 5942 delete next_vect; 5943 delete altomega; 5944 //Print("\n// Leaving the %d-th recursion with %d steps ",nlev, nwalks); 5945 //Print(" ** Overflow_Error? (%d)", Overflow_Error); 5946 5947 return (Gt); 5948 } 5949 else 5950 { 5951 //ivString(omega2, "tau'"); 5952 //Print("\n// tau' doesn't stay in the correct cone!!"); 5953 5954 #ifndef MSTDCC_FRACTAL 5955 //07.08.03 5956 //ivString(Xtau, "old Xtau"); 5957 intvec* Xtautmp = Mfpertvector(Gt, MivMatrixOrder(omtmp)); 5958 #ifdef TEST_OVERFLOW 5959 if(Overflow_Error == TRUE) 5960 Gt = NULL; return(Gt); 5961 #endif 5962 5963 if(MivSame(Xtau, Xtautmp) == 1) 5964 { 5965 //PrintS("\n// Update vectors are equal to the old vectors!!"); 5966 delete Xtautmp; 5967 goto FRACTAL_MSTDCC; 5968 } 5969 5970 Xtau = Xtautmp; 5971 Xtautmp = NULL; 5972 //ivString(Xtau, "new Xtau"); 5973 5974 for(i=nV-1; i>=0; i--) 5975 (*omega2)[i] = (*Xtau)[(nlev-1)*nV+i]; 5976 5977 //Print("\n// ring tau = %s;", rString(currRing)); 5978 rChangeCurrRing(oRing); 5979 G = idrMoveR(Gt, testring,currRing); 5980 5981 goto NEXT_VECTOR_FRACTAL; 5982 #endif 5983 5984 FRACTAL_MSTDCC: 5985 //Print("\n// apply BB-Alg in ring = %s;", rString(currRing)); 5986 to=clock(); 5987 G = MstdCC(Gt); 5988 xtextra=xtextra+clock()-to; 5989 5990 oRing = currRing; 5991 5992 // update the original target vector w.r.t. the current GB 5993 if(MivSame(Xivinput, Xivlp) == 1) 5994 if (rParameter(currRing) != NULL) 5995 DefRingParlp(); 5996 else 5997 VMrDefaultlp(); 5998 else 5999 if (rParameter(currRing) != NULL) 6000 DefRingPar(Xivinput); 6001 else 6002 VMrDefault(Xivinput); 6003 6004 testring = currRing; 6005 Gt = idrMoveR(G, oRing,currRing); 6006 6007 delete Xtau; 6008 Xtau = NewVectorlp(Gt); 6009 6010 rChangeCurrRing(oRing); 6011 G = idrMoveR(Gt, testring,currRing); 6012 6013 delete omega2; 6014 delete next_vect; 6015 delete altomega; 6016 /* 6017 Print("\n// Leaving the %d-th recursion with %d steps,", nlev,nwalks); 6018 Print(" ** Overflow_Error? (%d)", Overflow_Error); 6019 */ 6020 if(Overflow_Error == TRUE) 6021 nnflow ++; 6022 6023 Overflow_Error = FALSE; 6024 return(G); 6025 } 6026 } 6027 6028 for(i=nV-1; i>=0; i--) { 6029 (*altomega)[i] = (*omega)[i]; 6030 (*omega)[i] = (*next_vect)[i]; 6031 } 6032 delete next_vect; 6033 6034 to=clock(); 6035 /* Take the initial form of <G> w.r.t. omega */ 6036 Gomega = MwalkInitialForm(G, omega); 6037 xtif=xtif+clock()-to; 6038 6039 #ifndef BUCHBERGER_ALG 6040 if(isNolVector(omega) == 0) 6041 hilb_func = hFirstSeries(Gomega,NULL,NULL,omega,currRing); 6042 else 6043 hilb_func = hFirstSeries(Gomega,NULL,NULL,last_omega,currRing); 6044 #endif // BUCHBERGER_ALG 6045 6046 if (rParameter(currRing) != NULL) 6047 DefRingPar(omega); 6048 else 6049 VMrDefault(omega); 6050 6051 Gomega1 = idrMoveR(Gomega, oRing,currRing); 6052 6053 /* Maximal recursion depth, to compute a red. GB */ 6054 /* Fractal walk with the alternative recursion */ 6055 /* alternative recursion */ 6056 // if(nlev == nV || lengthpoly(Gomega1) == 0) 6057 if(nlev == Xnlev || lengthpoly(Gomega1) == 0) 6058 //if(nlev == nV) // blind recursion 6059 { 6060 /* 6061 if(Xnlev != nV) 6062 { 6063 Print("\n// ** Xnlev = %d", Xnlev); 6064 ivString(Xtau, "Xtau"); 6065 } 6066 */ 6067 to=clock(); 6068 #ifdef BUCHBERGER_ALG 6069 Gresult = MstdhomCC(Gomega1); 6070 #else 6071 Gresult =kStd(Gomega1,NULL,isHomog,NULL,hilb_func,0,NULL,omega); 6072 delete hilb_func; 6073 #endif // BUCHBERGER_ALG 6074 xtstd=xtstd+clock()-to; 6075 } 6076 else { 6077 rChangeCurrRing(oRing); 6078 Gomega1 = idrMoveR(Gomega1, oRing,currRing); 6079 Gresult = rec_fractal_call(idCopy(Gomega1),nlev+1,omega); 6080 } 6081 6082 //convert a Groebner basis from a ring to another ring, 6083 new_ring = currRing; 6084 6085 rChangeCurrRing(oRing); 6086 Gresult1 = idrMoveR(Gresult, new_ring,currRing); 6087 Gomega2 = idrMoveR(Gomega1, new_ring,currRing); 6088 6089 to=clock(); 6090 /* Lifting process */ 6091 F = MLifttwoIdeal(Gomega2, Gresult1, G); 6092 xtlift=xtlift+clock()-to; 6093 idDelete(&Gresult1); 6094 idDelete(&Gomega2); 6095 idDelete(&G); 6096 6097 rChangeCurrRing(new_ring); 6098 F1 = idrMoveR(F, oRing,currRing); 6099 6100 to=clock(); 6101 /* Interreduce G */ 6102 G = kInterRedCC(F1, NULL); 6103 xtred=xtred+clock()-to; 6104 idDelete(&F1); 6105 } 6106 } 6107 6108 6109 6110 /******************************************************************************* 6111 * The implementation of the fractal walk algorithm * 6112 * * 6113 * The main procedur Mfwalk calls the recursive Subroutine * 6114 * rec_fractal_call to compute the wanted Grï¿œbner basis. * 6115 * At the main procedur we compute the reduced Grï¿œbner basis w.r.t. a "fast" * 6116 * order, e.g. "dp" and a sequence of weight vectors which are row vectors * 6117 * of a matrix. This matrix defines the given monomial order, e.g. "lp" * 6118 *******************************************************************************/ 4656 6119 ideal Mfwalk(ideal G, intvec* ivstart, intvec* ivtarget) 4657 6120 { … … 4684 6147 for(i=IDELEMS(Gw)-1; i>=0; i--) 4685 6148 { 4686 if((Gw->m[i]!=NULL) / * len >=0 */4687 && (Gw->m[i]->next!=NULL) / * len >=1 */4688 && (Gw->m[i]->next->next!=NULL)) / * len >=2 */4689 { 4690 intvec* iv_dp = MivUnit(nV); // define (1,1,...,1)6149 if((Gw->m[i]!=NULL) // len >=0 6150 && (Gw->m[i]->next!=NULL) // len >=1 6151 && (Gw->m[i]->next->next!=NULL)) // len >=2 6152 { 6153 intvec* iv_dp = MivUnit(nV); // define (1,1,...,1) 4691 6154 intvec* Mdp; 4692 6155 … … 4780 6243 } 4781 6244 4782 /* Tran algorithm */ 4783 /* use kStd, if nP = 0, else call Ab_Rec_Pert (LastGB) */ 6245 ideal Mfrwalk(ideal G, intvec* ivstart, intvec* ivtarget,int weight_rad) 6246 { 6247 Set_Error(FALSE); 6248 Overflow_Error = FALSE; 6249 //Print("// pSetm_Error = (%d)", ErrorCheck()); 6250 //Print("\n// ring ro = %s;", rString(currRing)); 6251 6252 nnflow = 0; 6253 Xngleich = 0; 6254 Xcall = 0; 6255 xtif=0; xtstd=0; xtlift=0; xtred=0; xtnw=0; xtextra=0; 6256 xftinput = clock(); 6257 6258 ring oldRing = currRing; 6259 int i, nV = currRing->N; 6260 XivNull = new intvec(nV); 6261 Xivinput = ivtarget; 6262 ngleich = 0; 6263 to=clock(); 6264 ideal I = MstdCC(G); 6265 G = NULL; 6266 xftostd=clock()-to; 6267 Xsigma = ivstart; 6268 6269 Xnlev=nV; 6270 6271 #ifdef FIRST_STEP_FRACTAL 6272 ideal Gw = MwalkInitialForm(I, ivstart); 6273 for(i=IDELEMS(Gw)-1; i>=0; i--) 6274 { 6275 if((Gw->m[i]!=NULL) // len >=0 6276 && (Gw->m[i]->next!=NULL) // len >=1 6277 && (Gw->m[i]->next->next!=NULL)) // len >=2 6278 { 6279 intvec* iv_dp = MivUnit(nV); // define (1,1,...,1) 6280 intvec* Mdp; 6281 6282 if(MivSame(ivstart, iv_dp) != 1) 6283 { 6284 Mdp = MivWeightOrderdp(ivstart); 6285 } 6286 else 6287 { 6288 Mdp = MivMatrixOrderdp(nV); 6289 } 6290 Xsigma = Mfpertvector(I, Mdp); 6291 Overflow_Error = FALSE; 6292 6293 delete Mdp; 6294 delete iv_dp; 6295 break; 6296 } 6297 } 6298 idDelete(&Gw); 6299 #endif 6300 6301 ideal I1; 6302 intvec* Mlp; 6303 Xivlp = Mivlp(nV); 6304 6305 if(MivComp(ivtarget, Xivlp) != 1) 6306 { 6307 if (rParameter(currRing) != NULL) 6308 DefRingPar(ivtarget); 6309 else 6310 VMrDefault(ivtarget); 6311 6312 I1 = idrMoveR(I, oldRing,currRing); 6313 Mlp = MivWeightOrderlp(ivtarget); 6314 Xtau = Mfpertvector(I1, Mlp); 6315 } 6316 else 6317 { 6318 if (rParameter(currRing) != NULL) 6319 DefRingParlp(); 6320 else 6321 VMrDefaultlp(); 6322 6323 I1 = idrMoveR(I, oldRing,currRing); 6324 Mlp = MivMatrixOrderlp(nV); 6325 Xtau = Mfpertvector(I1, Mlp); 6326 } 6327 delete Mlp; 6328 Overflow_Error = FALSE; 6329 6330 //ivString(Xsigma, "Xsigma"); 6331 //ivString(Xtau, "Xtau"); 6332 6333 id_Delete(&I, oldRing); 6334 ring tRing = currRing; 6335 6336 if (rParameter(currRing) != NULL) 6337 DefRingPar(ivstart); 6338 else 6339 VMrDefault(ivstart); 6340 6341 I = idrMoveR(I1,tRing,currRing); 6342 to=clock(); 6343 ideal J = MstdCC(I); 6344 idDelete(&I); 6345 xftostd=xftostd+clock()-to; 6346 6347 ideal resF; 6348 ring helpRing = currRing; 6349 J = rec_r_fractal_call(J,1,ivtarget,weight_rad); 6350 rChangeCurrRing(oldRing); 6351 resF = idrMoveR(J, helpRing,currRing); 6352 idSkipZeroes(resF); 6353 6354 delete Xivlp; 6355 delete Xsigma; 6356 delete Xtau; 6357 delete XivNull; 6358 6359 #ifdef TIME_TEST 6360 TimeStringFractal(xftinput, xftostd, xtif, xtstd, xtextra, 6361 xtlift, xtred, xtnw); 6362 6363 6364 Print("\n// pSetm_Error = (%d)", ErrorCheck()); 6365 Print("\n// Overflow_Error? (%d)\n", Overflow_Error); 6366 Print("\n// the numbers of Overflow_Error (%d)", nnflow); 6367 #endif 6368 6369 return(resF); 6370 } 6371 6372 /******************************************************* 6373 * Tran's algorithm * 6374 * * 6375 * use kStd, if nP = 0, else call Ab_Rec_Pert (LastGB) * 6376 *******************************************************/ 4784 6377 ideal TranMImprovwalk(ideal G,intvec* curr_weight,intvec* target_tmp, int nP) 4785 6378 { 6379 #ifdef TIME_TEST 4786 6380 clock_t mtim = clock(); 6381 #endif 4787 6382 Set_Error(FALSE ); 4788 6383 Overflow_Error = FALSE; … … 4791 6386 4792 6387 clock_t tostd, tif=0, tstd=0, tlift=0, tred=0, tnw=0, textra=0; 6388 #ifdef TIME_TEST 4793 6389 clock_t tinput = clock(); 4794 6390 #endif 4795 6391 int nsteppert=0, i, nV = currRing->N, nwalk=0, npert_tmp=0; 4796 6392 int *npert=(int*)omAlloc(2*nV*sizeof(int)); 4797 6393 ideal Gomega, M,F, G1, Gomega1, Gomega2, M1, F1; 4798 ring endRing, newRing, oldRing, lpRing; 6394 //ring endRing; 6395 ring newRing, oldRing, lpRing; 4799 6396 intvec* next_weight; 4800 6397 intvec* ivNull = new intvec(nV); //define (0,...,0) 4801 6398 intvec* iv_dp = MivUnit(nV);// define (1,1,...,1) 4802 6399 intvec* iv_lp = Mivlp(nV); //define (1,0,...,0) 4803 ideal H0, H1, H2, Glp; 6400 ideal H0; 6401 //ideal H1; 6402 ideal H2, Glp; 4804 6403 int nGB, endwalks = 0, nwalkpert=0, npertstep=0; 4805 6404 intvec* Mlp = MivMatrixOrderlp(nV); 4806 6405 intvec* vector_tmp = new intvec(nV); 6406 #ifndef BUCHBERGER_ALG 4807 6407 intvec* hilb_func; 4808 6408 #endif 4809 6409 /* to avoid (1,0,...,0) as the target vector */ 4810 6410 intvec* last_omega = new intvec(nV); … … 4939 6539 nwalkpert++; 4940 6540 to=clock(); 4941 / * compute a next weight vector */6541 // compute a next weight vector 4942 6542 next_weight = MwalkNextWeightCC(curr_weight,target_weight, G); 4943 6543 tnw=tnw+clock()-to; … … 4947 6547 4948 6548 4949 /* check whether the computed intermediate weight vector i n4950 the correct cone is, sincesometimes it is very big e.g. s7, cyc7.4951 If it is NOT in the co ne, then one has directly to compute6549 /* check whether the computed intermediate weight vector is in 6550 the correct cone; sometimes it is very big e.g. s7, cyc7. 6551 If it is NOT in the correct cone, then compute directly 4952 6552 a reduced Groebner basis with respect to the lexicographic ordering 4953 6553 for the known Groebner basis that it is computed in the last step. … … 5005 6605 } 5006 6606 5007 /* check whether the computed Groebner basis a really Groebner basisis.5008 if no, we perturb the target vector with the maximal "perturbation"6607 /* check whether the computed Groebner basis is really a Groebner basis. 6608 If not, we perturb the target vector with the maximal "perturbation" 5009 6609 degree.*/ 5010 6610 if(MivComp(next_weight, target_weight) == 1 || … … 5144 6744 delete next_weight; 5145 6745 }//while 5146 6746 #ifdef TEST_OVERFLOW 5147 6747 BE_FINISH: 6748 #endif 5148 6749 rChangeCurrRing(XXRing); 5149 6750 G = idrMoveR(G, lpRing,currRing); … … 5168 6769 } 5169 6770 5170 5171 /* compute the reduced Grï¿œbner basis of an ideal <Go> w.r.t. lp */ 6771 /******************************************************* 6772 * Tran's algorithm with random element * 6773 * * 6774 * use kStd, if nP = 0, else call Ab_Rec_Pert (LastGB) * 6775 *******************************************************/ 6776 ideal TranMrImprovwalk(ideal G,intvec* curr_weight,intvec* target_tmp, int nP, int weight_rad, int pert_deg) 6777 { 6778 #ifdef TIME_TEST 6779 clock_t mtim = clock(); 6780 #endif 6781 Set_Error(FALSE ); 6782 Overflow_Error = FALSE; 6783 //Print("// pSetm_Error = (%d)", ErrorCheck()); 6784 //Print("\n// ring ro = %s;", rString(currRing)); 6785 6786 clock_t tostd, tif=0, tstd=0, tlift=0, tred=0, tnw=0, textra=0; 6787 #ifdef TIME_TEST 6788 clock_t tinput = clock(); 6789 #endif 6790 int nsteppert=0, i, k, nV = currRing->N, nwalk=0, npert_tmp=0; 6791 int *npert=(int*)omAlloc(2*nV*sizeof(int)); 6792 ideal Gomega, M,F, G1, Gomega1, Gomega2, M1, F1; 6793 //ring endRing; 6794 ring newRing, oldRing, lpRing; 6795 intvec* next_weight; 6796 intvec* ivNull = new intvec(nV); //define (0,...,0) 6797 intvec* iv_dp = MivUnit(nV);// define (1,1,...,1) 6798 intvec* iv_lp = Mivlp(nV); //define (1,0,...,0) 6799 ideal H0; 6800 //ideal H1; 6801 ideal H2, Glp; 6802 int weight_norm, nGB, endwalks = 0, nwalkpert=0, npertstep=0; 6803 intvec* Mlp = MivMatrixOrderlp(nV); 6804 intvec* vector_tmp = new intvec(nV); 6805 #ifndef BUCHBERGER_ALG 6806 intvec* hilb_func; 6807 #endif 6808 // to avoid (1,0,...,0) as the target vector 6809 intvec* last_omega = new intvec(nV); 6810 for(i=nV-1; i>0; i--) 6811 { 6812 (*last_omega)[i] = 1; 6813 } 6814 (*last_omega)[0] = 10000; 6815 6816 //intvec* extra_curr_weight = new intvec(nV); 6817 intvec* target_weight = new intvec(nV); 6818 for(i=nV-1; i>=0; i--) 6819 { 6820 (*target_weight)[i] = (*target_tmp)[i]; 6821 } 6822 ring XXRing = currRing; 6823 newRing = currRing; 6824 6825 to=clock(); 6826 // compute a red. GB w.r.t. the help ring 6827 if(MivComp(curr_weight, iv_dp) == 1) 6828 { 6829 //rOrdStr(currRing) = "dp" 6830 G = MstdCC(G); 6831 } 6832 else 6833 { 6834 //rOrdStr(currRing) = (a(.c_w..),lp,C) 6835 if (rParameter(currRing) != NULL) 6836 { 6837 DefRingPar(curr_weight); 6838 } 6839 else 6840 { 6841 VMrDefault(curr_weight); 6842 } 6843 G = idrMoveR(G, XXRing,currRing); 6844 G = MstdCC(G); 6845 } 6846 tostd=clock()-to; 6847 6848 #ifdef REPRESENTATION_OF_SIGMA 6849 ideal Gw = MwalkInitialForm(G, curr_weight); 6850 6851 if(islengthpoly2(Gw)==1) 6852 { 6853 intvec* MDp; 6854 if(MivComp(curr_weight, iv_dp) == 1) 6855 { 6856 MDp = MatrixOrderdp(nV); //MivWeightOrderlp(iv_dp); 6857 } 6858 else 6859 { 6860 MDp = MivWeightOrderlp(curr_weight); 6861 } 6862 curr_weight = RepresentationMatrix_Dp(G, MDp); 6863 6864 delete MDp; 6865 6866 ring exring = currRing; 6867 6868 if (rParameter(currRing) != NULL) 6869 { 6870 DefRingPar(curr_weight); 6871 } 6872 else 6873 { 6874 VMrDefault(curr_weight); 6875 } 6876 to=clock(); 6877 Gw = idrMoveR(G, exring,currRing); 6878 G = MstdCC(Gw); 6879 Gw = NULL; 6880 tostd=tostd+clock()-to; 6881 //ivString(curr_weight,"rep. sigma"); 6882 goto COMPUTE_NEW_VECTOR; 6883 } 6884 6885 idDelete(&Gw); 6886 delete iv_dp; 6887 #endif 6888 6889 6890 while(1) 6891 { 6892 to=clock(); 6893 // compute an initial form ideal of <G> w.r.t. "curr_vector" 6894 Gomega = MwalkInitialForm(G, curr_weight); 6895 tif=tif+clock()-to; 6896 6897 #ifndef BUCHBERGER_ALG 6898 if(isNolVector(curr_weight) == 0) 6899 { 6900 hilb_func = hFirstSeries(Gomega,NULL,NULL,curr_weight,currRing); 6901 } 6902 else 6903 { 6904 hilb_func = hFirstSeries(Gomega,NULL,NULL,last_omega,currRing); 6905 } 6906 #endif // BUCHBERGER_ALG 6907 6908 oldRing = currRing; 6909 6910 // define a new ring with ordering "(a(curr_weight),lp) 6911 if (rParameter(currRing) != NULL) 6912 { 6913 DefRingPar(curr_weight); 6914 } 6915 else 6916 { 6917 VMrDefault(curr_weight); 6918 } 6919 newRing = currRing; 6920 Gomega1 = idrMoveR(Gomega, oldRing,currRing); 6921 6922 to=clock(); 6923 // compute a reduced Groebner basis of <Gomega> w.r.t. "newRing" 6924 #ifdef BUCHBERGER_ALG 6925 M = MstdhomCC(Gomega1); 6926 #else 6927 M=kStd(Gomega1,NULL,isHomog,NULL,hilb_func,0,NULL,curr_weight); 6928 delete hilb_func; 6929 #endif 6930 tstd=tstd+clock()-to; 6931 6932 // change the ring to oldRing 6933 rChangeCurrRing(oldRing); 6934 M1 = idrMoveR(M, newRing,currRing); 6935 Gomega2 = idrMoveR(Gomega1, newRing,currRing); 6936 6937 to=clock(); 6938 // compute a representation of the generators of submod (M) with respect to those of mod (Gomega). 6939 // Gomega is a reduced Groebner basis w.r.t. the current ring 6940 F = MLifttwoIdeal(Gomega2, M1, G); 6941 tlift=tlift+clock()-to; 6942 6943 idDelete(&M1); 6944 idDelete(&Gomega2); 6945 idDelete(&G); 6946 6947 // change the ring to newRing 6948 rChangeCurrRing(newRing); 6949 F1 = idrMoveR(F, oldRing,currRing); 6950 6951 to=clock(); 6952 // reduce the Groebner basis <G> w.r.t. new ring 6953 G = kInterRedCC(F1, NULL); 6954 tred=tred+clock()-to; 6955 idDelete(&F1); 6956 6957 COMPUTE_NEW_VECTOR: 6958 newRing = currRing; 6959 nwalk++; 6960 nwalkpert++; 6961 to=clock(); 6962 // compute a next weight vector 6963 //next_weight = MwalkNextWeightCC(curr_weight,target_weight, G); 6964 next_weight = MWalkRandomNextWeight(G, curr_weight, target_weight, weight_rad, pert_deg); 6965 /* 6966 next_weight = MkInterRedNextWeight(curr_weight,target_weight,G); 6967 6968 if(MivComp(next_weight, target_weight) != 1) 6969 { 6970 // compute a perturbed next weight vector "next_weight1" 6971 intvec* next_weight1 = MkInterRedNextWeight(MPertVectors(G, MivMatrixOrder(curr_weight), pert_deg), target_weight, G); 6972 6973 // compare next_weight and next_weight1 6974 ideal G_test = MwalkInitialForm(G, next_weight); 6975 ideal G_test1 = MwalkInitialForm(G, next_weight1); 6976 if(IDELEMS(G_test1) <= IDELEMS(G_test)) 6977 { 6978 next_weight = ivCopy(next_weight1); 6979 } 6980 delete next_weight1; 6981 // compute a random next weight vector "next_weight2" 6982 intvec* next_weight22 = ivCopy(target_weight); 6983 // Print("\n// size of target_weight = %d", sizeof((*target_weight))); 6984 k = 0; 6985 6986 while(test_w_in_ConeCC(G, next_weight22) == 0 && k < 11) 6987 { 6988 k++; 6989 if(k>10) 6990 { 6991 break; 6992 } 6993 weight_norm = 0; 6994 while(weight_norm == 0) 6995 { 6996 for(i=nV-1; i>=0; i--) 6997 { 6998 // Print("\n// next_weight[%d] = %d", i, (*next_weight)[i]); 6999 (*next_weight22)[i] = rand() % 60000 - 30000; 7000 weight_norm = weight_norm + (*next_weight22)[i]*(*next_weight22)[i]; 7001 } 7002 weight_norm = 1 + floor(sqrt(weight_norm)); 7003 } 7004 for(i=nV-1; i>=0; i--) 7005 { 7006 if((*next_weight22)[i] < 0) 7007 { 7008 (*next_weight22)[i] = 1 + (*curr_weight)[i] + floor(weight_rad*(*next_weight22)[i]/weight_norm); 7009 } 7010 else 7011 { 7012 (*next_weight22)[i] = (*curr_weight)[i] + floor(weight_rad*(*next_weight22)[i]/weight_norm); 7013 } 7014 // Print("\n// next_weight22[%d] = %d", i, (*next_weight22)[i]); 7015 } 7016 } 7017 7018 if(test_w_in_ConeCC(G, next_weight22) == 1) 7019 { 7020 // compare next_weight and next_weight2 7021 // Print("\n// ZUFALL IM KEGEL"); 7022 intvec* next_weight2 = MkInterRedNextWeight(next_weight22, target_weight, G); 7023 7024 ideal G_test2 = MwalkInitialForm(G, next_weight2); 7025 if(IDELEMS(G_test2) <= IDELEMS(G_test)) 7026 { 7027 if(IDELEMS(G_test2) <= IDELEMS(G_test1)) 7028 { 7029 // Print("\n// ZUFALL BENUTZT!\n"); 7030 next_weight = ivCopy(next_weight2); 7031 } 7032 } 7033 idDelete(&G_test2); 7034 delete next_weight2; 7035 } 7036 delete next_weight22; 7037 idDelete(&G_test); 7038 idDelete(&G_test1); 7039 }*/ 7040 7041 tnw=tnw+clock()-to; 7042 #ifdef PRINT_VECTORS 7043 MivString(curr_weight, target_weight, next_weight); 7044 #endif 7045 7046 /* check whether the computed intermediate weight vector is in 7047 the correct cone; sometimes it is very big e.g. s7, cyc7. 7048 If it is NOT in the correct cone, then compute directly 7049 a reduced Groebner basis with respect to the lexicographic ordering 7050 for the known Groebner basis that it is computed in the last step. 7051 */ 7052 //if(test_w_in_ConeCC(G, next_weight) != 1) 7053 if(Overflow_Error == TRUE) 7054 { 7055 OMEGA_OVERFLOW_TRAN_NEW: 7056 //Print("\n// takes %d steps!", nwalk-1); 7057 //Print("\n//ring lastRing = %s;", rString(currRing)); 7058 #ifdef TEST_OVERFLOW 7059 goto BE_FINISH; 7060 #endif 7061 7062 #ifdef CHECK_IDEAL_MWALK 7063 idElements(G, "G"); 7064 //headidString(G, "G"); 7065 #endif 7066 7067 if(MivSame(target_tmp, iv_lp) == 1) 7068 { 7069 if (rParameter(currRing) != NULL) 7070 { 7071 DefRingParlp(); 7072 } 7073 else 7074 { 7075 VMrDefaultlp(); 7076 } 7077 } 7078 else 7079 { 7080 if (rParameter(currRing) != NULL) 7081 { 7082 DefRingPar(target_tmp); 7083 } 7084 else 7085 { 7086 VMrDefault(target_tmp); 7087 } 7088 } 7089 lpRing = currRing; 7090 G1 = idrMoveR(G, newRing,currRing); 7091 7092 to=clock(); 7093 // apply kStd or LastGB to compute a lex. red. Groebner basis of <G> 7094 if(nP == 0 || MivSame(target_tmp, iv_lp) == 0) 7095 { 7096 //Print("\n\n// calls \"std in ring r_%d = %s;", nwalk, rString(currRing)); 7097 G = MstdCC(G1);//no result for qnt1 7098 } 7099 else 7100 { 7101 rChangeCurrRing(newRing); 7102 G1 = idrMoveR(G1, lpRing,currRing); 7103 7104 //Print("\n\n// calls \"LastGB\" (%d) to compute a GB", nV-1); 7105 G = LastGB(G1, curr_weight, nV-1); //no result for kats7 7106 7107 rChangeCurrRing(lpRing); 7108 G = idrMoveR(G, newRing,currRing); 7109 } 7110 textra=clock()-to; 7111 npert[endwalks]=nwalk-npert_tmp; 7112 npert_tmp = nwalk; 7113 endwalks ++; 7114 break; 7115 } 7116 7117 // check whether the computed Groebner basis is really a Groebner basis. 7118 // If not, we perturb the target vector with the maximal "perturbation" degree. 7119 7120 if(MivComp(next_weight, target_weight) == 1 || MivComp(next_weight, curr_weight) == 1 ) 7121 { 7122 //Print("\n//ring r_%d = %s;", nwalk, rString(currRing)); 7123 7124 7125 //compute the number of perturbations and its step 7126 npert[endwalks]=nwalk-npert_tmp; 7127 npert_tmp = nwalk; 7128 7129 endwalks ++; 7130 7131 // it is very important if the walk only uses one step, e.g. Fate, liu 7132 if(endwalks == 1 && MivComp(next_weight, curr_weight) == 1) 7133 { 7134 rChangeCurrRing(XXRing); 7135 G = idrMoveR(G, newRing,currRing); 7136 goto FINISH; 7137 } 7138 H0 = id_Head(G,currRing); 7139 7140 if(MivSame(target_tmp, iv_lp) == 1) 7141 { 7142 if (rParameter(currRing) != NULL) 7143 { 7144 DefRingParlp(); 7145 } 7146 else 7147 { 7148 VMrDefaultlp(); 7149 } 7150 } 7151 else 7152 { 7153 if (rParameter(currRing) != NULL) 7154 { 7155 DefRingPar(target_tmp); 7156 } 7157 else 7158 { 7159 VMrDefault(target_tmp); 7160 } 7161 } 7162 lpRing = currRing; 7163 Glp = idrMoveR(G, newRing,currRing); 7164 H2 = idrMoveR(H0, newRing,currRing); 7165 7166 // Apply Lemma 2.2 in Collart et. al (1997) to check whether cone(k-1) is equal to cone(k) 7167 nGB = 1; 7168 for(i=IDELEMS(Glp)-1; i>=0; i--) 7169 { 7170 poly t; 7171 if((t=pSub(pHead(Glp->m[i]), pCopy(H2->m[i]))) != NULL) 7172 { 7173 pDelete(&t); 7174 idDelete(&H2);//5.5.02 7175 nGB = 0; //i.e. Glp is no reduced Groebner basis 7176 break; 7177 } 7178 pDelete(&t); 7179 } 7180 7181 idDelete(&H2);//5.5.02 7182 7183 if(nGB == 1) 7184 { 7185 G = Glp; 7186 Glp = NULL; 7187 break; 7188 } 7189 7190 // perturb the target weight vector, if the vector target_tmp stays in many cones 7191 poly p; 7192 BOOLEAN plength3 = FALSE; 7193 for(i=IDELEMS(Glp)-1; i>=0; i--) 7194 { 7195 p = MpolyInitialForm(Glp->m[i], target_tmp); 7196 if(p->next != NULL && 7197 p->next->next != NULL && 7198 p->next->next->next != NULL) 7199 { 7200 Overflow_Error = FALSE; 7201 7202 for(i=0; i<nV; i++) 7203 { 7204 (*vector_tmp)[i] = (*target_weight)[i]; 7205 } 7206 delete target_weight; 7207 target_weight = MPertVectors(Glp, Mlp, nV); 7208 7209 if(MivComp(vector_tmp, target_weight)==1) 7210 { 7211 //PrintS("\n// The old and new representaion vector are the same!!"); 7212 G = Glp; 7213 newRing = currRing; 7214 goto OMEGA_OVERFLOW_TRAN_NEW; 7215 } 7216 7217 if(Overflow_Error == TRUE) 7218 { 7219 rChangeCurrRing(newRing); 7220 G = idrMoveR(Glp, lpRing,currRing); 7221 goto OMEGA_OVERFLOW_TRAN_NEW; 7222 } 7223 7224 plength3 = TRUE; 7225 pDelete(&p); 7226 break; 7227 } 7228 pDelete(&p); 7229 } 7230 7231 if(plength3 == FALSE) 7232 { 7233 rChangeCurrRing(newRing); 7234 G = idrMoveR(Glp, lpRing,currRing); 7235 goto TRAN_LIFTING; 7236 } 7237 7238 7239 npertstep = nwalk; 7240 nwalkpert = 1; 7241 nsteppert ++; 7242 7243 /* 7244 Print("\n// Subroutine needs (%d) steps.", nwalk); 7245 idElements(Glp, "last G in walk:"); 7246 PrintS("\n// ****************************************"); 7247 Print("\n// Perturb the original target vector (%d): ", nsteppert); 7248 ivString(target_weight, "new target"); 7249 PrintS("\n// ****************************************\n"); 7250 */ 7251 rChangeCurrRing(newRing); 7252 G = idrMoveR(Glp, lpRing,currRing); 7253 7254 delete next_weight; 7255 7256 //Print("\n// ring rNEW = %s;", rString(currRing)); 7257 goto COMPUTE_NEW_VECTOR; 7258 } 7259 7260 TRAN_LIFTING: 7261 for(i=nV-1; i>=0; i--) 7262 { 7263 (*curr_weight)[i] = (*next_weight)[i]; 7264 } 7265 delete next_weight; 7266 } // end of while 7267 #ifdef TEST_OVERFLOW 7268 BE_FINISH: 7269 #endif 7270 rChangeCurrRing(XXRing); 7271 G = idrMoveR(G, lpRing,currRing); 7272 7273 FINISH: 7274 delete ivNull; 7275 delete next_weight; 7276 delete iv_lp; 7277 omFree(npert); 7278 7279 #ifdef TIME_TEST 7280 Print("\n// Computation took %d steps and %.2f sec", nwalk, ((double) (clock()-mtim)/1000000)); 7281 7282 TimeStringFractal(tinput, tostd, tif, tstd, textra, tlift, tred, tnw); 7283 7284 Print("\n// pSetm_Error = (%d)", ErrorCheck()); 7285 Print("\n// Overflow_Error? (%d)\n", Overflow_Error); 7286 #endif 7287 7288 return(G); 7289 } 7290 7291 7292 /***************************************************************** 7293 * compute the reduced Groebner basis of an ideal <Go> w.r.t. lp * 7294 *****************************************************************/ 5172 7295 static ideal Mpwalk_MAltwalk1(ideal Go, intvec* curr_weight, int tp_deg) 5173 7296 { 5174 7297 Overflow_Error = FALSE; 5175 BOOLEAN nOverflow_Error = FALSE;7298 // BOOLEAN nOverflow_Error = FALSE; 5176 7299 clock_t tproc=0; 5177 7300 clock_t tinput=clock(); … … 5180 7303 int tp_deg_tmp = tp_deg; 5181 7304 ideal Gomega, M, F, G, M1, F1, Gomega1, Gomega2, G1; 5182 ring endRing,newRing, oldRing, TargetRing;7305 ring newRing, oldRing, TargetRing; 5183 7306 intvec* next_weight; 5184 7307 intvec* ivNull = new intvec(nV); 5185 intvec* extra_curr_weight = new intvec(nV);5186 7308 5187 7309 ring YXXRing = currRing; … … 5191 7313 ideal ssG; 5192 7314 5193 / * perturb the target vector */7315 // perturb the target vector 5194 7316 while(1) 5195 7317 { … … 5197 7319 { 5198 7320 if (rParameter(currRing) != NULL) 7321 { 5199 7322 DefRingParlp(); 7323 } 5200 7324 else 7325 { 5201 7326 VMrDefaultlp(); 5202 7327 } 5203 7328 TargetRing = currRing; 5204 7329 ssG = idrMoveR(Go,YXXRing,currRing); … … 5206 7331 Overflow_Error = FALSE; 5207 7332 if(tp_deg != 1) 7333 { 5208 7334 target_weight = MPertVectors(ssG, iv_M_dpp, tp_deg); 7335 } 5209 7336 else 5210 7337 { … … 5213 7340 } 5214 7341 if(Overflow_Error == FALSE) 7342 { 5215 7343 break; 5216 7344 } 5217 7345 Overflow_Error = TRUE; 5218 7346 tp_deg --; … … 5221 7349 { 5222 7350 Overflow_Error = TRUE; 5223 nOverflow_Error = TRUE;7351 //nOverflow_Error = TRUE; 5224 7352 } 5225 7353 … … 5228 7356 5229 7357 delete iv_M_dpp; 7358 #ifndef BUCHBERGER_ALG 5230 7359 intvec* hilb_func; 5231 5232 / * to avoid (1,0,...,0) as the target vector */7360 #endif 7361 // to avoid (1,0,...,0) as the target vector 5233 7362 intvec* last_omega = new intvec(nV); 5234 7363 for(i=nV-1; i>0; i--) 7364 { 5235 7365 (*last_omega)[i] = 1; 7366 } 5236 7367 (*last_omega)[0] = 10000; 5237 7368 … … 5245 7376 5246 7377 if(nwalk==1) 7378 { 5247 7379 goto FIRST_STEP; 5248 7380 } 5249 7381 to=clock(); 5250 / * compute an initial form ideal of <G> w.r.t. "curr_vector" */7382 // compute an initial form ideal of <G> w.r.t. "curr_vector" 5251 7383 Gomega = MwalkInitialForm(G, curr_weight); 5252 7384 xtif=xtif+clock()-to; 5253 #if 0 5254 if(Overflow_Error == TRUE) 5255 { 5256 for(i=nV-1; i>=0; i--) 5257 (*curr_weight)[i] = (*extra_curr_weight)[i]; 5258 delete extra_curr_weight; 5259 goto LASTGB_ALT1; 5260 } 5261 #endif 7385 5262 7386 #ifndef BUCHBERGER_ALG 5263 7387 if(isNolVector(curr_weight) == 0) … … 5265 7389 else 5266 7390 hilb_func = hFirstSeries(Gomega,NULL,NULL,last_omega,currRing); 5267 #endif // BUCHBERGER_ALG7391 #endif 5268 7392 5269 7393 oldRing = currRing; 5270 7394 5271 / * define a new ring that its ordering is "(a(curr_weight),lp) */7395 // define a new ring that its ordering is "(a(curr_weight),lp) 5272 7396 if (rParameter(currRing) != NULL) 7397 { 5273 7398 DefRingPar(curr_weight); 7399 } 5274 7400 else 7401 { 5275 7402 VMrDefault(curr_weight); 5276 7403 } 5277 7404 newRing = currRing; 5278 7405 Gomega1 = idrMoveR(Gomega, oldRing,currRing); … … 5288 7415 5289 7416 to=clock(); 5290 / * compute a reduced Groebner basis of <Gomega> w.r.t. "newRing" */7417 // compute a reduced Groebner basis of <Gomega> w.r.t. "newRing" 5291 7418 #ifdef BUCHBERGER_ALG 5292 7419 M = MstdhomCC(Gomega1); … … 5297 7424 xtstd=xtstd+clock()-to; 5298 7425 5299 / * change the ring to oldRing */7426 // change the ring to oldRing 5300 7427 rChangeCurrRing(oldRing); 5301 7428 M1 = idrMoveR(M, newRing,currRing); … … 5303 7430 to=clock(); 5304 7431 5305 // if(endwalks == 1) PrintS("\n// Lifting is still working:"); 5306 5307 /* compute a reduced Groebner basis of <G> w.r.t. "newRing" by the 5308 lifting process */ 7432 // if(endwalks == 1){PrintS("\n// Lifting is still working:");} 7433 7434 // compute a reduced Groebner basis of <G> w.r.t. "newRing" by the lifting process 5309 7435 F = MLifttwoIdeal(Gomega2, M1, G); 5310 7436 xtlift=xtlift+clock()-to; … … 5314 7440 idDelete(&G); 5315 7441 5316 / * change the ring to newRing */7442 // change the ring to newRing 5317 7443 rChangeCurrRing(newRing); 5318 7444 F1 = idrMoveR(F, oldRing,currRing); 5319 7445 to=clock(); 5320 //if(endwalks == 1) PrintS("\n// InterRed is still working:");5321 / * reduce the Groebner basis <G> w.r.t. the new ring */7446 //if(endwalks == 1){ PrintS("\n// InterRed is still working:");} 7447 // reduce the Groebner basis <G> w.r.t. the new ring 5322 7448 G = kInterRedCC(F1, NULL); 5323 7449 xtred=xtred+clock()-to; … … 5330 7456 Overflow_Error=FALSE; 5331 7457 to=clock(); 5332 / * compute a next weight vector */7458 // compute a next weight vector 5333 7459 next_weight = MkInterRedNextWeight(curr_weight,target_weight, G); 5334 7460 xtnw=xtnw+clock()-to; … … 5340 7466 { 5341 7467 delete next_weight; 5342 5343 LASTGB_ALT1:5344 7468 if(tp_deg > 1){ 5345 nOverflow_Error = Overflow_Error;7469 //nOverflow_Error = Overflow_Error; 5346 7470 tproc = tproc+clock()-tinput; 5347 /* 5348 Print("\n// A subroutine takes %d steps and calls \"Mpwalk\" (1,%d):", 5349 nwalk, tp_deg-1); 5350 */ 7471 //Print("\n// A subroutine takes %d steps and calls \"Mpwalk\" (1,%d):", nwalk, tp_deg-1); 5351 7472 G1 = Mpwalk_MAltwalk1(G, curr_weight, tp_deg-1); 5352 7473 goto MPW_Finish; … … 5366 7487 } 5367 7488 if(MivComp(next_weight, target_weight) == 1) 7489 { 5368 7490 endwalks = 1; 5369 7491 } 5370 7492 for(i=nV-1; i>=0; i--) 5371 7493 { … … 5376 7498 }//while 5377 7499 5378 / * check wheather the pertubed target vector is correct */7500 // check whether the pertubed target vector is correct 5379 7501 5380 7502 //define and execute ring with lex. order 5381 7503 if (rParameter(currRing) != NULL) 7504 { 5382 7505 DefRingParlp(); 7506 } 5383 7507 else 7508 { 5384 7509 VMrDefaultlp(); 5385 7510 } 5386 7511 G1 = idrMoveR(G, newRing,currRing); 5387 7512 … … 5389 7514 { 5390 7515 PrintS("\n// The perturbed target vector doesn't STAY in the correct cone!!"); 5391 if(tp_deg == 1){ 7516 if(tp_deg == 1) 7517 { 5392 7518 //Print("\n// subroutine takes %d steps and applys \"std\"", nwalk); 5393 7519 to=clock(); … … 5398 7524 G2 = NULL; 5399 7525 } 5400 else { 5401 nOverflow_Error = Overflow_Error; 7526 else 7527 { 7528 //nOverflow_Error = Overflow_Error; 5402 7529 tproc = tproc+clock()-tinput; 5403 /* 5404 Print("\n// B subroutine takes %d steps and calls \"Mpwalk\" (1,%d) :", 5405 nwalk, tp_deg-1); 5406 */ 7530 // Print("\n// B subroutine takes %d steps and calls \"Mpwalk\" (1,%d) :", nwalk, tp_deg-1); 5407 7531 G1 = Mpwalk_MAltwalk1(G1, curr_weight, tp_deg-1); 5408 7532 } … … 5417 7541 delete target_weight; 5418 7542 5419 /* 5420 Print("\n// \"Mpwalk\" (1,%d) took %d steps and %.2f sec. Overflow_Error (%d)", tp_deg, 5421 nwalk, ((double) clock()-tinput)/1000000, nOverflow_Error); 5422 */ 7543 //Print("\n// \"Mpwalk\" (1,%d) took %d steps and %.2f sec. Overflow_Error (%d)", tp_deg, nwalk, ((double) clock()-tinput)/1000000, nOverflow_Error); 5423 7544 5424 7545 return(result); 5425 7546 } 5426 7547 5427 /* August 2003 */ 7548 /******************************************************************* 7549 * Implementation of the first alternative Groebner Walk Algorithm * 7550 *******************************************************************/ 5428 7551 ideal MAltwalk1(ideal Go, int op_deg, int tp_deg, intvec* curr_weight, 5429 7552 intvec* target_weight) … … 5431 7554 Set_Error(FALSE ); 5432 7555 Overflow_Error = FALSE; 7556 #ifdef TIME_TEST 5433 7557 BOOLEAN nOverflow_Error = FALSE; 7558 #endif 5434 7559 // Print("// pSetm_Error = (%d)", ErrorCheck()); 5435 7560 … … 5442 7567 int nwalk=0, endwalks=0; 5443 7568 int op_tmp = op_deg; 5444 ideal Gomega, M, F, G, G 1, Gomega1, Gomega2, M1, F1;5445 ring endRing, newRing, oldRing, TargetRing;7569 ideal Gomega, M, F, G, Gomega1, Gomega2, M1, F1; 7570 ring newRing, oldRing; 5446 7571 intvec* next_weight; 5447 7572 intvec* iv_M_dp; … … 5449 7574 intvec* iv_dp = MivUnit(nV);// define (1,1,...,1) 5450 7575 intvec* exivlp = Mivlp(nV); 5451 intvec* extra_curr_weight = new intvec(nV); 7576 //intvec* extra_curr_weight = new intvec(nV); 7577 #ifndef BUCHBERGER_ALG 5452 7578 intvec* hilb_func; 5453 7579 #endif 5454 7580 intvec* cw_tmp = curr_weight; 5455 7581 5456 / * to avoid (1,0,...,0) as the target vector */7582 // to avoid (1,0,...,0) as the target vector 5457 7583 intvec* last_omega = new intvec(nV); 5458 7584 for(i=nV-1; i>0; i--) 7585 { 5459 7586 (*last_omega)[i] = 1; 7587 } 5460 7588 (*last_omega)[0] = 10000; 5461 7589 … … 5470 7598 if(Overflow_Error == FALSE) 5471 7599 { 5472 if(MivComp(curr_weight, iv_dp) == 1) //rOrdStr(currRing) = "dp" 5473 if(op_tmp == op_deg) { 7600 if(MivComp(curr_weight, iv_dp) == 1) 7601 { 7602 //rOrdStr(currRing) = "dp" 7603 if(op_tmp == op_deg) 7604 { 5474 7605 G = MstdCC(Go); 5475 7606 if(op_deg != 1) 7607 { 5476 7608 iv_M_dp = MivMatrixOrderdp(nV); 7609 } 5477 7610 } 7611 } 5478 7612 } 5479 7613 else 5480 7614 { 5481 if(op_tmp == op_deg) { 7615 if(op_tmp == op_deg) 7616 { 5482 7617 //rOrdStr(currRing) = (a(...),lp,C) 5483 7618 if (rParameter(currRing) != NULL) 7619 { 5484 7620 DefRingPar(cw_tmp); 7621 } 5485 7622 else 7623 { 5486 7624 VMrDefault(cw_tmp); 5487 7625 } 5488 7626 G = idrMoveR(Go, XXRing,currRing); 5489 7627 G = MstdCC(G); … … 5494 7632 Overflow_Error = FALSE; 5495 7633 if(op_deg != 1) 7634 { 5496 7635 curr_weight = MPertVectors(G, iv_M_dp, op_deg); 5497 else { 7636 } 7637 else 7638 { 5498 7639 curr_weight = cw_tmp; 5499 7640 break; 5500 7641 } 5501 7642 if(Overflow_Error == FALSE) 7643 { 5502 7644 break; 5503 7645 } 5504 7646 Overflow_Error = TRUE; 5505 7647 op_deg --; … … 5520 7662 5521 7663 to = clock(); 5522 / * compute an initial form ideal of <G> w.r.t. "curr_vector" */7664 // compute an initial form ideal of <G> w.r.t. "curr_vector" 5523 7665 Gomega = MwalkInitialForm(G, curr_weight); 5524 7666 xtif=xtif+clock()-to; … … 5536 7678 #ifndef BUCHBERGER_ALG 5537 7679 if(isNolVector(curr_weight) == 0) 7680 { 5538 7681 hilb_func = hFirstSeries(Gomega,NULL,NULL,curr_weight,currRing); 7682 } 5539 7683 else 7684 { 5540 7685 hilb_func = hFirstSeries(Gomega,NULL,NULL,last_omega,currRing); 7686 } 5541 7687 #endif // BUCHBERGER_ALG 5542 7688 5543 7689 oldRing = currRing; 5544 7690 5545 / * define a new ring which ordering is "(a(curr_weight),lp) */7691 // define a new ring which ordering is "(a(curr_weight),lp) 5546 7692 if (rParameter(currRing) != NULL) 7693 { 5547 7694 DefRingPar(curr_weight); 7695 } 5548 7696 else 7697 { 5549 7698 VMrDefault(curr_weight); 5550 7699 } 5551 7700 newRing = currRing; 5552 7701 Gomega1 = idrMoveR(Gomega, oldRing,currRing); 5553 7702 5554 7703 to=clock(); 5555 / * compute a reduced Groebner basis of <Gomega> w.r.t. "newRing" */7704 // compute a reduced Groebner basis of <Gomega> w.r.t. "newRing" 5556 7705 #ifdef BUCHBERGER_ALG 5557 7706 M = MstdhomCC(Gomega1); … … 5562 7711 xtstd=xtstd+clock()-to; 5563 7712 5564 / * change the ring to oldRing */7713 // change the ring to oldRing 5565 7714 rChangeCurrRing(oldRing); 5566 7715 M1 = idrMoveR(M, newRing,currRing); … … 5568 7717 5569 7718 to=clock(); 5570 /* compute a reduced Groebner basis of <G> w.r.t. "newRing" by the 5571 lifting process */ 7719 // compute a reduced Groebner basis of <G> w.r.t. "newRing" by the lifting process 5572 7720 F = MLifttwoIdeal(Gomega2, M1, G); 5573 7721 xtlift=xtlift+clock()-to; … … 5577 7725 idDelete(&G); 5578 7726 5579 /* change the ring to newRing */7727 // change the ring to newRing 5580 7728 rChangeCurrRing(newRing); 5581 7729 F1 = idrMoveR(F, oldRing,currRing); 5582 7730 5583 7731 to=clock(); 5584 / * reduce the Groebner basis <G> w.r.t. new ring */7732 // reduce the Groebner basis <G> w.r.t. new ring 5585 7733 G = kInterRedCC(F1, NULL); 5586 7734 xtred=xtred+clock()-to; … … 5588 7736 5589 7737 if(endwalks == 1) 7738 { 5590 7739 break; 5591 7740 } 5592 7741 NEXT_VECTOR: 5593 7742 to=clock(); 5594 / * compute a next weight vector */7743 // compute a next weight vector 5595 7744 next_weight = MkInterRedNextWeight(curr_weight,target_weight, G); 5596 7745 xtnw=xtnw+clock()-to; … … 5604 7753 5605 7754 if (rParameter(currRing) != NULL) 7755 { 5606 7756 DefRingPar(target_weight); 7757 } 5607 7758 else 7759 { 5608 7760 VMrDefault(target_weight); 5609 7761 } 5610 7762 F1 = idrMoveR(G, newRing,currRing); 5611 7763 G = MstdCC(F1); … … 5630 7782 else 5631 7783 { 5632 MSTD_ALT1: 7784 // MSTD_ALT1: 7785 #ifdef TIME_TEST 5633 7786 nOverflow_Error = Overflow_Error; 7787 #endif 5634 7788 tproc = clock()-xftinput; 5635 /* 5636 Print("\n// main routine takes %d steps and calls \"Mpwalk\" (1,%d):", 5637 nwalk, tp_deg); 5638 */ 5639 // compute the red. GB of <G> w.r.t. the lex order by 5640 // the "recursive-modified" perturbation walk alg (1,tp_deg) 7789 7790 Print("\n// main routine takes %d steps and calls \"Mpwalk\" (1,%d):", nwalk, tp_deg); 7791 7792 // compute the red. GB of <G> w.r.t. the lex order by the "recursive-modified" perturbation walk alg (1,tp_deg) 5641 7793 G = Mpwalk_MAltwalk1(G, curr_weight, tp_deg); 5642 7794 delete next_weight; … … 5645 7797 } 5646 7798 5647 / * 06.11.01 NOT Changed, to free memory*/7799 //NOT Changed, to free memory 5648 7800 for(i=nV-1; i>=0; i--) 5649 7801 { … … 5660 7812 delete ivNull; 5661 7813 if(op_deg != 1 ) 7814 { 5662 7815 delete curr_weight; 5663 7816 } 5664 7817 delete exivlp; 5665 7818 #ifdef TIME_TEST … … 5676 7829 return(result); 5677 7830 } 5678 -
Singular/walk.h
r8659a9 re92c6a 55 55 ideal Mwalk(ideal G, intvec* curr_weight, intvec* target_weight); 56 56 57 // random walk algorithm to compute a Groebner basis 58 ideal Mrwalk(ideal Go, intvec* curr_weight, intvec* target_weight, int weight_rad, int pert_deg); 59 57 60 /* the perturbation walk algorithm */ 58 61 ideal Mpwalk(ideal G,int op,int tp,intvec* curr_weight,intvec* target_weight, int nP); … … 61 64 ideal Mfwalk(ideal G, intvec* ivstart, intvec* ivtarget); 62 65 66 /* The fractal walk algorithm with random element */ 67 ideal Mfrwalk(ideal G, intvec* ivstart, intvec* ivtarget,int weight_rad); 63 68 64 69 /* Implement Tran's idea */ 65 70 intvec* TranMPertVectorslp(ideal G); 66 71 ideal TranMImprovwalk(ideal Go, intvec* curr_weight,intvec* target_weight, int nP); 72 73 /* Implement Tran's idea with random element*/ 74 ideal TranMrImprovwalk(ideal G,intvec* curr_weight,intvec* target_tmp, int nP, int weight_rad, int pert_deg); 67 75 68 76 /* the first alternative algorithm */
Note: See TracChangeset
for help on using the changeset viewer.