Changeset 197a82 in git for kernel/gfan.cc


Ignore:
Timestamp:
Nov 17, 2009, 5:15:13 PM (14 years ago)
Author:
Martin Monerjan
Branches:
(u'spielwiese', 'd0474371d8c5d8068ab70bfb42719c97936b18a6')
Children:
53baa975e55756198c3de9aae898e792473ce1c2
Parents:
87353a6b9eb171d12e73ede9e24e8d8e68eaf111
Message:
removed reverse search related stuff
idrCopyR_NoSort in lprepareResult
Check whether input has no local ordering


git-svn-id: file:///usr/local/Singular/svn/trunk@12288 2c84dea3-7e68-4137-9b89-c4e89433aadc
File:
1 edited

Legend:

Unmodified
Added
Removed

  • kernel/gfan.cc

    r87353a6 r197a82  
    14111411//NOTE: Should be replaced by kNF or kNF2
    14121412//NOTE: Done
    1413 poly gcone::restOfDiv(poly const &f, ideal const &I)
    1414 {
    1415 //                      cout << "Entering restOfDiv" << endl;
    1416 //      poly p=f;
    1417 //                      //pWrite(p);                   
    1418 //      poly r=NULL;    //The zero polynomial
    1419 //      int ii;
    1420 //      bool divOccured;
    1421 //                     
    1422 //      while (p!=NULL)
    1423 //      {
    1424 //              ii=1;
    1425 //              divOccured=FALSE;
    1426 //                             
    1427 //              while( (ii<=IDELEMS(I) && (divOccured==FALSE) ))
    1428 //              {                                       
    1429 //                      if (pDivisibleBy(I->m[ii-1],p)) //does LM(I->m[ii]) divide LM(p) ?
    1430 //                      {                                               
    1431 //                              poly step1,step2,step3;
    1432 //                                              //cout << "dividing "; pWrite(pHead(p));cout << "by ";pWrite(pHead(I->m[ii-1])); cout << endl;
    1433 //                              step1 = pDivideM(pHead(p),pHead(I->m[ii-1]));
    1434 //                                              //cout << "LT(p)/LT(f_i)="; pWrite(step1); cout << endl;                               
    1435 //                              step2 = ppMult_qq(step1, I->m[ii-1]);                                           
    1436 //                              step3 = pSub(pCopy(p), step2);
    1437 //                                              //p=pSub(p,pMult( pDivide(pHead(p),pHead(I->m[ii])), I->m[ii]));                       
    1438 //                                              //pSetm(p);
    1439 //                              pSort(step3); //must be here, otherwise strange behaviour with many +o+o+o+o+ terms
    1440 //                              p=step3;
    1441 //                                              //pWrite(p);                                           
    1442 //                              divOccured=TRUE;
    1443 //                      }
    1444 //                      else
    1445 //                      {
    1446 //                              ii += 1;
    1447 //                      }//if (pLmDivisibleBy(I->m[ii],p,currRing))
    1448 //              }//while( (ii<IDELEMS(I) && (divOccured==FALSE) ))
    1449 //              if (divOccured==FALSE)
    1450 //              {
    1451 //                                      //cout << "TICK 5" << endl;
    1452 //                      r=pAdd(pCopy(r),pHead(p));
    1453 //                      pSetm(r);
    1454 //                      pSort(r);
    1455 //                                      //cout << "r="; pWrite(r); cout << endl;
    1456 //                                     
    1457 //                      if (pLength(p)!=1)
    1458 //                      {
    1459 //                              p=pSub(pCopy(p),pHead(p));      //Here it may occur that p=0 instead of p=NULL
    1460 //                      }
    1461 //                      else
    1462 //                      {
    1463 //                              p=NULL; //Hack to correct this situation                                               
    1464 //                      }                                       
    1465 //                                      //cout << "p="; pWrite(p);
    1466 //              }//if (divOccured==FALSE)
    1467 //      }//while (p!=0)
    1468 //      return r;
    1469 }//poly restOfDiv(poly const &f, ideal const &I)
     1413//NOTE: removed with r12286
    14701414               
    14711415/** \brief Compute \f$ f-f^{\mathcal{G}} \f$
    … …  
    14831427                res->m[ii]=kNF(G, NULL,H->m[ii],0,0);
    14841428                temp1=H->m[ii];
    1485                 temp2=res->m[ii];                               
     1429                temp2=res->m[ii];
    14861430                temp3=pSub(temp1, temp2);
    14871431                res->m[ii]=temp3;
    … …  
    14891433                //pSort(res->m[ii]);
    14901434                //pSetm(res->m[ii]);
    1491                 //cout << "res->m["<<ii<<"]=";pWrite(res->m[ii]);                                              
     1435                //cout << "res->m["<<ii<<"]=";pWrite(res->m[ii]);       
    14921436        }                       
    14931437        return res;
    … …  
    17171661 * If this is the case, then our facet is indeed a search facet and TRUE is retuned.
    17181662*/
    1719 bool gcone::isSearchFacet(gcone &gcTmp, facet *testfacet)
    1720 {                               
    1721         ring actRing=currRing;
    1722         facet *facetPtr=(facet*)gcTmp.facetPtr;                 
    1723         facet *fMin=new facet(*facetPtr);       //Pointer to the "minimal" facet
    1724                         //facet *fMin = new facet(tmpcone.facetPtr);
    1725                         //fMin=tmpcone.facetPtr;                //Initialise to first facet of tmpcone
    1726         facet *fAct;    //Ptr to alpha_i
    1727         facet *fCmp;    //Ptr to alpha_j
    1728         fAct = fMin;
    1729         fCmp = fMin->next;                             
    1730                        
    1731         rChangeCurrRing(this->rootRing);        //because we compare the monomials in the rootring                     
    1732         poly p=pInit();
    1733         poly q=pInit();
    1734         intvec *alpha_i = new intvec(this->numVars);                   
    1735         intvec *alpha_j = new intvec(this->numVars);
    1736         intvec *sigma = new intvec(this->numVars);
    1737         sigma=gcTmp.getIntPoint();
    1738                        
    1739         int *u=(int *)omAlloc((this->numVars+1)*sizeof(int));
    1740         int *v=(int *)omAlloc((this->numVars+1)*sizeof(int));
    1741         u[0]=0; v[0]=0;
    1742         int weight1,weight2;
    1743         while(fAct->next->next!=NULL)   //NOTE this is ugly. Can it be done without fCmp?
    1744         {
    1745                 /* Get alpha_i and alpha_{i+1} */
    1746                 alpha_i=fAct->getFacetNormal();
    1747                 alpha_j=fCmp->getFacetNormal();
    1748 #ifdef gfan_DEBUG
    1749                 alpha_i->show();
    1750                 alpha_j->show();
    1751 #endif
    1752                 /*Compute the dot product of sigma and alpha_{i,j}*/
    1753                 weight1=dotProduct(sigma,alpha_i);
    1754                 weight2=dotProduct(sigma,alpha_j);
    1755 #ifdef gfan_DEBUG
    1756                 cout << "weight1=" << weight1 << " " << "weight2=" << weight2 << endl;
    1757 #endif
    1758                 /*Adjust alpha_i and alpha_i+1 accordingly*/
    1759                 for(int ii=1;ii<=this->numVars;ii++)
    1760                 {                                       
    1761                         u[ii]=weight1*(*alpha_i)[ii-1];
    1762                         v[ii]=weight2*(*alpha_j)[ii-1];
    1763                 }                               
    1764                                
    1765                 /*Now p_weight and q_weight need to be compared as exponent vectors*/
    1766                 pSetCoeff0(p,nInit(1));
    1767                 pSetCoeff0(q,nInit(1));
    1768                 pSetExpV(p,u);
    1769                 pSetm(p);                       
    1770                 pSetExpV(q,v);
    1771                 pSetm(q);
    1772 #ifdef gfan_DEBUG                               
    1773                 pWrite(p);pWrite(q);
    1774 #endif 
    1775                                 /*We want to check whether x^p < x^q
    1776                 => want to check for return value 1 */
    1777                 if (pLmCmp(p,q)==1)     //i.e. x^q is smaller
    1778                 {
    1779                         fMin=fCmp;
    1780                         fAct=fMin;
    1781                         fCmp=fCmp->next;
    1782                 }
    1783                 else
    1784                 {
    1785                                         //fAct=fAct->next;
    1786                         if(fCmp->next!=NULL)
    1787                         {
    1788                                 fCmp=fCmp->next;
    1789                         }
    1790                         else
    1791                         {
    1792                                 fAct=fAct->next;
    1793                         }
    1794                 }
    1795                                 //fAct=fAct->next;
    1796         }//while(fAct.facetPtr->next!=NULL)
    1797         delete alpha_i,alpha_j,sigma;
    1798                        
    1799         /*If testfacet was minimal then fMin should still point there */
    1800                        
    1801                         //if(fMin->getFacetNormal()==ivNeg(testfacet.getFacetNormal()))                 
    1802 #ifdef gfan_DEBUG
    1803         cout << "Checking for parallelity" << endl <<" fMin is";
    1804         fMin->printNormal();
    1805         cout << "testfacet is ";
    1806         testfacet->printNormal();
    1807         cout << endl;
    1808 #endif
    1809         if (fMin==gcTmp.facetPtr)                       
    1810                         //if(areEqual(fMin->getFacetNormal(),ivNeg(testfacet.getFacetNormal())))
    1811                         //if (isParallel(fMin->getFacetNormal(),testfacet->getFacetNormal()))
    1812         {                               
    1813                 cout << "Parallel" << endl;
    1814                 rChangeCurrRing(actRing);
    1815                                 //delete alpha_min, test;
    1816                 return TRUE;
    1817         }
    1818         else
    1819         {
    1820                 cout << "Not parallel" << endl;
    1821                 rChangeCurrRing(actRing);
    1822                                 //delete alpha_min, test;
    1823                 return FALSE;
    1824         }
    1825 }//bool isSearchFacet
     1663//removed with r12286
    18261664               
    18271665/** \brief Check for equality of two intvecs
    … …  
    18431681/** \brief The reverse search algorithm
    18441682 */
    1845 void gcone::reverseSearch(gcone *gcAct) //no const possible here since we call gcAct->flip
    1846 {
    1847 //      facet *fAct=new facet();
    1848 //      fAct = gcAct->facetPtr;                 
    1849 //                     
    1850 //      while(fAct!=NULL)
    1851 //      {
    1852 //              cout << "======================"<< endl;
    1853 //              gcAct->flip(gcAct->gcBasis,gcAct->facetPtr);
    1854 //              gcone *gcTmp = new gcone(*gcAct);
    1855 //                              //idShow(gcTmp->gcBasis);
    1856 //              gcTmp->getConeNormals(gcTmp->gcBasis, TRUE);
    1857 // #ifdef gfan_DEBUG
    1858 //              facet *f = new facet();
    1859 //              f=gcTmp->facetPtr;
    1860 //              while(f!=NULL)
    1861 //              {
    1862 //                      f->printNormal();
    1863 //                      f=f->next;                                     
    1864 //              }
    1865 // #endif
    1866 //              gcTmp->showIntPoint();
    1867 //              /*recursive part goes gere*/
    1868 //              if (isSearchFacet(*gcTmp,(facet*)gcAct->facetPtr))
    1869 //              {
    1870 //                      gcAct->next=gcTmp;
    1871 //                      cout << "PING"<< endl;
    1872 //                      reverseSearch(gcTmp);
    1873 //              }
    1874 //              else
    1875 //              {
    1876 //                      delete gcTmp;
    1877 //                      /*NOTE remove fAct from linked list. It's no longer needed*/
    1878 //              }
    1879 //              /*recursion ends*/
    1880 //              fAct = fAct->next;             
    1881 //      }//while(fAct->next!=NULL)
    1882 }//reverseSearch
     1683//removed with r12286
    18831684               
    18841685/** \brief The new method of Markwig and Jensen
    … …  
    29072708                intvec iv=(gcAct->f2M(gcAct,gcAct->facetPtr));
    29082709                l->m[1].data=(void*)ivCopy(&iv);
    2909                 l->m[2].data=(void*)idCopy(gcAct->gcBasis);
     2710                l->m[2].data=(void*)idrCopyR_NoSort(gcAct->gcBasis,gcAct->getBaseRing());
    29102711                l->m[3].data=(void*)(gcAct->getBaseRing());
    29112712//              l->m[4].data=(void*)0;
    … …  
    29142715        }               
    29152716       
    2916 //      while( gcAct!=NULL )
    2917 //      {
    2918 //              l->m[0].data=(void*)gc->getUCN();
    2919 //              l->m[1].data=(intvec*)(&gcAct->f2M(gcAct,gcAct->facetPtr));
    2920 //              gcAct = gcAct->next;
    2921 //              lAdd((sleftv*)res,(sleftv*)res,(sleftv*)l);
    2922 //      }
    2923 //      char *v=lString(res);
    2924 //      cout << v;
    29252717        return res;
    29262718}
    … …  
    29632755lists gfan(ideal inputIdeal, int h)
    29642756{
    2965         lists lResList;
     2757        lists lResList; //this is the object we return
     2758       
     2759        if(rHasGlobalOrdering(currRing))
     2760        {       
    29662761        int numvar = pVariables;
    29672762        gfanHeuristic = h;
    … …  
    29742769        searchMethod method;
    29752770        method = noRevS;
    2976 //      method = reverseSearch;
     2771
    29772772       
    29782773#ifdef gfan_DEBUG
    … …  
    29822777        ring rootRing;                  // The ring associated to the target ordering
    29832778        ideal res;     
    2984         facet *fRoot;
    2985        
    2986         if (method==reverseSearch)
    2987         {
    2988        
    2989                 /* Construct a new ring which will serve as our root*/
    2990                 rootRing=rCopy0(currRing);
    2991                 rootRing->order[0]=ringorder_lp;
    2992                
    2993                 rComplete(rootRing);
    2994                 rChangeCurrRing(rootRing);
    2995                
    2996                 /* Fetch the inputIdeal into our rootRing */
    2997                 map theMap=(map)idMaxIdeal(1);  //evil hack!
    2998                 theMap->preimage=NULL;  //neccessary?
    2999                 ideal rootIdeal;
    3000                 rootIdeal=fast_map(inputIdeal,inputRing,(ideal)theMap, currRing);
    3001 #ifndef NDEBUG
    3002         #ifdef gfan_DEBUG
    3003                 cout << "Root ideal is " << endl;
    3004                 idShow(rootIdeal);
    3005                 cout << "The root ring is " << endl;
    3006                 rWrite(rootRing);
    3007                 cout << endl;
    3008         #endif
    3009 #endif
    3010                
    3011                 //gcone *gcRoot = new gcone();  //Instantiate the sink
    3012                 gcone *gcRoot = new gcone(rootRing,rootIdeal);
    3013                 gcone *gcAct;
    3014                 gcAct = gcRoot;
    3015                 gcAct->numVars=pVariables;
    3016                 gcAct->getGB(rootIdeal);        //sets gcone::gcBasis
    3017 #ifndef NDEBUG
    3018                 idShow(gcAct->gcBasis);
    3019 #endif
    3020                 gcAct->getConeNormals(gcAct->gcBasis);  //hopefully compute the normals
    3021                 //gcAct->flip(gcAct->gcBasis,gcAct->facetPtr); 
    3022                 /*Now it is time to compute the search facets, respectively start the reverse search.
    3023                 But since we are in the root all facets should be search facets. IS THIS TRUE?
    3024                 NOTE: Check for flippability is not very sophisticated
    3025                 */     
    3026                 //gcAct->reverseSearch(gcAct); 
    3027                 rChangeCurrRing(rootRing);
    3028                 res=gcRoot->gcBasis;   
    3029         }//if method==reverSearch
     2779        facet *fRoot;   
    30302780       
    30312781        if(method==noRevS)
    … …  
    30352785                gcAct = gcRoot;
    30362786                gcAct->numVars=pVariables;
    3037                 gcAct->getGB(inputIdeal);
     2787                gcAct->getGB(inputIdeal);               
     2788#ifndef NDEBUG
    30382789                cout << "GB of input ideal is:" << endl;
    3039 #ifndef NDEBUG
    30402790                idShow(gcAct->gcBasis);
    30412791#endif
    … …  
    30522802                //Below is a workaround, since gcAct->gcBasis gets deleted in noRevS
    30532803                res = inputIdeal;
    3054 //              intvec mRes=gcRoot->f2M(gcRoot,gcRoot->facetPtr);               
    30552804                lResList=lprepareResult(gcRoot,gcRoot->getCounter());
    30562805//              res=lResList;
    30572806        }
    30582807        dd_free_global_constants();
     2808        }//rHasGlobalOrdering
     2809        else
     2810        {
     2811                WerrorS("Ring has non-global ordering - terminating");
     2812                lResList->Init(1);
     2813                lResList->m[0].rtyp=INT_CMD;
     2814                int ires=0;
     2815                lResList->m[0].data=(void*)&ires;
     2816        }
    30592817       
    3060         /*As of now extra.cc expects gfan to return type ideal. Probably this will change in near future.
    3061         The return type will then be of type LIST_CMD
    3062         Assume gfan has finished, thus we have enumerated all the cones
    3063         Create an array of size of #cones. Let each entry in the array contain a pointer to the respective
    3064         Groebner Basis and merge this somehow with LIST_CMD
    3065         => Count the cones!
    3066         */
    3067         //rChangeCurrRing(rootRing);
    3068         //res=gcAct->gcBasis;
    3069         //res=gcRoot->gcBasis; 
    3070 //      return res;
    30712818        return lResList;
    3072         //return GBlist;
    3073 }
    3074 /*
    3075 Since gfan.cc is #included from extra.cc there must not be a int main(){} here
    3076 */
    3077 #endif
     2819}
     2820
     2821#endif
Note: See TracChangeset for help on using the changeset viewer.