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