source: git/kernel/f5gb.cc @ c3da59

spielwiese
Last change on this file since c3da59 was c3da59, checked in by Christian Eder, 15 years ago
still problems with rTag, old structures with new idea of reduction of GB after each iteration step => lots of optimizations need to be done git-svn-id: file:///usr/local/Singular/svn/trunk@11525 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 38.8 KB
Line 
1/****************************************
2*  Computer Algebra System SINGULAR     *
3****************************************/
4/* $Id: f5gb.cc,v 1.41 2009-03-08 20:20:21 ederc Exp $ */
5/*
6* ABSTRACT: f5gb interface
7*/
8#include "mod2.h"
9#ifdef HAVE_F5
10#include "kutil.h"
11#include "structs.h"
12#include "omalloc.h"
13#include "polys.h"
14#include "p_polys.h"
15#include "ideals.h"
16#include "febase.h"
17#include "kstd1.h"
18#include "khstd.h"
19#include "kbuckets.h"
20#include "weight.h"
21#include "intvec.h"
22#include "pInline1.h"
23#include "f5gb.h"
24#include "f5data.h"
25#include "f5lists.h"
26
27int reductionsToZero   =  0;
28
29/*
30====================================================================
31sorting ideals by decreasing total degree "left" and "right" are the
32pointer of the first and last polynomial in the considered ideal
33====================================================================
34*/
35void qsortDegree(poly* left, poly* right) {
36    poly* ptr1 = left;
37    poly* ptr2 = right;
38    poly p1,p2;
39    p2 = *(left + (right - left >> 1));
40    do {
41            while(pTotaldegree(*ptr1, currRing) < pTotaldegree(p2, currRing)) {
42                    ptr1++;
43            } 
44            while(pTotaldegree(*ptr2, currRing) > pTotaldegree(p2,currRing)) {
45                    ptr2--;
46            }
47            if(ptr1 > ptr2) {
48                    break;
49            }
50            p1    = *ptr1;
51            *ptr1 = *ptr2;
52            *ptr2 = p1;
53    } while(++ptr1 <= --ptr2);
54
55    if(left < ptr2) {
56            qsortDegree(left,ptr2);
57    }
58    if(ptr1 < right) {
59            qsortDegree(ptr1,right);
60    }
61}
62
63
64
65/*
66==================================================
67computes incrementally gbs of subsets of the input
68gb{f_m} -> gb{f_m,f_(m-1)} -> gb{f_m,...,f_1}
69==================================================
70*/
71LList* F5inc(int i, poly f_i, LList* gPrev, ideal gbPrev, poly ONE, LTagList* lTag, RList* rules, RTagList* rTag) {
72    //Print("in f5inc\n");           
73    //pWrite(rules->getFirst()->getRuleTerm());
74    int j;
75    //Print("%p\n",gPrev->getFirst());
76    //pWrite(gPrev->getFirst()->getPoly());
77    poly tempNF =   kNF(gbPrev,currQuotient,f_i);
78    f_i         =   tempNF;
79    //gPrev->insert(ONE,i,f_i);
80    gPrev->insert(ONE,gPrev->getLength()+1,f_i);
81    // tag the first element in gPrev of the current index for findReductor()
82    lTag->setFirstCurrentIdx(gPrev->getLast());
83    //Print("1st gPrev: ");
84    //pWrite(gPrev->getFirst()->getPoly());
85    //Print("2nd gPrev: ");
86    //pWrite(gPrev->getFirst()->getNext()->getPoly());
87    //pWrite(gPrev->getFirst()->getNext()->getPoly());   
88    CList* critPairs        =   new CList();
89    CNode* critPairsMinDeg  =   new CNode();   
90    // computation of critical pairs with checking of criterion 1 and criterion 2 and saving them
91    // in the list critPairs
92    criticalPair(gPrev, critPairs, lTag, rTag, rules);
93    static LList* sPolyList        =   new LList();
94    // labeled polynomials which have passed reduction() and have to be added to list gPrev
95    static LList* completed        =   new LList();
96    // the reduced labeled polynomials which are returned from subalgorithm reduction()
97    static LList* reducedLPolys     =   new LList();
98    // while there are critical pairs to be further checked and deleted/computed
99    while(NULL != critPairs->getFirst()->getData()) { 
100        // critPairs->getMinDeg() deletes the first elements of minimal degree from
101        // critPairs, thus the while loop is not infinite.
102        critPairsMinDeg =   critPairs->getMinDeg();
103        // adds all to be reduced S-polynomials in the list sPolyList and adds
104        // the corresponding rules to the list rules
105        // NOTE: inside there is a second check of criterion 2 if new rules are
106        // added
107        computeSPols(critPairsMinDeg,rTag,rules,sPolyList);
108       
109        // DEBUG STUFF FOR SPOLYLIST
110        LNode* temp     =   sPolyList->getFirst();
111        //while(NULL != temp && NULL != temp->getLPoly()) {
112            //Print("Spolylist element: ");
113            //pWrite(temp->getPoly());
114            //temp    =   temp->getNext();
115        //}
116        // reduction process of new S-polynomials and also adds new critical pairs to critPairs
117        reduction(sPolyList, critPairs, gPrev, rules, lTag, rTag, gbPrev);
118       
119        // DEBUG STUFF FOR GPREV
120        //temp    =   gPrev->getFirst();
121        //int number  =   1;
122        //Print("\n\n");
123        //while(NULL != temp) {
124        //    Print("%d.  ",number);
125        //    pWrite(temp->getPoly());
126        //    temp    =   temp->getNext();
127        //    number++;
128        //    Print("\n");
129        //}
130        //sleep(5);
131   
132    }
133    //Print("REDUCTION DONE\n");
134    //Print("%p\n",rules->getFirst());
135    //Print("%p\n",rTag->getFirst());
136    //if(rules->getFirst() != rTag->getFirst()) {
137        //Print("+++++++++++++++++++++++++++++++++++++NEW RULES+++++++++++++++++++++++++++++++++++++\n");
138        //rTag->insert(rules->getFirst());
139    //}
140    //else {
141        //Print("+++++++++++++++++++++++++++++++++++NO NEW RULES++++++++++++++++++++++++++++++++++++\n");
142    //}
143    lTag->insert(lTag->getFirstCurrentIdx());
144    //Print("INDEX: %d\n",tempTag->getIndex());
145    //pWrite(tempTag->getPoly());
146    //Print("COMPLETED FIRST IN F5INC: \n");
147    //Print("1st gPrev: ");
148    //pWrite(gPrev->getFirst()->getPoly());
149    //Print("2nd gPrev: ");
150    //pWrite(gPrev->getFirst()->getNext()->getPoly());
151    //Print("3rd gPrev: ");
152    //pWrite(gPrev->getFirst()->getNext()->getNext()->getPoly());
153    //delete sPolyList;
154    //critPairs->print();
155    delete critPairs;
156    //Print("IN F5INC\n");
157   
158    //gPrev->print();
159    return gPrev;
160}
161
162
163
164/*
165================================================================
166computes a list of critical pairs for the next reduction process
167first element in gPrev is always the newest element which must
168build critical pairs with all other elements in gPrev
169================================================================
170*/
171inline void criticalPair(LList* gPrev, CList* critPairs, LTagList* lTag, RTagList* rTag, RList* rules) {
172    //Print("IN CRITPAIRS\n");
173    // initialization for usage in pLcm()
174    number nOne         =   nInit(1);
175    LNode* newElement   =   gPrev->getLast();
176    LNode* temp         =   gPrev->getFirst();
177    poly u1             =   pOne();
178    poly u2             =   pOne();
179    poly lcm            =   pOne();
180    poly t              =   pHead(newElement->getPoly());
181    Rule* testedRule    =   rules->getFirst()->getRule();
182    // computation of critical pairs
183    while( gPrev->getLast() != temp) {
184        //pWrite( *(gPrev->getFirst()->getPoly()) );
185       // pWrite( *(l->getPoly()) );
186        pLcm(newElement->getPoly(), temp->getPoly(), lcm);
187        pSetCoeff(lcm,nOne);
188        // computing factors u2 for new labels
189        u1 = pDivide(lcm,t);
190        pSetCoeff(u1,nOne);
191        u2 = pDivide(lcm, pHead(temp->getPoly()));
192        pSetCoeff(u2,nOne);
193        //if(gPrev->getLast()->getIndex()==5) {
194            //Print("IN CRITPAIRS\n");
195        //    pWrite(u1);
196        //    Print("1st ELEMENT: ");
197        //    pWrite(newElement->getPoly());
198        //    Print("2nd ELEMENT: ");
199        //    pWrite(temp->getPoly());
200        //}
201        // testing both new labels by the F5 Criterion
202        if(!criterion1(gPrev,u1,newElement,lTag) && !criterion1(gPrev,u2,temp,lTag) && 
203           !criterion2(u1, newElement, rules, rTag) && !criterion2(u2, temp, rules, rTag)) {
204            // if they pass the test, add them to CList critPairs, having the LPoly with greater
205            // label as first element in the CPair
206            if(newElement->getIndex() == temp->getIndex() && 
207            -1 == pLmCmp(ppMult_qq(u1, newElement->getTerm()),ppMult_qq(u2, temp->getTerm()))) {
208                //Print("zweites groesser\n");
209                CPair* cp   =   new CPair(pDeg(ppMult_qq(u2,pHead(temp->getPoly()))), u2, 
210                                temp->getLPoly(), u1, newElement->getLPoly(), testedRule);                   
211                critPairs->insert(cp);
212            }
213            else {
214                //Print("erstes groesser\n");
215                CPair* cp   =   new CPair(pDeg(ppMult_qq(u2,pHead(temp->getPoly()))), u1, 
216                                newElement->getLPoly(), u2, temp->getLPoly(), testedRule);                   
217                critPairs->insert(cp);
218            }
219        }
220        else {
221        }
222       
223        //Print("\n\n");
224        temp    =   temp->getNext();
225    }
226    // for debugging
227    //if(NULL != critPairs) {
228        //critPairs->print();
229        //sleep(5);
230    //}
231}
232
233
234
235
236
237
238
239/*
240========================================
241Criterion 1, i.e. Faugere's F5 Criterion
242========================================
243*/
244inline bool criterion1(LList* gPrev, poly t, LNode* l, LTagList* lTag) {
245    // starts at the first element in gPrev with index = (index of l)-1, these tags are saved in lTag
246        int idx =   l->getIndex();
247    if(idx == 1) {
248        return false;
249    }
250    else {
251        LNode* testNode =   gPrev->getFirst();
252        // save the monom t1*label_term(l) as it is tested various times in the following
253        poly u1 = ppMult_qq(t,l->getTerm());
254        //Print("------------------------------IN CRITERION 1-----------------------------------------\n");
255        //Print("TESTED ELEMENT: ");
256        //pWrite(l->getPoly());
257        //pWrite(l->getTerm());
258        //pWrite(ppMult_qq(t,l->getTerm()));
259        //Print("%d\n\n",l->getIndex());
260        while(testNode->getIndex() < idx && NULL != testNode->getLPoly()) {
261            //pWrite(testNode->getPoly());
262            //Print("%d\n",testNode->getIndex());
263            if(pLmDivisibleByNoComp(testNode->getPoly(),u1)) {
264                //Print("Criterion 1 NOT passed!\n");
265                return true;
266            }
267            //pWrite(testNode->getNext()->getPoly());
268            testNode    =   testNode->getNext();
269        }
270        return false;
271    }
272}
273
274
275
276/*
277=====================================
278Criterion 2, i.e. Rewritten Criterion
279=====================================
280*/
281inline bool criterion2(poly t, LNode* l, RList* rules, RTagList* rTag) {
282    //Print("------------------------------IN CRITERION 2/1-----------------------------------------\n");
283        //Print("RULES: \n");
284        //RNode* tempR    =   rules->getFirst();
285        //while(NULL != tempR->getRule()) {
286            //Print("ADDRESS OF RNODE: %p\n",tempR);
287            //Print("ADDRESS OF RULE: %p\n",tempR->getRule());
288            //pWrite(tempR->getRuleTerm());
289            //Print("ADDRESS OF TERM: %p\n",tempR->getRuleTerm());
290            //Print("%d\n\n",tempR->getRuleIndex());
291            //tempR   =   tempR->getNext();
292        //}
293        //Print("TESTED ELEMENT: ");
294        //pWrite(l->getPoly());
295        //pWrite(l->getTerm());
296        //pWrite(ppMult_qq(t,l->getTerm()));
297        //Print("%d\n\n",l->getIndex());
298// start at the previously added element to gPrev, as all other elements will have the same index for sure
299    RNode* testNode =   new RNode();
300    /*
301    if(NULL == rTag->getFirst()->getRule()) {
302        testNode    =   rules->getFirst();
303    }
304    else {
305        if(l->getIndex() > rTag->getFirst()->getRuleIndex()) {
306            testNode    =   rules->getFirst();
307        }
308        else {
309        //Print("DEBUG\n");
310        //Print("L INDEX: %d\n",l->getIndex());
311            /*-------------------------------------
312             * TODO: WHEN INTERREDUCING THE GB THE
313             *       INDEX OF THE PREVIOUS ELEMENTS
314             *       GETS HIGHER!
315             *-----------------------------------
316            testNode    =   rules->getFirst();
317            //testNode    =   rTag->get(l->getIndex());
318            //Print("TESTNODE ADDRESS: %p\n",testNode);
319        }
320    }
321    */
322    testNode    =   rules->getFirst();
323        // save the monom t1*label_term(l) as it is tested various times in the following
324    poly u1 = ppMult_qq(t,l->getTerm());
325    // first element added to rTag was NULL, check for this
326    //Print("%p\n",testNode->getRule());
327    // NOTE: testNode is possibly NULL as rTag->get() returns NULL for elements of index <=1!
328    //Print("TESTNODE: %p\n",testNode);
329    //pWrite(testNode->getRuleTerm());
330    if(NULL != testNode && NULL != testNode->getRule()) {   
331        //pWrite(testNode->getRuleTerm());
332    }
333    if(NULL != testNode) {
334        if(testNode->getRule() == l->getRule()) {
335            //Print("%p\n%p\n",testNode->getRule(),l->getRule());
336            //Print("EQUAL\n");
337        }
338        else {
339            //Print("NOT EQUAL\n");
340        }
341    }
342    while(NULL != testNode && NULL != testNode->getRule() && testNode->getRule() != l->getRule() 
343          && l->getIndex() == testNode->getRuleIndex()) {
344        //Print("%p\n",testNode);
345        //pWrite(testNode->getRuleTerm());
346        //pWrite(t);
347        //pWrite(l->getTerm());
348        //pWrite(u1);
349        //Print("%d\n",testNode->getRuleIndex());
350        if(pLmDivisibleBy(testNode->getRuleTerm(),u1)) {
351            //Print("Criterion 2 NOT passed!\n");
352            pDelete(&u1);
353            return true;
354        }
355                testNode    =   testNode->getNext();
356    }
357    //delete testNode;
358    pDelete(&u1);
359    return false;
360}
361
362
363
364/*
365=================================================================================================================
366Criterion 2, i.e. Rewritten Criterion, for its second call in computeSPols(), with added lastRuleTested parameter
367=================================================================================================================
368*/
369inline bool criterion2(poly t, LPoly* l, RList* rules, Rule* testedRule) {
370    //Print("------------------------------IN CRITERION 2/2-----------------------------------------\n");
371    //Print("LAST RULE TESTED: %p",testedRule);
372    //Print("RULES: \n");
373        RNode* tempR    =   rules->getFirst();
374        //while(NULL != tempR->getRule()) {
375            //pWrite(tempR->getRuleTerm());
376            //Print("%d\n\n",tempR->getRuleIndex());
377            //tempR   =   tempR->getNext();
378        //}
379        //Print("TESTED ELEMENT: ");
380        //pWrite(l->getPoly());
381        //pWrite(l->getTerm());
382        //pWrite(ppMult_qq(t,l->getTerm()));
383        //Print("%d\n\n",l->getIndex());
384// start at the previously added element to gPrev, as all other elements will have the same index for sure
385        RNode* testNode =   rules->getFirst();
386    // save the monom t1*label_term(l) as it is tested various times in the following
387    poly u1 = ppMult_qq(t,l->getTerm());
388        // first element added to rTag was NULL, check for this
389        while(NULL != testNode->getRule() && testNode->getRule() != testedRule) {
390        //pWrite(testNode->getRuleTerm());
391        if(pLmDivisibleBy(testNode->getRuleTerm(),u1)) {
392            //Print("Criterion 2 NOT passed!\n");
393            pDelete(&u1);
394            return true;
395        }
396                testNode    =   testNode->getNext();
397    }
398    pDelete(&u1);
399    return false;
400}
401
402
403
404/*
405==================================
406Computation of S-Polynomials in F5
407==================================
408*/
409void computeSPols(CNode* first, RTagList* rTag, RList* rules, LList* sPolyList) { 
410    CNode* temp         =   first;
411    poly sp      =   pInit();
412    //Print("###############################IN SPOLS##############################\n");
413    while(NULL != temp->getData()) {
414        // only if a new rule was added since the last test in subalgorithm criticalPair()
415        //if(rules->getFirst() != rTag->getFirst()) {
416            if(!criterion2(temp->getT1(),temp->getAdLp1(),rules,temp->getTestedRule())) {
417                // the second component is tested only when it has the actual index, otherwise there is
418                // no new rule to test since the last test in subalgorithm criticalPair()
419                if(temp->getLp2Index() == temp->getLp1Index()) {
420                    if(!criterion2(temp->getT2(),temp->getAdLp2(),rules,temp->getTestedRule())) {
421                        // computation of S-polynomial
422                        sp      =   pSub(ppMult_qq(temp->getT1(),temp->getLp1Poly()),
423                                         ppMult_qq(temp->getT2(),temp->getLp2Poly()));
424                        //Print("BEGIN SPOLY1\n====================\n");
425                        pNorm(sp);
426                        //pWrite(sp);
427                        //Print("END SPOLY1\n====================\n");
428                        if(NULL == sp) {
429                            // as rules consist only of pointers we need to save the labeled
430                            // S-polynomial also of a zero S-polynomial
431                            //zeroList->insert(temp->getAdT1(),temp->getLp1Index(),&sp);
432                            // origin of rule can be set NULL as the labeled polynomial
433                            // will never be used again as it is zero => no problems with
434                            // further criterion2() tests and termination conditions
435                            //Print("~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ZERO REDUCTION~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~\n");
436                                                        reductionsToZero++;
437                        //Print("IN SPOLS 1\n");
438                            rules->insert(temp->getLp1Index(),ppMult_qq(temp->getT1(),temp->getLp1Term()));
439                            //Print("RULE ADDED: \n");
440                            //pWrite(rules->getFirst()->getRuleTerm());
441                            // as sp = NULL, delete it
442                            pDelete(&sp);
443                            //Print("HIER\n");
444                        }
445                        else {
446                        //Print("IN SPOLS 2\n");
447                            rules->insert(temp->getLp1Index(),ppMult_qq(temp->getT1(),temp->getLp1Term()));
448                            //Print("RULE ADDED: \n");
449                            //pWrite(rules->getFirst()->getRuleTerm()); 
450                            sPolyList->insertByLabel(ppMult_qq(temp->getT1(),temp->getLp1Term()),temp->getLp1Index(),sp,rules->getFirst()->getRule());
451                        }
452                        // data is saved in sPolyList or zero => delete sp
453                    }
454                }
455                else { // temp->getLp2Index() < temp->getLp1Index()
456                    // computation of S-polynomial
457                    sp      =   pSub(ppMult_qq(temp->getT1(),temp->getLp1Poly()),
458                                     ppMult_qq(temp->getT2(),temp->getLp2Poly()));
459                    //Print("BEGIN SPOLY2\n====================\n");
460                    pNorm(sp);
461                    //pWrite(sp);
462                    //Print("END SPOLY2\n====================\n");
463                    if(NULL == sp) {
464                        // zeroList->insert(temp->getAdT1(),temp->getLp1Index(),&sp);
465                        // origin of rule can be set NULL as the labeled polynomial
466                        // will never be used again as it is zero => no problems with
467                        // further criterion2() tests and termination conditions
468                            //Print("~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ZERO REDUCTION~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~\n");
469                        reductionsToZero++;
470                        //Print("IN SPOLS 3\n");
471                        rules->insert(temp->getLp1Index(),ppMult_qq(temp->getT1(),temp->getLp1Term()));
472                        //Print("RULE ADDED: \n");
473                        //pWrite(rules->getFirst()->getRuleTerm());
474                        // as sp = NULL, delete it
475                        pDelete(&sp);
476                    }
477                    else {
478                        //Print("IN SPOLS 4\n");
479                        rules->insert(temp->getLp1Index(),ppMult_qq(temp->getT1(),temp->getLp1Term()));
480                        //Print("RULE ADDED: \n");
481                        //pWrite(rules->getFirst()->getRuleTerm());
482                        //Print("%d\n",rules->getFirst()->getRuleIndex());
483                        //Print("%p\n",sPolyList->getFirst());
484                        sPolyList->insertByLabel(ppMult_qq(temp->getT1(),temp->getLp1Term()),temp->getLp1Index(),sp,rules->getFirst()->getRule());
485                    }
486                    // data is saved in sPolyList or zero => delete sp
487                }
488            }
489        //}
490        temp    =   temp->getNext();
491    }
492    // these critical pairs can be deleted now as they are either useless for further computations or
493    // already saved as an S-polynomial to be reduced in the following
494    delete first;   
495}
496
497
498
499/*
500========================================================================
501reduction including subalgorithm topReduction() using Faugere's criteria
502========================================================================
503*/
504void reduction(LList* sPolyList, CList* critPairs, LList* gPrev, RList* rules, LTagList* lTag, RTagList* rTag, 
505                 ideal gbPrev) { 
506    //Print("##########################################In REDUCTION!########################################\n");
507    // check if sPolyList has any elements
508    // NOTE: due to initialization sPolyList always has a default NULL element
509    while(sPolyList->getLength() > 0) {
510        //Print("SPOLYLIST LENGTH: %d\n",sPolyList->getLength());
511        if(sPolyList->getLength() > 1) {
512        //Print("%p\n",sPolyList->getFirst());
513        //Print("%p\n",sPolyList->getFirst()->getLPoly());
514        //Print("%p\n",sPolyList->getFirst()->getNext());
515        //Print("%p\n",sPolyList->getFirst()->getNext()->getLPoly());
516        //Print("%p\n",sPolyList->getFirst());
517        }
518        //if(gPrev->getLast()->getIndex() == 5) {
519            //sPolyList->print();
520            //sleep(5);
521        //}
522       
523        // temp is the first element in the sPolyList which should be reduced
524        // due to earlier sorting this is the element of minimal degree AND
525        // minimal label
526        LNode* temp =   sPolyList->getFirst();
527        //pWrite(temp->getPoly());
528        // delete the above first element from sPolyList, temp will be either reduced to
529        // zero or added to gPrev, but never come back to sPolyList
530        if(NULL != temp && NULL != temp->getLPoly()) {
531            sPolyList->setFirst(temp->getNext());
532        //Print("HIER\n");
533        }
534        else {
535            break;
536        }
537        //Print("HALLO %p\n",temp->getPoly());
538        //Print("%p\n",temp->getPoly());
539        //pWrite(temp->getPoly());
540        //idShow(gbPrev);
541        poly tempNF = kNF(gbPrev,currQuotient,temp->getPoly());
542        //Print("LENGTH: %d\n",sPolyList->getLength());
543        //pWrite(tempNF);
544        //pWrite(temp->getPoly());
545        if(NULL != tempNF) {
546            pNorm(tempNF);
547            // write the reduced polynomial in temp
548            temp->setPoly(tempNF);
549            // try further reductions of temp with polynomials in gPrev
550            // with label index = current label index: this is done such that there
551            // is no label corruption during the reduction process
552            topReduction(temp,sPolyList,gPrev,rules,lTag,rTag,gbPrev);
553       
554        }
555        if(NULL != temp->getPoly()) {
556            //CList* newCritPairs = new CList;
557            //Print("##################IN CRITPAIRS IN REDUCTION#####################\n");
558            criticalPair(gPrev,critPairs,lTag,rTag,rules);
559            //Print("+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++H I E R++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n");
560        }
561        else {
562            //delete temp;
563            LNode* tempLoop = gPrev->getFirst();
564            //Print("AUSGABE IN REDUCTION:\n");       
565            //while(NULL != tempLoop) {
566                //pWrite(tempLoop->getPoly());
567                //tempLoop = tempLoop->getNext();
568            //}
569            //sleep(10);
570        }
571       
572    }
573}   
574
575
576
577/*
578=====================================================================================
579top reduction in F5, i.e. reduction of a given S-polynomial by labeled polynomials of
580the same index whereas the labels are taken into account
581=====================================================================================
582*/
583void topReduction(LNode* l, LList* sPolyList, LList* gPrev, RList* rules, LTagList* lTag, RTagList* rTag, ideal gbPrev) {
584    //Print("##########################################In topREDUCTION!########################################\n");
585    // try l as long as there are reductors found by findReductor()
586    do {
587        LNode* gPrevRedCheck    =   new LNode(lTag->getFirstCurrentIdx());
588        //Print("gPrevCheckPOLY: ");
589        //pWrite(gPrevRedCheck->getPoly());
590        LNode* tempRed          =   new LNode();
591        //Print("TESTED POLYNOMIAL IN THE FOLLOWING: ");
592        //pWrite(l->getPoly());
593        //Print("HIER\n");
594        tempRed  =   findReductor(l,gPrevRedCheck,gPrev,rules,lTag,rTag);
595        //Print("TEMPRED NEXT: %p\n",tempRed->getNext());
596        //Print("--------------------------------HIER DEBUG 2----------------------------------\n");
597        // if a reductor for l is found and saved in tempRed
598        if(NULL != tempRed) {
599            // if label of reductor is greater than the label of l we have to built a new element
600            // and add it to sPolyList
601            //Print("REDUCTOR POLYNOMIAL: ");
602            //pWrite(tempRed->getPoly());
603            //Print("TEMPRED: %p\n",tempRed);
604            //Print("TERM: ");
605            //pWrite(tempRed->getTerm());
606            if(pLmCmp(tempRed->getTerm(),l->getTerm()) == 1) {
607                // needed sinc pSub destroys the arguments!
608                //Print("H----------I------------E--------------R\n");
609                poly temp_poly_l    =   pInit();
610                temp_poly_l         =   pCopy(l->getPoly());
611                //Print("POLYNOMIAL L: ");
612                //pWrite(l->getPoly());
613                //pWrite(temp_poly_l);
614                poly temp           =   pSub(tempRed->getPoly(),temp_poly_l);
615                //Print("POLYNOMIAL L: ");
616                //pWrite(l->getPoly());
617                //pWrite(temp_poly_l);
618                //Print("AFTER REDUCTION STEP: ");
619                //pWrite(temp);
620                //sleep(20);
621                //pWrite(temp);
622                if(NULL != temp) {
623                    pNorm(temp);
624                    tempRed->setPoly(temp);
625                    //pWrite(tempRed->getPoly());
626                    // for debugging
627                    //pWrite(tempRed->getPoly());
628                    //Print("RULE ADDED\n");
629                    rules->insert(tempRed->getIndex(),tempRed->getTerm());
630                    tempRed->getLPoly()->setRule(rules->getFirst()->getRule());
631                    //Print("%p\n",sPolyList->getFirst());
632                    //Print("%p\n",sPolyList->getFirst()->getLPoly());
633                    //Print("SPOLYLIST LENGTH: %d\n",sPolyList->getLength());
634                    //sPolyList->print();
635                   
636                    sPolyList->insertByLabel(tempRed);
637                    //sPolyList->print();
638                    //Print("SPOLYLIST LENGTH: %d\n",sPolyList->getLength());
639                    //Print("OH JE\n");
640                }
641                else {
642                    pDelete(&temp);
643                    reductionsToZero++;
644                    //Print("RULE ADDED\n");
645                    //rules->insert(tempRed->getIndex(),tempRed->getTerm());
646                    //pWrite(rules->getFirst()->getRuleTerm());
647                    //Print("DELETE TEMPRED\n");
648                    delete tempRed;
649                }
650            }
651            // label of reductor is smaller than the label of l, subtract reductor from l and delete the
652            // gPrevRedCheck pointer added to l during findReductor() as the head term of l changes
653            // after subtraction
654            else {
655                poly temp_poly_l    =   pInit();
656                temp_poly_l         =   pCopy(l->getPoly());
657                poly temp   =   pSub(temp_poly_l,tempRed->getPoly());
658                //Print("AFTER REDUCTION STEP: ");
659                //pWrite(temp);
660                if(NULL != temp) {
661                    pNorm(temp);
662                    poly tempNF =   kNF(gbPrev,currQuotient,temp); 
663                    pNorm(tempNF);
664                    //pWrite(tempNF);
665                    if(NULL == tempNF) {
666                        reductionsToZero++;
667                        pDelete(&tempNF);
668                        l->setPoly(NULL);
669                        break;
670                    }
671                    l->setPoly(tempNF);
672                   
673                    //pWrite(l->getPoly());
674                    gPrevRedCheck   =   lTag->getFirstCurrentIdx();
675                }
676                else {
677                    //Print("ZERO REDUCTION!\n");
678                    reductionsToZero++;
679                    pDelete(&temp);
680                    l->setPoly(NULL);
681                    //pWrite(gPrev->getLast()->getPoly());
682                    break;
683                }
684            }   
685        }
686        else {
687            if(NULL != l->getPoly()) {
688                pNorm(l->getPoly());
689                //Print("----------------------------------ADDED TO GPREV IN TOPREDUCTION:-------------------------------------- ");
690                //pWrite(l->getPoly());
691                //pWrite(l->getTerm());
692                //Print("INDEX: %d\n\n\n", l->getIndex());
693                //sleep(20);
694                gPrev->insert(l->getLPoly());
695                //Print("GPREV: \n");
696                LNode* tempLoop = gPrev->getFirst();
697                //while(NULL != tempLoop) {
698                    //Print("HERE\n");
699                    //pWrite(tempLoop->getPoly());
700                    //tempLoop = tempLoop->getNext();
701                //}
702            }
703            break;
704        }
705    } while(1);
706}
707
708
709/*
710=====================================================================
711subalgorithm to find a possible reductor for the labeled polynomial l
712=====================================================================
713*/
714LNode* findReductor(LNode* l, LNode* gPrevRedCheck, LList* gPrev, RList* rules, LTagList* lTag,RTagList* rTag) {
715    // allociation of memory for the possible reductor
716    //Print("IN FIND REDUCTOR\n");
717    poly u      =   pOne();
718    poly red    =   pOne();
719    number nOne =   nInit(1);
720    LNode* temp =   new LNode();
721    // head term of actual element such that we do not have to call pHead at each new reductor test
722    poly t      =   pHead(l->getPoly());
723    // if l was already checked use the information in gPrevRedCheck such
724    // that we can start searching for new reducers from this point and
725    // not from the first element of gPrev with the current index
726    temp    =   gPrevRedCheck;
727    //temp    =   lTag->getFirstCurrentIdx(); 
728    //Print("GPREVREDCHECK: %p\n",gPrevRedCheck);
729    //pWrite(gPrevRedCheck->getPoly());
730    // search for reductors until we are at the end of gPrev resp. at the
731    // end of the elements of the current index
732    while(NULL != temp && temp->getIndex() == l->getIndex()) {
733        //pWrite(temp->getPoly());
734        //Print("INDEX: %d\n",temp->getIndex());
735        // does the head of the element of gPrev divides the head of
736        // the to be reduced element?
737        //Print("-------------FOUND REDUCTORS----------------------\n");
738        //Print("\n");
739        //pWrite(temp->getPoly());
740        //pWrite(temp->getTerm());
741        //pWrite(t);
742        //Print("HALLO\n");
743        if(pLmDivisibleByNoComp(temp->getPoly(),t)) {
744        //Print("HALLO\n");
745            // get all the information needed for the following tests
746            // of the criteria
747            u   =   pDivide(t,pHead(temp->getPoly()));
748            pSetCoeff(u,nOne);
749            //Print("HIER FINDRED\n");
750            //pWrite(u);
751            //Print("\n");
752            red =   ppMult_qq(u,temp->getPoly());
753            pNorm(red);
754            //u   =   ppMult_qq(u,temp->getTerm());
755            //pSetCoeff(u,nOne);
756            // check if both have the same label
757        //Print("HALLO\n");
758            if(pLmCmp(u,l->getTerm()) != 0) {
759        //Print("HALLO\n");
760                // passing criterion2 ?
761                if(!criterion2(u,temp,rules,rTag)) {
762                    // passing criterion1 ?
763                    if(!criterion1(gPrev,u,temp,lTag)) {
764                            //Print("HIER DEBUG\n");
765                            gPrevRedCheck   =   temp->getNext();
766                            LNode* redNode  =   new LNode(ppMult_qq(u,temp->getTerm()),temp->getIndex(),red,NULL,NULL);
767                            return redNode;
768                    }
769                }
770            }
771        }
772        //Print("%p\n",temp->getNext());
773        //pWrite(temp->getPoly());
774        //Print("HALLO\n");
775        temp = temp->getNext();
776    }
777   
778//    delete temp;
779   //Print("1st gPrev: ");
780    //pWrite(gPrev->getFirst()->getPoly());
781    //Print("2nd gPrev: ");
782    //pWrite(gPrev->getFirst()->getNext()->getPoly());
783 return NULL;
784}
785
786
787
788/*
789==========================================================================
790MAIN:computes a gb of the ideal i in the ring r with our F5 implementation
791==========================================================================
792*/
793ideal F5main(ideal id, ring r) {
794    int i,j,k,l;
795    int gbLength;
796    // 1 polynomial for defining initial labels & further tests
797    poly ONE = pOne();
798    poly pOne   =   pOne();
799    number nOne =   nInit(1);
800    // tag the first element of index i-1 for criterion 1
801    LTagList* lTag  =   new LTagList();
802    //Print("LTAG BEGINNING: %p\n",lTag);
803   
804    // first element in rTag is first element of rules which is NULL RNode,
805    // this must be done due to possible later improvements
806    RList* rules    =   new RList();
807    //Print("RULES FIRST: %p\n",rules->getFirst());
808    //Print("RULES FIRST DATA: %p\n",rules->getFirst()->getRule());
809    RTagList* rTag  =   new RTagList(rules->getFirst());
810    i = 1;
811    for(j=0; j<IDELEMS(id); j++) {
812        if(NULL != id->m[j]) { 
813            if(pComparePolys(id->m[j],ONE)) {
814                Print("One Polynomial in Input => Computations stopped\n");
815                ideal idNew = idInit(1,1);
816                idNew->m[0] = ONE;
817                return(idNew);
818            }   
819        }
820    } 
821    ideal idNew     =   kInterRed(id); 
822    id              =   idNew;
823    qsortDegree(&id->m[0],&id->m[IDELEMS(id)-1]);
824    idShow(id);
825    LList* gPrev    =   new LList(ONE, i, id->m[0]);
826    //idShow(id);
827    //Print("%p\n",id->m[0]);
828    //pWrite(id->m[0]);
829    //Print("%p\n",gPrev->getFirst()->getPoly());
830    //pWrite(gPrev->getFirst()->getPoly());
831
832    lTag->insert(gPrev->getFirst());
833    lTag->setFirstCurrentIdx(gPrev->getFirst());
834    // computing the groebner basis of the elements of index < actual index
835    gbLength    =   gPrev->getLength();
836    //Print("Laenge der bisherigen Groebner Basis: %d\n",gPrev->getLength());
837    ideal gbPrev    =   idInit(gbLength,1);
838    // initializing the groebner basis of elements of index < actual index
839    gbPrev->m[0]    =   gPrev->getFirst()->getPoly();
840    //idShow(gbPrev);
841    //idShow(currQuotient);
842
843    for(i=2; i<=IDELEMS(id); i++) {
844        LNode* gPrevTag =   gPrev->getLast();
845        //Print("Last POlynomial in GPREV: ");
846        //Print("%p\n",gPrevTag);   
847        //pWrite(gPrevTag->getPoly());
848        gPrev   =   F5inc(i, id->m[i-1], gPrev, gbPrev, ONE, lTag, rules, rTag);
849       
850        // DEBUGGING STUFF
851        LNode* temp    =   gPrev->getFirst();
852    // computing new groebner basis gbPrev
853        if(gPrev->getLength() > gbLength) {
854            ideal gbAdd =   idInit(gPrev->getLength()-gbLength,1);
855            LNode*  temp =   gPrevTag;
856        //Print("%p\n",gPrevTag);   
857        //Print("%p\n",gPrev->getLast());   
858        //pWrite(temp->getPoly());
859        //Print("LENGTH OF GPREV LIST: %d\n",gPrev->getLength());
860        //Print("%d\n",gbLength);
861        //Print("%d\n",gPrev->getLength()-gbLength-1);
862            for(j=0;j<=gPrev->getLength()-gbLength-1;j++) {
863                //Print("YES\n");
864                temp        =   temp->getNext();
865                if(temp) {
866                    //Print("%p\n",temp);
867                    //pWrite(temp->getPoly());
868                    //Print("%p\n",temp->getNext());
869                }
870                gbAdd->m[j] =   temp->getPoly();
871                //pWrite(temp->getPoly());
872            }
873            //Print("HIER AUCH\n");
874            gbPrev          =   idAdd(gbPrev,gbAdd);
875           
876           
877            // interreduction stuff
878           
879            if(i<IDELEMS(id)) {
880                ideal tempId    =   kInterRed(gbPrev);
881                gbPrev          =   tempId;
882                qsortDegree(&gbPrev->m[0],&gbPrev->m[IDELEMS(gbPrev)-1]);
883                delete gPrev;
884                //Print("RULES FIRST NOW1: %p\n",rules->getFirst());
885                //Print("HIER\n");
886                delete rules;
887                //delete rTag;
888                //Print("HIER AUCH\n");
889                //Print("%p\n",rules->getFirst());
890                gPrev    =   new LList(pOne,1,gbPrev->m[0]);
891                gPrev->insert(pOne,1,gbPrev->m[1]);
892                poly tempPoly = pInit();
893                pSetCoeff(tempPoly,nOne);
894                pLcm(pHead(gbPrev->m[0]),pHead(gbPrev->m[1]),tempPoly);
895                rules    =   new RList();
896               
897                //Print("%p\n",rules->getFirst());
898                //pWrite(tempPoly);
899                rules->insert(2,tempPoly);
900                rTag    =   new RTagList(rules->getFirst());
901                //Print("%p\n",rules->getFirst());
902                //Print("%p\n",rules->getFirst()->getNext()->getNext());
903                //Print("HIERLALA\n");
904            //pWrite(rules->getFirst()->getRuleTerm());
905           // Print("RULES FIRST NOW2: %p\n",rules->getFirst());
906                for(k=2; k<IDELEMS(gbPrev); k++) {
907                    gPrev->insert(pOne,k+1,gbPrev->m[k]);
908                    for(l=0; l<k; l++) {
909                        //pWrite(gbPrev->m[k]);
910                        //pWrite(gbPrev->m[l]);
911                        pLcm(pHead(gbPrev->m[k]),pHead(gbPrev->m[l]),tempPoly);
912                        pSetCoeff(tempPoly,nOne);
913                        rules->insert(k+1,tempPoly);
914                    }
915                    rTag->insert(rules->getFirst());
916                }
917            }
918           
919            gbLength    =   gPrev->getLength(); 
920
921        }
922        //gPrev->print();
923        //int anzahl  =   1;
924        //while(NULL != temp) {
925        //    Print("%d. Element: ",anzahl);
926        //    pWrite(temp->getPoly());
927        //    Print("\n");
928        //    temp    =   temp->getNext();
929        //    anzahl++;
930        //}
931        //sleep(5);
932        //Print("GROEBNER BASIS:\n====================================================\n");
933        //idShow(gbPrev);
934        //Print("===================================================\n");
935        //Print("JA\n");
936    } 
937        idShow(gbPrev);
938    Print("\n\nNumber of zero-reductions:  %d\n",reductionsToZero);
939    //LNode* temp    =   gPrev->getFirst();
940    //while(NULL != temp) {
941    //    pWrite(temp->getPoly());
942    //    temp    =   temp->getNext();
943    // }
944    //gPrev->print();
945    //delete lTag;
946    //delete rTag;
947    //delete gPrev;
948    return(gbPrev);
949
950
951}
952
953#endif
Note: See TracBrowser for help on using the repository browser.