Changeset 8d7ca0 in git
- Timestamp:
- Jul 2, 2015, 7:33:23 PM (8 years ago)
- Branches:
- (u'spielwiese', '8e0ad00ce244dfd0756200662572aef8402f13d5')
- Children:
- 654a230847a8318a2a297f1520550f19f34a4401
- Parents:
- bbfc4a825e46a6bf0bdd30896d5189f16c4e4d99
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
Singular/walk.cc
rbbfc4a8 r8d7ca0 16 16 17 17 //#define TEST_OVERFLOW 18 //#define CHECK_IDEAL 19 //#define CHECK_IDEAL_MWALK 18 19 #define CHECK_IDEAL_MWALK //to print intermediate results 20 20 21 21 //#define NEXT_VECTORS_CC 22 //#define PRINT_VECTORS //to print vectors (sigma, tau, omega) 22 #define PRINT_VECTORS //to print weight vectors 23 23 24 24 #define INVEPS_SMALL_IN_FRACTAL //to choose the small invers of epsilon … … 27 27 28 28 #define FIRST_STEP_FRACTAL // to define the first step of the fractal 29 #define MSTDCC_FRACTAL // apply Buchberger alg to compute a red GB, if 30 // tau doesn't stay in the correct cone 29 #define MSTDCC_FRACTAL // apply Buchberger alg to compute a red GB, if tau doesn't stay in the correct cone 31 30 32 31 //#define TIME_TEST // print the used time of each subroutine … … 122 121 ************************************/ 123 122 // unused 124 #if 0 123 /* 125 124 static void initSSpecialCC (ideal F, ideal Q, ideal P,kStrategy strat) 126 125 { … … 270 269 #endif 271 270 } 272 #endif 271 */ 273 272 274 273 /***************** … … 279 278 int j; 280 279 kStrategy strat = new skStrategy; 281 282 // if (TEST_OPT_PROT) 283 // { 284 // writeTime("start InterRed:"); 285 // mflush(); 286 // } 287 //strat->syzComp = 0; 280 /* 281 if (TEST_OPT_PROT) 282 { 283 writeTime("start InterRed:"); 284 mflush(); 285 } 286 strat->syzComp = 0; 287 */ 288 288 strat->kHEdgeFound = (currRing->ppNoether) != NULL; 289 289 strat->kNoether=pCopy((currRing->ppNoether)); … … 348 348 strat->fromQ = NULL; 349 349 } 350 // if (TEST_OPT_PROT) 351 // { 352 // writeTime("end Interred:"); 353 // mflush(); 354 // } 350 /* 351 if (TEST_OPT_PROT) 352 { 353 writeTime("end Interred:"); 354 mflush(); 355 } 356 */ 355 357 ideal shdl=strat->Shdl; 356 358 idSkipZeroes(shdl); … … 360 362 } 361 363 362 //unused 363 #if 0 364 #ifdef TIME_TEST 364 365 static void TimeString(clock_t tinput, clock_t tostd, clock_t tif,clock_t tstd, 365 366 clock_t tlf,clock_t tred, clock_t tnw, int step) … … 399 400 ((((double) xtextra)/1000000)/totm)*100); 400 401 } 401 #endif 402 403 //unused 404 #if 0 402 405 403 static void TimeStringFractal(clock_t tinput, clock_t tostd, clock_t tif,clock_t tstd, 406 404 clock_t textra, clock_t tlf,clock_t tred, clock_t tnw) … … 431 429 #endif 432 430 433 //#ifdef CHECK_IDEAL_MWALK431 #ifdef CHECK_IDEAL_MWALK 434 432 static void idString(ideal L, const char* st) 435 433 { … … 443 441 Print(" %s;", pString(L->m[nL-1])); 444 442 } 445 //#endif446 443 #endif 444 /* 447 445 #if defined(CHECK_IDEAL_MWALK) || defined(ENDWALKS) 448 446 static void headidString(ideal L, char* st) … … 498 496 } 499 497 #endif 500 498 */ 501 499 502 500 static void ivString(intvec* iv, const char* ch) … … 512 510 } 513 511 514 //unused 515 //#if 0 512 #ifdef PRINT_VECTORS 516 513 static void MivString(intvec* iva, intvec* ivb, intvec* ivc) 517 514 { … … 536 533 Print("%d)", (*ivc)[nV]); 537 534 } 538 //#endif535 #endif 539 536 540 537 /******************************************************************** … … 1040 1037 * print the max total degree and the max coefficient of G * 1041 1038 *****************************************************************************/ 1042 #if 0 1039 /* 1043 1040 static void checkComplexity(ideal G, char* cG) 1044 1041 { … … 1081 1078 PrintLn(); 1082 1079 } 1083 #endif 1080 */ 1084 1081 1085 1082 /***************************************************************************** … … 1103 1100 intvec* v_null = new intvec(nV); 1104 1101 1105 1106 1102 // Check that the perturbed degree is valid 1107 1103 if(pdeg > nV || pdeg <= 0) … … 1117 1113 } 1118 1114 mpz_t *pert_vector = (mpz_t*)omAlloc(nV*sizeof(mpz_t)); 1119 //mpz_t *pert_vector1 = (mpz_t*)omAlloc(nV*sizeof(mpz_t));1115 mpz_t *pert_vector1 = (mpz_t*)omAlloc(nV*sizeof(mpz_t)); 1120 1116 1121 1117 for(i=0; i<nV; i++) 1122 1118 { 1123 1119 mpz_init_set_si(pert_vector[i], (*ivtarget)[i]); 1124 //mpz_init_set_si(pert_vector1[i], (*ivtarget)[i]);1120 mpz_init_set_si(pert_vector1[i], (*ivtarget)[i]); 1125 1121 } 1126 1122 // Calculate max1 = Max(A2)+Max(A3)+...+Max(Apdeg), … … 1202 1198 } 1203 1199 } 1200 1201 // 2147483647 is max. integer representation in SINGULAR 1202 mpz_t sing_int; 1203 mpz_init_set_ui(sing_int, 2147483647); 1204 1205 mpz_t check_int; 1206 mpz_init_set_ui(check_int, 100000); 1207 1204 1208 mpz_t ztemp; 1205 1209 mpz_init(ztemp); … … 1221 1225 } 1222 1226 1223 intvec *pert_vector1= new intvec(nV);1224 j = 0;1225 1227 for(i=0; i<nV; i++) 1226 1228 { 1227 (* pert_vector1)[i] = mpz_get_si(pert_vector[i]); 1228 (* pert_vector1)[i] = 0.1*(* pert_vector1)[i]; 1229 (* pert_vector1)[i] = floor((* pert_vector1)[i] + 0.5); 1230 if((* pert_vector1)[i] == 0) 1231 { 1232 j++; 1233 } 1234 } 1235 if(j > nV - 1) 1236 { 1237 // Print("\n// MPertVectors: geaenderter vector gleich Null! \n"); 1238 delete pert_vector1; 1239 goto CHECK_OVERFLOW; 1240 } 1241 1242 // check that the perturbed weight vector lies in the Groebner cone 1243 if(test_w_in_ConeCC(G,pert_vector1) != 0) 1244 { 1245 // Print("\n// MPertVectors: geaenderter vector liegt in Groebnerkegel! \n"); 1229 if(mpz_cmp(pert_vector[i], check_int)>=0) 1230 { 1231 for(j=0; j<nV; j++) 1232 { 1233 mpz_fdiv_q_ui(pert_vector1[j], pert_vector[j], 100); 1234 } 1235 } 1236 } 1237 1238 intvec* result = new intvec(nV); 1239 1240 int ntrue=0; 1241 1242 for(i=0; i<nV; i++) 1243 { 1244 (*result)[i] = mpz_get_si(pert_vector1[i]); 1245 if(mpz_cmp(pert_vector1[i], sing_int)>=0) 1246 { 1247 ntrue++; 1248 } 1249 } 1250 if(ntrue > 0 || test_w_in_ConeCC(G,result)==0) 1251 { 1252 ntrue=0; 1246 1253 for(i=0; i<nV; i++) 1247 1254 { 1248 mpz_set_si(pert_vector[i], (*pert_vector1)[i]); 1249 } 1250 } 1251 else 1252 { 1253 //Print("\n// MpertVectors: geaenderter vector liegt nicht in Groebnerkegel! \n"); 1254 } 1255 delete pert_vector1; 1256 1257 CHECK_OVERFLOW: 1258 intvec* result = new intvec(nV); 1259 1260 /* 2147483647 is max. integer representation in SINGULAR */ 1261 mpz_t sing_int; 1262 mpz_init_set_ui(sing_int, 2147483647); 1263 1264 int ntrue=0; 1265 for(i=0; i<nV; i++) 1266 { 1267 (*result)[i] = mpz_get_si(pert_vector[i]); 1268 if(mpz_cmp(pert_vector[i], sing_int)>=0) 1269 { 1270 ntrue++; 1271 if(Overflow_Error == FALSE) 1272 { 1273 Overflow_Error = TRUE; 1274 PrintS("\n// ** OVERFLOW in \"MPertvectors\": "); 1275 mpz_out_str( stdout, 10, pert_vector[i]); 1276 PrintS(" is greater than 2147483647 (max. integer representation)"); 1277 Print("\n// So vector[%d] := %d is wrong!!", i+1, (*result)[i]); 1278 } 1279 } 1280 } 1281 1282 if(Overflow_Error == TRUE) 1283 { 1284 ivString(result, "pert_vector"); 1285 Print("\n// %d element(s) of it is overflow!!", ntrue); 1255 (*result)[i] = mpz_get_si(pert_vector[i]); 1256 if(mpz_cmp(pert_vector[i], sing_int)>=0) 1257 { 1258 ntrue++; 1259 if(Overflow_Error == FALSE) 1260 { 1261 Overflow_Error = TRUE; 1262 PrintS("\n// ** OVERFLOW in \"MPertvectors\": "); 1263 mpz_out_str( stdout, 10, pert_vector[i]); 1264 PrintS(" is greater than 2147483647 (max. integer representation)"); 1265 Print("\n// So vector[%d] := %d is wrong!!", i+1, (*result)[i]); 1266 } 1267 } 1268 } 1269 1270 if(Overflow_Error == TRUE) 1271 { 1272 ivString(result, "pert_vector"); 1273 Print("\n// %d element(s) of it is overflow!!", ntrue); 1274 } 1286 1275 } 1287 1276 1288 1277 mpz_clear(ztemp); 1289 1278 mpz_clear(sing_int); 1279 mpz_clear(check_int); 1290 1280 omFree(pert_vector); 1291 //omFree(pert_vector1);1281 omFree(pert_vector1); 1292 1282 mpz_clear(tot_deg); 1293 1283 mpz_clear(maxdeg); … … 1491 1481 1492 1482 //unused 1493 #if 0 1483 /* 1494 1484 static intvec* MatrixOrderdp(int nV) 1495 1485 { … … 1507 1497 return(ivM); 1508 1498 } 1509 #endif 1499 */ 1510 1500 1511 1501 intvec* MivUnit(int nV) … … 1584 1574 mpz_cdiv_q_ui(inveps, inveps, nV); 1585 1575 } 1586 // PrintS("\n// choose the \"small\" inverse epsilon!");1576 // choose the small inverse epsilon 1587 1577 #endif 1588 1578 … … 1618 1608 1619 1609 for(j=0; j<nV; j++) 1620 1610 { 1621 1611 mpz_init_set(pert_vector[i*nV+j],ivtemp[j]); 1622 1623 } 1624 1625 / * 2147483647 is max. integer representation in SINGULAR */1612 } 1613 } 1614 1615 // 2147483647 is max. integer representation in SINGULAR 1626 1616 mpz_t sing_int; 1627 1617 mpz_init_set_ui(sing_int, 2147483647); … … 1647 1637 (* result)[i] = mpz_get_si(pert_vector[i]); 1648 1638 } 1649 /* 1650 j = 0; 1651 for(i=0; i<niv; i++) 1652 { 1653 (* result1)[i] = mpz_get_si(pert_vector[i]); 1654 (* result1)[i] = 0.1*(* result1)[i]; 1655 (* result1)[i] = floor((* result1)[i] + 0.5); 1656 if((* result1)[i] == 0) 1657 { 1658 j++; 1659 } 1660 } 1661 if(j > niv - 1) 1662 { 1663 // Print("\n// MfPertwalk: geaenderter vector gleich Null! \n"); 1664 delete result1; 1665 goto CHECK_OVERFLOW; 1666 } 1667 /* 1668 // check that the perturbed weight vector lies in the Groebner cone 1669 Print("\n========================================\n//** MfPertvector: test in Cone.\n"); 1670 if(test_w_in_ConeCC(G,result1) != 0) 1671 { 1672 // Print("\n// MfPertwalk: geaenderter vector liegt in Groebnerkegel! \n"); 1673 delete result; 1674 result = result1; 1675 for(i=0; i<nV; i++) 1676 { 1677 mpz_set_si(pert_vector[i], (*result1)[i]); 1678 } 1679 } 1680 else 1681 { 1682 delete result1; 1683 // Print("\n// Mfpertwalk: geaenderter vector liegt nicht in Groebnerkegel! \n"); 1684 } 1685 Print("\n ========================================\n");*/ 1639 1686 1640 CHECK_OVERFLOW: 1687 1641 … … 1721 1675 while(p!=NULL) 1722 1676 { 1723 p_Setm(p,currRing); pIter(p); 1677 p_Setm(p,currRing); 1678 pIter(p); 1724 1679 } 1725 1680 } … … 1804 1759 1805 1760 //unused 1806 #if 0 1761 /* 1807 1762 static void checkidealCC(ideal G, char* Ch) 1808 1763 { … … 1830 1785 PrintLn(); 1831 1786 } 1832 #endif 1787 */ 1833 1788 1834 1789 //unused 1835 #if 0 1790 /* 1836 1791 static void HeadidString(ideal L, char* st) 1837 1792 { … … 1845 1800 Print(" %s;\n", pString(pHead(L->m[nL]))); 1846 1801 } 1847 #endif 1848 1802 1803 */ 1849 1804 static inline int MivComp(intvec* iva, intvec* ivb) 1850 1805 { … … 1919 1874 * with respect to an ideal <G>. * 1920 1875 **********************************************************************/ 1876 /* 1921 1877 static intvec* MwalkNextWeightCC(intvec* curr_weight, intvec* target_weight, 1922 1878 ideal G) … … 2040 1996 if(mpz_cmp(t_nenner, t_null) == 0) 2041 1997 { 2042 2043 Print("\n//MwalkNextWeightCC: t_nenner ist Null!");2044 1998 #ifndef SING_NDEBUG 1999 Print("\n//MwalkNextWeightCC: t_nenner=0\n"); 2000 #endif 2045 2001 delete diff_weight; 2046 2002 diff_weight = ivCopy(curr_weight);//take memory … … 2171 2127 2172 2128 for(j=0; j<nRing; j++) 2173 2129 { 2174 2130 if(mpz_cmp(vec[j], sing_int_half) >= 0) 2175 2131 { 2176 2132 goto REDUCTION; 2177 2178 2133 } 2134 } 2179 2135 checkRed = 1; 2180 2136 for (j=0; j<nRing; j++) … … 2271 2227 return diff_weight; 2272 2228 } 2229 */ 2230 /********************************************************************** 2231 * Compute a next weight vector between curr_weight and target_weight * 2232 * with respect to an ideal <G>. * 2233 **********************************************************************/ 2234 static intvec* MwalkNextWeightCC(intvec* curr_weight, intvec* target_weight, 2235 ideal G) 2236 { 2237 BOOLEAN nError = Overflow_Error; 2238 Overflow_Error = FALSE; 2239 2240 assume(currRing != NULL && curr_weight != NULL && 2241 target_weight != NULL && G != NULL); 2242 2243 int nRing = currRing->N; 2244 int checkRed, j, nG = IDELEMS(G); 2245 intvec* ivtemp; 2246 2247 mpz_t t_zaehler, t_nenner; 2248 mpz_init(t_zaehler); 2249 mpz_init(t_nenner); 2250 2251 mpz_t s_zaehler, s_nenner, temp, MwWd; 2252 mpz_init(s_zaehler); 2253 mpz_init(s_nenner); 2254 mpz_init(temp); 2255 mpz_init(MwWd); 2256 2257 mpz_t sing_int; 2258 mpz_init(sing_int); 2259 mpz_set_si(sing_int, 2147483647); 2260 2261 mpz_t sing_int_half; 2262 mpz_init(sing_int_half); 2263 mpz_set_si(sing_int_half, 3*(1073741824/2)); 2264 2265 mpz_t deg_w0_p1, deg_d0_p1; 2266 mpz_init(deg_w0_p1); 2267 mpz_init(deg_d0_p1); 2268 2269 mpz_t sztn, sntz; 2270 mpz_init(sztn); 2271 mpz_init(sntz); 2272 2273 mpz_t t_null; 2274 mpz_init(t_null); 2275 2276 mpz_t ggt; 2277 mpz_init(ggt); 2278 2279 mpz_t dcw; 2280 mpz_init(dcw); 2281 2282 int gcd_tmp; 2283 //intvec* diff_weight = MivSub(target_weight, curr_weight); 2284 2285 intvec* diff_weight1 = new intvec(nRing); //MivSub(target_weight, curr_weight); 2286 poly g; 2287 2288 // reduce the size of the entries of the current weight vector 2289 if(TEST_OPT_REDSB) 2290 { 2291 for (j=0; j<nRing; j++) 2292 { 2293 (*diff_weight1)[j] = (*curr_weight)[j]; 2294 } 2295 while(MivAbsMax(diff_weight1)>10000 && test_w_in_ConeCC(G,diff_weight1)==1) 2296 { 2297 for(j=0; j<nRing; j++) 2298 { 2299 (*curr_weight)[j] = (*diff_weight1)[j]; 2300 } 2301 for(j=0; j<nRing; j++) 2302 { 2303 (*diff_weight1)[j] = floor(0.1*(*diff_weight1)[j] + 0.5); 2304 } 2305 } 2306 2307 if(MivAbsMax(curr_weight)>100000) 2308 { 2309 for(j=0; j<nRing; j++) 2310 { 2311 (*diff_weight1)[j] = (*curr_weight)[j]; 2312 } 2313 j = 0; 2314 while(test_w_in_ConeCC(G,diff_weight1)==1 && MivAbsMax(diff_weight1)>1000) 2315 { 2316 (*curr_weight)[j] = (*diff_weight1)[j]; 2317 j = MivAbsMaxArg(diff_weight1); 2318 (*diff_weight1)[j] = floor(0.1*(*diff_weight1)[j] + 0.5); 2319 } 2320 } 2321 2322 } 2323 intvec* diff_weight = MivSub(target_weight, curr_weight); 2324 2325 // compute a suitable next weight vector 2326 for (j=0; j<nG; j++) 2327 { 2328 g = G->m[j]; 2329 if (g != NULL) 2330 { 2331 ivtemp = MExpPol(g); 2332 mpz_set_si(deg_w0_p1, MivDotProduct(ivtemp, curr_weight)); 2333 mpz_set_si(deg_d0_p1, MivDotProduct(ivtemp, diff_weight)); 2334 delete ivtemp; 2335 2336 pIter(g); 2337 while (g != NULL) 2338 { 2339 ivtemp = MExpPol(g); 2340 mpz_set_si(MwWd, MivDotProduct(ivtemp, curr_weight)); 2341 mpz_sub(s_zaehler, deg_w0_p1, MwWd); 2342 if(mpz_cmp(s_zaehler, t_null) != 0) 2343 { 2344 mpz_set_si(MwWd, MivDotProduct(ivtemp, diff_weight)); 2345 mpz_sub(s_nenner, MwWd, deg_d0_p1); 2346 // check for 0 < s <= 1 2347 if( (mpz_cmp(s_zaehler,t_null) > 0 && 2348 mpz_cmp(s_nenner, s_zaehler)>=0) || 2349 (mpz_cmp(s_zaehler, t_null) < 0 && 2350 mpz_cmp(s_nenner, s_zaehler)<=0)) 2351 { 2352 // make both positive 2353 if (mpz_cmp(s_zaehler, t_null) < 0) 2354 { 2355 mpz_neg(s_zaehler, s_zaehler); 2356 mpz_neg(s_nenner, s_nenner); 2357 } 2358 2359 //compute a simple fraction of s 2360 cancel(s_zaehler, s_nenner); 2361 2362 if(mpz_cmp(t_nenner, t_null) != 0) 2363 { 2364 mpz_mul(sztn, s_zaehler, t_nenner); 2365 mpz_mul(sntz, s_nenner, t_zaehler); 2366 2367 if(mpz_cmp(sztn,sntz) < 0) 2368 { 2369 mpz_add(t_nenner, t_null, s_nenner); 2370 mpz_add(t_zaehler,t_null, s_zaehler); 2371 } 2372 } 2373 else 2374 { 2375 mpz_add(t_nenner, t_null, s_nenner); 2376 mpz_add(t_zaehler,t_null, s_zaehler); 2377 } 2378 } 2379 } 2380 pIter(g); 2381 delete ivtemp; 2382 } 2383 } 2384 } 2385 //Print("\n// Alloc Size = %d \n", nRing*sizeof(mpz_t)); 2386 mpz_t *vec=(mpz_t*)omAlloc(nRing*sizeof(mpz_t)); 2387 2388 2389 // there is no 0<t<1 and define the next weight vector that is equal 2390 // to the current weight vector 2391 if(mpz_cmp(t_nenner, t_null) == 0) 2392 { 2393 #ifndef SING_NDEBUG 2394 Print("\n//MwalkNextWeightCC: t_nenner=0\n"); 2395 #endif 2396 delete diff_weight; 2397 diff_weight = ivCopy(curr_weight);//take memory 2398 goto FINISH; 2399 } 2400 2401 // define the target vector as the next weight vector, if t = 1 2402 if(mpz_cmp_si(t_nenner, 1)==0 && mpz_cmp_si(t_zaehler,1)==0) 2403 { 2404 delete diff_weight; 2405 diff_weight = ivCopy(target_weight); //this takes memory 2406 goto FINISH; 2407 } 2408 2409 //checkRed = 0; 2410 2411 SIMPLIFY_GCD: 2412 2413 // simplify the vectors curr_weight and diff_weight (C-int) 2414 gcd_tmp = (*curr_weight)[0]; 2415 2416 for (j=1; j<nRing; j++) 2417 { 2418 gcd_tmp = gcd(gcd_tmp, (*curr_weight)[j]); 2419 if(gcd_tmp == 1) 2420 { 2421 break; 2422 } 2423 } 2424 if(gcd_tmp != 1) 2425 { 2426 for (j=0; j<nRing; j++) 2427 { 2428 gcd_tmp = gcd(gcd_tmp, (*diff_weight)[j]); 2429 if(gcd_tmp == 1) 2430 { 2431 break; 2432 } 2433 } 2434 } 2435 if(gcd_tmp != 1) 2436 { 2437 for (j=0; j<nRing; j++) 2438 { 2439 (*curr_weight)[j] = (*curr_weight)[j]/gcd_tmp; 2440 (*diff_weight)[j] = (*diff_weight)[j]/gcd_tmp; 2441 } 2442 } 2443 if(checkRed > 0) 2444 { 2445 for (j=0; j<nRing; j++) 2446 { 2447 mpz_set_si(vec[j], (*diff_weight)[j]); 2448 } 2449 goto TEST_OVERFLOW; 2450 } 2451 2452 #ifdef NEXT_VECTORS_CC 2453 Print("\n// gcd of the weight vectors (current and target) = %d", gcd_tmp); 2454 ivString(curr_weight, "new cw"); 2455 ivString(diff_weight, "new dw"); 2456 2457 PrintS("\n// t_zaehler: "); mpz_out_str( stdout, 10, t_zaehler); 2458 PrintS(", t_nenner: "); mpz_out_str( stdout, 10, t_nenner); 2459 #endif 2460 2461 // construct a new weight vector and check whether vec[j] is overflow, i.e. vec[j] > 2^31. 2462 // If vec[j] doesn't overflow, define a weight vector. Otherwise, report that overflow 2463 // appears. In the second case, test whether the the correctness of the new vector plays 2464 // an important role 2465 2466 for (j=0; j<nRing; j++) 2467 { 2468 mpz_set_si(dcw, (*curr_weight)[j]); 2469 mpz_mul(s_nenner, t_nenner, dcw); 2470 2471 if( (*diff_weight)[j]>0) 2472 { 2473 mpz_mul_ui(s_zaehler, t_zaehler, (*diff_weight)[j]); 2474 } 2475 else 2476 { 2477 mpz_mul_ui(s_zaehler, t_zaehler, -(*diff_weight)[j]); 2478 mpz_neg(s_zaehler, s_zaehler); 2479 } 2480 mpz_add(sntz, s_nenner, s_zaehler); 2481 mpz_init_set(vec[j], sntz); 2482 2483 #ifdef NEXT_VECTORS_CC 2484 Print("\n// j = %d ==> ", j); 2485 PrintS("("); 2486 mpz_out_str( stdout, 10, t_nenner); 2487 Print(" * %d)", (*curr_weight)[j]); 2488 Print(" + ("); mpz_out_str( stdout, 10, t_zaehler); 2489 Print(" * %d) = ", (*diff_weight)[j]); 2490 mpz_out_str( stdout, 10, s_nenner); 2491 PrintS(" + "); 2492 mpz_out_str( stdout, 10, s_zaehler); 2493 PrintS(" = "); mpz_out_str( stdout, 10, sntz); 2494 Print(" ==> vector[%d]: ", j); mpz_out_str(stdout, 10, vec[j]); 2495 #endif 2496 2497 if(j==0) 2498 { 2499 mpz_set(ggt, sntz); 2500 } 2501 else 2502 { 2503 if(mpz_cmp_si(ggt,1) != 0) 2504 { 2505 mpz_gcd(ggt, ggt, sntz); 2506 } 2507 } 2508 } 2509 // reduce the vector with the gcd 2510 if(mpz_cmp_si(ggt,1) != 0) 2511 { 2512 for (j=0; j<nRing; j++) 2513 { 2514 mpz_divexact(vec[j], vec[j], ggt); 2515 } 2516 } 2517 #ifdef NEXT_VECTORS_CC 2518 PrintS("\n// gcd of elements of the vector: "); 2519 mpz_out_str( stdout, 10, ggt); 2520 #endif 2521 2522 for (j=0; j<nRing; j++) 2523 { 2524 (*diff_weight)[j] = mpz_get_si(vec[j]); 2525 } 2526 2527 TEST_OVERFLOW: 2528 2529 for (j=0; j<nRing; j++) 2530 { 2531 if(mpz_cmp(vec[j], sing_int)>=0) 2532 { 2533 if(Overflow_Error == FALSE) 2534 { 2535 Overflow_Error = TRUE; 2536 PrintS("\n// ** OVERFLOW in \"MwalkNextWeightCC\": "); 2537 mpz_out_str( stdout, 10, vec[j]); 2538 PrintS(" is greater than 2147483647 (max. integer representation)\n"); 2539 Print("// So vector[%d] := %d is wrong!!\n",j+1, vec[j]);// vec[j] is mpz_t 2540 } 2541 } 2542 } 2543 2544 FINISH: 2545 delete diff_weight1; 2546 mpz_clear(t_zaehler); 2547 mpz_clear(t_nenner); 2548 mpz_clear(s_zaehler); 2549 mpz_clear(s_nenner); 2550 mpz_clear(sntz); 2551 mpz_clear(sztn); 2552 mpz_clear(temp); 2553 mpz_clear(MwWd); 2554 mpz_clear(deg_w0_p1); 2555 mpz_clear(deg_d0_p1); 2556 mpz_clear(ggt); 2557 omFree(vec); 2558 mpz_clear(sing_int_half); 2559 mpz_clear(sing_int); 2560 mpz_clear(dcw); 2561 mpz_clear(t_null); 2562 2563 if(Overflow_Error == FALSE) 2564 { 2565 Overflow_Error = nError; 2566 } 2567 rComplete(currRing); 2568 for(j=0; j<IDELEMS(G); j++) 2569 { 2570 poly p=G->m[j]; 2571 while(p!=NULL) 2572 { 2573 p_Setm(p,currRing); 2574 pIter(p); 2575 } 2576 } 2577 return diff_weight; 2578 } 2579 2273 2580 2274 2581 /********************************************************************** … … 2956 3263 int i,j,N = IDELEMS(Gomega); 2957 3264 poly p,lm,factor1,factor2; 2958 //PrintS("\n//** idCopy\n"); 3265 2959 3266 ideal Go = idCopy(G); 2960 3267 2961 //PrintS("\n//** jetzt for-Loop!\n");2962 2963 3268 // check whether leading monomials of G and Gomega coincide 2964 3269 // and return NULL if not 2965 3270 for(i=0; i<N; i++) 2966 3271 { 2967 p = pCopy(Gomega->m[i]); 2968 lm = pCopy(pHead(G->m[i])); 2969 if(!pIsConstant(pSub(p,lm))) 2970 { 2971 //pDelete(&p); 2972 //pDelete(&lm); 3272 if(!pIsConstant(pSub(pCopy(Gomega->m[i]),pCopy(pHead(G->m[i]))))) 3273 { 2973 3274 idDelete(&Go); 2974 3275 return NULL; 2975 3276 } 2976 //pDelete(&p);2977 //pDelete(&lm);2978 3277 } 2979 3278 for(i=0; i<N; i++) … … 3007 3306 } 3008 3307 } 3009 3010 //PrintS("\n//** jetzt Delete!\n"); 3011 //pDelete(&p); 3012 //pDelete(&factor); 3013 //pDelete(&lm); 3308 3014 3309 if(middle == TRUE) 3015 3310 { 3016 //PrintS("\n//** middle TRUE!\n");3017 3311 return Go; 3018 3312 } 3019 //PrintS("\n//** middle FALSE!\n");3020 3313 idDelete(&Go); 3021 3314 return NULL; … … 3151 3444 { 3152 3445 Print("\n// ring r%d_%d = %s;\n", tp_deg, nwalk, rString(currRing)); 3446 /* 3153 3447 idElements(Gomega, "Gw"); 3154 3448 headidString(Gomega, "Gw"); 3449 */ 3155 3450 } 3156 3451 #endif … … 3302 3597 3303 3598 /********************************************************** 3304 * check whether a polynomial of G has least 3monomials *3599 * check whether a polynomial of G has least 4 monomials * 3305 3600 **********************************************************/ 3306 3601 static int lengthpoly(ideal G) … … 3309 3604 for(i=IDELEMS(G)-1; i>=0; i--) 3310 3605 { 3311 #if 03312 if(pLength(G->m[i])>2)3313 {3314 return 1;3315 }3316 #else3317 3606 if((G->m[i]!=NULL) /* len >=0 */ 3318 3607 && (G->m[i]->next!=NULL) /* len >=1 */ 3319 3608 && (G->m[i]->next->next!=NULL) /* len >=2 */ 3320 3609 && (G->m[i]->next->next->next!=NULL) /* len >=3 */ 3321 //&& (G->m[i]->next->next->next->next!=NULL) /* len >=4 */ 3322 ) 3323 { 3324 return 1; 3325 } 3326 #endif 3610 && (G->m[i]->next->next->next->next!=NULL) /* len >=4*/ ) 3611 { 3612 return 1; 3613 } 3327 3614 } 3328 3615 return 0; … … 3414 3701 for(i=nG-1; i>=0; i--) 3415 3702 { 3416 #if 0 3703 /* 3417 3704 poly t; 3418 3705 if((t=pSub(pCopy(H0->m[i]), pCopy(H1->m[i]))) != NULL) … … 3422 3709 } 3423 3710 pDelete(&t); 3424 #else 3711 */ 3425 3712 if(!pEqualPolys(H0->m[i],H1->m[i])) 3426 3713 { 3427 3714 return 0; 3428 3715 } 3429 #endif3430 3716 } 3431 3717 return 1; … … 3436 3722 * find the maximal total degree of polynomials in G * 3437 3723 *****************************************************/ 3438 #if 0 3724 /* 3439 3725 static int Trandegreebound(ideal G) 3440 3726 { … … 3457 3743 return result; 3458 3744 } 3459 #endif 3745 */ 3460 3746 3461 3747 //unused … … 3804 4090 * basis or n times, where n is the numbers of variables. * 3805 4091 *****************************************************************************/ 3806 3807 //unused3808 #if 03809 static int testnegintvec(intvec* v)3810 {3811 int n = v->length();3812 int i;3813 for(i=0; i<n; i++)3814 {3815 if((*v)[i]<0)3816 {3817 return(1);3818 }3819 }3820 return(0);3821 }3822 #endif3823 4092 3824 4093 // npwinc = 0, if curr_weight doesn't stay in the correct Groebner cone … … 4032 4301 //nOverflow_Error = Overflow_Error; 4033 4302 tproc=tproc+clock()-tinput; 4034 /*4035 4036 4037 */4303 4304 Print("\n// takes %d steps and calls \"Rec_LastGB\" (%d):", 4305 nwalk, tp_deg+1); 4306 4038 4307 G = Rec_LastGB(G,curr_weight, orig_target_weight, tp_deg+1,nnwinC); 4039 4308 newRing = currRing; … … 4069 4338 { 4070 4339 // nOverflow_Error = Overflow_Error; 4071 //Print("\n// takes %d steps and calls \"Rec_LastGB (%d):", tp_deg+1);4340 Print("\n// takes %d steps and calls \"Rec_LastGB (%d):", tp_deg+1); 4072 4341 tproc=tproc+clock()-tinput; 4073 4342 F1 = Rec_LastGB(F1,curr_weight, orig_target_weight, tp_deg+1,nnwinC); … … 4123 4392 Overflow_Error=nError; 4124 4393 } 4125 // Print("\n// \"Rec_LastGB\" (%d) took %d steps and %.2f sec.Overflow_Error (%d)", tp_deg, nwalk, ((double) tproc)/1000000, nOverflow_Error); 4394 #ifdef TIME_TEST 4395 Print("\n// \"Rec_LastGB\" (%d) took %d steps and %.2f sec.Overflow_Error (%d)", tp_deg, nwalk, ((double) tproc)/1000000, nOverflow_Error); 4396 #endif 4126 4397 return(result); 4127 4398 } … … 4145 4416 //BOOLEAN nOverflow_Error = FALSE; 4146 4417 //Print("// pSetm_Error = (%d)", ErrorCheck()); 4147 4418 #ifdef TIME_TEST 4148 4419 xtif=0; xtstd=0; xtlift=0; xtred=0; xtnw=0; xtextra=0; 4149 4420 xftinput = clock(); 4150 4421 clock_t tostd, tproc; 4151 4422 #endif 4152 4423 nstep = 0; 4153 4424 int i, nV = currRing->N; … … 4160 4431 intvec* ivNull = new intvec(nV); 4161 4432 intvec* next_weight; 4162 #if 0 4163 intvec* extra_curr_weight = new intvec(nV); 4164 #endif 4433 //intvec* extra_curr_weight = new intvec(nV); 4165 4434 //intvec* hilb_func; 4166 4435 intvec* exivlp = Mivlp(nV); 4167 4168 4436 ring XXRing = currRing; 4169 4437 4170 4438 //Print("\n// ring r_input = %s;", rString(currRing)); 4439 #ifdef TIME_TEST 4171 4440 to = clock(); 4441 #endif 4172 4442 /* compute the reduced Groebner basis of the given ideal w.r.t. 4173 4443 a "fast" monomial order, e.g. degree reverse lex. order (dp) */ 4174 4444 G = MstdCC(Go); 4445 #ifdef TIME_TEST 4175 4446 tostd=clock()-to; 4176 4447 4177 /*4178 4448 Print("\n// Computation of the first std took = %.2f sec", 4179 4449 ((double) tostd)/1000000); 4180 */ 4450 #endif 4181 4451 if(currRing->order[0] == ringorder_a) 4182 4452 { … … 4187 4457 nwalk ++; 4188 4458 nstep ++; 4459 #ifdef TIME_TEST 4189 4460 to = clock(); 4461 #endif 4190 4462 /* compute an initial form ideal of <G> w.r.t. "curr_vector" */ 4191 4463 Gomega = MwalkInitialForm(G, curr_weight); 4464 #ifdef TIME_TEST 4192 4465 xtif=xtif+clock()-to; 4193 #if 0 4466 #endif 4467 /* 4194 4468 if(Overflow_Error == TRUE) 4195 4469 { … … 4199 4473 goto LAST_GB_ALT2; 4200 4474 } 4201 #endif 4475 */ 4202 4476 oldRing = currRing; 4203 4477 … … 4213 4487 newRing = currRing; 4214 4488 Gomega1 = idrMoveR(Gomega, oldRing,currRing); 4489 #ifdef TIME_TEST 4215 4490 to = clock(); 4491 #endif 4216 4492 /* compute a reduced Groebner basis of <Gomega> w.r.t. "newRing" */ 4217 4493 M = MstdhomCC(Gomega1); 4494 #ifdef TIME_TEST 4218 4495 xtstd=xtstd+clock()-to; 4496 #endif 4219 4497 /* change the ring to oldRing */ 4220 4498 rChangeCurrRing(oldRing); 4221 4499 M1 = idrMoveR(M, newRing,currRing); 4222 4500 Gomega2 = idrMoveR(Gomega1, newRing,currRing); 4223 4501 #ifdef TIME_TEST 4224 4502 to = clock(); 4503 #endif 4225 4504 /* compute the reduced Groebner basis of <G> w.r.t. "newRing" 4226 4505 by the liftig process */ 4227 4506 F = MLifttwoIdeal(Gomega2, M1, G); 4507 #ifdef TIME_TEST 4228 4508 xtlift=xtlift+clock()-to; 4509 #endif 4229 4510 idDelete(&M1); 4230 4511 idDelete(&Gomega2); … … 4234 4515 rChangeCurrRing(newRing); 4235 4516 F1 = idrMoveR(F, oldRing,currRing); 4236 4517 #ifdef TIME_TEST 4237 4518 to = clock(); 4519 #endif 4238 4520 /* reduce the Groebner basis <G> w.r.t. newRing */ 4239 4521 G = kInterRedCC(F1, NULL); 4522 #ifdef TIME_TEST 4240 4523 xtred=xtred+clock()-to; 4524 #endif 4241 4525 idDelete(&F1); 4242 4526 … … 4245 4529 4246 4530 NEXT_VECTOR: 4531 #ifdef TIME_TEST 4247 4532 to = clock(); 4533 #endif 4248 4534 /* compute a next weight vector */ 4249 4535 next_weight = MkInterRedNextWeight(curr_weight,target_weight, G); 4536 #ifdef TIME_TEST 4250 4537 xtnw=xtnw+clock()-to; 4538 #endif 4251 4539 #ifdef PRINT_VECTORS 4252 4540 MivString(curr_weight, target_weight, next_weight); … … 4292 4580 // LAST_GB_ALT2: 4293 4581 //nOverflow_Error = Overflow_Error; 4582 #ifdef TIME_TEST 4294 4583 tproc = clock()-xftinput; 4584 #endif 4295 4585 //Print("\n// takes %d steps and calls the recursion of level 2:", nwalk); 4296 4586 /* call the changed perturbation walk algorithm with degree 2 */ … … 4319 4609 4320 4610 #ifdef TIME_TEST 4321 // Print("\n// \"Main procedure\" took %d steps dnd %.2f sec. Overflow_Error (%d)", nwalk, ((double) tproc)/1000000, nOverflow_Error); 4611 Print("\n// \"Main procedure\" took %d steps dnd %.2f sec. Overflow_Error (%d)", 4612 nwalk, ((double) tproc)/1000000, nOverflow_Error); 4322 4613 4323 4614 TimeStringFractal(xftinput, tostd, xtif, xtstd, xtextra,xtlift, xtred,xtnw); … … 4355 4646 * compute a next weight vector * 4356 4647 ********************************/ 4357 static intvec* MWalkRandomNextWeight(ideal G, intvec* curr_weight,4358 intvec* target_weight,int weight_rad, int pert_deg)4648 static intvec* MWalkRandomNextWeight(ideal G, intvec* orig_M, intvec* target_weight, 4649 int weight_rad, int pert_deg) 4359 4650 { 4360 assume(currRing != NULL && curr_weight!= NULL &&4651 assume(currRing != NULL && orig_M != NULL && 4361 4652 target_weight != NULL && G->m[0] != NULL); 4362 4653 4654 //BOOLEAN nError = Overflow_Error; 4655 Overflow_Error = FALSE; 4656 4657 BOOLEAN found_random_weight = FALSE; 4658 int i,nV = currRing->N; 4659 intvec* curr_weight = new intvec(nV); 4660 4661 for(i=0; i<nV; i++) 4662 { 4663 (*curr_weight)[i] = (*orig_M)[i]; 4664 } 4665 4666 int k=0,weight_norm; 4667 intvec* next_weight; 4363 4668 intvec* next_weight1 = MkInterRedNextWeight(curr_weight,target_weight,G); 4364 if(weight_rad == 0)4365 {4366 return(next_weight1);4367 }4368 4369 int i,weight_norm,nV = currRing->N;4370 4669 intvec* next_weight2 = new intvec(nV); 4371 4670 intvec* next_weight22 = new intvec(nV); 4372 4671 intvec* result = new intvec(nV); 4373 4374 ideal G_test, G_test 2;4375 4376 // computea random next weight vector "next_weight2"4377 while(1)4672 intvec* curr_weight1; 4673 ideal G_test, G_test1, G_test2; 4674 4675 //try to find a random next weight vector "next_weight2" 4676 if(weight_rad > 0){ while(k<10) 4378 4677 { 4379 4678 weight_norm = 0; … … 4401 4700 } 4402 4701 4403 if(test_w_in_ConeCC(G, next_weight2) == 1) 4404 { 4405 PrintS("\n Hier 1\n"); 4406 G_test2 = MwalkInitialForm(G, next_weight2); 4407 PrintS("\n Hier 2\n"); 4408 if(maxlengthpoly(G_test2) <= 1) 4702 if(test_w_in_ConeCC(G,next_weight2) == 1) 4703 { 4704 if(maxlengthpoly(MwalkInitialForm(G,next_weight2))<2) 4409 4705 { 4410 4706 next_weight2 = MkInterRedNextWeight(next_weight2,target_weight,G); 4411 4707 } 4412 idDelete(&G_test2); 4413 4708 /* 4414 4709 if(MivAbsMax(next_weight2)>1147483647) 4415 4710 { … … 4428 4723 delete next_weight22; 4429 4724 } 4725 */ 4726 G_test2 = MwalkInitialForm(G, next_weight2); 4727 found_random_weight = TRUE; 4430 4728 break; 4431 4729 } 4432 weight_rad++; 4433 } 4434 4435 // compute "usual" next weight vector 4436 intvec* next_weight = MwalkNextWeightCC(curr_weight,target_weight, G); 4437 G_test = MwalkInitialForm(G, next_weight); 4438 G_test2 = MwalkInitialForm(G, next_weight2); 4439 4730 k++; 4731 }} 4732 Print("\n MWalkRandomNextWeight: compute perurbation...\n"); 4733 // compute "perturbed" next weight vector 4734 if(pert_deg > 1) 4735 { 4736 curr_weight1 = MPertVectors(G,orig_M,pert_deg); 4737 next_weight = MkInterRedNextWeight(curr_weight1,target_weight,G); 4738 delete curr_weight1; 4739 } 4740 else 4741 { 4742 next_weight = MkInterRedNextWeight(curr_weight,target_weight,G); 4743 } 4744 if(MivSame(curr_weight,next_weight)==1 || Overflow_Error == TRUE) 4745 { 4746 Overflow_Error = FALSE; 4747 delete next_weight; 4748 next_weight = MkInterRedNextWeight(curr_weight,target_weight,G); 4749 } 4750 G_test=MwalkInitialForm(G,next_weight); 4751 G_test1=MwalkInitialForm(G,next_weight1); 4752 Print("\n MWalkRandomNextWeight: finished...\n"); 4440 4753 // compare next weights 4441 4754 if(Overflow_Error == FALSE) 4442 4755 { 4443 i deal G_test1 = MwalkInitialForm(G, next_weight1);4444 if(G_test1->m[0] != NULL && maxlengthpoly(G_test1) < maxlengthpoly(G_test))4445 {4446 if(G_test 2->m[0] != NULL && maxlengthpoly(G_test2) < maxlengthpoly(G_test1))4447 { 4448 for(i=0; i<nV; i++)4756 if(found_random_weight == TRUE) 4757 { 4758 // random next weight vector found 4759 if(G_test1->m[0] != NULL && maxlengthpoly(G_test1) < maxlengthpoly(G_test)) 4760 { 4761 if(G_test2->m[0] != NULL && maxlengthpoly(G_test2) < maxlengthpoly(G_test1)) 4449 4762 { 4450 (*result)[i] = (*next_weight2)[i]; 4763 for(i=0; i<nV; i++) 4764 { 4765 (*result)[i] = (*next_weight2)[i]; 4766 } 4451 4767 } 4768 else 4769 { 4770 for(i=0; i<nV; i++) 4771 { 4772 (*result)[i] = (*next_weight1)[i]; 4773 } 4774 } 4452 4775 } 4453 4776 else 4454 4777 { 4455 for(i=0; i<nV; i++) 4778 if(G_test2->m[0] != NULL && maxlengthpoly(G_test2) < maxlengthpoly(G_test)) 4779 { 4780 for(i=0; i<nV; i++) 4781 { 4782 (*result)[i] = (*next_weight2)[i]; 4783 } 4784 } 4785 else 4786 { 4787 for(i=0; i<nV; i++) 4788 { 4789 (*result)[i] = (*next_weight)[i]; 4790 } 4791 } 4792 } 4793 } 4794 else 4795 { 4796 // no random next weight vector found 4797 if(G_test1->m[0] != NULL && maxlengthpoly(G_test1) < maxlengthpoly(G_test)) 4798 { 4799 for(i=0; i<nV; i++) 4456 4800 { 4457 4801 (*result)[i] = (*next_weight1)[i]; 4458 }4459 }4460 }4461 else4462 {4463 if(G_test2->m[0] != NULL && maxlengthpoly(G_test2) < maxlengthpoly(G_test))4464 {4465 for(i=0; i<nV; i++)4466 {4467 (*result)[i] = (*next_weight2)[i];4468 4802 } 4469 4803 } … … 4476 4810 } 4477 4811 } 4478 idDelete(&G_test1);4479 4812 } 4480 4813 else 4481 4814 { 4482 4815 Overflow_Error = FALSE; 4483 if(G_test2->m[0] != NULL && maxlengthpoly(G_test2) < maxlengthpoly(G_test)) 4484 { 4485 for(i=1; i<nV; i++) 4486 { 4487 (*result)[i] = (*next_weight2)[i]; 4816 if(found_random_weight == TRUE) 4817 { 4818 if(G_test2->m[0] != NULL && maxlengthpoly(G_test2) < maxlengthpoly(G_test)) 4819 { 4820 for(i=1; i<nV; i++) 4821 { 4822 (*result)[i] = (*next_weight2)[i]; 4823 } 4824 } 4825 else 4826 { 4827 for(i=0; i<nV; i++) 4828 { 4829 (*result)[i] = (*next_weight)[i]; 4830 } 4488 4831 } 4489 4832 } … … 4496 4839 } 4497 4840 } 4498 4841 // delete curr_weight1; 4842 delete next_weight; 4843 delete next_weight2; 4499 4844 idDelete(&G_test); 4500 idDelete(&G_test2); 4501 if(test_w_in_ConeCC(G, result) == 1) 4502 { 4503 delete next_weight2; 4504 delete next_weight; 4845 idDelete(&G_test1); 4846 if(found_random_weight == TRUE) 4847 { 4848 idDelete(&G_test2); 4849 } 4850 if(test_w_in_ConeCC(G, result) == 1 && MivSame(curr_weight,result)==0) 4851 { 4852 delete curr_weight; 4505 4853 delete next_weight1; 4506 4854 return result; … … 4508 4856 else 4509 4857 { 4858 delete curr_weight; 4510 4859 delete result; 4511 delete next_weight2; 4512 delete next_weight1; 4513 return next_weight; 4860 return next_weight1; 4514 4861 } 4515 4862 } … … 4519 4866 * The procedur REC_GB_Mwalk computes a GB for <G> w.r.t. the weight order * 4520 4867 * otw, where G is a reduced GB w.r.t. the weight order cw. * 4521 * The new procedur Mwalk calls REC_GB.*4868 * The new procedure Mwalk calls REC_GB. * 4522 4869 ***************************************************************************/ 4523 4870 static ideal REC_GB_Mwalk(ideal G, intvec* curr_weight, intvec* orig_target_weight, … … 4822 5169 4823 5170 ideal Gomega, M, F, Gomega1, Gomega2, M1, F1, G; 4824 //ideal G1; 4825 //ring endRing; 5171 4826 5172 ring newRing, oldRing; 4827 5173 intvec* ivNull = new intvec(nV); … … 4865 5211 the recursive changed perturbation walk alg. */ 4866 5212 tim = clock(); 4867 /* 5213 #ifdef CHECK_IDEAL_MWALK 4868 5214 Print("\n// **** Groebnerwalk took %d steps and ", nwalk); 4869 5215 PrintS("\n// **** call the rec. Pert. Walk to compute a red GB of:"); 4870 id Elements(Gomega, "G_omega");4871 */ 5216 idString(Gomega, "Gomega"); 5217 #endif 4872 5218 4873 5219 if(MivSame(exivlp, target_weight)==1) … … 4875 5221 else 4876 5222 goto NORMAL_GW; 4877 /* 5223 #ifdef TIME_TEST 4878 5224 Print("\n// time for the last std(Gw) = %.2f sec", 4879 5225 ((double) (clock()-tim)/1000000)); 4880 PrintS("\n// ***************************************************\n"); 4881 */ 5226 #endif 5227 /* 4882 5228 #ifdef CHECK_IDEAL_MWALK 4883 5229 idElements(Gomega, "G_omega"); … … 4886 5232 //headidString(M, "M"); 4887 5233 #endif 5234 */ 4888 5235 to = clock(); 4889 5236 F = MLifttwoIdeal(Gomega, M, G); … … 5045 5392 } 5046 5393 5047 5048 5394 /******************************* 5049 5395 * THE GROEBNER WALK ALGORITHM * … … 5103 5449 #endif 5104 5450 rComplete(currRing); 5105 //#ifdef CHECK_IDEAL_MWALK5451 #ifdef CHECK_IDEAL_MWALK 5106 5452 if(printout > 2) 5107 5453 { 5108 5454 idString(Go,"//** Mwalk: Go"); 5109 5455 } 5110 //#endif5456 #endif 5111 5457 5112 5458 if(target_M->length() == nV) … … 5133 5479 to = clock(); 5134 5480 #endif 5135 5136 5481 ideal G = MstdCC(idrMoveR(Go,baseRing,currRing)); 5137 5138 5482 #ifdef TIME_TEST 5139 5483 tostd = clock()-to; … … 5147 5491 nwalk ++; 5148 5492 nstep ++; 5149 5493 //compute an initial form ideal of <G> w.r.t. "curr_vector" 5150 5494 #ifdef TIME_TEST 5151 5495 to = clock(); 5152 5496 #endif 5153 // compute an initial form ideal of <G> w.r.t. "curr_vector"5154 5497 Gomega = MwalkInitialForm(G, curr_weight); 5155 5498 #ifdef TIME_TEST … … 5157 5500 #endif 5158 5501 5159 //#ifdef CHECK_IDEAL_MWALK5502 #ifdef CHECK_IDEAL_MWALK 5160 5503 if(printout > 1) 5161 5504 { 5162 5505 idString(Gomega,"//** Mwalk: Gomega"); 5163 5506 } 5164 //#endif5507 #endif 5165 5508 5166 5509 if(reduction == 0) 5167 5510 { 5168 5511 FF = middleOfCone(G,Gomega); 5169 if( 5512 if(FF != NULL) 5170 5513 { 5171 5514 idDelete(&G); … … 5229 5572 #endif 5230 5573 idSkipZeroes(M); 5231 //#ifdef CHECK_IDEAL_MWALK5574 #ifdef CHECK_IDEAL_MWALK 5232 5575 if(printout > 2) 5233 5576 { 5234 5577 idString(M, "//** Mwalk: M"); 5235 5578 } 5236 //#endif5579 #endif 5237 5580 //change the ring to baseRing 5238 5581 rChangeCurrRing(baseRing); … … 5247 5590 // where Gomega is a reduced Groebner basis w.r.t. the current ring 5248 5591 F = MLifttwoIdeal(Gomega2, M1, G); 5249 5250 5592 #ifdef TIME_TEST 5251 5593 tlift = tlift + clock() - to; 5252 5594 #endif 5253 //#ifdef CHECK_IDEAL_MWALK5595 #ifdef CHECK_IDEAL_MWALK 5254 5596 if(printout > 2) 5255 5597 { 5256 5598 idString(F, "//** Mwalk: F"); 5257 5599 } 5258 //#endif5600 #endif 5259 5601 idDelete(&Gomega2); 5260 5602 idDelete(&M1); … … 5265 5607 idSkipZeroes(G); 5266 5608 5267 //#ifdef CHECK_IDEAL_MWALK5609 #ifdef CHECK_IDEAL_MWALK 5268 5610 if(printout > 2) 5269 5611 { 5270 5612 idString(G, "//** Mwalk: G"); 5271 5613 } 5272 //#endif5614 #endif 5273 5615 5274 5616 rChangeCurrRing(targetRing); … … 5286 5628 baseRing = currRing; 5287 5629 5630 NEXT_VECTOR: 5631 #ifdef TIME_TEST 5632 to = clock(); 5633 #endif 5634 intvec* next_weight = MwalkNextWeightCC(curr_weight,target_weight,G); 5635 #ifdef TIME_TEST 5636 tnw = tnw + clock() - to; 5637 #endif 5638 #ifdef PRINT_VECTORS 5639 if(printout > 0) 5640 { 5641 MivString(curr_weight, target_weight, next_weight); 5642 } 5643 #endif 5644 if(MivComp(target_weight,curr_weight) == 1)// || endwalks == TRUE) 5645 { 5288 5646 /* 5289 #ifdef TIME_TEST 5290 to = clock(); 5291 #endif 5292 5293 #ifdef TIME_TEST 5294 tstd = tstd + clock() - to; 5295 #endif 5296 */ 5297 5298 5299 #ifdef TIME_TEST 5300 to = clock(); 5301 #endif 5302 NEXT_VECTOR: 5303 intvec* next_weight = MwalkNextWeightCC(curr_weight,target_weight,G); 5304 #ifdef TIME_TEST 5305 tnw = tnw + clock() - to; 5306 #endif 5307 //#ifdef PRINT_VECTORS 5308 if(printout > 0) 5309 { 5310 MivString(curr_weight, target_weight, next_weight); 5311 } 5312 //#endif 5313 if(MivComp(target_weight,curr_weight) == 1)// || endwalks == TRUE) 5314 {/* 5315 //#ifdef CHECK_IDEAL_MWALK 5647 #ifdef CHECK_IDEAL_MWALK 5316 5648 if(printout > 0) 5317 5649 { 5318 5650 PrintS("\n//** Mwalk: entering last cone.\n"); 5319 5651 } 5320 //#endif5652 #endif 5321 5653 5322 5654 Gomega = MwalkInitialForm(G, curr_weight); // compute an initial form ideal of <G> w.r.t. "curr_vector" … … 5332 5664 Gomega1 = idrMoveR(Gomega, baseRing,currRing); 5333 5665 idDelete(&Gomega); 5334 //#ifdef CHECK_IDEAL_MWALK5666 #ifdef CHECK_IDEAL_MWALK 5335 5667 if(printout > 1) 5336 5668 { … … 5338 5670 } 5339 5671 PrintS("\n //** Mwalk: kStd(Gomega)"); 5340 //#endif5672 #endif 5341 5673 M = kStd(Gomega1,NULL,testHomog,NULL,NULL,0,0,NULL); 5342 //#ifdef CHECK_IDEAL_MWALK5674 #ifdef CHECK_IDEAL_MWALK 5343 5675 if(printout > 1) 5344 5676 { 5345 5677 idString(M,"//** Mwalk: M"); 5346 5678 } 5347 //#endif5679 #endif 5348 5680 rChangeCurrRing(baseRing); 5349 5681 M1 = idrMoveR(M, newRing,currRing); … … 5353 5685 //PrintS("\n //** Mwalk: MLifttwoIdeal"); 5354 5686 F = MLifttwoIdeal(Gomega2, M1, G); 5355 //#ifdef CHECK_IDEAL_MWALK5687 #ifdef CHECK_IDEAL_MWALK 5356 5688 if(printout > 2) 5357 5689 { 5358 5690 idString(F,"//** Mwalk: F"); 5359 5691 } 5360 //#endif5692 #endif 5361 5693 idDelete(&Gomega2); 5362 5694 idDelete(&M1); … … 5365 5697 idDelete(&F); 5366 5698 baseRing = currRing; 5367 si_opt_1 = save1; //set original options, e. g. option(RedSB)5368 5699 idSkipZeroes(G); 5369 5700 #ifdef TIME_TEST … … 5377 5708 #endif 5378 5709 idSkipZeroes(G); 5379 delete next_weight; */ 5710 delete next_weight; 5711 */ 5380 5712 break; 5381 5713 } … … 5404 5736 #endif 5405 5737 Print("\n//** Mwalk: Groebner Walk took %d steps.\n", nstep); 5738 si_opt_1 = save1; //set original options 5406 5739 return(result); 5407 5740 } … … 5417 5750 si_opt_1 &= (~Sy_bit(OPT_REDTAIL)); // not tail reductions 5418 5751 } 5752 5419 5753 Set_Error(FALSE); 5420 5754 Overflow_Error = FALSE; … … 5427 5761 #endif 5428 5762 nstep=0; 5429 int i, polylength,nwalk;5763 int i,nwalk;//polylength; 5430 5764 int nV = currRing->N; 5765 5766 //check that weight radius is valid 5767 if(weight_rad < 0) 5768 { 5769 Werror("Invalid radius.\n"); 5770 return NULL; 5771 } 5772 5773 //check that perturbation degree is valid 5774 if(pert_deg > nV || pert_deg < 1) 5775 { 5776 Werror("Invalid perturbation degree.\n"); 5777 return NULL; 5778 } 5431 5779 5432 5780 ideal Gomega, M, F,FF, Gomega1, Gomega2, M1; … … 5435 5783 ring baseRing = currRing; 5436 5784 ring XXRing = currRing; 5785 intvec* iv_M; 5437 5786 intvec* ivNull = new intvec(nV); 5438 5787 intvec* curr_weight = new intvec(nV); … … 5445 5794 (*target_weight)[i] = (*target_M)[i]; 5446 5795 } 5796 5447 5797 #ifndef BUCHBERGER_ALG 5448 5798 intvec* hilb_func; … … 5456 5806 #endif 5457 5807 rComplete(currRing); 5808 5809 if(target_M->length() == nV) 5810 { 5811 targetRing = VMrDefault(target_weight); // define the target ring 5812 } 5813 else 5814 { 5815 targetRing = VMatrDefault(target_M); 5816 } 5817 if(orig_M->length() == nV) 5818 { 5819 newRing = VMrDefault(curr_weight); // define a new ring with ordering "(a(curr_weight),lp) 5820 } 5821 else 5822 { 5823 newRing = VMatrDefault(orig_M); 5824 } 5825 rChangeCurrRing(newRing); 5458 5826 #ifdef TIME_TEST 5459 5827 to = clock(); 5460 5828 #endif 5461 5462 if(target_M->length() == nV)5463 {5464 // define the target ring5465 targetRing = VMrDefault(target_weight);5466 }5467 else5468 {5469 targetRing = VMatrDefault(target_M);5470 }5471 if(orig_M->length() == nV)5472 {5473 newRing = VMrDefault(curr_weight); // define a new ring with ordering "(a(curr_weight),lp)5474 }5475 else5476 {5477 newRing = VMatrDefault(orig_M);5478 }5479 rChangeCurrRing(newRing);5480 5829 ideal G = MstdCC(idrMoveR(Go,baseRing,currRing)); 5830 #ifdef TIME_TEST 5831 tostd = clock()-to; 5832 #endif 5481 5833 baseRing = currRing; 5482 #ifdef TIME_TEST5483 tostd = clock()-to;5484 #endif5485 5486 5834 nwalk = 0; 5487 5835 … … 5498 5846 nwalk ++; 5499 5847 nstep ++; 5848 #ifdef CHECK_IDEAL_MWALK 5500 5849 if(printout > 1) 5501 5850 { 5502 5851 idString(Gomega,"//** Mrwalk: Gomega"); 5503 5852 } 5853 #endif 5504 5854 if(reduction == 0) 5505 5855 { … … 5510 5860 G = idCopy(FF); 5511 5861 idDelete(&FF); 5512 5513 5862 goto NEXT_VECTOR; 5514 5863 } … … 5553 5902 to = clock(); 5554 5903 #endif 5555 #ifndef 5904 #ifndef BUCHBERGER_ALG 5556 5905 M=kStd(Gomega1,NULL,isHomog,NULL,hilb_func,0,NULL,curr_weight); 5557 5906 delete hilb_func; … … 5563 5912 #endif 5564 5913 idSkipZeroes(M); 5565 //#ifdef CHECK_IDEAL_MWALK5914 #ifdef CHECK_IDEAL_MWALK 5566 5915 if(printout > 2) 5567 5916 { 5568 5917 idString(M, "//** Mrwalk: M"); 5569 5918 } 5570 //#endif5919 #endif 5571 5920 //change the ring to baseRing 5572 5921 rChangeCurrRing(baseRing); … … 5584 5933 tlift = tlift + clock() - to; 5585 5934 #endif 5586 //#ifdef CHECK_IDEAL_MWALK5935 #ifdef CHECK_IDEAL_MWALK 5587 5936 if(printout > 2) 5588 5937 { 5589 idString(F, 5590 } 5591 //#endif5938 idString(F,"//** Mrwalk: F"); 5939 } 5940 #endif 5592 5941 idDelete(&Gomega2); 5593 5942 idDelete(&M1); … … 5601 5950 #endif 5602 5951 idSkipZeroes(G); 5603 //#ifdef CHECK_IDEAL_MWALK5952 #ifdef CHECK_IDEAL_MWALK 5604 5953 if(printout > 2) 5605 5954 { 5606 idString(G, 5607 } 5608 //#endif5955 idString(G,"//** Mrwalk: G"); 5956 } 5957 #endif 5609 5958 5610 5959 rChangeCurrRing(targetRing); … … 5639 5988 #endif 5640 5989 5641 // polylength= 1 if there is a polynomial in Gomega with at least 3 monomials and 0 otherwise5642 polylength = lengthpoly(Gomega);5643 if( polylength> 0)5990 //lengthpoly(Gomega) = 1 if there is a polynomial in Gomega with at least 3 monomials and 0 otherwise 5991 //polylength = lengthpoly(Gomega); 5992 if(lengthpoly(Gomega) > 0) 5644 5993 { 5645 5994 //there is a polynomial in Gomega with at least 3 monomials, 5646 5995 //low-dimensional facet of the cone 5647 5996 delete next_weight; 5648 #ifdef TIME_TEST 5649 to = clock(); 5650 #endif 5651 next_weight = MWalkRandomNextWeight(G, curr_weight, target_weight, weight_rad, pert_deg); 5652 #ifdef TIME_TEST 5653 tnw = tnw + clock() - to; 5997 if(target_M->length() == nV) 5998 { 5999 iv_M = MivMatrixOrder(curr_weight); 6000 } 6001 else 6002 { 6003 iv_M = MivMatrixOrderRefine(curr_weight,target_M); 6004 } 6005 #ifdef TIME_TEST 6006 to = clock(); 6007 #endif 6008 next_weight = MWalkRandomNextWeight(G, iv_M, target_weight, weight_rad, pert_deg); 6009 #ifdef TIME_TEST 6010 tnw = tnw + clock() - to; 5654 6011 #endif 5655 6012 idDelete(&Gomega); 5656 6013 #ifdef TIME_TEST 5657 to = clock();6014 to = clock(); 5658 6015 #endif 5659 6016 Gomega = MwalkInitialForm(G, next_weight); 5660 6017 #ifdef TIME_TEST 5661 tif = tif + clock()-to; //time for computing initial form ideal 5662 #endif 6018 tif = tif + clock()-to; //time for computing initial form ideal 6019 #endif 6020 delete iv_M; 5663 6021 } 5664 6022 … … 5671 6029 } 5672 6030 5673 //#ifdef PRINT_VECTORS6031 #ifdef PRINT_VECTORS 5674 6032 if(printout > 0) 5675 6033 { 5676 6034 MivString(curr_weight, target_weight, next_weight); 5677 6035 } 5678 //#endif6036 #endif 5679 6037 5680 6038 for(i=nV-1; i>=0; i--) … … 5688 6046 ideal result = idrMoveR(G,baseRing,currRing); 5689 6047 idDelete(&G); 5690 si_opt_1 = save1; //set original options, e. g. option(RedSB)5691 6048 delete ivNull; 5692 6049 #ifndef BUCHBERGER_ALG … … 5699 6056 //Print("\n// Overflow_Error? (%d)\n", Overflow_Error); 5700 6057 #endif 6058 si_opt_1 = save1; //set original options 5701 6059 return(result); 5702 6060 } 5703 5704 //unused5705 #if 05706 ideal Mwalk_tst(ideal Go, intvec* curr_weight, intvec* target_weight)5707 {5708 //clock_t tinput=clock();5709 //idString(Go,"Ginp");5710 int i, nV = currRing->N;5711 int nwalk=0, endwalks=0;5712 5713 ideal Gomega, M, F, Gomega1, Gomega2, M1, F1, G;5714 // ideal G1; ring endRing;5715 ring newRing, oldRing;5716 intvec* ivNull = new intvec(nV);5717 ring XXRing = currRing;5718 5719 intvec* tmp_weight = new intvec(nV);5720 for(i=nV-1; i>=0; i--)5721 {5722 (*tmp_weight)[i] = (*curr_weight)[i];5723 }5724 /* the monomial ordering of this current ring would be "dp" */5725 G = MstdCC(Go);5726 #ifndef BUCHBERGER_ALG5727 intvec* hilb_func;5728 #endif5729 /* to avoid (1,0,...,0) as the target vector */5730 intvec* last_omega = new intvec(nV);5731 for(i=nV-1; i>0; i--)5732 (*last_omega)[i] = 1;5733 (*last_omega)[0] = 10000;5734 5735 while(1)5736 {5737 nwalk ++;5738 //Print("\n// Entering the %d-th step:", nwalk);5739 //Print("\n// ring r[%d] = %s;", nwalk, rString(currRing));5740 idString(G,"G");5741 /* compute an initial form ideal of <G> w.r.t. "curr_vector" */5742 Gomega = MwalkInitialForm(G, curr_weight);5743 //ivString(curr_weight, "omega");5744 idString(Gomega,"Gw");5745 5746 #ifndef BUCHBERGER_ALG5747 if(isNolVector(curr_weight) == 0)5748 hilb_func = hFirstSeries(Gomega,NULL,NULL,curr_weight,currRing);5749 else5750 hilb_func = hFirstSeries(Gomega,NULL,NULL,last_omega,currRing);5751 #endif // BUCHBERGER_ALG5752 5753 5754 oldRing = currRing;5755 5756 /* define a new ring that its ordering is "(a(curr_weight),lp) */5757 VMrDefault(curr_weight);5758 newRing = currRing;5759 5760 Gomega1 = idrMoveR(Gomega, oldRing,currRing);5761 5762 /* compute a reduced Groebner basis of <Gomega> w.r.t. "newRing" */5763 #ifdef BUCHBERGER_ALG5764 M = MstdhomCC(Gomega1);5765 #else5766 M=kStd(Gomega1,NULL,isHomog,NULL,hilb_func,0,NULL,curr_weight);5767 delete hilb_func;5768 #endif // BUCHBERGER_ALG5769 5770 idString(M,"M");5771 5772 /* change the ring to oldRing */5773 rChangeCurrRing(oldRing);5774 M1 = idrMoveR(M, newRing,currRing);5775 Gomega2 = idrMoveR(Gomega1, newRing,currRing);5776 5777 /* compute a representation of the generators of submod (M)5778 with respect to those of mod (Gomega).5779 Gomega is a reduced Groebner basis w.r.t. the current ring */5780 F = MLifttwoIdeal(Gomega2, M1, G);5781 idDelete(&M1);5782 idDelete(&Gomega2);5783 idDelete(&G);5784 idString(F,"F");5785 5786 /* change the ring to newRing */5787 rChangeCurrRing(newRing);5788 F1 = idrMoveR(F, oldRing,currRing);5789 5790 /* reduce the Groebner basis <G> w.r.t. new ring */5791 G = kInterRedCC(F1, NULL);5792 //idSkipZeroes(G);//done by kInterRed5793 idDelete(&F1);5794 idString(G,"G");5795 if(endwalks == 1)5796 break;5797 5798 /* compute a next weight vector */5799 intvec* next_weight = MkInterRedNextWeight(curr_weight,target_weight,G);5800 #ifdef PRINT_VECTORS5801 MivString(curr_weight, target_weight, next_weight);5802 #endif5803 5804 if(MivComp(next_weight, ivNull) == 1)5805 {5806 delete next_weight;5807 break;5808 }5809 if(MivComp(next_weight, target_weight) == 1)5810 endwalks = 1;5811 5812 for(i=nV-1; i>=0; i--)5813 (*tmp_weight)[i] = (*curr_weight)[i];5814 5815 /* 06.11.01 to free the memory: NOT Changed!!*/5816 for(i=nV-1; i>=0; i--)5817 (*curr_weight)[i] = (*next_weight)[i];5818 delete next_weight;5819 }5820 rChangeCurrRing(XXRing);5821 G = idrMoveR(G, newRing,currRing);5822 5823 delete tmp_weight;5824 delete ivNull;5825 PrintLn();5826 return(G);5827 }5828 #endif5829 6061 5830 6062 /**************************************************************/ … … 5840 6072 2) the changed perturbation walk algorithm with a decreased degree. 5841 6073 */ 5842 // use kStd, if nP = 0, else call LastGB6074 // if nP = 0 use kStd, else call LastGB 5843 6075 ideal Mpwalk(ideal Go, int op_deg, int tp_deg,intvec* curr_weight, 5844 6076 intvec* target_weight, int nP, int reduction, int printout) … … 5848 6080 { 5849 6081 si_opt_1 &= (~Sy_bit(OPT_REDSB)); // no reduced Groebner basis 5850 //si_opt_1 &= (~Sy_bit(OPT_REDTAIL)); // not tail reductions 5851 //si_opt_1|=(Sy_bit(OPT_REDTAIL)|Sy_bit(OPT_REDSB)); 6082 si_opt_1 &= (~Sy_bit(OPT_REDTAIL)); // not tail reductions 5852 6083 } 5853 6084 Set_Error(FALSE ); 5854 6085 Overflow_Error = FALSE; 5855 6086 //Print("// pSetm_Error = (%d)", ErrorCheck()); 5856 6087 #ifdef TIME_TEST 5857 6088 clock_t tinput, tostd, tif=0, tstd=0, tlift=0, tred=0, tnw=0; 5858 6089 xtextra=0; … … 5861 6092 5862 6093 clock_t tim; 5863 6094 #endif 5864 6095 nstep = 0; 5865 6096 int i, ntwC=1, ntestw=1, nV = currRing->N; 6097 6098 //check that perturbation degree is valid 6099 if(op_deg < 1 || tp_deg < 1 || op_deg > nV || tp_deg > nV) 6100 { 6101 Werror("Invalid perturbation degree.\n"); 6102 return NULL; 6103 } 6104 5866 6105 BOOLEAN endwalks = FALSE; 5867 5868 6106 ideal Gomega, M, F, FF, G, Gomega1, Gomega2, M1,F1,Eresult,ssG; 5869 6107 ring newRing, oldRing, TargetRing; … … 5887 6125 5888 6126 ring XXRing = currRing; 5889 6127 #ifdef TIME_TEST 5890 6128 to = clock(); 6129 #endif 5891 6130 // perturbs the original vector 5892 6131 if(MivComp(curr_weight, iv_dp) == 1) //rOrdStr(currRing) := "dp" 5893 6132 { 5894 6133 G = MstdCC(Go); 6134 #ifdef TIME_TEST 5895 6135 tostd = clock()-to; 6136 #endif 5896 6137 if(op_deg != 1){ 5897 6138 iv_M_dp = MivMatrixOrderdp(nV); … … 5910 6151 G = idrMoveR(Go, XXRing,currRing); 5911 6152 G = MstdCC(G); 6153 #ifdef TIME_TEST 5912 6154 tostd = clock()-to; 6155 #endif 5913 6156 if(op_deg != 1){ 5914 6157 iv_M_dp = MivMatrixOrder(curr_weight); … … 5946 6189 G = idrMoveR(ssG, TargetRing,currRing); 5947 6190 } 5948 if(printout > 0) 5949 { 5950 Print("\n//** Mpwalk: Perturbation Walk of degree (%d,%d):",op_deg,tp_deg); 5951 ivString(curr_weight, "//** Mpwalk: new current weight"); 5952 ivString(target_weight, "//** Mpwalk: new target weight"); 5953 } 6191 if(printout > 0) 6192 { 6193 Print("\n//** Mpwalk: Perturbation Walk of degree (%d,%d):",op_deg,tp_deg); 6194 #ifdef PRINT_VECTORS 6195 ivString(curr_weight, "//** Mpwalk: new current weight"); 6196 ivString(target_weight, "//** Mpwalk: new target weight"); 6197 #endif 6198 } 5954 6199 while(1) 5955 6200 { 5956 6201 nstep ++; 6202 #ifdef TIME_TEST 5957 6203 to = clock(); 6204 #endif 5958 6205 // compute an initial form ideal of <G> w.r.t. the weight vector 5959 6206 // "curr_weight" 5960 6207 Gomega = MwalkInitialForm(G, curr_weight); 5961 //#ifdef CHECK_IDEAL_MWALK 6208 #ifdef TIME_TEST 6209 tif = tif + clock()-to; 6210 #endif 6211 #ifdef CHECK_IDEAL_MWALK 5962 6212 if(printout > 1) 5963 6213 { 5964 6214 idString(Gomega,"//** Mpwalk: Gomega"); 5965 6215 } 5966 //#endif6216 #endif 5967 6217 if(reduction == 0 && nstep > 1) 5968 6218 { … … 5978 6228 idDelete(&FF); 5979 6229 next_weight = MwalkNextWeightCC(curr_weight,target_weight,G); 6230 #ifdef PRINT_VECTORS 5980 6231 if(printout > 0) 5981 6232 { 5982 6233 MivString(curr_weight, target_weight, next_weight); 5983 6234 } 6235 #endif 5984 6236 } 5985 6237 else … … 5992 6244 } 5993 6245 Gomega = MwalkInitialForm(G, curr_weight); 6246 #ifdef CHECK_IDEAL_MWALK 5994 6247 if(printout > 1) 5995 6248 { 5996 6249 idString(Gomega,"//** Mpwalk: Gomega"); 5997 6250 } 6251 #endif 5998 6252 } 5999 6253 */ … … 6012 6266 { 6013 6267 Print("\n// ring r%d = %s;\n", nstep, rString(currRing)); 6268 /* 6014 6269 idElements(G, "G"); 6015 6270 headidString(G, "G"); 6016 } 6017 #endif 6018 6019 tif = tif + clock()-to; 6271 */ 6272 } 6273 #endif 6020 6274 6021 6275 #ifndef BUCHBERGER_ALG … … 6041 6295 { 6042 6296 Print("\n// ring r%d = %s;\n", nstep, rString(currRing)); 6297 /* 6043 6298 idElements(Gomega1, "Gw"); 6044 6299 headidString(Gomega1, "headGw"); 6300 */ 6045 6301 PrintS("\n// compute a rGB of Gw:\n"); 6046 6047 6302 #ifndef BUCHBERGER_ALG 6048 6303 ivString(hilb_func, "w"); … … 6050 6305 } 6051 6306 #endif 6052 6307 #ifdef TIME_TEST 6053 6308 tim = clock(); 6054 6309 to = clock(); 6310 #endif 6055 6311 // compute a reduced Groebner basis of <Gomega> w.r.t. "newRing" 6056 6312 #ifdef BUCHBERGER_ALG … … 6060 6316 delete hilb_func; 6061 6317 #endif 6062 //#ifdef CHECK_IDEAL_MWALK 6063 if(printout > 2) 6064 { 6065 idString(M,"//** Mpwalk: M"); 6066 } 6067 //#endif 6068 6069 if(endwalks == TRUE){ 6318 6319 if(endwalks == TRUE) 6320 { 6321 #ifdef TIME_TEST 6070 6322 xtstd = xtstd+clock()-to; 6323 #endif 6071 6324 #ifdef ENDWALKS 6072 6325 Print("\n// time for the last std(Gw) = %.2f sec\n", … … 6075 6328 } 6076 6329 else 6330 { 6331 #ifdef TIME_TEST 6077 6332 tstd=tstd+clock()-to; 6078 6333 #endif 6334 } 6335 #ifdef CHECK_IDEAL_MWALK 6336 if(printout > 2) 6337 { 6338 idString(M,"//** Mpwalk: M"); 6339 } 6340 #endif 6079 6341 // change the ring to oldRing 6080 6342 rChangeCurrRing(oldRing); 6081 6343 M1 = idrMoveR(M, newRing,currRing); 6082 6344 Gomega2 = idrMoveR(Gomega1, newRing,currRing); 6083 6345 #ifdef TIME_TEST 6084 6346 to=clock(); 6347 #endif 6085 6348 /* compute a representation of the generators of submod (M) 6086 6349 with respect to those of mod (Gomega). 6087 6350 Gomega is a reduced Groebner basis w.r.t. the current ring */ 6088 6351 F = MLifttwoIdeal(Gomega2, M1, G); 6352 #ifdef TIME_TEST 6089 6353 if(endwalks == FALSE) 6090 6354 tlift = tlift+clock()-to; 6091 6355 else 6092 6356 xtlift=clock()-to; 6093 6094 //#ifdef CHECK_IDEAL_MWALK6357 #endif 6358 #ifdef CHECK_IDEAL_MWALK 6095 6359 if(printout > 2) 6096 6360 { 6097 6361 idString(F,"//** Mpwalk: F"); 6098 6362 } 6099 //#endif6363 #endif 6100 6364 6101 6365 idDelete(&M1); … … 6116 6380 PrintS("\n //** Mpwalk: reduce the Groebner basis.\n"); 6117 6381 } 6382 #ifdef TIME_TEST 6118 6383 to=clock(); 6384 #endif 6119 6385 G = kInterRedCC(F1, NULL); 6386 #ifdef TIME_TEST 6120 6387 if(endwalks == FALSE) 6121 6388 tred = tred+clock()-to; 6122 6389 else 6123 6390 xtred=clock()-to; 6391 #endif 6124 6392 idDelete(&F1); 6125 6393 } … … 6128 6396 6129 6397 NEXT_VECTOR: 6398 #ifdef TIME_TEST 6130 6399 to=clock(); 6400 #endif 6131 6401 // compute a next weight vector 6132 6402 next_weight = MkInterRedNextWeight(curr_weight,target_weight, G); 6403 #ifdef TIME_TEST 6133 6404 tnw=tnw+clock()-to; 6134 //#ifdef PRINT_VECTORS 6405 #endif 6406 #ifdef PRINT_VECTORS 6135 6407 if(printout > 0) 6136 6408 { 6137 6409 MivString(curr_weight, target_weight, next_weight); 6138 6410 } 6139 //#endif6411 #endif 6140 6412 6141 6413 if(Overflow_Error == TRUE) … … 6179 6451 TargetRing=currRing; 6180 6452 F1 = idrMoveR(G, newRing,currRing); 6181 #ifdef CHECK_IDEAL 6453 /* 6454 #ifdef CHECK_IDEAL_MWALK 6182 6455 headidString(G, "G"); 6183 6456 #endif 6184 6457 */ 6185 6458 6186 6459 // check whether the pertubed target vector stays in the correct cone … … 6265 6538 Overflow_Error = FALSE; 6266 6539 //Print("// pSetm_Error = (%d)", ErrorCheck()); 6267 6540 #ifdef TIME_TEST 6268 6541 clock_t tinput, tostd, tif=0, tstd=0, tlift=0, tred=0, tnw=0; 6269 6542 xtextra=0; … … 6272 6545 6273 6546 clock_t tim; 6274 6547 #endif 6275 6548 nstep = 0; 6276 int i, ntwC=1, ntestw=1, polylength, nV = currRing->N; 6549 int i, ntwC=1, ntestw=1, nV = currRing->N; //polylength 6550 6551 //check that weight radius is valid 6552 if(weight_rad < 0) 6553 { 6554 Werror("Invalid radius.\n"); 6555 return NULL; 6556 } 6557 6558 //check that perturbation degree is valid 6559 if(op_deg < 1 || tp_deg < 1 || op_deg > nV || tp_deg > nV) 6560 { 6561 Werror("Invalid perturbation degree.\n"); 6562 return NULL; 6563 } 6564 6277 6565 BOOLEAN endwalks = FALSE; 6278 6566 6279 6567 ideal Gomega, M, F, FF, G, Gomega1, Gomega2, M1,F1,Eresult,ssG; 6280 6568 ring newRing, oldRing, TargetRing; 6569 intvec* iv_M; 6281 6570 intvec* iv_M_dp; 6282 6571 intvec* iv_M_lp; … … 6306 6595 ring XXRing = currRing; 6307 6596 6308 to = clock();6309 6597 // perturbs the original vector 6310 6598 if(orig_M->length() == nV) … … 6312 6600 if(MivComp(curr_weight, iv_dp) == 1) //rOrdStr(currRing) := "dp" 6313 6601 { 6602 #ifdef TIME_TEST 6603 to = clock(); 6604 #endif 6314 6605 G = MstdCC(Go); 6606 #ifdef TIME_TEST 6315 6607 tostd = clock()-to; 6608 #endif 6316 6609 if(op_deg != 1) 6317 6610 { … … 6329 6622 6330 6623 G = idrMoveR(Go, XXRing,currRing); 6624 #ifdef TIME_TEST 6625 to = clock(); 6626 #endif 6331 6627 G = MstdCC(G); 6628 #ifdef TIME_TEST 6332 6629 tostd = clock()-to; 6630 #endif 6333 6631 if(op_deg != 1) 6334 6632 { … … 6342 6640 rChangeCurrRing(VMatrDefault(orig_M)); 6343 6641 G = idrMoveR(Go, XXRing,currRing); 6642 #ifdef TIME_TEST 6643 to = clock(); 6644 #endif 6344 6645 G = MstdCC(G); 6646 #ifdef TIME_TEST 6345 6647 tostd = clock()-to; 6648 #endif 6346 6649 if(op_deg != 1) 6347 6650 { … … 6411 6714 { 6412 6715 nstep ++; 6413 to = clock(); 6414 // compute an initial form ideal of G w.r.t. the weight vector "curr_weight" 6415 /* Gomega = MwalkInitialForm(G, curr_weight); 6416 polylength = lengthpoly(Gomega); 6417 */ 6418 //#ifdef CHECK_IDEAL_MWALK 6716 #ifdef CHECK_IDEAL_MWALK 6419 6717 if(printout > 1) 6420 6718 { 6421 6719 idString(Gomega,"//** Mprwalk: Gomega"); 6422 6720 } 6423 //#endif6721 #endif 6424 6722 6425 6723 if(reduction == 0 && nstep > 1) … … 6436 6734 idDelete(&FF); 6437 6735 next_weight = MwalkNextWeightCC(curr_weight,target_weight,G); 6736 #ifdef PRINT_VECTORS 6438 6737 if(printout > 0) 6439 6738 { 6440 6739 MivString(curr_weight, target_weight, next_weight); 6441 6740 } 6741 #endif 6442 6742 } 6443 6743 else … … 6450 6750 } 6451 6751 Gomega = MwalkInitialForm(G, curr_weight); 6752 #ifdef CHECK_IDEAL_MWALK 6452 6753 if(printout > 1) 6453 6754 { 6454 6755 idString(Gomega,"//** Mprwalk: Gomega"); 6455 6756 } 6757 #endif 6456 6758 } 6457 6759 */ … … 6470 6772 { 6471 6773 Print("\n// ring r%d = %s;\n", nstep, rString(currRing)); 6774 /* 6472 6775 idElements(G, "G"); 6473 6776 headidString(G, "G"); 6777 */ 6474 6778 } 6475 6779 #endif … … 6502 6806 { 6503 6807 Print("\n// ring r%d = %s;\n", nstep, rString(currRing)); 6808 /* 6504 6809 idElements(Gomega1, "Gw"); 6505 6810 headidString(Gomega1, "headGw"); 6811 */ 6506 6812 PrintS("\n// compute a rGB of Gw:\n"); 6507 6813 … … 6511 6817 } 6512 6818 #endif 6513 6819 #ifdef TIME_TEST 6514 6820 tim = clock(); 6515 6821 to = clock(); 6822 #endif 6516 6823 // compute a reduced Groebner basis of <Gomega> w.r.t. "newRing" 6517 6824 #ifdef BUCHBERGER_ALG … … 6521 6828 delete hilb_func; 6522 6829 #endif 6523 //#ifdef CHECK_IDEAL_MWALK6830 #ifdef CHECK_IDEAL_MWALK 6524 6831 if(printout > 2) 6525 6832 { 6526 6833 idString(M,"//** Mprwalk: M"); 6527 6834 } 6528 //#endif6529 6835 #endif 6836 #ifdef TIME_TEST 6530 6837 if(endwalks == TRUE) 6531 6838 { … … 6538 6845 else 6539 6846 tstd=tstd+clock()-to; 6540 6847 #endif 6541 6848 /* change the ring to oldRing */ 6542 6849 rChangeCurrRing(oldRing); 6543 6850 M1 = idrMoveR(M, newRing,currRing); 6544 6851 Gomega2 = idrMoveR(Gomega1, newRing,currRing); 6545 6852 #ifdef TIME_TEST 6546 6853 to=clock(); 6854 #endif 6547 6855 /* compute a representation of the generators of submod (M) 6548 6856 with respect to those of mod (Gomega). 6549 6857 Gomega is a reduced Groebner basis w.r.t. the current ring */ 6550 6858 F = MLifttwoIdeal(Gomega2, M1, G); 6859 #ifdef TIME_TEST 6551 6860 if(endwalks == FALSE) 6552 6861 tlift = tlift+clock()-to; 6553 6862 else 6554 6863 xtlift=clock()-to; 6555 6556 //#ifdef CHECK_IDEAL_MWALK6864 #endif 6865 #ifdef CHECK_IDEAL_MWALK 6557 6866 if(printout > 2) 6558 6867 { 6559 6868 idString(F,"//** Mprwalk: F"); 6560 6869 } 6561 //#endif6870 #endif 6562 6871 6563 6872 idDelete(&M1); … … 6578 6887 PrintS("\n //** Mprwalk: reduce the Groebner basis.\n"); 6579 6888 } 6889 #ifdef TIME_TEST 6580 6890 to=clock(); 6891 #endif 6581 6892 G = kInterRedCC(F1, NULL); 6893 #ifdef TIME_TEST 6582 6894 if(endwalks == FALSE) 6583 6895 tred = tred+clock()-to; 6584 6896 else 6585 6897 xtred=clock()-to; 6898 #endif 6586 6899 idDelete(&F1); 6587 6900 } … … 6590 6903 break; 6591 6904 6592 Print("\n Next weight");6593 6905 NEXT_VECTOR: 6594 6906 #ifdef TIME_TEST … … 6603 6915 to = clock(); 6604 6916 #endif 6605 Gomega = MwalkInitialForm(G, next_weight); // compute an initial form ideal of <G> w.r.t. "curr_vector" 6917 // compute an initial form ideal of <G> w.r.t. "next_vector" 6918 Gomega = MwalkInitialForm(G, next_weight); 6606 6919 #ifdef TIME_TEST 6607 6920 tif = tif + clock()-to; //time for computing initial form ideal 6608 6921 #endif 6609 6922 6610 //polylength = 1 if there is a polynomial in Gomega with at least 3 monomials and 0 otherwise 6611 polylength = lengthpoly(Gomega); 6612 if(polylength > 0) 6613 { 6614 Print("\n there is a polynomial in Gomega with at least 3 monomials"); 6615 //low-dimensional facet of the cone 6923 //lengthpoly(Gomega) = 1 if there is a polynomial in Gomega with at least 3 monomials and 0 otherwise 6924 if(lengthpoly(Gomega) > 0) 6925 { 6926 Print("\n there is a polynomial in Gomega with at least 3 monomials,\n"); 6927 // low-dimensional facet of the cone 6616 6928 delete next_weight; 6929 if(target_M->length() == nV) 6930 { 6931 iv_M = MivMatrixOrder(curr_weight); 6932 } 6933 else 6934 { 6935 iv_M = MivMatrixOrderRefine(curr_weight,target_M); 6936 } 6617 6937 #ifdef TIME_TEST 6618 6938 to = clock(); 6619 6939 #endif 6620 next_weight = MWalkRandomNextWeight(G, curr_weight, target_weight, weight_rad, op_deg);6940 next_weight = MWalkRandomNextWeight(G, iv_M, target_weight, weight_rad, op_deg); 6621 6941 #ifdef TIME_TEST 6622 6942 tnw = tnw + clock() - to; … … 6630 6950 tif = tif + clock()-to; //time for computing initial form ideal 6631 6951 #endif 6952 Print("delete\n"); 6953 delete iv_M; 6632 6954 } 6633 6955 … … 6646 6968 tnw=tnw+clock()-to; 6647 6969 */ 6648 //#ifdef PRINT_VECTORS6970 #ifdef PRINT_VECTORS 6649 6971 if(printout > 0) 6650 6972 { 6651 6973 MivString(curr_weight, target_weight, next_weight); 6652 6974 } 6653 //#endif6975 #endif 6654 6976 6655 6977 if(Overflow_Error == TRUE) … … 6674 6996 6675 6997 delete next_weight; 6676 }// while6998 }// end of while-loop 6677 6999 6678 7000 if(tp_deg != 1) … … 6698 7020 TargetRing=currRing; 6699 7021 F1 = idrMoveR(G, newRing,currRing); 6700 #ifdef CHECK_IDEAL6701 headidString(G, "G");6702 #endif6703 7022 6704 7023 // check whether the pertubed target vector stays in the correct cone … … 6711 7030 if(ntestw != 1 && printout > 2) 6712 7031 { 7032 #ifdef PRINT_VECTORS 6713 7033 ivString(pert_target_vector, "tau"); 7034 #endif 6714 7035 PrintS("\n// **Mprwalk: perturbed target vector doesn't stay in cone."); 6715 7036 Print("\n// ring r%d = %s;\n", nstep, rString(currRing)); … … 6717 7038 } 6718 7039 // LastGB is "better" than the kStd subroutine 7040 #ifdef TIME_TEST 6719 7041 to=clock(); 7042 #endif 6720 7043 ideal eF1; 6721 7044 if(nP == 0 || tp_deg == 1 || MivSame(orig_target, exivlp) != 1 || target_M->length() != nV) … … 6739 7062 F2=NULL; 6740 7063 } 7064 #ifdef TIME_TEST 6741 7065 xtextra=clock()-to; 7066 #endif 6742 7067 ring exTargetRing = currRing; 6743 7068 … … 6817 7142 intvec* hilb_func; 6818 7143 #endif 6819 //intvec* extXtau;7144 //intvec* extXtau; 6820 7145 intvec* next_vect; 6821 7146 intvec* omega2 = new intvec(nV); … … 6866 7191 nwalks ++; 6867 7192 NEXT_VECTOR_FRACTAL: 7193 #ifdef TIME_TEST 6868 7194 to=clock(); 7195 #endif 6869 7196 // determine the next border 6870 7197 next_vect = MkInterRedNextWeight(omega,omega2,G); 7198 #ifdef TIME_TEST 6871 7199 xtnw=xtnw+clock()-to; 6872 7200 #endif 6873 7201 oRing = currRing; 6874 7202 … … 6928 7256 6929 7257 delete next_vect; 7258 #ifdef TIME_TEST 6930 7259 to=clock(); 6931 7260 #endif 6932 7261 // to avoid the value of Overflow_Error that occur in Mfpertvector 6933 7262 Overflow_Error = FALSE; 6934 7263 next_vect = MkInterRedNextWeight(omega,omega2,G); 7264 #ifdef TIME_TEST 6935 7265 xtnw=xtnw+clock()-to; 7266 #endif 6936 7267 }// end of (if MivComp(next_vect, omega2) == 1) 6937 7268 6938 //#ifdef PRINT_VECTORS7269 #ifdef PRINT_VECTORS 6939 7270 if(printout > 0) 6940 7271 { 6941 7272 MivString(omega, omega2, next_vect); 6942 7273 } 6943 //#endif7274 #endif 6944 7275 6945 7276 // check whether the the computed vector is in the correct cone. … … 6969 7300 rString(currRing)); 6970 7301 } 7302 #ifdef TIME_TEST 6971 7303 to=clock(); 7304 #endif 6972 7305 Gt = idrMoveR(G, oRing,currRing); 6973 7306 G1 = MstdCC(Gt); 7307 #ifdef TIME_TEST 6974 7308 xtextra=xtextra+clock()-to; 7309 #endif 6975 7310 Gt = NULL; 6976 7311 … … 6988 7323 return (G1); 6989 7324 } 6990 6991 7325 6992 7326 /* If the perturbed target vector stays in the correct cone, … … 7076 7410 rString(currRing)); 7077 7411 } 7412 #ifdef TIME_TEST 7078 7413 to=clock(); 7414 #endif 7079 7415 G = MstdCC(Gt); 7416 #ifdef TIME_TEST 7080 7417 xtextra=xtextra+clock()-to; 7081 7418 #endif 7082 7419 oRing = currRing; 7083 7420 … … 7140 7477 } 7141 7478 delete next_vect; 7142 7479 #ifdef TIME_TEST 7143 7480 to=clock(); 7481 #endif 7144 7482 // Take the initial form of <G> w.r.t. omega 7145 7483 Gomega = MwalkInitialForm(G, omega); 7484 #ifdef TIME_TEST 7146 7485 xtif=xtif+clock()-to; 7486 #endif 7487 #ifdef CHECK_IDEAL_MWALK 7147 7488 if(printout > 1) 7148 7489 { 7149 7490 idString(Gomega,"//** rec_fractal_call: Gomega"); 7150 7491 } 7151 7492 #endif 7152 7493 if(reduction == 0) 7153 7494 { … … 7194 7535 Print("\n//** rec_fractal_call: Maximal recursion depth.\n"); 7195 7536 } 7196 7537 #ifdef TIME_TEST 7197 7538 to=clock(); 7539 #endif 7198 7540 #ifdef BUCHBERGER_ALG 7199 7541 Gresult = MstdhomCC(Gomega1); … … 7202 7544 delete hilb_func; 7203 7545 #endif 7546 #ifdef TIME_TEST 7204 7547 xtstd=xtstd+clock()-to; 7548 #endif 7205 7549 } 7206 7550 else … … 7210 7554 Gresult = rec_fractal_call(idCopy(Gomega1),nlev+1,omega,reduction,printout); 7211 7555 } 7556 #ifdef CHECK_IDEAL_MWALK 7212 7557 if(printout > 2) 7213 7558 { 7214 7559 idString(Gresult,"//** rec_fractal_call: M"); 7215 7560 } 7561 #endif 7216 7562 //convert a Groebner basis from a ring to another ring 7217 7563 new_ring = currRing; … … 7220 7566 Gresult1 = idrMoveR(Gresult, new_ring,currRing); 7221 7567 Gomega2 = idrMoveR(Gomega1, new_ring,currRing); 7222 7568 #ifdef TIME_TEST 7223 7569 to=clock(); 7570 #endif 7224 7571 // Lifting process 7225 7572 F = MLifttwoIdeal(Gomega2, Gresult1, G); 7573 #ifdef TIME_TEST 7226 7574 xtlift=xtlift+clock()-to; 7575 #endif 7576 #ifdef CHECK_IDEAL_MWALK 7227 7577 if(printout > 2) 7228 7578 { 7229 7579 idString(F,"//** rec_fractal_call: F"); 7230 7580 } 7581 #endif 7231 7582 idDelete(&Gresult1); 7232 7583 idDelete(&Gomega2); … … 7234 7585 7235 7586 rChangeCurrRing(new_ring); 7236 //F1 = idrMoveR(F, oRing,currRing);7237 7587 G = idrMoveR(F,oRing,currRing); 7588 /* 7589 F1 = idrMoveR(F, oRing,currRing); 7590 #ifdef TIME_TEST 7238 7591 to=clock(); 7592 #endif 7239 7593 // Interreduce G 7240 // G = kInterRedCC(F1, NULL); 7594 G = kInterRedCC(F1, NULL); 7595 #ifdef TIME_TEST 7241 7596 xtred=xtred+clock()-to; 7242 //idDelete(&F1); 7597 #endif 7598 idDelete(&F1); 7599 */ 7243 7600 } 7244 7601 } … … 7253 7610 //Print("\n\n// Entering the %d-th recursion:", nlev); 7254 7611 7255 int i, polylength, nV = currRing->N;7612 int nwalks = 0,i,nV=currRing->N;//polylength 7256 7613 ring new_ring, testring; 7257 7614 //ring extoRing; 7258 7615 ideal Gomega, Gomega1, Gomega2, F, FF, F1, Gresult, Gresult1, G1, Gt; 7259 int nwalks = 0;7260 7616 intvec* Mwlp; 7261 7617 #ifndef BUCHBERGER_ALG … … 7264 7620 // intvec* extXtau; 7265 7621 intvec* next_vect; 7622 intvec* iv_M; 7266 7623 intvec* omega2 = new intvec(nV); 7267 7624 intvec* omtmp = new intvec(nV); … … 7313 7670 nwalks ++; 7314 7671 NEXT_VECTOR_FRACTAL: 7672 #ifdef TIME_TEST 7315 7673 to=clock(); 7674 #endif 7316 7675 /* determine the next border */ 7317 7676 next_vect = MkInterRedNextWeight(omega,omega2,G); 7318 if(polylength > 0 && G->m[0] != NULL) 7677 #ifdef TIME_TEST 7678 xtnw=xtnw+clock()-to; 7679 #endif 7680 if(lengthpoly(MwalkInitialForm(G, next_vect)) > 0 && G->m[0] != NULL) 7319 7681 { 7320 7682 if(printout > 0) 7321 7683 { 7322 PrintS("\n**// rec_r_fractal_call: there is a polynomial in Gomega with at least 3 monomials , low-dimensional facet of the cone.\n");7684 PrintS("\n**// rec_r_fractal_call: there is a polynomial in Gomega with at least 3 monomials.\n"); 7323 7685 } 7324 7686 delete next_vect; 7325 next_vect = MWalkRandomNextWeight(G,omega,omega2,weight_rad,nlev); 7687 iv_M = MivMatrixOrder(omega); 7688 #ifdef TIME_TEST 7689 to=clock(); 7690 #endif 7691 next_vect = MWalkRandomNextWeight(G,iv_M,omega2,weight_rad,nlev); 7692 #ifdef TIME_TEST 7693 xtnw=xtnw+clock()-to; 7694 #endif 7326 7695 if(isNegNolVector(next_vect) == 1) 7327 7696 { 7328 7697 delete next_vect; 7698 #ifdef TIME_TEST 7699 to=clock(); 7700 #endif 7329 7701 next_vect = MkInterRedNextWeight(omega,omega2,G); 7330 } 7331 } 7332 xtnw=xtnw+clock()-to; 7333 7702 #ifdef TIME_TEST 7703 xtnw=xtnw+clock()-to; 7704 #endif 7705 } 7706 } 7334 7707 oRing = currRing; 7335 7708 … … 7386 7759 delete Mwlp; 7387 7760 7388 for(i=nV-1; i>=0; i--) { 7761 for(i=nV-1; i>=0; i--) 7762 { 7389 7763 (*omega2)[i] = (*Xtau)[nV+i]; 7390 7764 (*omega)[i] = (*Xsigma)[nV+i]; … … 7392 7766 7393 7767 delete next_vect; 7768 7769 //to avoid the value of Overflow_Error that occur in Mfpertvector 7770 Overflow_Error = FALSE; 7771 #ifdef TIME_TEST 7394 7772 to=clock(); 7395 7396 /* 7397 to avoid the value of Overflow_Error that occur in 7398 Mfpertvector 7399 */ 7400 Overflow_Error = FALSE; 7401 7773 #endif 7402 7774 next_vect = MkInterRedNextWeight(omega,omega2,G); 7403 if(G->m[0] != NULL && polylength > 0) 7775 #ifdef TIME_TEST 7776 xtnw=xtnw+clock()-to; 7777 #endif 7778 if(lengthpoly(MwalkInitialForm(G, next_vect)) > 0 && G->m[0] != NULL) 7404 7779 { 7405 7406 PrintS("//** rec_r_fractal_call: there is a polynomial in Gomega with at least 3 monomials, low-dimensional facet of the cone"); 7407 7780 // there is a polynomial in Gomega with at least 3 monomials 7781 iv_M = MivMatrixOrder(omega); 7408 7782 delete next_vect; 7409 next_vect = MWalkRandomNextWeight(G,omega,omega2,weight_rad,nlev); 7783 #ifdef TIME_TEST 7784 to=clock(); 7785 #endif 7786 next_vect = MWalkRandomNextWeight(G,iv_M,omega2,weight_rad,nlev); 7787 #ifdef TIME_TEST 7788 xtnw=xtnw+clock()-to; 7789 #endif 7790 delete iv_M; 7410 7791 if(isNegNolVector(next_vect) == 1) 7411 7792 { 7412 7793 delete next_vect; 7794 #ifdef TIME_TEST 7795 to=clock(); 7796 #endif 7413 7797 next_vect = MkInterRedNextWeight(omega,omega2,G); 7798 #ifdef TIME_TEST 7799 xtnw=xtnw+clock()-to; 7800 #endif 7414 7801 } 7415 7802 } 7416 xtnw=xtnw+clock()-to; 7417 } 7418 //#ifdef PRINT_VECTORS 7803 } 7804 #ifdef PRINT_VECTORS 7419 7805 if(printout > 0) 7420 7806 { 7421 7807 MivString(omega, omega2, next_vect); 7422 7808 } 7423 //#endif7424 7425 /*check whether the the computed vector is in the correct cone7809 #endif 7810 7811 /* check whether the the computed vector is in the correct cone 7426 7812 If no, the reduced GB of an omega-homogeneous ideal will be 7427 computed by Buchberger algorithm and stop this recursion step */7428 //if(test_w_in_ConeCC(G, next_vect) != 1) //e.g. Example s7, cyc6 7429 if(Overflow_Error == TRUE || test_w_in_ConeCC(G,next_vect) != 1) 7813 computed by Buchberger algorithm and stop this recursion step 7814 */ 7815 if(Overflow_Error == TRUE || test_w_in_ConeCC(G,next_vect) != 1)//e.g. Example s7, cyc6 7430 7816 { 7431 7817 delete next_vect; … … 7455 7841 rString(currRing)); 7456 7842 } 7843 Gt = idrMoveR(G, oRing,currRing); 7844 #ifdef TIME_TEST 7457 7845 to=clock(); 7458 Gt = idrMoveR(G, oRing,currRing); 7846 #endif 7459 7847 G1 = MstdCC(Gt); 7848 #ifdef TIME_TEST 7460 7849 xtextra=xtextra+clock()-to; 7850 #endif 7461 7851 Gt = NULL; 7462 7852 … … 7519 7909 7520 7910 #ifndef MSTDCC_FRACTAL 7521 //ivString(Xtau, "old Xtau"); 7911 #ifdef PRINT_VECTORS 7912 if(printout > 0) 7913 { 7914 ivString(Xtau, "old Xtau"); 7915 } 7916 #endif 7522 7917 intvec* Xtautmp; 7523 7918 if(ivtarget->length() == nV) … … 7543 7938 Xtau = Xtautmp; 7544 7939 Xtautmp = NULL; 7545 //ivString(Xtau, "new Xtau"); 7940 #ifdef PRINT_VECTORS 7941 if(printout > 0) 7942 { 7943 ivString(Xtau, "new Xtau"); 7944 } 7945 #endif 7546 7946 7547 7947 for(i=nV-1; i>=0; i--) … … 7561 7961 rString(currRing)); 7562 7962 } 7963 #ifdef TIME_TEST 7563 7964 to=clock(); 7965 #endif 7564 7966 G = MstdCC(Gt); 7967 #ifdef TIME_TEST 7565 7968 xtextra=xtextra+clock()-to; 7566 7969 #endif 7567 7970 oRing = currRing; 7568 7971 … … 7626 8029 } 7627 8030 delete next_vect; 7628 8031 #ifdef TIME_TEST 7629 8032 to=clock(); 8033 #endif 7630 8034 // Take the initial form of <G> w.r.t. omega 7631 8035 Gomega = MwalkInitialForm(G, omega); 8036 #ifdef TIME_TEST 7632 8037 xtif=xtif+clock()-to; 8038 #endif 7633 8039 //polylength = 1 if there is a polynomial in Gomega with at least 3 monomials and 0 otherwise 7634 polylength = lengthpoly(Gomega); 8040 //polylength = lengthpoly(Gomega); 8041 #ifdef CHECK_IDEAL_MWALK 7635 8042 if(printout > 1) 7636 8043 { 7637 8044 idString(Gomega,"//** rec_r_fractal_call: Gomega"); 7638 8045 } 7639 8046 #endif 7640 8047 if(reduction == 0) 7641 8048 { … … 7677 8084 if(nlev == Xnlev || lengthpoly(Gomega1) == 0) 7678 8085 { 8086 #ifdef TIME_TEST 7679 8087 to=clock(); 8088 #endif 7680 8089 #ifdef BUCHBERGER_ALG 7681 8090 Gresult = MstdhomCC(Gomega1); … … 7684 8093 delete hilb_func; 7685 8094 #endif 8095 #ifdef TIME_TEST 7686 8096 xtstd=xtstd+clock()-to; 8097 #endif 7687 8098 } 7688 8099 else … … 7692 8103 Gresult = rec_r_fractal_call(idCopy(Gomega1),nlev+1,omega,weight_rad,reduction,printout); 7693 8104 } 8105 #ifdef CHECK_IDEAL_MWALK 7694 8106 if(printout > 2) 7695 8107 { 7696 8108 idString(Gresult,"//** rec_r_fractal_call: M"); 7697 8109 } 8110 #endif 7698 8111 //convert a Groebner basis from a ring to another ring 7699 8112 new_ring = currRing; … … 7702 8115 Gresult1 = idrMoveR(Gresult, new_ring,currRing); 7703 8116 Gomega2 = idrMoveR(Gomega1, new_ring,currRing); 7704 8117 #ifdef TIME_TEST 7705 8118 to=clock(); 8119 #endif 7706 8120 // Lifting process 7707 8121 F = MLifttwoIdeal(Gomega2, Gresult1, G); 8122 #ifdef TIME_TEST 7708 8123 xtlift=xtlift+clock()-to; 7709 8124 #endif 8125 #ifdef CHECK_IDEAL_MWALK 7710 8126 if(printout > 2) 7711 8127 { 7712 8128 idString(F,"//** rec_r_fractal_call: F"); 7713 8129 } 7714 8130 #endif 7715 8131 idDelete(&Gresult1); 7716 8132 idDelete(&Gomega2); … … 7720 8136 //F1 = idrMoveR(F, oRing,currRing); 7721 8137 G = idrMoveR(F,oRing,currRing); 8138 /* 8139 #ifdef TIME_TEST 7722 8140 to=clock(); 8141 #endif 7723 8142 // Interreduce G 7724 //G = kInterRedCC(F1, NULL); 8143 G = kInterRedCC(F1, NULL); 8144 #ifdef TIME_TEST 7725 8145 xtred=xtred+clock()-to; 7726 //idDelete(&F1); 8146 #endif 8147 idDelete(&F1); 8148 */ 7727 8149 } 7728 8150 } … … 7755 8177 Xngleich = 0; 7756 8178 Xcall = 0; 8179 #ifdef TIME_TEST 7757 8180 xtif=0; xtstd=0; xtlift=0; xtred=0; xtnw=0; xtextra=0; 7758 8181 xftinput = clock(); 7759 8182 #endif 7760 8183 ring oldRing = currRing; 7761 8184 int i, nV = currRing->N; … … 7763 8186 Xivinput = ivtarget; 7764 8187 ngleich = 0; 8188 #ifdef TIME_TEST 7765 8189 to=clock(); 8190 #endif 7766 8191 ideal I = MstdCC(G); 7767 8192 G = NULL; 8193 #ifdef TIME_TEST 7768 8194 xftostd=clock()-to; 8195 #endif 7769 8196 Xsigma = ivstart; 7770 8197 … … 7861 8288 7862 8289 I = idrMoveR(I1,tRing,currRing); 8290 #ifdef TIME_TEST 7863 8291 to=clock(); 8292 #endif 7864 8293 ideal J = MstdCC(I); 7865 8294 idDelete(&I); 8295 #ifdef TIME_TEST 7866 8296 xftostd=xftostd+clock()-to; 7867 8297 #endif 7868 8298 ideal resF; 7869 8299 ring helpRing = currRing; … … 7907 8337 { 7908 8338 BITSET save1 = si_opt_1; // save current options 8339 //check that weight radius is valid 8340 if(weight_rad < 0) 8341 { 8342 Werror("Invalid radius.\n"); 8343 return NULL; 8344 } 7909 8345 if(reduction == 0) 7910 8346 { 7911 8347 si_opt_1 &= (~Sy_bit(OPT_REDSB)); // no reduced Groebner basis 7912 //si_opt_1 &= (~Sy_bit(OPT_REDTAIL)); // not tail reductions8348 si_opt_1 &= (~Sy_bit(OPT_REDTAIL)); // not tail reductions 7913 8349 } 7914 8350 Set_Error(FALSE); … … 7920 8356 Xngleich = 0; 7921 8357 Xcall = 0; 8358 #ifdef TIME_TEST 7922 8359 xtif=0; xtstd=0; xtlift=0; xtred=0; xtnw=0; xtextra=0; 7923 8360 xftinput = clock(); 7924 8361 #endif 7925 8362 ring oldRing = currRing; 7926 8363 int i, nV = currRing->N; … … 7928 8365 Xivinput = ivtarget; 7929 8366 ngleich = 0; 8367 #ifdef TIME_TEST 7930 8368 to=clock(); 8369 #endif 7931 8370 ideal I = MstdCC(G); 7932 8371 G = NULL; 8372 #ifdef TIME_TEST 7933 8373 xftostd=clock()-to; 8374 #endif 7934 8375 Xsigma = ivstart; 7935 8376 … … 8026 8467 8027 8468 I = idrMoveR(I1,tRing,currRing); 8469 #ifdef TIME_TEST 8028 8470 to=clock(); 8471 #endif 8029 8472 ideal J = MstdCC(I); 8030 8473 idDelete(&I); 8474 #ifdef TIME_TEST 8031 8475 xftostd=xftostd+clock()-to; 8032 8476 #endif 8033 8477 ideal resF; 8034 8478 ring helpRing = currRing; … … 8251 8695 goto BE_FINISH; 8252 8696 #endif 8253 8697 /* 8254 8698 #ifdef CHECK_IDEAL_MWALK 8255 8699 idElements(G, "G"); 8256 8700 //headidString(G, "G"); 8257 8701 #endif 8258 8702 */ 8259 8703 if(MivSame(target_tmp, iv_lp) == 1) 8260 8704 if (rParameter(currRing) != NULL) … … 8990 9434 clock_t tinput=clock(); 8991 9435 int i, nV = currRing->N; 9436 9437 //check that perturbation degree is valid 9438 if(tp_deg < 1 || tp_deg > nV) 9439 { 9440 Werror("Invalid perturbation degree.\n"); 9441 return NULL; 9442 } 9443 8992 9444 int nwalk=0, endwalks=0, ntestwinC=1; 8993 9445 int tp_deg_tmp = tp_deg; … … 9099 9551 { 9100 9552 Print("\n// it is %d-th step!!", nwalk); 9101 id Elements(Gomega1, "Gw");9553 idString(Gomega1, "Gw"); 9102 9554 PrintS("\n// compute a rGB of Gw:"); 9103 9555 }
Note: See TracChangeset
for help on using the changeset viewer.