[762407] | 1 | #include "config.h" |
---|
[599326] | 2 | #include <kernel/mod2.h> |
---|
| 3 | #include <kernel/F4.h> |
---|
| 4 | #include <kernel/tgb_internal.h> |
---|
| 5 | #include <kernel/tgbgauss.h> |
---|
[43bf98] | 6 | /**************************************** |
---|
| 7 | * Computer Algebra System SINGULAR * |
---|
| 8 | ****************************************/ |
---|
| 9 | /* |
---|
| 10 | * ABSTRACT: F4 implementation |
---|
| 11 | */ |
---|
[9cbb7a3] | 12 | static int posInPolys (poly* p, int pn, poly qe,slimgb_alg* c) |
---|
[b75d13] | 13 | { |
---|
| 14 | if(pn==0) return 0; |
---|
| 15 | |
---|
| 16 | int length=pn-1; |
---|
| 17 | int i; |
---|
| 18 | int an = 0; |
---|
| 19 | int en= length; |
---|
| 20 | |
---|
| 21 | if (pLmCmp(qe,p[en])==1) |
---|
| 22 | return length+1; |
---|
| 23 | |
---|
| 24 | while(1) |
---|
| 25 | { |
---|
| 26 | //if (an >= en-1) |
---|
| 27 | if(en-1<=an) |
---|
| 28 | { |
---|
| 29 | if (pLmCmp(p[an],qe)==1) return an; |
---|
| 30 | return en; |
---|
| 31 | } |
---|
| 32 | i=(an+en) / 2; |
---|
| 33 | if (pLmCmp(p[i],qe)==1) |
---|
| 34 | en=i; |
---|
| 35 | else an=i; |
---|
| 36 | } |
---|
| 37 | } |
---|
| 38 | |
---|
[9cbb7a3] | 39 | static tgb_sparse_matrix* build_sparse_matrix(poly* p,int p_index,poly* done, int done_index, slimgb_alg* c){ |
---|
[b75d13] | 40 | tgb_sparse_matrix* t=new tgb_sparse_matrix(p_index,done_index,c->r); |
---|
| 41 | int i, pos; |
---|
| 42 | // Print("\n 0:%s\n",pString(done[done_index-1])); |
---|
| 43 | //Print("\n 1:%s\n",pString(done[done_index-2])); |
---|
| 44 | // assume((!(pLmEqual(done[done_index-1],done[done_index-2])))); |
---|
| 45 | #ifdef TGB_DEGUG |
---|
| 46 | for(i=0;i<done_index;i++) |
---|
| 47 | { |
---|
| 48 | int j; |
---|
| 49 | for(j=0;j<i;j++) |
---|
| 50 | { |
---|
| 51 | assume((!(pLmEqual(done[i],done[j])))); |
---|
| 52 | } |
---|
| 53 | } |
---|
| 54 | #endif |
---|
| 55 | for(i=0;i<p_index;i++) |
---|
| 56 | { |
---|
| 57 | // Print("%i ter Eintrag:%s\n",i,pString(p[i])); |
---|
| 58 | mac_poly m=NULL; |
---|
| 59 | mac_poly* set_this=&m; |
---|
| 60 | poly p_i=p[i]; |
---|
| 61 | while(p_i) |
---|
| 62 | { |
---|
| 63 | |
---|
| 64 | int v=-1; |
---|
| 65 | pos=posInPolys (done, done_index, p_i,c); |
---|
| 66 | if((done_index>pos)&&(pLmEqual(p_i,done[pos]))) |
---|
| 67 | v=pos; |
---|
| 68 | if((pos>0) &&(pLmEqual(p_i,done[pos-1]))) |
---|
| 69 | v=pos-1; |
---|
| 70 | assume(v!=-1); |
---|
| 71 | //v is ascending ordered, we need descending order |
---|
| 72 | v=done_index-1-v; |
---|
| 73 | (*set_this)=new mac_poly_r(); |
---|
| 74 | (*set_this)->exp=v; |
---|
| 75 | |
---|
| 76 | (*set_this)->coef=nCopy(p_i->coef); |
---|
| 77 | set_this=&(*set_this)->next; |
---|
| 78 | p_i=p_i->next; |
---|
| 79 | |
---|
| 80 | } |
---|
| 81 | init_with_mac_poly(t,i,m); |
---|
| 82 | } |
---|
| 83 | return t; |
---|
| 84 | } |
---|
[9cbb7a3] | 85 | static tgb_matrix* build_matrix(poly* p,int p_index,poly* done, int done_index, slimgb_alg* c){ |
---|
[b75d13] | 86 | tgb_matrix* t=new tgb_matrix(p_index,done_index); |
---|
| 87 | int i, pos; |
---|
| 88 | // Print("\n 0:%s\n",pString(done[done_index-1])); |
---|
| 89 | //Print("\n 1:%s\n",pString(done[done_index-2])); |
---|
| 90 | assume((!(pLmEqual(done[done_index-1],done[done_index-2])))); |
---|
| 91 | #ifdef TGB_DEGUG |
---|
| 92 | for(i=0;i<done_index;i++) |
---|
| 93 | { |
---|
| 94 | int j; |
---|
| 95 | for(j=0;j<i;j++) |
---|
| 96 | { |
---|
| 97 | assume((!(pLmEqual(done[i],done[j])))); |
---|
| 98 | } |
---|
| 99 | } |
---|
| 100 | #endif |
---|
| 101 | for(i=0;i<p_index;i++) |
---|
| 102 | { |
---|
| 103 | // Print("%i ter Eintrag:%s\n",i,pString(p[i])); |
---|
| 104 | poly p_i=p[i]; |
---|
| 105 | while(p_i) |
---|
| 106 | { |
---|
| 107 | |
---|
| 108 | int v=-1; |
---|
| 109 | pos=posInPolys (done, done_index, p_i,c); |
---|
| 110 | if((done_index>pos)&&(pLmEqual(p_i,done[pos]))) |
---|
| 111 | v=pos; |
---|
| 112 | if((pos>0) &&(pLmEqual(p_i,done[pos-1]))) |
---|
| 113 | v=pos-1; |
---|
| 114 | assume(v!=-1); |
---|
| 115 | //v is ascending ordered, we need descending order |
---|
| 116 | v=done_index-1-v; |
---|
| 117 | number nt=t->get(i,v); |
---|
| 118 | nDelete(&nt); |
---|
| 119 | t->set(i,v,nCopy(p_i->coef)); |
---|
| 120 | p_i=p_i->next; |
---|
| 121 | } |
---|
| 122 | } |
---|
| 123 | return t; |
---|
| 124 | } |
---|
| 125 | |
---|
[9cbb7a3] | 126 | static int retranslate(poly* m,tgb_sparse_matrix* mat,poly* done, slimgb_alg* c){ |
---|
[b75d13] | 127 | int i; |
---|
| 128 | int m_index=0; |
---|
| 129 | for(i=0;i<mat->get_rows();i++) |
---|
| 130 | { |
---|
| 131 | if(mat->zero_row(i)) |
---|
| 132 | { |
---|
| 133 | mat->free_row(i); |
---|
| 134 | continue; |
---|
| 135 | } |
---|
| 136 | m[m_index++]= free_row_to_poly(mat, i, done, mat->get_columns()); |
---|
| 137 | |
---|
| 138 | } |
---|
| 139 | |
---|
| 140 | delete mat; |
---|
| 141 | return m_index; |
---|
| 142 | |
---|
| 143 | } |
---|
| 144 | //!returns m_index and destroys mat |
---|
[9cbb7a3] | 145 | static int retranslate(poly* m,tgb_matrix* mat,poly* done, slimgb_alg* c){ |
---|
[b75d13] | 146 | int i; |
---|
| 147 | int m_index=0; |
---|
| 148 | for(i=0;i<mat->get_rows();i++) |
---|
| 149 | { |
---|
| 150 | if(mat->zero_row(i)) |
---|
| 151 | { |
---|
| 152 | mat->free_row(i); |
---|
| 153 | continue; |
---|
| 154 | } |
---|
| 155 | |
---|
| 156 | m[m_index]=pInit(); |
---|
| 157 | int v=mat->min_col_not_zero_in_row(i); |
---|
| 158 | //v=done_index-1-pos; => pos=done_index-1-v=mat->get_columns()-1-v |
---|
| 159 | int pos=mat->get_columns()-1-v; |
---|
| 160 | int* ev=(int*) omalloc((c->r->N+1)*sizeof(int)); |
---|
| 161 | pGetExpV(done[pos],ev); |
---|
| 162 | pSetExpV(m[m_index],ev); |
---|
| 163 | omfree(ev); |
---|
| 164 | |
---|
| 165 | poly p=m[m_index]; |
---|
| 166 | pSetCoeff(p,mat->get(i,v)); |
---|
| 167 | while((v=mat->next_col_not_zero(i,v))!=mat->get_columns()) |
---|
| 168 | { |
---|
| 169 | poly pn=pInit(); |
---|
| 170 | |
---|
| 171 | //v=done_index-1-pos; => pos=done_index-1-v=mat->get_columns()-1-v |
---|
| 172 | pos=mat->get_columns()-1-v; |
---|
| 173 | ev=(int*) omalloc((c->r->N+1)*sizeof(int)); |
---|
| 174 | pGetExpV(done[pos],ev); |
---|
| 175 | pSetExpV(pn,ev); |
---|
| 176 | omfree(ev); |
---|
| 177 | pSetCoeff(pn,mat->get(i,v)); |
---|
| 178 | p->next=pn; |
---|
| 179 | p=pn; |
---|
| 180 | } |
---|
| 181 | p->next=NULL; |
---|
| 182 | mat->free_row(i, FALSE); |
---|
| 183 | m_index++; |
---|
| 184 | } |
---|
| 185 | |
---|
| 186 | delete mat; |
---|
| 187 | return m_index; |
---|
| 188 | |
---|
| 189 | } |
---|
| 190 | static int pLmCmp_func(const void* ap1, const void* ap2){ |
---|
| 191 | poly p1,p2; |
---|
| 192 | p1=*((poly*) ap1); |
---|
| 193 | p2=*((poly*)ap2); |
---|
| 194 | |
---|
| 195 | return pLmCmp(p1,p2); |
---|
| 196 | } |
---|
| 197 | static int monom_poly_crit(const void* ap1, const void* ap2){ |
---|
| 198 | monom_poly* p1; |
---|
| 199 | monom_poly* p2; |
---|
| 200 | p1=((monom_poly*) ap1); |
---|
| 201 | p2=((monom_poly*)ap2); |
---|
| 202 | if(((unsigned long) p1->f)>((unsigned long) p2->f)) return 1; |
---|
| 203 | if(((unsigned long) p1->f)<((unsigned long)p2->f)) return -1; |
---|
| 204 | |
---|
| 205 | return pLmCmp(p1->m,p2->m); |
---|
| 206 | |
---|
| 207 | } |
---|
[9cbb7a3] | 208 | static int posInMonomPolys (monom_poly* p, int pn, monom_poly & qe,slimgb_alg* c) |
---|
[b75d13] | 209 | { |
---|
| 210 | if(pn==0) return 0; |
---|
| 211 | |
---|
| 212 | int length=pn-1; |
---|
| 213 | int i; |
---|
| 214 | int an = 0; |
---|
| 215 | int en= length; |
---|
| 216 | |
---|
| 217 | if (monom_poly_crit(&qe,&p[en])==1) |
---|
| 218 | return length+1; |
---|
| 219 | |
---|
| 220 | while(1) |
---|
| 221 | { |
---|
| 222 | //if (an >= en-1) |
---|
| 223 | if(en-1<=an) |
---|
| 224 | { |
---|
| 225 | if (monom_poly_crit(&p[an],&qe)==1) return an; |
---|
| 226 | return en; |
---|
| 227 | } |
---|
| 228 | i=(an+en) / 2; |
---|
| 229 | if (monom_poly_crit(&p[i],&qe)==1) |
---|
| 230 | en=i; |
---|
| 231 | else an=i; |
---|
| 232 | } |
---|
| 233 | } |
---|
[9cbb7a3] | 234 | static void simplify(monom_poly& h, slimgb_alg* c){ |
---|
[b75d13] | 235 | mp_array_list* F=c->F; |
---|
| 236 | poly_array_list* F_minus=c->F_minus; |
---|
| 237 | while(F) |
---|
| 238 | { |
---|
| 239 | assume(F_minus!=NULL); |
---|
| 240 | int i; |
---|
| 241 | int posm=posInMonomPolys (F->mp, F->size, h,c); |
---|
| 242 | #ifdef TGB_DEBUG |
---|
| 243 | #endif |
---|
| 244 | //for(i=0;i<F->size;i++) |
---|
| 245 | for(i=si_min(posm,F->size-1);i>=0;i--) |
---|
| 246 | { |
---|
| 247 | // if((i>=posm)&&(F->mp[i].f!=h.f)) break; |
---|
| 248 | if((i<posm)&&(F->mp[i].f!=h.f)) break; |
---|
| 249 | if ((h.f==F->mp[i].f) &&(p_LmDivisibleBy(F->mp[i].m,h.m,c->r))) |
---|
| 250 | { |
---|
| 251 | |
---|
| 252 | // Print("found"); |
---|
| 253 | |
---|
| 254 | //according to the algorithm you should test (!(pIsConstant(F[i].m))) |
---|
| 255 | //but I think this is only because of bad formulation |
---|
| 256 | int j; |
---|
| 257 | |
---|
| 258 | poly lm=pLmInit(h.f); |
---|
| 259 | pSetCoeff(lm,nInit(1)); |
---|
| 260 | |
---|
| 261 | lm=pMult_mm(lm,F->mp[i].m); |
---|
| 262 | |
---|
| 263 | int pos; |
---|
| 264 | j=-1; |
---|
| 265 | pos=posInPolys (F_minus->p, F_minus->size, lm,c); |
---|
| 266 | if((F_minus->size>pos)&&(pLmEqual(lm,F_minus->p[pos]))) |
---|
| 267 | j=pos; |
---|
| 268 | if((pos>0) &&(pLmEqual(lm,F_minus->p[pos-1]))) |
---|
| 269 | j=pos-1; |
---|
| 270 | assume(j!=-1); |
---|
| 271 | // if(j==-1) Print("\n jAltert \n"); |
---|
| 272 | // for(j=0;j<F_minus->size;j++) |
---|
| 273 | // { |
---|
| 274 | // if (pLmEqual(F_minus->p[j],lm)) |
---|
| 275 | // break; |
---|
| 276 | // } |
---|
| 277 | // assume(j<F_minus->size); |
---|
| 278 | pDelete(&lm); |
---|
| 279 | if(j>=0) |
---|
| 280 | { |
---|
| 281 | pExpVectorSub(h.m,F->mp[i].m); |
---|
| 282 | h.f=F_minus->p[j]; |
---|
| 283 | } |
---|
| 284 | break; |
---|
| 285 | } |
---|
| 286 | } |
---|
| 287 | F=F->next; |
---|
| 288 | F_minus=F_minus->next; |
---|
| 289 | } |
---|
| 290 | assume(F_minus==NULL); |
---|
| 291 | } |
---|
| 292 | |
---|
[9cbb7a3] | 293 | void go_on_F4 (slimgb_alg* c){ |
---|
[b75d13] | 294 | //set limit of 1000 for multireductions, at the moment for |
---|
| 295 | //programming reasons |
---|
| 296 | int max_par; |
---|
| 297 | if (c->is_homog) |
---|
| 298 | max_par=PAR_N_F4; |
---|
| 299 | else |
---|
| 300 | max_par=200; |
---|
| 301 | int done_size=max_par*2; |
---|
| 302 | poly* done=(poly*) omalloc(done_size*sizeof(poly)); |
---|
| 303 | int done_index=0; //done_index must always be smaller than done_size |
---|
| 304 | int chosen_size=max_par*2; |
---|
| 305 | monom_poly* chosen=(monom_poly*) omalloc(chosen_size*sizeof(monom_poly)); |
---|
| 306 | int chosen_index=0; |
---|
| 307 | // monom_poly* vgl=(monom_poly*) omalloc(chosen_size*sizeof(monom_poly)); |
---|
| 308 | int i=0; |
---|
| 309 | c->average_length=0; |
---|
| 310 | for(i=0;i<c->n;i++){ |
---|
| 311 | c->average_length+=c->lengths[i]; |
---|
| 312 | } |
---|
| 313 | c->average_length=c->average_length/c->n; |
---|
| 314 | i=0; |
---|
| 315 | poly* p; |
---|
| 316 | int nfs=0; |
---|
| 317 | int curr_deg=-1; |
---|
| 318 | |
---|
| 319 | //choose pairs and preprocess symbolically |
---|
| 320 | while(chosen_index<max_par) |
---|
| 321 | { |
---|
| 322 | |
---|
| 323 | // sorted_pair_node* s=c->apairs[c->pair_top]; |
---|
| 324 | sorted_pair_node* s; |
---|
| 325 | if(c->pair_top>=0) |
---|
| 326 | s=top_pair(c);//here is actually chain criterium done |
---|
| 327 | else |
---|
| 328 | s=NULL; |
---|
| 329 | if (!s) break; |
---|
| 330 | nfs++; |
---|
| 331 | if(curr_deg>=0) |
---|
| 332 | { |
---|
| 333 | if (s->deg >curr_deg) break; |
---|
| 334 | } |
---|
| 335 | |
---|
| 336 | else curr_deg=s->deg; |
---|
| 337 | quick_pop_pair(c); |
---|
| 338 | if(s->i>=0) |
---|
| 339 | { |
---|
| 340 | //replace_pair(s->i,s->j,c); |
---|
| 341 | if(s->i==s->j) |
---|
| 342 | { |
---|
| 343 | free_sorted_pair_node(s,c->r); |
---|
| 344 | continue; |
---|
| 345 | } |
---|
| 346 | } |
---|
| 347 | assume(s->i>=0); |
---|
| 348 | monom_poly h; |
---|
| 349 | BOOLEAN only_free=FALSE; |
---|
| 350 | if(s->i>=0) |
---|
| 351 | { |
---|
| 352 | |
---|
| 353 | poly lcm=pOne(); |
---|
| 354 | |
---|
| 355 | pLcm(c->S->m[s->i], c->S->m[s->j], lcm); |
---|
| 356 | pSetm(lcm); |
---|
| 357 | assume(lcm!=NULL); |
---|
| 358 | poly factor1=pCopy(lcm); |
---|
| 359 | pExpVectorSub(factor1,c->S->m[s->i]); |
---|
| 360 | poly factor2=pCopy(lcm); |
---|
| 361 | pExpVectorSub(factor2,c->S->m[s->j]); |
---|
| 362 | |
---|
| 363 | if(done_index>=done_size) |
---|
| 364 | { |
---|
| 365 | done_size+=max_par; |
---|
| 366 | done=(poly*) omrealloc(done,done_size*sizeof(poly)); |
---|
| 367 | } |
---|
| 368 | done[done_index++]=lcm; |
---|
| 369 | if(chosen_index+1>=chosen_size) |
---|
| 370 | { |
---|
| 371 | //max_par must be greater equal 2 |
---|
| 372 | chosen_size+=si_max(max_par,2); |
---|
| 373 | chosen=(monom_poly*) omrealloc(chosen,chosen_size*sizeof(monom_poly)); |
---|
| 374 | } |
---|
| 375 | h.f=c->S->m[s->i]; |
---|
| 376 | h.m=factor1; |
---|
| 377 | chosen[chosen_index++]=h; |
---|
| 378 | h.f=c->S->m[s->j]; |
---|
| 379 | h.m=factor2; |
---|
| 380 | chosen[chosen_index++]=h; |
---|
| 381 | } |
---|
| 382 | else |
---|
| 383 | { |
---|
| 384 | if(chosen_index>=chosen_size) |
---|
| 385 | { |
---|
| 386 | chosen_size+=max_par; |
---|
| 387 | chosen=(monom_poly*) omrealloc(chosen,chosen_size*sizeof(monom_poly)); |
---|
| 388 | } |
---|
| 389 | h.f=s->lcm_of_lm; |
---|
| 390 | h.m=pOne(); |
---|
| 391 | //pSetm(h.m); done by pOne |
---|
| 392 | chosen[chosen_index++]=h; |
---|
| 393 | //must carefull remember to destroy such a h; |
---|
| 394 | poly_list_node* next=c->to_destroy; |
---|
| 395 | |
---|
| 396 | c->to_destroy=(poly_list_node*) omalloc(sizeof(poly_list_node)); |
---|
| 397 | c->to_destroy->p=h.f; |
---|
| 398 | c->to_destroy->next=next; |
---|
| 399 | only_free=TRUE; |
---|
| 400 | } |
---|
| 401 | |
---|
| 402 | if(s->i>=0) |
---|
| 403 | now_t_rep(s->j,s->i,c); |
---|
| 404 | |
---|
| 405 | if(!(only_free)) |
---|
| 406 | free_sorted_pair_node(s,c->r); |
---|
| 407 | else |
---|
| 408 | omfree(s); |
---|
| 409 | |
---|
| 410 | |
---|
| 411 | |
---|
| 412 | } |
---|
| 413 | c->normal_forms+=nfs; |
---|
| 414 | if (TEST_OPT_PROT) |
---|
| 415 | Print("M[%i, ",nfs); |
---|
| 416 | //next Step, simplify all pairs |
---|
| 417 | for(i=0;i<chosen_index;i++) |
---|
| 418 | { |
---|
| 419 | // PrintS("simplifying "); |
---|
| 420 | //Print("from %s",pString(chosen[i].m)); |
---|
| 421 | //Print(", %s", pString(chosen[i].f)); |
---|
| 422 | simplify(chosen[i],c); |
---|
| 423 | //Print(" to %s",pString(chosen[i].m)); |
---|
| 424 | //Print(", %s \n", pString(chosen[i].f)); |
---|
| 425 | } |
---|
| 426 | |
---|
| 427 | //next Step remove duplicate entries |
---|
| 428 | qsort(chosen,chosen_index,sizeof(monom_poly),monom_poly_crit); |
---|
| 429 | int pos=0; |
---|
| 430 | for(i=1;i<chosen_index;i++) |
---|
| 431 | { |
---|
| 432 | if(((chosen[i].f!=chosen[pos].f))||(!(pLmEqual(chosen[i].m,chosen[pos].m)))) |
---|
| 433 | chosen[++pos]=chosen[i]; |
---|
| 434 | else pDelete(&(chosen[i].m)); |
---|
| 435 | } |
---|
| 436 | if(chosen_index>0) |
---|
| 437 | chosen_index=pos+1; |
---|
| 438 | //next step process polys |
---|
| 439 | int p_size=2*chosen_index; |
---|
| 440 | int p_index=0; |
---|
| 441 | p=(poly*) omalloc(p_size*sizeof(poly)); |
---|
| 442 | for(p_index=0;p_index<chosen_index;p_index++) |
---|
| 443 | { |
---|
| 444 | p[p_index]=ppMult_mm(chosen[p_index].f,chosen[p_index].m); |
---|
| 445 | } |
---|
| 446 | |
---|
| 447 | qsort(done, done_index,sizeof(poly),pLmCmp_func); |
---|
| 448 | pos=0; |
---|
| 449 | //Print("Done_index:%i",done_index); |
---|
| 450 | if(done_index>0) |
---|
| 451 | { |
---|
| 452 | pTest(done[0]); |
---|
| 453 | } |
---|
| 454 | for(i=1;i<done_index;i++) |
---|
| 455 | { |
---|
| 456 | pTest(done[i]); |
---|
| 457 | if((!(pLmEqual(done[i],done[pos])))) |
---|
| 458 | done[++pos]=done[i]; |
---|
| 459 | else pDelete(&(done[i])); |
---|
| 460 | } |
---|
| 461 | if(done_index>0) |
---|
| 462 | done_index=pos+1; |
---|
| 463 | #ifdef TGB_DEBUG |
---|
| 464 | for(i=0;i<done_index;i++) |
---|
| 465 | { |
---|
| 466 | int j; |
---|
| 467 | for(j=0;j<i;j++) |
---|
| 468 | { |
---|
| 469 | assume((!pLmEqual(done[j],done[i]))); |
---|
| 470 | } |
---|
| 471 | } |
---|
| 472 | #endif |
---|
| 473 | #ifdef TGB_DEBUG |
---|
| 474 | for(i=0;i<done_index;i++) |
---|
| 475 | { |
---|
| 476 | pTest(done[i]); |
---|
| 477 | } |
---|
| 478 | #endif |
---|
| 479 | poly* m; |
---|
| 480 | int m_index=0; |
---|
| 481 | int m_size=0; |
---|
| 482 | for(i=0;i<p_index;i++) |
---|
| 483 | { |
---|
| 484 | m_size+=pLength(p[i]); |
---|
| 485 | } |
---|
| 486 | m=(poly*) omalloc(m_size*sizeof(poly)); |
---|
| 487 | //q=(poly*) omalloc(m_size*sizeof(poly)); |
---|
| 488 | |
---|
| 489 | |
---|
| 490 | |
---|
| 491 | for(i=0;i<p_index;i++) |
---|
| 492 | { |
---|
| 493 | |
---|
| 494 | poly p_i=p[i]; |
---|
| 495 | |
---|
| 496 | pTest(p[i]); |
---|
| 497 | |
---|
| 498 | while(p_i) |
---|
| 499 | { |
---|
| 500 | |
---|
| 501 | m[m_index]=pLmInit(p_i); |
---|
| 502 | |
---|
| 503 | pSetCoeff(m[m_index],nInit(1)); |
---|
| 504 | |
---|
| 505 | p_i=p_i->next; |
---|
| 506 | |
---|
| 507 | m_index++; |
---|
| 508 | |
---|
| 509 | } |
---|
| 510 | } |
---|
| 511 | |
---|
| 512 | int q_size=m_index; |
---|
| 513 | poly* q=(poly*) omalloc(q_size*sizeof(poly)); |
---|
| 514 | int q_index=0; |
---|
| 515 | //next Step additional reductors |
---|
| 516 | |
---|
| 517 | while(m_index>0) |
---|
| 518 | { |
---|
| 519 | #ifdef TGB_DEBUG |
---|
| 520 | |
---|
| 521 | for(i=0;i<done_index;i++) |
---|
| 522 | { |
---|
| 523 | pTest(done[i]); |
---|
| 524 | } |
---|
| 525 | #endif |
---|
| 526 | |
---|
| 527 | qsort(m, m_index,sizeof(poly),pLmCmp_func); |
---|
| 528 | |
---|
| 529 | |
---|
| 530 | pos=0; |
---|
| 531 | #ifdef TGB_DEBUG |
---|
| 532 | |
---|
| 533 | for(i=0;i<done_index;i++) |
---|
| 534 | { |
---|
| 535 | |
---|
| 536 | pTest(done[i]); |
---|
| 537 | } |
---|
| 538 | #endif |
---|
| 539 | for(i=1;i<m_index;i++) |
---|
| 540 | { |
---|
| 541 | pTest(m[i]); |
---|
| 542 | pTest(m[pos]); |
---|
| 543 | if((!(pLmEqual(m[i],m[pos])))) |
---|
| 544 | m[++pos]=m[i]; |
---|
| 545 | else pDelete(&(m[i])); |
---|
| 546 | } |
---|
| 547 | if(m_index>1) |
---|
| 548 | m_index=pos+1; |
---|
| 549 | if(done_size<done_index+m_index) |
---|
| 550 | { |
---|
| 551 | done_size=done_index+2*m_index; |
---|
| 552 | done=(poly*) omrealloc(done,done_size*sizeof(poly)); |
---|
| 553 | } |
---|
| 554 | if(chosen_size<chosen_index+m_index) |
---|
| 555 | { |
---|
| 556 | chosen_size=chosen_index+2*m_index; |
---|
| 557 | chosen=(monom_poly*) omrealloc(chosen,chosen_size*sizeof(monom_poly)); |
---|
| 558 | } |
---|
| 559 | q_index=0; |
---|
| 560 | if(q_size<m_index) |
---|
| 561 | { |
---|
| 562 | q_size=2*m_index; |
---|
| 563 | omfree(q); |
---|
| 564 | q=(poly*) omalloc(q_size*sizeof(poly)); |
---|
| 565 | } |
---|
| 566 | |
---|
| 567 | for(i=0;i<m_index;i++) |
---|
| 568 | { |
---|
| 569 | |
---|
| 570 | #ifdef TGB_DEBUG |
---|
| 571 | { |
---|
| 572 | int my_i; |
---|
| 573 | for(my_i=0;my_i<done_index;my_i++) |
---|
| 574 | { |
---|
| 575 | int my_j; |
---|
| 576 | for(my_j=0;my_j<my_i;my_j++) |
---|
| 577 | { |
---|
| 578 | assume((!pLmEqual(done[my_j],done[my_i]))); |
---|
| 579 | } |
---|
| 580 | } |
---|
| 581 | } |
---|
| 582 | #endif |
---|
| 583 | BOOLEAN in_done=FALSE; |
---|
| 584 | pTest(m[i]); |
---|
| 585 | |
---|
| 586 | pos=posInPolys (done, done_index, m[i],c); |
---|
| 587 | #ifdef TGB_DEBUG |
---|
| 588 | { |
---|
| 589 | int my_i; |
---|
| 590 | for (my_i=0;my_i<done_index;my_i++) |
---|
| 591 | pTest(done[my_i]); |
---|
| 592 | } |
---|
| 593 | #endif |
---|
| 594 | if(((done_index>pos)&&(pLmEqual(m[i],done[pos]))) ||(pos>0) &&(pLmEqual(m[i],done[pos-1]))) |
---|
| 595 | in_done=TRUE; |
---|
| 596 | if (!(in_done)) |
---|
| 597 | { |
---|
| 598 | |
---|
| 599 | int S_pos=kFindDivisibleByInS_easy(c->strat,m[i], pGetShortExpVector(m[i])); |
---|
| 600 | if(S_pos>=0) |
---|
| 601 | { |
---|
| 602 | monom_poly h; |
---|
| 603 | h.f=c->strat->S[S_pos]; |
---|
| 604 | h.m=pOne(); |
---|
| 605 | int* ev=(int*) omalloc((c->r->N+1)*sizeof(int)); |
---|
| 606 | pGetExpV(m[i],ev); |
---|
| 607 | pSetExpV(h.m,ev); |
---|
| 608 | omfree(ev); |
---|
| 609 | pExpVectorSub(h.m,c->strat->S[S_pos]); |
---|
| 610 | simplify(h,c); |
---|
| 611 | q[q_index]=ppMult_mm(h.f,h.m); |
---|
| 612 | chosen[chosen_index++]=h; |
---|
| 613 | q_index++; |
---|
| 614 | } |
---|
| 615 | pTest(m[i]); |
---|
| 616 | #ifdef TGB_DEBUG |
---|
| 617 | { |
---|
| 618 | int my_i; |
---|
| 619 | for (my_i=0;my_i<done_index;my_i++) |
---|
| 620 | pTest(done[my_i]); |
---|
| 621 | } |
---|
| 622 | #endif |
---|
| 623 | memmove(&(done[pos+1]),&(done[pos]), (done_index-pos)*sizeof(poly)); |
---|
| 624 | done[pos]=m[i]; |
---|
| 625 | done_index++; |
---|
| 626 | #ifdef TGB_DEBUG |
---|
| 627 | { |
---|
| 628 | int my_i; |
---|
| 629 | for (my_i=0;my_i<done_index;my_i++) |
---|
| 630 | { |
---|
| 631 | // Print("Position %i pos %i size %i\n",my_i,pos,done_index); |
---|
| 632 | pTest(done[my_i]); |
---|
| 633 | } |
---|
| 634 | } |
---|
| 635 | #endif |
---|
| 636 | } |
---|
| 637 | else |
---|
| 638 | pDelete(&m[i]); |
---|
| 639 | #ifdef TGB_DEBUG |
---|
| 640 | { |
---|
| 641 | int my_i; |
---|
| 642 | for(my_i=0;my_i<done_index;my_i++) |
---|
| 643 | { |
---|
| 644 | pTest(done[my_i]); |
---|
| 645 | int my_j; |
---|
| 646 | for(my_j=0;my_j<my_i;my_j++) |
---|
| 647 | { |
---|
| 648 | assume((!pLmEqual(done[my_j],done[my_i]))); |
---|
| 649 | } |
---|
| 650 | } |
---|
| 651 | } |
---|
| 652 | #endif |
---|
| 653 | } |
---|
| 654 | if(p_size<p_index+q_index) |
---|
| 655 | { |
---|
| 656 | p_size=p_index+2*q_index; |
---|
| 657 | p=(poly*) omrealloc(p,p_size*sizeof(poly)); |
---|
| 658 | } |
---|
| 659 | for (i=0;i<q_index;i++) |
---|
| 660 | p[p_index++]=q[i]; |
---|
| 661 | m_index=0; |
---|
| 662 | int sum=0; |
---|
| 663 | for(i=0;i<p_index;i++) |
---|
| 664 | { |
---|
| 665 | sum+=pLength(p[i])-1; |
---|
| 666 | } |
---|
| 667 | if (m_size<sum) |
---|
| 668 | { |
---|
| 669 | omfree(m); |
---|
| 670 | m=(poly*) omalloc(sum*sizeof(poly)); |
---|
| 671 | } |
---|
| 672 | m_size=sum; |
---|
| 673 | for(i=0;i<q_index;i++) |
---|
| 674 | { |
---|
| 675 | poly p_i=q[i]->next; |
---|
| 676 | while(p_i) |
---|
| 677 | { |
---|
| 678 | m[m_index]=pLmInit(p_i); |
---|
| 679 | pSetCoeff(m[m_index],nInit(1)); |
---|
| 680 | p_i=p_i->next; |
---|
| 681 | m_index++; |
---|
| 682 | } |
---|
| 683 | } |
---|
| 684 | //qsort(done, done_index,sizeof(poly),pLmCmp_func); |
---|
| 685 | #ifdef TGB_DEBUG |
---|
| 686 | for(i=0;i<done_index;i++) |
---|
| 687 | { |
---|
| 688 | pTest(done[i]); |
---|
| 689 | } |
---|
| 690 | #endif |
---|
| 691 | } |
---|
| 692 | #ifdef TGB_DEBUG |
---|
| 693 | for(i=0;i<done_index;i++) |
---|
| 694 | { |
---|
| 695 | int j; |
---|
| 696 | for(j=0;j<i;j++) |
---|
| 697 | { |
---|
| 698 | assume((!pLmEqual(done[j],done[i]))); |
---|
| 699 | } |
---|
| 700 | } |
---|
| 701 | #endif |
---|
| 702 | omfree(m); |
---|
| 703 | omfree(q); |
---|
| 704 | if (TEST_OPT_PROT) |
---|
| 705 | Print("Mat[%i x %i], ",chosen_index, done_index); |
---|
| 706 | //next Step build matrix |
---|
| 707 | #ifdef TGB_DEBUG |
---|
| 708 | for(i=0;i<done_index;i++) |
---|
| 709 | { |
---|
| 710 | int j; |
---|
| 711 | for(j=0;j<i;j++) |
---|
| 712 | { |
---|
| 713 | assume((!pLmEqual(done[j],done[i]))); |
---|
| 714 | } |
---|
| 715 | } |
---|
| 716 | #endif |
---|
| 717 | assume(p_index==chosen_index); |
---|
| 718 | |
---|
| 719 | // tgb_matrix* mat=build_matrix(p,p_index,done, done_index,c); |
---|
| 720 | |
---|
| 721 | // simple_gauss2(mat); |
---|
| 722 | tgb_sparse_matrix* mat=build_sparse_matrix(p,p_index,done, done_index,c); |
---|
| 723 | simple_gauss(mat,c); |
---|
| 724 | m_size=mat->get_rows(); |
---|
| 725 | m=(poly*) omalloc(m_size*sizeof(poly)); |
---|
| 726 | m_index=retranslate(m,mat,done,c); |
---|
| 727 | |
---|
| 728 | mat=NULL; |
---|
| 729 | for(i=0;i<done_index;i++) |
---|
| 730 | pDelete(&done[i]); |
---|
| 731 | omfree(done); |
---|
| 732 | done=NULL; |
---|
| 733 | //next Step addElements to basis |
---|
| 734 | int F_plus_size=m_index; |
---|
| 735 | poly* F_plus=(poly*)omalloc(F_plus_size*sizeof(poly)); |
---|
| 736 | int F_plus_index=0; |
---|
| 737 | int F_minus_size=m_index; |
---|
| 738 | poly* F_minus=(poly*) omalloc(F_minus_size*sizeof(poly)); |
---|
| 739 | int F_minus_index=0; |
---|
| 740 | |
---|
| 741 | //better algorithm replace p by its monoms, qsort,delete duplicates and binary search for testing if monom is contained in array |
---|
| 742 | qsort(p, p_index,sizeof(poly),pLmCmp_func); |
---|
| 743 | for(i=0;i<p_index;i++) |
---|
| 744 | pDelete(&p[i]->next); |
---|
| 745 | for(i=m_index-1;i>=0;--i) |
---|
| 746 | { |
---|
| 747 | int j; |
---|
| 748 | int pos=posInPolys (p,p_index, m[i],c); |
---|
| 749 | BOOLEAN minus=FALSE; |
---|
| 750 | if(((p_index>pos)&&(pLmEqual(m[i],p[pos]))) ||(pos>0) &&(pLmEqual(m[i],p[pos-1]))) |
---|
| 751 | minus=TRUE; |
---|
| 752 | |
---|
| 753 | if(minus) |
---|
| 754 | { |
---|
| 755 | F_minus[F_minus_index++]=m[i]; |
---|
| 756 | m[i]=NULL; |
---|
| 757 | } |
---|
| 758 | else |
---|
| 759 | { |
---|
| 760 | F_plus[F_plus_index++]=m[i]; |
---|
| 761 | m[i]=NULL; |
---|
| 762 | } |
---|
| 763 | } |
---|
| 764 | |
---|
| 765 | // for(i=m_index-1;i>=0;--i) |
---|
| 766 | // { |
---|
| 767 | // int j; |
---|
| 768 | // BOOLEAN minus=FALSE; |
---|
| 769 | // for(j=0;j<p_index;j++) |
---|
| 770 | // if (pLmEqual(p[j],m[i])) |
---|
| 771 | // { |
---|
| 772 | // minus=TRUE; |
---|
| 773 | // break; |
---|
| 774 | // } |
---|
| 775 | // if(minus) |
---|
| 776 | // { |
---|
| 777 | // F_minus[F_minus_index++]=m[i]; |
---|
| 778 | // m[i]=NULL; |
---|
| 779 | // } |
---|
| 780 | // else |
---|
| 781 | // { |
---|
| 782 | // F_plus[F_plus_index++]=m[i]; |
---|
| 783 | // m[i]=NULL; |
---|
| 784 | // } |
---|
| 785 | // } |
---|
| 786 | //in this order F_minus will be automatically ascending sorted |
---|
| 787 | //to make this sure for foreign gauss |
---|
| 788 | //uncomment the following line |
---|
| 789 | // qsort(F_minus, F_minus_index,sizeof(poly),pLmCmp_func); |
---|
| 790 | assume((F_plus_index+F_minus_index)==m_index); |
---|
| 791 | if (TEST_OPT_PROT) |
---|
| 792 | Print("%i]", F_plus_index); |
---|
| 793 | for(i=0;i<p_index;i++) |
---|
| 794 | pDelete(&p[i]); |
---|
| 795 | omfree(p); |
---|
| 796 | p=NULL; |
---|
| 797 | omfree(m); |
---|
| 798 | m=NULL; |
---|
| 799 | |
---|
| 800 | //the F_minus list must be cleared separately at the end |
---|
| 801 | mp_array_list** F_i; |
---|
| 802 | poly_array_list** F_m_i; |
---|
| 803 | F_i=&(c->F); |
---|
| 804 | F_m_i=&(c->F_minus); |
---|
| 805 | |
---|
| 806 | while((*F_i)!=NULL) |
---|
| 807 | { |
---|
| 808 | assume((*F_m_i)!=NULL); |
---|
| 809 | F_i=(&((*F_i)->next)); |
---|
| 810 | F_m_i=(&((*F_m_i)->next)); |
---|
| 811 | } |
---|
| 812 | assume((*F_m_i)==NULL); |
---|
| 813 | //should resize the array to save memory |
---|
| 814 | //F and F_minus |
---|
| 815 | qsort(chosen,chosen_index,sizeof(monom_poly),monom_poly_crit);//important for simplify |
---|
| 816 | (*F_m_i)=(poly_array_list*) omalloc(sizeof(poly_array_list)); |
---|
| 817 | (*F_m_i)->size=F_minus_index; |
---|
| 818 | (*F_m_i)->p=F_minus; |
---|
| 819 | (*F_m_i)->next=NULL; |
---|
| 820 | (*F_i)=(mp_array_list*) omalloc(sizeof(poly_array_list)); |
---|
| 821 | (*F_i)->size=chosen_index; |
---|
| 822 | (*F_i)->mp=chosen; |
---|
| 823 | (*F_i)->next=NULL; |
---|
| 824 | |
---|
| 825 | if(F_plus_index>0) |
---|
| 826 | { |
---|
| 827 | int j; |
---|
| 828 | int* ibuf=(int*) omalloc(F_plus_index*sizeof(int)); |
---|
| 829 | sorted_pair_node*** sbuf=(sorted_pair_node***) omalloc(F_plus_index*sizeof(sorted_pair_node**)); |
---|
| 830 | |
---|
| 831 | for(j=0;j<F_plus_index;j++) |
---|
| 832 | { |
---|
| 833 | int len; |
---|
| 834 | poly p=F_plus[j]; |
---|
| 835 | |
---|
| 836 | // delete buf[j]; |
---|
| 837 | //remember to free res here |
---|
| 838 | // p=redTailShort(p, c->strat); |
---|
[0f2374] | 839 | sbuf[j]=add_to_basis_ideal_quotient(p,c,ibuf+j); |
---|
[b75d13] | 840 | |
---|
| 841 | } |
---|
| 842 | int sum=0; |
---|
| 843 | for(j=0;j<F_plus_index;j++) |
---|
| 844 | { |
---|
| 845 | sum+=ibuf[j]; |
---|
| 846 | } |
---|
| 847 | sorted_pair_node** big_sbuf=(sorted_pair_node**) omalloc(sum*sizeof(sorted_pair_node*)); |
---|
| 848 | int partsum=0; |
---|
| 849 | for(j=0;j<F_plus_index;j++) |
---|
| 850 | { |
---|
| 851 | memmove(big_sbuf+partsum, sbuf[j],ibuf[j]*sizeof(sorted_pair_node*)); |
---|
| 852 | omfree(sbuf[j]); |
---|
| 853 | partsum+=ibuf[j]; |
---|
| 854 | } |
---|
| 855 | |
---|
| 856 | qsort(big_sbuf,sum,sizeof(sorted_pair_node*),tgb_pair_better_gen2); |
---|
| 857 | c->apairs=spn_merge(c->apairs,c->pair_top+1,big_sbuf,sum,c); |
---|
| 858 | c->pair_top+=sum; |
---|
| 859 | clean_top_of_pair_list(c); |
---|
| 860 | omfree(big_sbuf); |
---|
| 861 | omfree(sbuf); |
---|
| 862 | omfree(ibuf); |
---|
| 863 | } |
---|
| 864 | omfree(F_plus); |
---|
| 865 | #ifdef TGB_DEBUG |
---|
| 866 | int z; |
---|
| 867 | for(z=1;z<=c->pair_top;z++) |
---|
| 868 | { |
---|
| 869 | assume(pair_better(c->apairs[z],c->apairs[z-1],c)); |
---|
| 870 | } |
---|
| 871 | #endif |
---|
| 872 | |
---|
| 873 | return; |
---|
| 874 | } |
---|