source: git/kernel/f5gb.cc @ bb02ea

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