source: git/kernel/f5gb.cc @ d92713

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