source: git/kernel/f5gb.cc @ 416ea2

spielwiese
Last change on this file since 416ea2 was 416ea2, checked in by Christian Eder, 15 years ago
new try on reduction procedures, no longer list of completed elements, adding a "last" node to LLists for more flexibility using the lists git-svn-id: file:///usr/local/Singular/svn/trunk@11401 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 20.4 KB
Line 
1/****************************************
2*  Computer Algebra System SINGULAR     *
3****************************************/
4/* $Id: f5gb.cc,v 1.29 2009-02-16 14:23:42 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    int j;
73    Print("%p\n",gPrev->getFirst());
74    pWrite(gPrev->getFirst()->getPoly());
75    gPrev->insert(ONE,i,f_i);
76    Print("1st gPrev: ");
77    pWrite(gPrev->getFirst()->getPoly());
78    Print("2nd gPrev: ");
79    pWrite(gPrev->getFirst()->getNext()->getPoly());
80    //pWrite(gPrev->getFirst()->getNext()->getPoly());   
81    CList* critPairs        =   new CList();
82    CNode* critPairsMinDeg  =   new CNode();   
83    // computation of critical pairs with checking of criterion 1 and criterion 2
84    critPairs               =   criticalPair(gPrev, critPairs, lTag, rTag, rules);
85    static LList* sPolyList        =   new LList();
86    // labeled polynomials which have passed reduction() and have to be added to list gPrev
87    static LList* completed        =   new LList();
88    Print("_____________________________________________________________________________%p\n",completed->getFirst()->getLPoly());
89    // the reduced labeled polynomials which are returned from subalgorithm reduction()
90    static LList* reducedLPolys     =   new LList();
91    // while there are critical pairs to be further checked and deleted/computed
92    while(NULL != critPairs->getFirst()->getData()) { 
93        // critPairs->getMinDeg() deletes the first elements of minimal degree from
94        // critPairs, thus the while loop is not infinite.
95        critPairs->print();
96        critPairsMinDeg =   critPairs->getMinDeg();
97        // adds all to be reduced S-polynomials in the list sPolyList and adds
98        // the corresponding rules to the list rules
99        // NOTE: inside there is a second check of criterion 2 if new rules are
100        // added
101        computeSPols(critPairsMinDeg,rTag,rules,sPolyList);
102       
103        // DEBUG STUFF FOR SPOLYLIST
104        LNode* temp     =   sPolyList->getFirst();
105        while(NULL != temp && NULL != temp->getLPoly()) {
106            Print("Spolylist element: ");
107            pWrite(temp->getPoly());
108            temp    =   temp->getNext();
109        } 
110        //while(NULL != sPolyList->getFirst()->getLPoly()) {
111            reduction(sPolyList, completed, gPrev, rules, lTag, rTag, gbPrev);
112        //}
113    }
114    Print("HIER123\n");
115    Print("%p\n",rules->getFirst());
116    Print("%p\n",rTag->getFirst());
117    if(rules->getFirst() != rTag->getFirst()) {
118        Print("+++++++++++++++++++++++++++++++++++++NEW RULES+++++++++++++++++++++++++++++++++++++\n");
119        rTag->insert(rules->getFirst());
120    }
121    else {
122        Print("+++++++++++++++++++++++++++++++++++NO NEW RULES++++++++++++++++++++++++++++++++++++\n");
123    }
124    lTag->insert(gPrev->getFirst());
125    Print("COMPLETED FIRST IN F5INC: \n");
126    return gPrev;
127}
128
129
130
131/*
132================================================================
133computes a list of critical pairs for the next reduction process
134first element in gPrev is always the newest element which must
135build critical pairs with all other elements in gPrev
136================================================================
137*/
138CList* criticalPair(LList* gPrev, CList* critPairs, LTagList* lTag, RTagList* rTag, RList* rules) {
139    // initialization for usage in pLcm()
140    number nOne         =   nInit(1);
141    LNode* first        =   gPrev->getFirst();
142    LNode* l            =   first->getNext();
143    poly* u1            =   new poly(pInit());
144    poly* u2            =   new poly(pInit());
145    poly* lcm           =   new poly(pInit());
146    poly t              =   pHead(first->getPoly());
147    // computation of critical pairs
148    while( NULL != l) {
149        //pWrite( *(gPrev->getFirst()->getPoly()) );
150        //pWrite( *(l->getPoly()) );
151        pLcm(first->getPoly(), l->getPoly(), *lcm);
152        pSetCoeff(*lcm,nOne);
153        // computing factors u2 for new labels
154        pWrite(t);
155        *u1 = pDivide(*lcm,t);
156        pWrite(*u1);
157        pSetCoeff(*u1,nOne);
158        pWrite(pHead(l->getPoly()));
159        *u2 = pDivide(*lcm, pHead(l->getPoly()));
160        pSetCoeff(*u2,nOne);
161        Print("IN CRITPAIRS\n");
162        pWrite(*u2);
163        // testing both new labels by the F5 Criterion
164        if(!criterion1(u1, first, lTag) && !criterion1(u2, l, lTag) && 
165           !criterion2(u1, first, rules, rTag) && !criterion2(u2, l, rules, rTag)) {
166            // if they pass the test, add them to CList critPairs, having the LPoly with greater
167            // label as first element in the CPair
168            if(first->getIndex() == l->getIndex() && 
169            pMult(*u1, first->getTerm()) < pMult(*u2, l->getTerm())) {
170                CPair* cp   =   new CPair(pDeg(ppMult_qq(*u2,pHead(l->getPoly()))), *u2, 
171                                l->getLPoly(), *u1, first->getLPoly());                   
172                    critPairs->insert(cp);
173            }
174            else {
175                CPair* cp   =   new CPair(pDeg(ppMult_qq(*u2,pHead(l->getPoly()))), *u1, 
176                                first->getLPoly(), *u2, l->getLPoly());                   
177                    critPairs->insert(cp);
178            }
179        }
180        else {
181        }
182       
183        // for debugging
184        if(NULL != critPairs) {
185            critPairs->print(); 
186        }
187        l   =   l->getNext();
188    }
189    return critPairs;
190}
191
192
193
194/*
195========================================
196Criterion 1, i.e. Faugere's F5 Criterion
197========================================
198*/
199bool criterion1(poly* t, LNode* l, LTagList* lTag) {
200    // starts at the first element in gPrev with index = (index of l)-1, these tags are saved in lTag
201        LNode* testNode =   lTag->get(l->getIndex());
202        // save the monom t1*label_term(l) as it is tested various times in the following
203    poly u1 = ppMult_qq(*t,l->getTerm());
204    while(NULL != testNode && NULL != testNode->getLPoly()) {
205        if(pLmDivisibleByNoComp(testNode->getPoly(),u1)) {
206            Print("Criterion 1 NOT passed!\n");
207            return true;
208        }
209        //pWrite(testNode->getNext()->getPoly());
210                testNode    =   testNode->getNext();
211        Print("ADDRESS OF TEST NODE: %p\n",testNode);
212    Print("Hier\n");
213    }
214    Print("HIER DRIN CRITERION 1\n");
215       
216    return false;
217}
218
219
220
221/*
222=====================================
223Criterion 2, i.e. Rewritten Criterion
224=====================================
225*/
226bool criterion2(poly* t, LNode* l, RList* rules, RTagList* rTag) {
227        Print("HIER DRIN CRITERION 2:=========================\n");
228    // start at the previously added element to gPrev, as all other elements will have the same index for sure
229        Print("%d\n",l->getIndex());
230    RNode* testNode =   new RNode();
231    if(NULL == rTag->getFirst()->getRule()) {
232        testNode    =   rules->getFirst();
233    }
234    else {
235        Print("%d\n",l->getIndex());
236        Print("%d\n",rTag->getFirst()->getRuleIndex());
237pWrite(rTag->getFirst()->getRuleTerm());
238        if(l->getIndex() > rTag->getFirst()->getRuleIndex()) {
239            testNode    =   rules->getFirst();
240            Print("TEST NODE: %p\n",testNode);
241        }
242        else {
243            testNode    =   rTag->get(l->getIndex());
244        }
245    }
246        // save the monom t1*label_term(l) as it is tested various times in the following
247    poly u1 = ppMult_qq(*t,l->getTerm());
248        pWrite(u1);
249    // first element added to rTag was NULL, check for this
250    pWrite(l->getPoly());
251    //Print("%p\n",testNode->getRule());
252    Print("HIER !!!!\n");
253    // NOTE: testNode is possibly NULL as rTag->get() returns NULL for elements of index <=1!
254    while(NULL != testNode && NULL != testNode->getRule() && testNode->getRule() != l->getRule() 
255          && l->getIndex() == testNode->getRuleIndex()) {
256                pWrite(ppMult_qq(*t,l->getTerm()));
257                pWrite(testNode->getRuleTerm());
258                if(pLmDivisibleBy(ppMult_qq(*t,l->getTerm()),testNode->getRuleTerm())) {
259            Print("Criterion 2 NOT passed!\n");
260            return true;
261        }
262                testNode    =   testNode->getNext();
263    }
264    return false;
265}
266
267
268
269/*
270=================================================================================================================
271Criterion 2, i.e. Rewritten Criterion, for its second call in computeSPols(), with added lastRuleTested parameter
272=================================================================================================================
273*/
274bool criterion2(poly* t, LPoly* l, RList* rules, Rule* lastRuleTested) {
275    // start at the previously added element to gPrev, as all other elements will have the same index for sure
276        RNode* testNode =   rules->getFirst();
277    // save the monom t1*label_term(l) as it is tested various times in the following
278    poly u1 = ppMult_qq(*t,l->getTerm());
279        // first element added to rTag was NULL, check for this
280        while(NULL != testNode->getRule() && l->getRule() != lastRuleTested) {
281        if(pLmDivisibleBy(testNode->getRuleTerm(),ppMult_qq(*t,l->getTerm()))) {
282            Print("Criterion 2 NOT passed!\n");
283            return true;
284        }
285                testNode    =   testNode->getNext();
286    }
287    return false;
288}
289
290
291
292/*
293==================================
294Computation of S-Polynomials in F5
295==================================
296*/
297void computeSPols(CNode* first, RTagList* rTag, RList* rules, LList* sPolyList) { 
298    CNode* temp         =   first;
299    poly sp      =   pInit();
300    Print("###############################IN SPOLS##############################\n");
301    while(NULL != temp->getData()) {
302        // only if a new rule was added since the last test in subalgorithm criticalPair()
303        //if(rules->getFirst() != rTag->getFirst()) {
304            Print("IN SPOLS => NEW RULES AVAILABLE\n");
305            if(!criterion2(temp->getAdT1(),temp->getAdLp1(),rules,temp->getLastRuleTested())) {
306                // the second component is tested only when it has the actual index, otherwise there is
307                // no new rule to test since the last test in subalgorithm criticalPair()
308                if(temp->getLp2Index() == temp->getLp1Index()) {
309                    if(!criterion2(temp->getAdT2(),temp->getAdLp2(),rules,temp->getLastRuleTested())) {
310                        // computation of S-polynomial
311                        sp      =   pSub(ppMult_qq(temp->getT1(),temp->getLp1Poly()),
312                                         ppMult_qq(temp->getT2(),temp->getLp2Poly()));
313                        Print("BEGIN SPOLY1\n====================\n");
314                        pWrite(sp);
315                        Print("END SPOLY1\n====================\n");
316                        if(NULL == sp) {
317                            // as rules consist only of pointers we need to save the labeled
318                            // S-polynomial also of a zero S-polynomial
319                            //zeroList->insert(temp->getAdT1(),temp->getLp1Index(),&sp);
320                            // origin of rule can be set NULL as the labeled polynomial
321                            // will never be used again as it is zero => no problems with
322                            // further criterion2() tests and termination conditions
323                            Print("~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ZERO REDUCTION~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~\n");
324                                                        reductionsToZero++;
325                            rules->insert(temp->getLp1Index(),temp->getT1());
326                            // as sp = NULL, delete it
327                            delete(&sp);
328                        }
329                        else {
330                            rules->insert(temp->getLp1Index(),temp->getT1());
331                            sPolyList->insert(temp->getT1(),temp->getLp1Index(),sp,rules->getFirst()->getRule());
332                        }
333                        // data is saved in sPolyList or zero => delete sp
334                    }
335                }
336                else { // temp->getLp2Index() < temp->getLp1Index()
337                    // computation of S-polynomial
338                    sp      =   pSub(ppMult_qq(temp->getT1(),temp->getLp1Poly()),
339                                     ppMult_qq(temp->getT2(),temp->getLp2Poly()));
340                    Print("BEGIN SPOLY2\n====================\n");
341                    pWrite(sp);
342                    Print("END SPOLY2\n====================\n");
343                    if(NULL == sp) {
344                        // zeroList->insert(temp->getAdT1(),temp->getLp1Index(),&sp);
345                        // origin of rule can be set NULL as the labeled polynomial
346                        // will never be used again as it is zero => no problems with
347                        // further criterion2() tests and termination conditions
348                            Print("~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ZERO REDUCTION~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~\n");
349                        reductionsToZero++;
350                        rules->insert(temp->getLp1Index(),temp->getT1());
351                        // as sp = NULL, delete it
352                        delete(&sp);
353                    }
354                    else {
355                        rules->insert(temp->getLp1Index(),temp->getT1());
356                        sPolyList->insert(temp->getT1(),temp->getLp1Index(),sp,rules->getFirst()->getRule());
357                    }
358                    // data is saved in sPolyList or zero => delete sp
359                }
360            }
361        //}
362        temp    =   temp->getNext();
363    }
364    // these critical pairs can be deleted now as they are either useless for further computations or
365    // already saved as an S-polynomial to be reduced in the following
366    //pDelete(&sp);
367    delete first;   
368}
369
370
371
372/*
373========================================================================
374reduction including subalgorithm topReduction() using Faugere's criteria
375========================================================================
376*/
377void reduction(LList* sPolyList, LList* completed, LList* gPrev, RList* rules, LTagList* lTag, RTagList* rTag, 
378                 ideal gbPrev) { 
379    Print("##########################################In REDUCTION!########################################\n");
380     // LNode* saving the last element in gPrev before the completed elements from the reduction process are added
381    // this is needed to get the new critical pairs computed
382    LNode* gPrevTag = gPrev->getLast();
383    //for debugging
384    pWrite(gPrevTag->getPoly());
385    // check if sPolyList has any elements
386    // NOTE: due to initialization sPolyList always has a default NULL element
387    while(sPolyList->getLength() > 1) {
388        // temp is the first element in the sPolyList which should be reduced
389        // due to earlier sorting this is the element of minimal degree AND
390        // minimal label
391        LNode* temp =   sPolyList->getFirst();
392        // delete the above first element from sPolyList, temp will be either reduced to
393        // zero or added to gPrev, but never come back to sPolyList
394        sPolyList->setFirst(temp->getNext());
395        pWrite(temp->getPoly());
396        poly tempNF = kNF(gbPrev,currQuotient,temp->getPoly());
397        Print("LENGTH: %d\n",sPolyList->getLength());
398        pWrite(tempNF);
399        pWrite(temp->getPoly());
400        if(NULL != tempNF) {
401            // write the reduced polynomial in temp
402            temp->setPoly(tempNF);
403            // try further reductions of temp with polynomials in gPrev
404            // with label index = current label index: this is done such that there
405            // is no label corruption during the reduction process
406            topReduction(temp,completed,gPrev,rules,lTag,rTag);
407        }
408    }
409}   
410
411
412
413/*
414=====================================================================================
415top reduction in F5, i.e. reduction of a given S-polynomial by labeled polynomials of
416the same index whereas the labels are taken into account
417=====================================================================================
418*/
419void topReduction(LNode* l, LList* completed, LList* gPrev, RList* rules, LTagList* lTag, RTagList* rTag) {
420    Print("##########################################In topREDUCTION!########################################\n");
421
422}
423
424
425/*
426=====================================================================
427subalgorithm to find a possible reductor for the labeled polynomial l
428=====================================================================
429*/
430LNode* findReductor(LNode* l,LList* completed,LList* gPrev, RList* rules, LTagList* lTag,RTagList* rTag,
431                    LNode* gPrevRedCheck) {
432   Print("HIER DU NULL!\n");
433    return NULL;
434}
435
436
437
438/*
439==========================================================================
440MAIN:computes a gb of the ideal i in the ring r with our F5 implementation
441==========================================================================
442*/
443ideal F5main(ideal id, ring r) {
444    int i,j;
445    int gbLength;
446    // 1 polynomial for defining initial labels & further tests
447    poly ONE = pOne();
448    // tag the first element of index i-1 for criterion 1
449    LTagList* lTag  =   new LTagList();
450    Print("LTAG BEGINNING: %p\n",lTag->getFirst());
451   
452    // first element in rTag is first element of rules which is NULL RNode,
453    // this must be done due to possible later improvements
454    RList* rules    =   new RList();
455    RTagList* rTag  =   new RTagList(rules->getFirst());
456    i = 1;
457    for(j=0; j<IDELEMS(id); j++) {
458        if(NULL != id->m[j]) { 
459            if(pComparePolys(id->m[j],ONE)) {
460                Print("One Polynomial in Input => Computations stopped\n");
461                ideal idNew = idInit(1,1);
462                idNew->m[0] = ONE;
463                return(idNew);
464            }   
465        }
466    } 
467    LList* gPrev    =   new LList(ONE, i, id->m[0]);
468    Print("%p\n",id->m[0]);
469    pWrite(id->m[0]);
470    Print("%p\n",gPrev->getFirst()->getPoly());
471    pWrite(gPrev->getFirst()->getPoly());
472
473    lTag->insert(gPrev->getFirst());
474    // computing the groebner basis of the elements of index < actual index
475    gbLength    =   gPrev->getLength();
476    Print("Laenge der bisherigen Groebner Basis: %d\n",gPrev->getLength());
477    ideal gbPrev    =   idInit(gbLength,1);
478    // initializing the groebner basis of elements of index < actual index
479    gbPrev->m[0]    =   gPrev->getFirst()->getPoly();
480    idShow(gbPrev);
481    idShow(currQuotient);
482
483    for(i=2; i<=IDELEMS(id); i++) {
484        gPrev   =   F5inc(i, id->m[i-1], gPrev, gbPrev, ONE, lTag, rules, rTag);
485        // comuting new groebner basis gbPrev
486        ideal gbAdd =   idInit(gPrev->getLength()-gbLength,1);
487        LNode*  temp    =   gPrev->getFirst();
488        for(j=0;j<=gPrev->getLength()-gbLength-1;j++) {
489            gbAdd->m[j] =   temp->getPoly();
490            temp        =   temp->getNext();
491        }
492        gbPrev  =   idAdd(gbPrev,gbAdd);
493        Print("GROEBNER BASIS:\n====================================================\n");
494        idShow(gbPrev);
495        Print("===================================================\n");
496
497        Print("JA\n");
498    } 
499    Print("\n\nNumber of zero-reductions:  %d\n",reductionsToZero);
500    return(gbPrev);
501
502
503}
504
505#endif
Note: See TracBrowser for help on using the repository browser.