Changeset 0d278cd in git


Ignore:
Timestamp:
Nov 17, 2009, 3:53:24 PM (14 years ago)
Author:
Martin Monerjan
Branches:
(u'spielwiese', 'fe61d9c35bf7c61f2b6cbf1b56e25e2f08d536cc')
Children:
87353a6b9eb171d12e73ede9e24e8d8e68eaf111
Parents:
26300d06f8088e506f0c3c44b085718c58333635
Message:
gcone::getCounter
gcone::getBaseRing
Changed coputation of diffs of expvects in getConeNormals
Cleanup of no longer used stuff
lprepareResult


git-svn-id: file:///usr/local/Singular/svn/trunk@12286 2c84dea3-7e68-4137-9b89-c4e89433aadc
Location:
kernel
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • kernel/gfan.cc

    r26300d r0d278cd  
    386386        //dd_FreeMatrix(this->ddFacets);
    387387}                       
    388                
     388       
     389int gcone::getCounter()
     390{
     391        return this->counter;
     392}
     393       
    389394/** \brief Set the interior point of a cone */
    390395void gcone::setIntPoint(intvec *iv)
    … …  
    549554        int pCompCount;                 // # of terms in a poly
    550555        poly aktpoly;
    551         int numvar = pVariables;        // # of variables in a polynomial (or ring?)
    552         int leadexp[numvar];            // dirty hack of exp.vects
    553         int aktexp[numvar];
     556        int numvar = pVariables;        // # of variables in currRing
     557//      int leadexp[numvar];            // dirty hack of exp.vects
     558//      int aktexp[numvar];
    554559        int cols,rows;                  // will contain the dimensions of the ineq matrix - deprecated by
    555560        dd_rowrange ddrows;
    … …  
    562567        // End of var declaration
    563568#ifdef gfan_DEBUG
    564 //                      cout << "The Groebner basis has " << lengthGB << " elements" << endl;
    565 //                      cout << "The current ring has " << numvar << " variables" << endl;
     569//      cout << "The Groebner basis has " << lengthGB << " elements" << endl;
     570//      cout << "The current ring has " << numvar << " variables" << endl;
    566571#endif
    567572        cols = numvar;
    568573               
    569                         //Compute the # inequalities i.e. rows of the matrix
     574        //Compute the # inequalities i.e. rows of the matrix
    570575        rows=0; //Initialization
    571576        for (int ii=0;ii<IDELEMS(I);ii++)
    … …  
    575580        }
    576581#ifdef gfan_DEBUG
    577 //                      cout << "rows=" << rows << endl;
    578 //                      cout << "Will create a " << rows << " x " << cols << " matrix to store inequalities" << endl;
     582//      cout << "rows=" << rows << endl;
     583//      cout << "Will create a " << rows << " x " << cols << " matrix to store inequalities" << endl;
    579584#endif
    580585        dd_rowrange aktmatrixrow=0;     // needed to store the diffs of the expvects in the rows of ddineq
    … …  
    594599#endif
    595600               
    596                 int *v=(int *)omAlloc((numvar+1)*sizeof(int));
    597                 pGetExpV(aktpoly,v);    //find the exp.vect in v[1],...,v[n], use pNext(p)
    598                                
    599                 //Store leadexp for aktpoly
    600                 for (int kk=0;kk<numvar;kk++)
    601                 {
    602                         leadexp[kk]=v[kk+1];
    603                         //Since we need to know the difference of leadexp with the other expvects we do nothing here
    604                         //but compute the diff below
    605                 }
    606                
    607                                
    608                 while (pNext(aktpoly)!=NULL) //move to next term until NULL
     601//              int *v=(int *)omAlloc((numvar+1)*sizeof(int));
     602//              pGetExpV(aktpoly,v);    //find the exp.vect in v[1],...,v[n], use pNext(p)
     603//                             
     604//              //Store leadexp for aktpoly
     605//              for (int kk=0;kk<numvar;kk++)
     606//              {
     607//                      leadexp[kk]=v[kk+1];
     608//                      //Since we need to know the difference of leadexp with the other expvects we do nothing here
     609//                      //but compute the diff below
     610//              }
     611//             
     612//                             
     613//              while (pNext(aktpoly)!=NULL) //move to next term until NULL
     614//              {
     615//                      aktpoly=pNext(aktpoly);
     616//                      pSetm(aktpoly);         //doesn't seem to help anything
     617//                      pGetExpV(aktpoly,v);
     618//                     
     619//                      for (int kk=0;kk<numvar;kk++)
     620//                      {
     621//                              aktexp[kk]=v[kk+1];
     622//                                              //ineq[aktmatrixrow][kk]=leadexp[kk]-aktexp[kk];        //dito
     623//                              dd_set_si(ddineq->matrix[(dd_rowrange)aktmatrixrow][kk+1],leadexp[kk]-aktexp[kk]); //because of the 1st col being const 0
     624//                      }
     625//                      aktmatrixrow=aktmatrixrow+1;
     626//              }//while
     627//              omFree(v);
     628               
     629                //simpler version of the above
     630                int *leadexpv=(int*)omAlloc((numvar+1)*sizeof(int));
     631                int *tailexpv=(int*)omAlloc((numvar+1)*sizeof(int));
     632                pGetExpV(aktpoly,leadexpv);
     633                while(pNext(aktpoly)!=NULL)
    609634                {
    610635                        aktpoly=pNext(aktpoly);
    611                         pSetm(aktpoly);         //doesn't seem to help anything
    612                         pGetExpV(aktpoly,v);
    613                        
    614                         for (int kk=0;kk<numvar;kk++)
     636                        pGetExpV(aktpoly,tailexpv);
     637                        for(int kk=1;kk<=numvar;kk++)
    615638                        {
    616                                 aktexp[kk]=v[kk+1];
    617                                                 //ineq[aktmatrixrow][kk]=leadexp[kk]-aktexp[kk];        //dito
    618                                 dd_set_si(ddineq->matrix[(dd_rowrange)aktmatrixrow][kk+1],leadexp[kk]-aktexp[kk]); //because of the 1st col being const 0
     639                                dd_set_si(ddineq->matrix[(dd_rowrange)aktmatrixrow][kk],leadexpv[kk]-tailexpv[kk]);
    619640                        }
    620                         aktmatrixrow=aktmatrixrow+1;
    621                 }//while
    622                 omFree(v);
     641                        aktmatrixrow += 1;
     642                }
     643                omFree(tailexpv);
     644                omFree(leadexpv);
     645               
    623646        } //for
    624647               
    … …  
    10571080        computeInv(gb,initialForm,*fNormal);
    10581081        //NOTE The code below went into gcone::computeInv
    1059 //      for (int ii=0;ii<IDELEMS(gb);ii++)
    1060 //      {
    1061 //              aktpoly = (poly)gb->m[ii];                                                             
    1062 //              int *v=(int *)omAlloc((this->numVars+1)*sizeof(int));
    1063 //              int *leadExpV=(int *)omAlloc((this->numVars+1)*sizeof(int));
    1064 //              pGetExpV(aktpoly,leadExpV);     //find the leading exponent in leadExpV[1],...,leadExpV[n], use pNext(p)
    1065 //              initialFormElement[ii]=pHead(aktpoly);
    1066 //                             
    1067 //              while(pNext(aktpoly)!=NULL)     //*loop trough terms and check for parallelity*/
    1068 //              {
    1069 //                      aktpoly=pNext(aktpoly); //next term
    1070 //                      pSetm(aktpoly);
    1071 //                      pGetExpV(aktpoly,v);           
    1072 //                      /* Convert (int)v into (intvec)check */                 
    1073 //                      for (int jj=0;jj<this->numVars;jj++)
    1074 //                      {
    1075 //                                              //cout << "v["<<jj+1<<"]="<<v[jj+1]<<endl;
    1076 //                                              //cout << "leadExpV["<<jj+1<<"]="<<leadExpV[jj+1]<<endl;
    1077 //                              (*check)[jj]=v[jj+1]-leadExpV[jj+1];
    1078 //                      }
    1079 // #ifdef gfan_DEBUG
    1080 // //                                   cout << "check=";                       
    1081 // //                                   check->show();
    1082 // //                                   cout << endl;
    1083 // #endif
    1084 //                      if (isParallel(*check,*fNormal)) //pass *check when
    1085 //                      {
    1086 // //                           cout << "Parallel vector found, adding to initialFormElement" << endl;                 
    1087 //                              initialFormElement[ii] = pAdd(pCopy(initialFormElement[ii]),(poly)pHead(aktpoly));
    1088 //                      }                                               
    1089 //              }//while
    1090 // #ifdef gfan_DEBUG
    1091 // //                           cout << "Initial Form=";                               
    1092 // //                           pWrite(initialFormElement[ii]);
    1093 // //                           cout << "---" << endl;
    1094 // #endif
    1095 //              //*Now initialFormElement must be added to (ideal)initialForm */
    1096 //              initialForm->m[ii]=initialFormElement[ii];
    1097 //              omFree(leadExpV);
    1098 //              omFree(v);
    1099 //      }//for                 
    11001082#ifdef gfan_DEBUG
    11011083/*      cout << "Initial ideal is: " << endl;
    … …  
    11421124#ifndef NDEBUG                 
    11431125#ifdef gfan_DEBUG
    1144                         cout << "ina=";
    1145                         idShow(ina); cout << endl;
     1126        cout << "ina=";
     1127        idShow(ina); cout << endl;
    11461128#endif
    11471129#endif
    11481130        ideal H;
    1149                         //H=kStd(ina,NULL,isHomog,NULL);        //we know it is homogeneous
     1131        //H=kStd(ina,NULL,isHomog,NULL);        //we know it is homogeneous
    11501132        H=kStd(ina,NULL,testHomog,NULL);        //This is \mathcal(G)_{>_-\alpha}(in_v(I))
    11511133        idSkipZeroes(H);
    … …  
    11641146#ifndef NDEBUG
    11651147#ifdef gfan_DEBUG
    1166                         cout << "srcRing_H = ";
    1167                         idShow(srcRing_H); cout << endl;
     1148        cout << "srcRing_H = ";
     1149        idShow(srcRing_H); cout << endl;
    11681150#endif
    11691151#endif
    … …  
    11711153#ifndef NDEBUG
    11721154#ifdef gfan_DEBUG
    1173                         cout << "srcRing_HH = ";
    1174                         idShow(srcRing_HH); cout << endl;
     1155        cout << "srcRing_HH = ";
     1156        idShow(srcRing_HH); cout << endl;
    11751157#endif
    11761158#endif
    … …  
    12281210                                markingsAreCorrect=TRUE; //everything is fine
    12291211#ifdef gfan_DEBUG
    1230 //                                              cout << "correct markings" << endl;
     1212//                              cout << "correct markings" << endl;
    12311213#endif
    12321214                        }//if (pHead(aktpoly)==pHead(H->m[jj])
    1233 //                      delete src_ExpV;
    1234 //                      delete dst_ExpV;
    12351215                        omFree(src_ExpV);
    12361216                        omFree(dst_ExpV);
    … …  
    12721252                        aktrow +=1;
    12731253                }
    1274 //              delete v;
    1275 //              delete leadExpV;
    12761254                omFree(v);
    12771255                omFree(leadExpV);
    12781256        }//for (int ii=0;ii<IDELEMS(srcRing_HH);ii++)
    12791257        /*Now we add the constraint for the standard simplex*/
    1280         /*NOTE:Might actually work without the standard simplex*/
     1258        /*NOTE:Might actually work without the standard simplex
     1259        * No, it won't. MM
     1260        */
    12811261        dd_set_si(intPointMatrix->matrix[aktrow][0],-1);
    12821262        for (int jj=1;jj<=this->numVars;jj++)
    … …  
    13261306        rComplete(dstRing);                                     
    13271307        rChangeCurrRing(dstRing);
    1328                         //rDelete(tmpRing);
     1308        //rDelete(tmpRing);
    13291309        delete iv_weight;
    13301310//#ifdef gfan_DEBUG
    … …  
    13451325#endif
    13461326        ideal tmpI;
    1347                         //NOTE Any of the two variants of tmpI={idrCopy(),dstRing_I} does the trick
    1348                         //tmpI = idrCopyR(this->inputIdeal,this->baseRing);
     1327        //NOTE Any of the two variants of tmpI={idrCopy(),dstRing_I} does the trick
     1328        //tmpI = idrCopyR(this->inputIdeal,this->baseRing);
    13491329        tmpI = dstRing_I;
    13501330        dstRing_I=kStd(tmpI,NULL,testHomog,NULL);
    … …  
    13991379                        }
    14001380#ifdef gfan_DEBUG
    1401 //                                      cout << "check=";                       
    1402 //                                      check->show();
    1403 //                                      cout << endl;
     1381//                      cout << "check=";                       
     1382//                      check->show();
     1383//                      cout << endl;
    14041384#endif
    14051385                        if (isParallel(*check,fNormal)) //pass *check when
    … …  
    14101390                }//while
    14111391#ifdef gfan_DEBUG
    1412 //                              cout << "Initial Form=";                               
    1413 //                              pWrite(initialFormElement[ii]);
    1414 //                              cout << "---" << endl;
     1392//              cout << "Initial Form=";                               
     1393//              pWrite(initialFormElement[ii]);
     1394//              cout << "---" << endl;
    14151395#endif
    14161396                /*Now initialFormElement must be added to (ideal)initialForm */
    … …  
    14301410 */
    14311411//NOTE: Should be replaced by kNF or kNF2
     1412//NOTE: Done
    14321413poly gcone::restOfDiv(poly const &f, ideal const &I)
    14331414{
    14341415//                      cout << "Entering restOfDiv" << endl;
    1435         poly p=f;
    1436                         //pWrite(p);                   
    1437         poly r=NULL;    //The zero polynomial
    1438         int ii;
    1439         bool divOccured;
    1440                        
    1441         while (p!=NULL)
    1442         {
    1443                 ii=1;
    1444                 divOccured=FALSE;
    1445                                
    1446                 while( (ii<=IDELEMS(I) && (divOccured==FALSE) ))
    1447                 {                                       
    1448                         if (pDivisibleBy(I->m[ii-1],p)) //does LM(I->m[ii]) divide LM(p) ?
    1449                         {                                               
    1450                                 poly step1,step2,step3;
    1451                                                 //cout << "dividing "; pWrite(pHead(p));cout << "by ";pWrite(pHead(I->m[ii-1])); cout << endl;
    1452                                 step1 = pDivideM(pHead(p),pHead(I->m[ii-1]));
    1453                                                 //cout << "LT(p)/LT(f_i)="; pWrite(step1); cout << endl;                               
    1454                                 step2 = ppMult_qq(step1, I->m[ii-1]);                                           
    1455                                 step3 = pSub(pCopy(p), step2);
    1456                                                 //p=pSub(p,pMult( pDivide(pHead(p),pHead(I->m[ii])), I->m[ii]));                       
    1457                                                 //pSetm(p);
    1458                                 pSort(step3); //must be here, otherwise strange behaviour with many +o+o+o+o+ terms
    1459                                 p=step3;
    1460                                                 //pWrite(p);                                           
    1461                                 divOccured=TRUE;
    1462                         }
    1463                         else
    1464                         {
    1465                                 ii += 1;
    1466                         }//if (pLmDivisibleBy(I->m[ii],p,currRing))
    1467                 }//while( (ii<IDELEMS(I) && (divOccured==FALSE) ))
    1468                 if (divOccured==FALSE)
    1469                 {
    1470                                         //cout << "TICK 5" << endl;
    1471                         r=pAdd(pCopy(r),pHead(p));
    1472                         pSetm(r);
    1473                         pSort(r);
    1474                                         //cout << "r="; pWrite(r); cout << endl;
    1475                                        
    1476                         if (pLength(p)!=1)
    1477                         {
    1478                                 p=pSub(pCopy(p),pHead(p));      //Here it may occur that p=0 instead of p=NULL
    1479                         }
    1480                         else
    1481                         {
    1482                                 p=NULL; //Hack to correct this situation                                               
    1483                         }                                       
    1484                                         //cout << "p="; pWrite(p);
    1485                 }//if (divOccured==FALSE)
    1486         }//while (p!=0)
    1487         return r;
     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;
    14881469}//poly restOfDiv(poly const &f, ideal const &I)
    14891470               
    14901471/** \brief Compute \f$ f-f^{\mathcal{G}} \f$
    1491         */
     1472*/
    14921473//NOTE: use kNF or kNF2 instead of restOfDivision
    14931474ideal gcone::ffG(ideal const &H, ideal const &G)
    … …  
    15051486                temp3=pSub(temp1, temp2);
    15061487                res->m[ii]=temp3;
    1507                                 //res->m[ii]=pSub(temp1,temp2); //buggy
    1508                                 //pSort(res->m[ii]);
    1509                                 //pSetm(res->m[ii]);
    1510                                 //cout << "res->m["<<ii<<"]=";pWrite(res->m[ii]);                                               
     1488                //res->m[ii]=pSub(temp1,temp2); //buggy
     1489                //pSort(res->m[ii]);
     1490                //pSetm(res->m[ii]);
     1491                //cout << "res->m["<<ii<<"]=";pWrite(res->m[ii]);                                               
    15111492        }                       
    15121493        return res;
    … …  
    18641845void gcone::reverseSearch(gcone *gcAct) //no const possible here since we call gcAct->flip
    18651846{
    1866         facet *fAct=new facet();
    1867         fAct = gcAct->facetPtr;                 
    1868                        
    1869         while(fAct!=NULL)
    1870         {
    1871                 cout << "======================"<< endl;
    1872                 gcAct->flip(gcAct->gcBasis,gcAct->facetPtr);
    1873                 gcone *gcTmp = new gcone(*gcAct);
    1874                                 //idShow(gcTmp->gcBasis);
    1875                 gcTmp->getConeNormals(gcTmp->gcBasis, TRUE);
    1876 #ifdef gfan_DEBUG
    1877                 facet *f = new facet();
    1878                 f=gcTmp->facetPtr;
    1879                 while(f!=NULL)
    1880                 {
    1881                         f->printNormal();
    1882                         f=f->next;                                     
    1883                 }
    1884 #endif
    1885                 gcTmp->showIntPoint();
    1886                 /*recursive part goes gere*/
    1887                 if (isSearchFacet(*gcTmp,(facet*)gcAct->facetPtr))
    1888                 {
    1889                         gcAct->next=gcTmp;
    1890                         cout << "PING"<< endl;
    1891                         reverseSearch(gcTmp);
    1892                 }
    1893                 else
    1894                 {
    1895                         delete gcTmp;
    1896                         /*NOTE remove fAct from linked list. It's no longer needed*/
    1897                 }
    1898                 /*recursion ends*/
    1899                 fAct = fAct->next;             
    1900         }//while(fAct->next!=NULL)
     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)
    19011882}//reverseSearch
    19021883               
    … …  
    28962877//      }
    28972878// }
    2898 
    2899 void gcone::lprepareResult(lists l, gcone &gc)
     2879ring gcone::getBaseRing()
     2880{
     2881        return this->baseRing;
     2882}
     2883/** \brief Gather the output
     2884* List of lists
     2885*\param n the number of cones
     2886*/
     2887lists lprepareResult(gcone *gc, int n)
    29002888{
    29012889        gcone *gcAct;
    2902         gcAct = &gc;
    2903         matrix mFacetNormal=mpNew(counter,gc.numVars);
     2890        gcAct = gc;     
     2891        facet *fAct;
     2892        fAct = gc->facetPtr;
    29042893       
    2905         facet *fAct;
    2906         fAct = gc.facetPtr;
    2907         while( gcAct!=NULL )
    2908         {
    2909                 l->m[0].data=(int*)gc.getUCN();
    2910                 l->m[1].data=(intvec*)(&gcAct->f2M(gcAct,gcAct->facetPtr));
     2894        lists res=(lists)omAllocBin(slists_bin);
     2895        res->Init(n);
     2896        for(int ii=0;ii<n;ii++)
     2897        {
     2898                res->m[ii].rtyp=LIST_CMD;
     2899                lists l=(lists)omAllocBin(slists_bin);
     2900                l->Init(4);
     2901                l->m[0].rtyp=INT_CMD;
     2902                l->m[1].rtyp=INTVEC_CMD;
     2903                l->m[2].rtyp=IDEAL_CMD;
     2904                l->m[3].rtyp=RING_CMD;
     2905//              l->m[4].rtyp=LIST_CMD;
     2906                l->m[0].data=(void*)gcAct->getUCN();
     2907                intvec iv=(gcAct->f2M(gcAct,gcAct->facetPtr));
     2908                l->m[1].data=(void*)ivCopy(&iv);
     2909                l->m[2].data=(void*)idCopy(gcAct->gcBasis);
     2910                l->m[3].data=(void*)(gcAct->getBaseRing());
     2911//              l->m[4].data=(void*)0;
     2912                res->m[ii].data=(void*)l;
    29112913                gcAct = gcAct->next;
    2912         }
     2914        }               
     2915       
     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;
     2925        return res;
    29132926}
    29142927/** \brief Write facets of a cone into a matrix
    … …  
    29272940                codim=2;
    29282941        }
    2929         intvec mRes=new intvec(counter,gc->numVars,0);//nrows==counter, ncols==numVars
     2942        /** mrRes is a matrix containing the facet normals AS ROWS*/
     2943        intvec *mRes=new intvec(this->numFacets,gc->numVars,0);//nrows==numFacets, ncols==numVars
    29302944        intvec *fNormal = new intvec(this->numVars);
    2931         for(int ii=0;ii<gc->numFacets*(this->numVars);ii++)
     2945        int ii=0;
     2946//      for(int ii=0;ii<gc->numFacets*(this->numVars);ii++)
     2947        while(fAct!=NULL && ii<gc->numFacets*(this->numVars))
    29322948        {
    29332949                fNormal=fAct->getFacetNormal();
    29342950                for(int jj=0;jj<this->numVars ;jj++ )
    29352951                {                       
    2936                         mRes[ii]=(*fNormal)[jj];
     2952                        (*mRes)[ii]=(*fNormal)[jj];
    29372953                        ii++;
    29382954                }
    29392955                fAct = fAct->next;
    29402956        }
    2941         return mRes;
     2957        return *mRes;
    29422958}
    29432959
    … …  
    29472963lists gfan(ideal inputIdeal, int h)
    29482964{
    2949         lists lResList=(lists)omAllocBin(slists_bin);
    2950         lResList->Init(5);
    2951         lResList->m[0].rtyp=INT_CMD;
    2952         lResList->m[0].data=(int*)255;
    2953         lResList->m[1].rtyp=INTVEC_CMD;
    2954         lResList->m[1].data=(intvec*)NULL;
    2955         lResList->m[2].rtyp=IDEAL_CMD;
    2956         lResList->m[2].data=(ideal*)NULL;
    2957         lResList->m[3].rtyp=RING_CMD;
    2958         lResList->m[3].data=(ring*)NULL;
    2959         lResList->m[4].rtyp=LIST_CMD;
    2960         lResList->m[4].data=(lists)NULL;
     2965        lists lResList;
    29612966        int numvar = pVariables;
    29622967        gfanHeuristic = h;
    … …  
    30473052                //Below is a workaround, since gcAct->gcBasis gets deleted in noRevS
    30483053                res = inputIdeal;
    3049 //              intvec mRes=gcRoot->f2M(gcRoot,gcRoot->facetPtr);
    3050                 gcRoot->lprepareResult(lResList,*gcRoot);
     3054//              intvec mRes=gcRoot->f2M(gcRoot,gcRoot->facetPtr);               
     3055                lResList=lprepareResult(gcRoot,gcRoot->getCounter());
     3056//              res=lResList;
    30513057        }
    30523058        dd_free_global_constants();

  • kernel/gfan.h

    r26300d r0d278cd  
    2424// ideal gfan(ideal inputIdeal, int heuristic);
    2525lists gfan(ideal inputIdeal, int heuristic);
     26
    2627//int dotProduct(intvec a, intvec b);
    2728//bool isParallel(intvec a, intvec b);
    … …  
    123124                intvec *ivIntPt;        //an interior point of the cone
    124125                int UCN;                //unique number of the cone
    125                 static int counter;
     126                static int counter;
    126127               
    127128        public:
    … …  
    153154                gcone(const gcone& gc, const facet &f);
    154155                ~gcone();
     156                int getCounter();
     157                ring getBaseRing();
    155158                void setIntPoint(intvec *iv);
    156159                intvec *getIntPoint();
    … …  
    188191                void writeConeToFile(gcone const &gc, bool usingIntPoints=FALSE);
    189192                void readConeFromFile(int gcNum, gcone *gc);
    190                 void lprepareResult(lists l, gcone &gc);
     193               
    191194//              static void gcone::idPrint(ideal &I);
    192195                intvec f2M(gcone *gc, facet *f);
    193196                friend class facet;     
    194197};
    195 
     198lists lprepareResult(gcone *gc, int n);
    196199#endif
    197200#endif
Note: See TracChangeset for help on using the changeset viewer.