source: git/kernel/f5gb.cc @ 7bd779

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