Changeset 3b45cf in git
- Timestamp:
- May 5, 2009, 5:15:02 PM (14 years ago)
- Branches:
- (u'spielwiese', '0d6b7fcd9813a1ca1ed4220cfa2b104b97a0a003')
- Children:
- 8ca3c766834164238a40705d1b906a5f5e003d3f
- Parents:
- def568e9169e67d58842356d8b9718f9a255b706
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
kernel/gfan.cc
rdef568 r3b45cf 2 2 Compute the Groebner fan of an ideal 3 3 $Author: monerjan $ 4 $Date: 2009-05-0 4 15:14:38$5 $Header: /exports/cvsroot-2/cvsroot/kernel/gfan.cc,v 1.4 4 2009-05-04 15:14:38monerjan Exp $6 $Id: gfan.cc,v 1.4 4 2009-05-04 15:14:38monerjan Exp $4 $Date: 2009-05-05 15:15:02 $ 5 $Header: /exports/cvsroot-2/cvsroot/kernel/gfan.cc,v 1.45 2009-05-05 15:15:02 monerjan Exp $ 6 $Id: gfan.cc,v 1.45 2009-05-05 15:15:02 monerjan Exp $ 7 7 */ 8 8 … … 75 75 76 76 77 public: 77 public: 78 //bool isFlippable; //flippable facet? Want to have cone->isflippable.facet[i] 79 bool isIncoming; //Is the facet incoming or outgoing? 80 facet *next; //Pointer to next facet 81 78 82 /** The default constructor. Do I need a constructor of type facet(intvec)? */ 79 83 facet() … … 125 129 } 126 130 127 bool isFlippable; //flippable facet? Want to have cone->isflippable.facet[i] 128 bool isIncoming; //Is the facet incoming or outgoing? 129 facet *next; //Pointer to next facet 131 /*bool isFlippable(intvec &load) 132 { 133 bool res=TRUE; 134 int jj; 135 for (int jj = 0; jj<load.length(); jj++) 136 { 137 intvec *ivCanonical = new intvec(load.length()); 138 (*ivCanonical)[jj]=1; 139 if (ivMult(&load,ivCanonical)<0) 140 { 141 res=FALSE; 142 break; 143 } 144 } 145 return res; 146 147 /*while (dotProduct(load,ivCanonical)>=0) 148 { 149 if (jj!=this->numVars) 150 { 151 intvec *ivCanonical = new intvec(this->numVars); 152 (*ivCanonical)[jj]=1; 153 res=TRUE; 154 jj += 1; 155 } 156 } 157 if (jj==this->numVars) 158 { 159 delete ivCanonical; 160 return FALSE; 161 } 162 else 163 { 164 delete ivCanonical; 165 return TRUE; 166 }*/ 167 }//bool isFlippable(facet &f) 168 130 169 131 170 friend class gcone; //Bad style … … 403 442 }//for (int jj = 1; jj <ddcols; jj++) 404 443 /*Quick'n'dirty hack for flippability*/ 405 bool isFlippable ;444 bool isFlippable=FALSE; 406 445 //NOTE BUG HERE 407 446 /* The flippability check is wrong! … … 409 448 REWRITE THIS 410 449 */ 411 for (int jj = 0; jj<this->numVars; jj++)450 /*for (int jj = 0; jj<this->numVars; jj++) 412 451 { 413 452 intvec *ivCanonical = new intvec(this->numVars); … … 423 462 delete ivCanonical; 424 463 }//for (int jj = 0; jj<this->numVars; jj++) 464 */ 465 for (int jj = 0; jj<load->length(); jj++) 466 { 467 intvec *ivCanonical = new intvec(load->length()); 468 (*ivCanonical)[jj]=1; 469 if (dotProduct(load,ivCanonical)<0) 470 { 471 isFlippable=TRUE; 472 break; //URGHS 473 } 474 } 425 475 if (isFlippable==FALSE) 426 476 { … … 1034 1084 * If this is the case, then our facet is indeed a search facet and TRUE is retuned. 1035 1085 */ 1036 bool isSearchFacet(gcone & tmpcone, facet &testfacet)1086 bool isSearchFacet(gcone &gcTmp, facet &testfacet) 1037 1087 { 1038 1088 ring actRing=currRing; 1039 facet *facetPtr=(facet*) tmpcone.facetPtr;1089 facet *facetPtr=(facet*)gcTmp.facetPtr; 1040 1090 facet *fMin=new facet(*facetPtr); //Pointer to the "minimal" facet 1041 1091 //facet *fMin = new facet(tmpcone.facetPtr); 1042 1092 //fMin=tmpcone.facetPtr; //Initialise to first facet of tmpcone 1043 facet *fAct; 1044 //fMin = testfacet; 1045 //fMin->next=testfacet->next; 1093 facet *fAct; //Ptr to alpha_i 1094 facet *fCmp; //Ptr to alpha_j 1046 1095 fAct = fMin; 1047 1048 cout << endl<< fMin->next << endl; 1049 cout << tmpcone.facetPtr->next << endl; 1050 tmpcone.facetPtr->printNormal(); 1051 fAct=tmpcone.facetPtr->next; 1052 fAct->printNormal(); 1053 /*fMin->printNormal(); 1054 fMin->next->printNormal();*/ 1096 fCmp = fMin->next; 1097 1098 //cout << endl<< fMin->next << endl; 1099 //cout << gcTmp.facetPtr->next << endl; 1100 //gcTmp.facetPtr->printNormal(); 1101 //fAct=gcTmp.facetPtr->next; 1102 //fAct->printNormal(); 1055 1103 1056 1104 rChangeCurrRing(this->rootRing); //because we compare the monomials in the rootring 1057 poly p= NULL;1058 poly q= NULL;1059 intvec * p_weight= new intvec(this->numVars);1060 intvec * q_weight= new intvec(this->numVars);1105 poly p=pInit(); 1106 poly q=pInit(); 1107 intvec *alpha_i = new intvec(this->numVars); 1108 intvec *alpha_j = new intvec(this->numVars); 1061 1109 intvec *sigma = new intvec(this->numVars); 1062 sigma= tmpcone.getIntPoint();1063 cout << "pling" << endl;1110 sigma=gcTmp.getIntPoint(); 1111 1064 1112 int *u=(int *)omAlloc((this->numVars+1)*sizeof(int)); 1065 1113 int *v=(int *)omAlloc((this->numVars+1)*sizeof(int)); 1114 u[0]=0; v[0]=0; 1066 1115 int weight1,weight2; 1067 while(fAct->next !=NULL) //ersetzen durch fAct1116 while(fAct->next->next!=NULL) //NOTE this is ugly. Can it be done without fCmp? 1068 1117 { 1069 1118 /* Get alpha_i and alpha_{i+1} */ 1070 p_weight=fMin->getFacetNormal(); 1071 q_weight=fMin->next->getFacetNormal(); 1072 p_weight->show(); 1073 q_weight->show(); 1119 alpha_i=fAct->getFacetNormal(); 1120 alpha_j=fCmp->getFacetNormal(); 1121 #ifdef gfan_DEBUG 1122 alpha_i->show(); 1123 alpha_j->show(); 1124 #endif 1074 1125 /*Compute the dot product of sigma and alpha_{i,j}*/ 1075 weight1=dotProduct(sigma,p_weight); 1076 weight2=dotProduct(sigma,q_weight); 1126 weight1=dotProduct(sigma,alpha_i); 1127 weight2=dotProduct(sigma,alpha_j); 1128 #ifdef gfan_DEBUG 1077 1129 cout << "weight1=" << weight1 << " " << "weight2=" << weight2 << endl; 1130 #endif 1078 1131 /*Adjust alpha_i and alpha_i+1 accordingly*/ 1079 1132 for(int ii=1;ii<=this->numVars;ii++) 1080 { 1081 /*p_weight[ii]=weight1*(*p_weight)[ii]; 1082 q_weight[ii]=weight2*(*q_weight)[ii];*/ 1083 u[ii]=weight1*(*p_weight)[ii]; 1084 v[ii]=weight2*(*q_weight)[ii]; 1085 } 1086 cout << "PLONK" << endl; 1133 { 1134 u[ii]=weight1*(*alpha_i)[ii-1]; 1135 v[ii]=weight2*(*alpha_j)[ii-1]; 1136 } 1087 1137 1088 1138 /*Now p_weight and q_weight need to be compared as exponent vectors*/ 1089 pSetExpV(p,u); 1090 pSetExpV(q,v); 1091 cout << "AARGH" << endl; 1139 pSetCoeff0(p,nInit(1)); 1140 pSetCoeff0(q,nInit(1)); 1141 pSetExpV(p,u); 1142 pSetm(p); 1143 pSetExpV(q,v); 1144 pSetm(q); 1145 #ifdef gfan_DEBUG 1146 pWrite(p);pWrite(q); 1147 #endif 1092 1148 /*We want to check whether x^p < x^q 1093 1149 => want to check for return value 1 */ 1094 1150 if (pLmCmp(p,q)==1) //i.e. x^q is smaller 1095 1151 { 1096 fMin=f Min->next;1152 fMin=fCmp; 1097 1153 fAct=fMin; 1098 1154 } 1099 1155 else 1100 1156 { 1101 fAct=fAct->next; 1157 //fAct=fAct->next; 1158 if(fCmp->next!=NULL) 1159 { 1160 fCmp=fCmp->next; 1161 } 1162 else 1163 { 1164 fAct=fAct->next; 1165 } 1102 1166 } 1103 fAct=fAct->next;1104 }//while( tmpcone.facetPtr->next!=NULL)1167 //fAct=fAct->next; 1168 }//while(fAct.facetPtr->next!=NULL) 1105 1169 1106 1170 /*If testfacet was minimal then fMin should still point there */ 1107 //if (isParallel(fMin->getFacetNormal,testfacet.getFacetNormal)) 1108 if (fMin==tmpcone.facetPtr) 1171 intvec *alpha_min = new intvec(this->numVars); 1172 alpha_min=fMin->getFacetNormal(); 1173 intvec *test = new intvec(this->numVars); 1174 test=testfacet.getFacetNormal(); 1175 if (isParallel(alpha_min,test)) 1176 //if (fMin==gcTmp.facetPtr) 1109 1177 { 1110 1178 rChangeCurrRing(actRing); … … 1156 1224 }//while(fAct->next!=NULL) 1157 1225 }//reverseSearch 1226 friend class facet; 1158 1227 };//class gcone 1159 1228
Note: See TracChangeset
for help on using the changeset viewer.