Changeset 0cf2ccd in git for kernel/tgb.cc
- Timestamp:
- Feb 14, 2007, 2:55:47 PM (17 years ago)
- Branches:
- (u'spielwiese', '5b153614cbc72bfa198d75b1e9e33dab2645d9fe')
- Children:
- b3daab03c4cc77a25533037671204cd7b2ba146f
- Parents:
- 9a26cb0ddc5a028ba02a4f2087848c04e7ea6035
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
kernel/tgb.cc
r9a26cb r0cf2ccd 5 5 * Computer Algebra System SINGULAR * 6 6 ****************************************/ 7 /* $Id: tgb.cc,v 1.13 0 2007-02-13 14:58:52bricken Exp $ */7 /* $Id: tgb.cc,v 1.131 2007-02-14 13:55:47 bricken Exp $ */ 8 8 /* 9 9 * ABSTRACT: slimgb and F4 implementation … … 2085 2085 PrintS("StopGauss\n"); 2086 2086 } 2087 static void write_poly_to_row(number* row, poly h, poly*terms, int tn ){2087 static void write_poly_to_row(number* row, poly h, poly*terms, int tn, ring r){ 2088 2088 //poly* base=row; 2089 2089 while(h!=NULL){ 2090 2090 //Print("h:%i\n",h); 2091 number coef=p_GetCoeff(h, c->r);2091 number coef=p_GetCoeff(h,r); 2092 2092 poly* ptr_to_h=(poly*) bsearch(&h,terms,tn,sizeof(poly),terms_sort_crit); 2093 2093 assume(ptr_to_h!=NULL); … … 2098 2098 } 2099 2099 } 2100 static void linalg_step_modp(poly *p, poly* p_out, int&pn, poly* terms,int tn, slimgb_alg* c){ 2100 static poly row_to_poly(number* row, poly* terms, int tn, ring r){ 2101 poly h=NULL; 2102 int j; 2103 for(j=tn-1;j>=0;j--){ 2104 if (!(npIsZero(row[j]))){ 2105 poly t=terms[j]; 2106 t=p_LmInit(t,r); 2107 p_SetCoeff(t,row[j],r); 2108 pNext(t)=h; 2109 h=t; 2110 } 2111 2112 } 2113 return h; 2114 } 2115 static void linalg_step_modp(poly *p, poly* p_out, int& pn, poly* terms,int tn, slimgb_alg* c){ 2101 2116 static int export_n=0; 2102 2117 assume(terms[tn-1]!=NULL); … … 2113 2128 poly h=p[i]; 2114 2129 //int base=tn*i; 2115 write_poly_to_row(number_array+tn*i,h,terms,tn );2130 write_poly_to_row(number_array+tn*i,h,terms,tn,c->r); 2116 2131 2117 2132 } … … 2129 2144 int base=tn*i; 2130 2145 number* row=number_array+base; 2131 for(j=tn-1;j>=0;j--){ 2132 if (!(npIsZero(row[j]))){ 2133 poly t=terms[j]; 2134 t=p_LmInit(t,c->r); 2135 p_SetCoeff(t,row[j],c->r); 2136 pNext(t)=h; 2137 h=t; 2138 } 2139 2140 } 2146 h=row_to_poly(row,terms,tn,c->r); 2147 2141 2148 if (h!=NULL){ 2142 2149 p_out[p_pos++]=h; … … 2200 2207 branches=NULL; 2201 2208 branches_len=0; 2209 2202 2210 } 2203 2211 NoroCacheNode* setNode(int branch, NoroCacheNode* node){ … … 2247 2255 public: 2248 2256 number* array; 2249 int start;2257 int begin; 2250 2258 int end; 2259 DenseRow(){ 2260 array=NULL; 2261 } 2262 ~DenseRow(){ 2263 omfree(array); 2264 } 2251 2265 }; 2266 2252 2267 class DataNoroCacheNode:public NoroCacheNode{ 2253 2268 public: … … 2261 2276 value_poly=p; 2262 2277 row=NULL; 2278 term_index=-1; 2279 } 2280 ~DataNoroCacheNode(){ 2281 //p_Delete(&value_poly,currRing); 2282 if (row) delete row; 2263 2283 } 2264 2284 }; … … 2266 2286 public: 2267 2287 poly temp_term; 2288 void evaluatePlaceHolder(number* row,std::vector<NoroPlaceHolder>& place_holders); 2268 2289 void collectIrreducibleMonomials( std::vector<DataNoroCacheNode*>& res); 2269 2290 void collectIrreducibleMonomials(int level, NoroCacheNode* node, std::vector<DataNoroCacheNode*>& res); 2291 void evaluateRows(); 2292 void evaluateRows(int level, NoroCacheNode* node); 2270 2293 #ifdef NORO_RED_ARRAY_RESERVER 2271 2294 int reserved; … … 2276 2299 //assume(impl.find(p_Copy(term,currRing))==impl.end()); 2277 2300 //assume(len==pLength(nf)); 2301 assume(npIsOne(p_GetCoeff(term,currRing))); 2278 2302 if (term==nf){ 2279 2303 term=p_Copy(term,currRing); 2280 2304 2281 2305 ressources.push_back(term); 2306 nIrreducibleMonomials++; 2282 2307 return treeInsertBackLink(term); 2283 2308 2284 2309 } else{ 2285 2310 2286 2311 if (nf){ 2287 2312 //nf=p_Copy(nf,currRing); 2313 assume(p_LmCmp(nf,term,currRing)==-1); 2288 2314 ressources.push_back(nf); 2289 2315 } … … 2295 2321 } 2296 2322 DataNoroCacheNode* insertAndTransferOwnerShip(poly t, ring r){ 2323 nIrreducibleMonomials++; 2324 ressources.push_back(t); 2297 2325 return treeInsertBackLink(t); 2298 2326 } … … 2302 2330 #ifdef NORO_RED_ARRAY_RESERVER 2303 2331 reserved=0; 2332 2304 2333 recursionPolyBuffer=(poly*)omalloc(1000000*sizeof(poly)); 2305 2334 #endif 2335 nIrreducibleMonomials=0; 2306 2336 temp_term=pOne(); 2307 2337 } … … 2327 2357 #endif 2328 2358 } 2359 2360 int nIrreducibleMonomials; 2329 2361 protected: 2330 2362 DataNoroCacheNode* treeInsert(poly term,poly nf,int len){ … … 2344 2376 parent=parent->getOrInsertBranch(p_GetExp(term,i,currRing)); 2345 2377 } 2346 return (DataNoroCacheNode*) parent->setNode(p_GetExp(term,nvars,currRing),new DataNoroCacheNode( NULL,backLinkCode));2378 return (DataNoroCacheNode*) parent->setNode(p_GetExp(term,nvars,currRing),new DataNoroCacheNode(term,backLinkCode)); 2347 2379 } 2348 2380 … … 2353 2385 //cache_map impl; 2354 2386 NoroCacheNode root; 2387 2355 2388 }; 2389 void NoroCache::evaluateRows(){ 2390 //after that can evaluate placeholders 2391 int i; 2392 for(i=0;i<root.branches_len;i++){ 2393 evaluateRows(1,root.branches[i]); 2394 } 2395 } 2396 void NoroCache::evaluateRows(int level, NoroCacheNode* node){ 2397 assume(level>=0); 2398 if (node==NULL) return; 2399 if (level<pVariables){ 2400 int i,sum; 2401 for(i=0;i<node->branches_len;i++){ 2402 evaluateRows(level+1,node->branches[i]); 2403 } 2404 } else { 2405 DataNoroCacheNode* dn=(DataNoroCacheNode*) node; 2406 if (dn->value_len!=backLinkCode){ 2407 dn->row=new DenseRow(); 2408 dn->row->array=(number*) omalloc(nIrreducibleMonomials*sizeof(number)); 2409 dn->row->begin=0; 2410 dn->row->end=nIrreducibleMonomials; 2411 number* row= dn->row->array; 2412 int i; 2413 for(i=0;i<nIrreducibleMonomials;i++){ 2414 row[i]=npInit(0); 2415 } 2416 poly p=dn->value_poly; 2417 i=0; 2418 while(p){ 2419 DataNoroCacheNode* ref=getCacheReference(p); 2420 2421 int idx=ref->term_index; 2422 if (i==0){ 2423 dn->row->begin=idx; 2424 } 2425 dn->row->end=idx+1; 2426 i++; 2427 assume(idx>=0); 2428 //Print("idx:%d\n",idx); 2429 row[idx]=p_GetCoeff(p,currRing); 2430 pIter(p); 2431 } 2432 //p_Delete(&dn->value_poly,currRing); 2433 } 2434 } 2435 } 2436 void NoroCache::evaluatePlaceHolder(number* row,std::vector<NoroPlaceHolder>& place_holders){ 2437 int i; 2438 int s=place_holders.size(); 2439 for(i=0;i<s;i++){ 2440 DataNoroCacheNode* ref=place_holders[i].ref; 2441 number coef=place_holders[i].coef; 2442 if (ref->value_len==backLinkCode){ 2443 row[ref->term_index]=npAdd(row[ref->term_index],coef); 2444 } else { 2445 DenseRow* ref_row=ref->row; 2446 int j; 2447 if (!(npIsOne(coef))){ 2448 for(j=ref_row->begin;j<ref_row->end;j++){ 2449 row[j]=npAdd(row[j],npMult(coef,ref_row->array[j])); 2450 } 2451 } else{ 2452 for(j=ref_row->begin;j<ref_row->end;j++) 2453 row[j]=npAdd(row[j],ref_row->array[j]); 2454 } 2455 } 2456 } 2457 2458 } 2356 2459 void NoroCache::collectIrreducibleMonomials( std::vector<DataNoroCacheNode*>& res){ 2357 2460 int i; … … 2417 2520 MonRedRes noro_red_mon(poly t, BOOLEAN force_unique, NoroCache* cache,slimgb_alg* c){ 2418 2521 MonRedRes res_holder; 2419 2522 2420 2523 //wrp(t); 2421 2524 res_holder.changed=TRUE; … … 2427 2530 res_holder.len=1; 2428 2531 2429 2430 } 2532 2533 } 2534 res_holder.coef=p_GetCoeff(t,c->r); 2431 2535 res_holder.p=ref->value_poly; 2432 2536 res_holder.ref=ref; 2433 2537 res_holder.onlyBorrowed=TRUE; 2434 2538 res_holder.changed=TRUE; 2539 p_Delete(&t,c->r); 2435 2540 return res_holder; 2436 2541 } 2437 2542 } else{ 2438 BOOLEAN succ;2439 poly cache_lookup=cache->lookup(t,succ, res_holder.len);//don't own this yet2440 if (succ){2441 if (cache_lookup==t){2543 BOOLEAN succ; 2544 poly cache_lookup=cache->lookup(t,succ, res_holder.len);//don't own this yet 2545 if (succ){ 2546 if (cache_lookup==t){ 2442 2547 //know they are equal 2443 2548 //res_holder.len=1; 2444 2549 2550 res_holder.changed=FALSE; 2551 res_holder.p=t; 2552 res_holder.coef=npInit(1); 2553 2554 res_holder.onlyBorrowed=FALSE; 2555 return res_holder; 2556 } 2557 //poly t_new=ppMult_nn(cache_lookup,pGetCoeff(t)); 2558 res_holder.coef=p_GetCoeff(t,c->r); 2559 p_Delete(&t,c->r); 2560 //res_holder.len=1; 2561 //res_holder.changed=TRUE; 2562 res_holder.p=cache_lookup; 2563 2564 res_holder.onlyBorrowed=TRUE; 2565 return res_holder; 2566 //return t_new; 2567 } 2568 } 2569 2570 unsigned long sev=p_GetShortExpVector(t,currRing); 2571 int i=kFindDivisibleByInS_easy(c->strat,t,sev); 2572 if (i>=0){ 2573 number coef_bak=p_GetCoeff(t,c->r); 2574 2575 p_SetCoeff(t,npInit(1),c->r); 2576 assume(npIsOne(p_GetCoeff(c->strat->S[i],c->r))); 2577 number coefstrat=p_GetCoeff(c->strat->S[i],c->r); 2578 2579 //poly t_copy_mon=p_Copy(t,c->r); 2580 poly exp_diff=cache->temp_term; 2581 p_ExpVectorDiff(exp_diff,t,c->strat->S[i],c->r); 2582 p_SetCoeff(exp_diff,npNeg(npInvers(coefstrat)),c->r); 2583 p_Setm(exp_diff,c->r); 2584 assume(c->strat->S[i]!=NULL); 2585 //poly t_to_del=t; 2586 poly res; 2587 res=pp_Mult_mm(pNext(c->strat->S[i]),exp_diff,c->r); 2588 2589 res_holder.len=c->strat->lenS[i]-1; 2590 res=noro_red_non_unique(res,res_holder.len,cache,c); 2591 2592 DataNoroCacheNode* ref=cache->insert(t,res,res_holder.len); 2593 p_Delete(&t,c->r); 2594 //p_Delete(&t_copy_mon,c->r); 2595 //res=pMult_nn(res,coef_bak); 2596 2597 res_holder.changed=TRUE; 2598 res_holder.p=res; 2599 res_holder.coef=coef_bak; 2600 res_holder.onlyBorrowed=TRUE; 2601 res_holder.ref=ref; 2602 return res_holder; 2603 2604 } else { 2605 number coef_bak=p_GetCoeff(t,c->r); 2606 number one=npInit(1); 2607 p_SetCoeff(t,one,c->r); 2608 res_holder.len=1; 2609 if (!(force_unique)){ 2610 res_holder.ref=cache->insert(t,t,res_holder.len); 2611 p_SetCoeff(t,coef_bak,c->r); 2612 //return t; 2613 2614 //we need distinction 2445 2615 res_holder.changed=FALSE; 2446 2616 res_holder.p=t; 2617 2447 2618 res_holder.coef=npInit(1); 2448 2619 res_holder.onlyBorrowed=FALSE; 2449 2620 return res_holder; 2450 } 2451 //poly t_new=ppMult_nn(cache_lookup,pGetCoeff(t)); 2452 res_holder.coef=p_GetCoeff(t,c->r); 2453 p_Delete(&t,c->r); 2454 //res_holder.len=1; 2455 //res_holder.changed=TRUE; 2456 res_holder.p=cache_lookup; 2457 2458 res_holder.onlyBorrowed=TRUE; 2459 return res_holder; 2460 //return t_new; 2461 } 2462 } 2463 2464 unsigned long sev=p_GetShortExpVector(t,currRing); 2465 int i=kFindDivisibleByInS_easy(c->strat,t,sev); 2466 if (i>=0){ 2467 number coef_bak=p_GetCoeff(t,c->r); 2468 2469 p_SetCoeff(t,npInit(1),c->r); 2470 assume(npIsOne(p_GetCoeff(c->strat->S[i],c->r))); 2471 number coefstrat=p_GetCoeff(c->strat->S[i],c->r); 2472 2473 //poly t_copy_mon=p_Copy(t,c->r); 2474 poly exp_diff=cache->temp_term; 2475 p_ExpVectorDiff(exp_diff,t,c->strat->S[i],c->r); 2476 p_SetCoeff(exp_diff,npNeg(npInvers(coefstrat)),c->r); 2477 p_Setm(exp_diff,c->r); 2478 assume(c->strat->S[i]!=NULL); 2479 //poly t_to_del=t; 2480 poly res; 2481 res=pp_Mult_mm(pNext(c->strat->S[i]),exp_diff,c->r); 2482 2483 res_holder.len=c->strat->lenS[i]-1; 2484 res=noro_red_non_unique(res,res_holder.len,cache,c); 2485 2486 DataNoroCacheNode* ref=cache->insert(t,res,pLength(res)); 2487 p_Delete(&t,c->r); 2488 //p_Delete(&t_copy_mon,c->r); 2489 //res=pMult_nn(res,coef_bak); 2490 2491 res_holder.changed=TRUE; 2492 res_holder.p=res; 2621 } else { 2622 res_holder.ref=cache->insertAndTransferOwnerShip(t,c->r); 2493 2623 res_holder.coef=coef_bak; 2494 2624 res_holder.onlyBorrowed=TRUE; 2495 res_holder.ref=ref;2496 return res_holder;2497 2498 } else {2499 number coef_bak=p_GetCoeff(t,c->r);2500 number one=npInit(1);2501 p_SetCoeff(t,one,c->r);2502 res_holder.len=1;2503 if (!(force_unique)){2504 res_holder.ref=cache->insert(t,t,res_holder.len);2505 p_SetCoeff(t,coef_bak,c->r);2506 //return t;2507 2508 //we need distinction2509 2625 res_holder.changed=FALSE; 2510 2626 res_holder.p=t; 2511 2512 res_holder.coef=npInit(1);2513 res_holder.onlyBorrowed=FALSE;2514 2627 return res_holder; 2515 } else { 2516 res_holder.ref=cache->insertAndTransferOwnerShip(t,c->r); 2517 res_holder.coef=coef_bak; 2518 res_holder.onlyBorrowed=TRUE; 2519 res_holder.changed=FALSE; 2520 res_holder.p=t; 2521 return res_holder; 2522 } 2523 } 2524 2628 } 2629 } 2630 2525 2631 } 2526 2632 poly tree_add(poly* a,int begin, int end,ring r){ … … 2557 2663 pIter(p); 2558 2664 pNext(t)=NULL; 2559 2665 #ifndef NDEBUG 2666 number coef_debug=p_GetCoeff(t,currRing); 2667 #endif 2560 2668 MonRedRes red=noro_red_mon(t,FALSE,cache,c); 2561 2669 if ((!(red.changed))&&(!(red.onlyBorrowed))){ 2562 2670 unchanged_size++; 2671 assume(npIsOne(red.coef)); 2672 assume(p_GetCoeff(red.p,currRing)==coef_debug); 2563 2673 if (unchanged_head){ 2564 2674 pNext(unchanged_tail)=red.p; … … 2570 2680 } else{ 2571 2681 assume(red.len==pLength(red.p)); 2572 kBucket_Add_q(bucket,red.p,&red.len); 2682 if (red.onlyBorrowed){ 2683 t=pp_Mult_nn(red.p,red.coef,currRing); 2684 } else { 2685 t=p_Mult_nn(red.p,red.coef,currRing); 2686 } 2687 kBucket_Add_q(bucket,t,&red.len); 2573 2688 } 2574 2689 } … … 2601 2716 h.ref=red.ref; 2602 2717 h.coef=red.coef; 2718 assume(!((h.ref->value_poly==NULL) &&(h.ref->value_len!=0))); 2603 2719 if (h.ref->value_poly) 2604 2720 res.push_back(h); … … 2695 2811 } 2696 2812 void noro_step(poly*p,int &pn,slimgb_alg* c){ 2697 2813 //Print("Input rows %d\n",pn); 2698 2814 int j; 2699 2815 … … 2716 2832 cache.collectIrreducibleMonomials(irr_nodes); 2717 2833 //now can build up terms array 2834 //Print("historic irred Mon%d\n",cache.nIrreducibleMonomials); 2718 2835 int n=irr_nodes.size();//cache.countIrreducibleMonomials(); 2719 2836 cache.nIrreducibleMonomials=n; 2837 if (TEST_OPT_PROT) 2838 Print("Irred Mon:%d\n",n); 2720 2839 TermNoroDataNode* term_nodes=(TermNoroDataNode*) omalloc(n*sizeof(TermNoroDataNode)); 2721 2840 … … 2735 2854 terms[j]=term_nodes[j].t; 2736 2855 } 2737 2738 2856 if (TEST_OPT_PROT) 2857 Print("Evaluate Rows \n"); 2858 cache.evaluateRows(); 2859 number* number_array=(number*) omalloc(n*pn*sizeof(number)); 2860 number zero=npInit(0); 2861 if (TEST_OPT_PROT) 2862 Print("Evaluate Place Holders\n"); 2863 for(j=0;j<pn;j++){ 2864 int i; 2865 number* row=number_array+n*j; 2866 for(i=0;i<n;i++){ 2867 row[i]=zero; 2868 } 2869 cache.evaluatePlaceHolder(row,place_holders[j]); 2870 } 2871 if (TEST_OPT_PROT){ 2872 Print("Input rows %d\n",pn); 2873 } 2874 static int export_n=0; 2875 //export_mat(number_array,pn,n,"mat%i.py",++export_n); 2876 simplest_gauss_modp(number_array,pn,n); 2877 2878 int p_pos=0; 2879 for(j=0;j<pn;j++){ 2880 poly h=row_to_poly(number_array+j*n,terms,n,c->r); 2881 if(h!=NULL){ 2882 p[p_pos++]=h; 2883 } 2884 } 2885 pn=p_pos; 2739 2886 omfree(terms); 2740 2887 omfree(term_nodes); 2888 omfree(number_array); 2741 2889 //don't forget the rank 2742 2890 … … 2748 2896 //programming reasons 2749 2897 #ifdef USE_NORO 2750 const BOOLEAN use_noro=((!(c->nc))&&(rField_is_Zp(c->r))); 2898 //Print("module rank%d\n",c->S->rank); 2899 const BOOLEAN use_noro=((!(c->nc))&&(c->S->rank<=1)&&(rField_is_Zp(c->r))); 2751 2900 #else 2752 2901 const BOOLEAN use_noro=FALSE; … … 2833 2982 if (use_noro){ 2834 2983 int pn=i; 2984 if (pn==0) return; 2835 2985 noro_step(p,pn,c); 2836 2986 if (TEST_OPT_PROT){
Note: See TracChangeset
for help on using the changeset viewer.