1 | |
---|
2 | |
---|
3 | |
---|
4 | #include <kernel/mod2.h> |
---|
5 | #include <kernel/GBEngine/F4.h> |
---|
6 | #include <kernel/GBEngine/tgb_internal.h> |
---|
7 | #include <kernel/GBEngine/tgbgauss.h> |
---|
8 | /**************************************** |
---|
9 | * Computer Algebra System SINGULAR * |
---|
10 | ****************************************/ |
---|
11 | /* |
---|
12 | * ABSTRACT: F4 implementation |
---|
13 | */ |
---|
14 | static int posInPolys (poly* p, int pn, poly qe,slimgb_alg* c) |
---|
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 | |
---|
41 | static tgb_sparse_matrix* build_sparse_matrix(poly* p,int p_index,poly* done, int done_index, slimgb_alg* c){ |
---|
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 | } |
---|
87 | static tgb_matrix* build_matrix(poly* p,int p_index,poly* done, int done_index, slimgb_alg* c){ |
---|
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 | |
---|
128 | static int retranslate(poly* m,tgb_sparse_matrix* mat,poly* done, slimgb_alg* c){ |
---|
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 |
---|
147 | static int retranslate(poly* m,tgb_matrix* mat,poly* done, slimgb_alg* c){ |
---|
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 | } |
---|
210 | static int posInMonomPolys (monom_poly* p, int pn, monom_poly & qe,slimgb_alg* c) |
---|
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 | } |
---|
236 | static void simplify(monom_poly& h, slimgb_alg* c){ |
---|
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 | |
---|
295 | void go_on_F4 (slimgb_alg* c){ |
---|
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); |
---|
841 | sbuf[j]=add_to_basis_ideal_quotient(p,c,ibuf+j); |
---|
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 | } |
---|