1 // Copyright (c) 2014-2017 John Abbott and Anna M. Bigatti 2 // Authors: 2014-2017 John Abbott, Anna M. Bigatti 3 4 // This file is part of the source of CoCoALib, the CoCoA Library. 5 6 // CoCoALib is free software: you can redistribute it and/or modify 7 // it under the terms of the GNU General Public License as published by 8 // the Free Software Foundation, either version 3 of the License, or 9 // (at your option) any later version. 10 11 // CoCoALib is distributed in the hope that it will be useful, 12 // but WITHOUT ANY WARRANTY; without even the implied warranty of 13 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 14 // GNU General Public License for more details. 15 16 // You should have received a copy of the GNU General Public License 17 // along with CoCoALib. If not, see <http://www.gnu.org/licenses/>. 18 19 20 #include "CoCoA/SparsePolyOps-implicit.H" 21 22 #include "CoCoA/DenseMatrix.H" // for elim (weights) and ImplicitDirectOrd2 23 #include "CoCoA/MatrixForOrdering.H" // for elim (weights) and ImplicitDirectOrd2 24 #include "CoCoA/MatrixOps.H" // for ImplicitDirectOrd2 25 #include "CoCoA/MatrixView.H" // for ImplicitDirectOrd2 26 #include "CoCoA/NumTheory-prime.H" 27 #include "CoCoA/PPMonoid.H" 28 #include "CoCoA/PPMonoidHom.H" // for PPMonoidHom 29 #include "CoCoA/PPOrdering.H" 30 #include "CoCoA/PPOrdering.H" // for elim (weights) and ImplicitDirectOrd2 31 #include "CoCoA/QBGenerator.H" 32 #include "CoCoA/QuotientRing.H" // for NewZZmod 33 #include "CoCoA/RingDistrMPolyInlFpPP.H" 34 #include "CoCoA/RingDistrMPolyInlPP.H" 35 #include "CoCoA/RingHom.H" 36 #include "CoCoA/RingQQ.H" // for elim (weights) 37 #include "CoCoA/RingZZ.H" // for ord matrix (ImplicitDirectOrd2) 38 #include "CoCoA/SmallFpImpl.H" 39 #include "CoCoA/SparsePolyIter.H" 40 #include "CoCoA/SparsePolyOps-RingElem.H" 41 //#include "CoCoA/SparsePolyRing.H" 42 #include "CoCoA/TmpGOperations.H" // for ComputeElimFirst 43 #include "CoCoA/VectorOps.H" // (tmp) for printing 44 #include "CoCoA/apply.H" 45 #include "CoCoA/ideal.H" // for NF: checking result 46 #include "CoCoA/interrupt.H" // for CheckForInterrupt 47 #include "CoCoA/random.H" 48 #include "CoCoA/time.H" // (tmp) for timings 49 #include "CoCoA/verbose.H" // for VerboseLog 50 51 #include <map> 52 using std::map; 53 #include <string> 54 using std::string; 55 //#include <vector> 56 using std::vector; 57 58 namespace CoCoA 59 { 60 61 // ImplicitDirect 62 63 namespace // anonymous for file local fns 64 { 65 ComputeImage(const PPMonoidElem & t,const vector<RingElem> & L)66 RingElem ComputeImage(const PPMonoidElem& t, const vector<RingElem>& L) 67 { 68 RingElem ans = one(owner(L[0])); 69 vector<long> exp; exponents(exp, t); 70 for (int i=0; i < len(exp); ++i) 71 ans *= power(L[i],exp[i]); 72 return ans; 73 } 74 reduce(RingElem & NewRed,RingElem & NewPolyx,map<PPMonoidElem,int> & FindReducerIndex,const vector<RingElem> & reducer,const vector<RingElem> & polyx)75 void reduce(RingElem& NewRed, RingElem& NewPolyx, map<PPMonoidElem, int>& FindReducerIndex, const vector<RingElem>& reducer, const vector<RingElem>& polyx) 76 { 77 // RingHom embed1 = CoeffEmbeddingHom(owner(NewRed)); 78 // RingHom embed2 = CoeffEmbeddingHom(owner(NewPolyx)); 79 while (true) 80 { 81 if (IsZero(NewRed)) return; 82 const PPMonoidElem PP = LPP(NewRed); 83 const int i = -1 + FindReducerIndex[PP]; // remove shift of 1 84 if (i == -1) return; 85 const RingElem c = LC(NewRed)/LC(reducer[i]); 86 // NewRed -= embed1(c)*reducer[i]; 87 // NewPolyx -= embed2(c)*polyx[i]; 88 RingElem red = reducer[i]; 89 RingElem pol = polyx[i]; 90 SparsePolyRingPtr(owner(red))->myMulByCoeff(raw(red), raw(c)); 91 SparsePolyRingPtr(owner(pol))->myMulByCoeff(raw(pol), raw(c)); 92 NewRed -= red; 93 NewPolyx -= pol; 94 } 95 } 96 97 // this allows the computation in QQ[...] 98 // in fact it is propably better to compute in Fp[...] and lift NewPolyRingForImplicit(const ring & K,const string & name,long n)99 SparsePolyRing NewPolyRingForImplicit(const ring& K, const string& name, long n) 100 { 101 if (IsQQ(K)) 102 return NewPolyRing_DMPI(K, SymbolRange(name,1,n)); 103 return NewPolyRing_DMPII(K, SymbolRange(name,1,n)); 104 //this should be done properly: check if K is SmallFpImpl 105 } 106 NewPolyRingForImplicit(const ring & K,long n,PPOrdering ord)107 SparsePolyRing NewPolyRingForImplicit(const ring& K, long n, PPOrdering ord) 108 { 109 if (IsQQ(K)) 110 return NewPolyRing_DMPI(K, NewSymbols(n), ord); 111 return NewPolyRing_DMPII(K, NewSymbols(n), ord); 112 //this should be done properly: check if K is SmallFpImpl 113 } 114 NewPolyRingForImplicit(const ring & K,const string & name,PPOrdering ord)115 SparsePolyRing NewPolyRingForImplicit(const ring& K, const string& name, PPOrdering ord) 116 { 117 if (IsQQ(K)) 118 return NewPolyRing_DMPI(K, SymbolRange(name,1,NumIndets(ord)), ord); 119 return NewPolyRing_DMPII(K, SymbolRange(name,1,NumIndets(ord)), ord); 120 //this should be done properly: check if K is SmallFpImpl 121 } 122 123 } // end of anonymous namespace /////////////////////////////////// 124 125 ImplicitDirect(const std::vector<RingElem> & ParamDescrOrig)126 RingElem ImplicitDirect(const std::vector<RingElem>& ParamDescrOrig) 127 { 128 if (ParamDescrOrig.empty() || len(ParamDescrOrig) != 1+NumIndets(owner(ParamDescrOrig[0]))) 129 CoCoA_THROW_ERROR(ERR::BadArg, "ImplicitDirectLPP"); 130 const ring& Porig = owner(ParamDescrOrig[0]); 131 const int n = len(ParamDescrOrig); 132 133 ring Kt = NewPolyRingForImplicit(CoeffRing(Porig), "t", NumIndets(Porig)); 134 RingHom phi = PolyAlgebraHom(Porig, Kt, indets(Kt)); 135 vector<RingElem> ParamDescr = apply(phi, ParamDescrOrig); 136 // Originally the code was like this (ie. compute in the originally given ring) 137 // if (ParamDescr.empty() || len(ParamDescr) != 1+NumIndets(owner(ParamDescr[0]))) 138 // CoCoA_THROW_ERROR(ERR::BadArg, "ImplicitDirect"); 139 // ring P = owner(ParamDescr[0]); 140 // const int n = len(ParamDescr); 141 ring Kx = NewPolyRingForImplicit(CoeffRing(Porig), "x", n); 142 vector<RingElem> reducer; 143 vector<RingElem> polyx; 144 map<PPMonoidElem,int> FindReducerIndex; 145 146 reducer.push_back(one(Kt)); 147 polyx.push_back(one(Kx)); 148 FindReducerIndex[LPP(reducer[0])] = 1; // deliberately add 1 149 QBGenerator QBG(PPM(Kx)); 150 PPMonoidElem PP1 = QBG.myCorners().front(); 151 QBG.myCornerPPIntoQB(PP1); 152 while (true) 153 { 154 const PPMonoidElem PP = QBG.myCorners().front(); 155 QBG.myCornerPPIntoQB(PP); 156 RingElem NewReducer = ComputeImage(PP, ParamDescr); 157 RingElem NewPolyx = monomial(Kx, PP); 158 reduce(NewReducer, NewPolyx, FindReducerIndex, reducer, polyx); 159 if (IsZero(NewReducer)) return NewPolyx; 160 reducer.push_back(NewReducer); 161 polyx.push_back(NewPolyx); 162 FindReducerIndex[LPP(NewReducer)] = len(reducer); 163 } 164 } 165 166 ComputeLPP(const PPMonoidElem & t,const vector<PPMonoidElem> & VecLPP)167 PPMonoidElem ComputeLPP(const PPMonoidElem& t, const vector<PPMonoidElem>& VecLPP) 168 { 169 const PPMonoid& PPM = owner(VecLPP[0]); 170 PPMonoidElem ans = one(PPM); 171 const int n = len(VecLPP); 172 CoCoA_ASSERT(n == NumIndets(owner(t))); 173 vector<long> exp(n); exponents(exp,t); 174 for (int i=0;i<n;++i) 175 ans *= power(VecLPP[i],exp[i]); 176 return ans; 177 } 178 179 ImplicitDirectLPP(const std::vector<RingElem> & ParamDescrOrig)180 RingElem ImplicitDirectLPP(const std::vector<RingElem>& ParamDescrOrig) 181 { 182 if (ParamDescrOrig.empty() || len(ParamDescrOrig) != 1+NumIndets(owner(ParamDescrOrig[0]))) 183 CoCoA_THROW_ERROR(ERR::BadArg, "ImplicitDirectLPP"); 184 const ring& Porig = owner(ParamDescrOrig[0]); 185 const int n = len(ParamDescrOrig); 186 187 ring Kt = NewPolyRingForImplicit(CoeffRing(Porig), "t", NumIndets(Porig)); 188 RingHom phi = PolyAlgebraHom(Porig, Kt, indets(Kt)); 189 vector<RingElem> ParamDescr = apply(phi, ParamDescrOrig); 190 ring Kx = NewPolyRingForImplicit(CoeffRing(Porig), "x",n); 191 vector<PPMonoidElem> VecLPP; VecLPP.reserve(len(ParamDescr)); for(int i=0;i<len(ParamDescr);++i)VecLPP.push_back(LPP(ParamDescr[i])); 192 vector<RingElem> reducer; 193 vector<RingElem> polyx; 194 map<PPMonoidElem,int> FindReducerIndex; 195 196 polyx.push_back(one(Kx)); 197 reducer.push_back(one(Kt)); 198 FindReducerIndex[LPP(reducer[0])] = 1; // deliberately add 1 199 QBGenerator QBG(PPM(Kx)); 200 PPMonoidElem PP1 = QBG.myCorners().front(); 201 QBG.myCornerPPIntoQB(PP1); 202 while (true) 203 { 204 using std::list; 205 // pick corner elem giving smallest LPP in image 206 list<PPMonoidElem>::const_iterator BestIter = QBG.myCorners().begin(); 207 PPMonoidElem BestLPP = ComputeLPP(*BestIter, VecLPP); 208 209 for(list<PPMonoidElem>::const_iterator it=QBG.myCorners().begin(); it != QBG.myCorners().end(); ++it) 210 // for(int i=1;i<len(QBG.myCorners());++i) 211 { 212 const PPMonoidElem tmp = ComputeLPP(*it, VecLPP); 213 if (tmp < BestLPP) { BestLPP = tmp; BestIter = it; } 214 } 215 const PPMonoidElem PP = *BestIter; 216 QBG.myCornerPPIntoQB(PP); 217 RingElem NewReducer = ComputeImage(PP,ParamDescr); 218 RingElem NewPolyx = monomial(Kx, PP); 219 reduce(NewReducer, NewPolyx, FindReducerIndex, reducer, polyx); 220 if (IsZero(NewReducer)) return NewPolyx; 221 reducer.push_back(NewReducer); 222 polyx.push_back(NewPolyx); 223 FindReducerIndex[LPP(NewReducer)] = len(reducer); 224 } 225 } 226 227 228 ChooseIndet(const PPMonoidElem & t,const vector<RingElem> & ParamDescr,std::map<PPMonoidElem,int> & FindPP,const std::vector<RingElem> & polyx)229 int ChooseIndet(const PPMonoidElem& t, const vector<RingElem>& ParamDescr, std::map<PPMonoidElem, int>& FindPP, const std::vector<RingElem>& polyx) 230 { 231 vector<long> exp; exponents(exp, t); 232 const int n = len(exp); // NumIndets(PPM) 233 PPMonoid PPMt = owner(t); 234 int BestI = -1; 235 long BestCost = 0; // should never be used 236 for (int i=0; i<n; ++i) 237 { 238 if (exp[i] == 0) continue; 239 // const PPMonoidElem& x_i = indet(PPMt,i); 240 // const PPMonoidElem tfactor = t/x_i; 241 const PPMonoidElem tfactor = t/indet(PPMt,i); 242 const int j = FindPP[tfactor]; 243 const long EstCost = NumTerms(ParamDescr[i])*NumTerms(polyx[j]); 244 if (BestCost == 0 || BestCost > EstCost) 245 { 246 BestCost = EstCost; 247 BestI = i; 248 } 249 } 250 return BestI; 251 } 252 253 // void ComputeImage(RingElem& NewReducer, RingElem& NewPolyx, const PPMonoidElem& t, const vector<RingElem>& L, const std::map<PPMonoidElem, int>& FindFactorIndex, const std::vector<RingElem>& polyx) 254 // { 255 // /// RingElem ans = one(owner(L[0])); 256 // vector<long> exp; exponents(exp, t); 257 // const int n = len(exp); // NumIndets(PPM) 258 // PPMonoid PPM = owner(t); 259 // int BestI = -1; 260 // long BestCost = 0; // should never be used 261 // for (int i=0; i<n; ++i) 262 // { 263 // if (exp[i] == 0) continue; 264 // const int j = FindFactorIndex(t/indet(PPM,i)); 265 // const long EstCost = NumTerms(L[i])*NumTerms(polyx[j]); 266 // if (BestI >=0 && EstCost > BestCost) continue; 267 // BestCost = EstCost; 268 // BestI = i; 269 // BestFactor = t/indet(PPM,i); 270 // BestPoly = j; 271 // } 272 // const ring& P = owner(L[0]); 273 // NewPolyx = L[BestI]*polyx[BestPoly]; 274 // NewReducer = reducer[FindReducerIndex[BestFactor]]*indet(PP,BestI); 275 // } 276 277 // PPMonoidElem ComputeLPP(const PPMonoidElem& t, const vector<PPMonoidElem>& VecLPP) 278 // { 279 // const PPMonoid& PPM = owner(VecLPP[0]); 280 // PPMonoidElem ans = one(PPM); 281 // const int n = len(VecLPP); 282 // CoCoA_ASSERT(n == NumIndets(owner(t))); 283 // vector<long> exp(n); exponents(exp,t); 284 // for (int i=0;i<n;++i) 285 // ans *= power(VecLPP[i],exp[i]); 286 // return ans; 287 // } 288 289 RingElem ImplicitDirectLPP2Homog(const std::vector<RingElem>& ParamDescrOrig); 290 RingElem ImplicitDirectLPP2HomogBis(const std::vector<RingElem>& ParamDescrOrig); 291 ImplicitDirectLPP2(const std::vector<RingElem> & ParamDescrOrig)292 RingElem ImplicitDirectLPP2(const std::vector<RingElem>& ParamDescrOrig) 293 { 294 if (ParamDescrOrig.empty() || len(ParamDescrOrig) != 1+NumIndets(owner(ParamDescrOrig[0]))) 295 CoCoA_THROW_ERROR(ERR::BadArg, "ImplicitDirectLPP"); 296 if (IsHomog(ParamDescrOrig)) 297 { 298 long i; 299 long d = deg(ParamDescrOrig[0]); 300 for (i=1; i<len(ParamDescrOrig); ++i) 301 if (d != deg(ParamDescrOrig[i])) break; 302 if (i==len(ParamDescrOrig)) return ImplicitDirectLPP2HomogBis(ParamDescrOrig); 303 } 304 const ring& Porig = owner(ParamDescrOrig[0]); 305 const int n = len(ParamDescrOrig); 306 307 ring Kt = NewPolyRingForImplicit(CoeffRing(Porig), "t", NumIndets(Porig)); 308 RingHom phi = PolyAlgebraHom(Porig, Kt, indets(Kt)); 309 vector<RingElem> ParamDescr = apply(phi, ParamDescrOrig); 310 // ring Kx = NewPolyRing(CoeffRing(P), SymbolRange("x",1,n)); 311 ring Kx = NewPolyRingForImplicit(CoeffRing(Porig), "x",n); 312 vector<PPMonoidElem> VecLPP; VecLPP.reserve(len(ParamDescr)); for(int i=0;i<len(ParamDescr);++i)VecLPP.push_back(LPP(ParamDescr[i])); 313 vector<RingElem> reducer; 314 vector<RingElem> polyx; 315 map<PPMonoidElem,int> FindReducerIndex; 316 map<PPMonoidElem,int> FindFactorIndex; 317 318 polyx.push_back(one(Kx)); 319 FindFactorIndex[LPP(polyx[0])] = 0; 320 reducer.push_back(one(Kt)); 321 FindReducerIndex[LPP(reducer[0])] = 1; // deliberately add 1 322 QBGenerator QBG(PPM(Kx)); 323 PPMonoidElem PP1 = QBG.myCorners().front(); 324 QBG.myCornerPPIntoQB(PP1); 325 while (true) 326 { 327 using std::list; 328 // pick corner elem giving smallest LPP in image 329 list<PPMonoidElem>::const_iterator BestIter = QBG.myCorners().begin(); 330 PPMonoidElem BestLPP = ComputeLPP(*BestIter, VecLPP); 331 332 for(list<PPMonoidElem>::const_iterator it=QBG.myCorners().begin(); it != QBG.myCorners().end(); ++it) 333 // for(int i=1;i<len(QBG.myCorners());++i) 334 { 335 const PPMonoidElem tmp = ComputeLPP(*it, VecLPP); 336 if (tmp < BestLPP) { BestLPP = tmp; BestIter = it; } 337 } 338 const PPMonoidElem PP = *BestIter; 339 QBG.myCornerPPIntoQB(PP); 340 const int var = ChooseIndet(PP, ParamDescr, FindFactorIndex, polyx); 341 const int ReducerIndex = FindFactorIndex[PP/indet(PPM(Kx),var)]; 342 RingElem NewPolyx = indet(Kx,var)*polyx[ReducerIndex]; 343 RingElem NewReducer = reducer[ReducerIndex]*ParamDescr[var]; 344 ////// ComputeImage(NewReducer, NewPolyx, PP,ParamDescr, FindFactorIndex, polyx); 345 reduce(NewReducer, NewPolyx, FindReducerIndex, reducer, polyx); 346 if (IsZero(NewReducer)) return NewPolyx; 347 reducer.push_back(NewReducer); 348 polyx.push_back(NewPolyx); 349 FindReducerIndex[LPP(NewReducer)] = len(reducer); 350 FindFactorIndex[PP] = len(reducer)-1; 351 } 352 } 353 354 RingElem ImplicitDirectOrd2HomogBis(const std::vector<RingElem>& ParamDescrOrig); 355 ImplicitDirectOrd2(const std::vector<RingElem> & ParamDescrOrig)356 RingElem ImplicitDirectOrd2(const std::vector<RingElem>& ParamDescrOrig) 357 { 358 if (ParamDescrOrig.empty() || len(ParamDescrOrig) != 1+NumIndets(owner(ParamDescrOrig[0]))) 359 CoCoA_THROW_ERROR(ERR::BadArg, "ImplicitDirectOrd2"); 360 if (IsHomog(ParamDescrOrig)) 361 { 362 long i; 363 long d = deg(ParamDescrOrig[0]); 364 for (i=1; i<len(ParamDescrOrig); ++i) 365 if (d != deg(ParamDescrOrig[i])) break; 366 if (i==len(ParamDescrOrig)) return ImplicitDirectOrd2HomogBis(ParamDescrOrig); 367 } 368 const ring& Porig = owner(ParamDescrOrig[0]); 369 const int n = len(ParamDescrOrig); 370 371 ring Kt = NewPolyRingForImplicit(CoeffRing(Porig), "t", NumIndets(Porig)); 372 RingHom phi = PolyAlgebraHom(Porig, Kt, indets(Kt)); 373 vector<RingElem> ParamDescr = apply(phi, ParamDescrOrig); 374 375 // odering compatible with LPP in Kt 376 vector<long> V(n,0); 377 vector<vector<long> > VV(n,V); 378 for (long i=0; i<n; ++i) exponents(VV[i], LPP(ParamDescr[i])); 379 380 matrix M = MakeTermOrd(StdDegRevLexMat(n-1) * transpose(NewDenseMat(RingZZ(),VV))); 381 // std::cout << M << std::endl; 382 383 /// we should always check that there are no constants in ParamDescr 384 ring Kx = NewPolyRingForImplicit(CoeffRing(Porig), "x", NewMatrixOrdering(M, 1)); 385 386 vector<RingElem> reducer; 387 vector<RingElem> polyx; 388 map<PPMonoidElem,int> FindReducerIndex; 389 map<PPMonoidElem,int> FindFactorIndex; 390 391 polyx.push_back(one(Kx)); 392 FindFactorIndex[LPP(polyx[0])] = 0; 393 reducer.push_back(one(Kt)); 394 FindReducerIndex[LPP(reducer[0])] = 1; // deliberately add 1 395 QBGenerator QBG(PPM(Kx)); 396 PPMonoidElem PP1 = QBG.myCorners().front(); 397 QBG.myCornerPPIntoQB(PP1); 398 while (true) 399 { 400 using std::list; 401 // pick corner elem giving smallest LPP in image 402 const PPMonoidElem PP = QBG.myCorners().front(); // special ord 403 QBG.myCornerPPIntoQB(PP); 404 const int var = ChooseIndet(PP, ParamDescr, FindFactorIndex, polyx); 405 const int ReducerIndex = FindFactorIndex[PP/indet(PPM(Kx),var)]; 406 RingElem NewPolyx = indet(Kx,var)*polyx[ReducerIndex]; 407 RingElem NewReducer = reducer[ReducerIndex]*ParamDescr[var]; 408 //// ComputeImage(NewReducer, NewPolyx, PP,ParamDescr, FindFactorIndex, polyx); 409 reduce(NewReducer, NewPolyx, FindReducerIndex, reducer, polyx); 410 if (IsZero(NewReducer)) return NewPolyx; 411 reducer.push_back(NewReducer); 412 polyx.push_back(NewPolyx); 413 FindReducerIndex[LPP(NewReducer)] = len(reducer); 414 FindFactorIndex[PP] = len(reducer)-1; 415 } 416 } 417 418 ImplicitDirectLPP2Homog(const std::vector<RingElem> & ParamDescrOrig)419 RingElem ImplicitDirectLPP2Homog(const std::vector<RingElem>& ParamDescrOrig) 420 { 421 using std::list; 422 VerboseLog VERBOSE("ImplicitDirectLPP2Homog"); 423 VERBOSE(90) << "start" << std::endl; 424 const ring& Porig = owner(ParamDescrOrig[0]); 425 const int n = len(ParamDescrOrig); 426 427 ring Kt = NewPolyRingForImplicit(CoeffRing(Porig), "t", NumIndets(Porig)); 428 RingHom phi = PolyAlgebraHom(Porig, Kt, indets(Kt)); 429 vector<RingElem> ParamDescr = apply(phi, ParamDescrOrig); 430 ring Kx = NewPolyRingForImplicit(CoeffRing(Porig), "x",n); 431 vector<PPMonoidElem> VecLPP; VecLPP.reserve(len(ParamDescr)); for(int i=0;i<len(ParamDescr);++i)VecLPP.push_back(LPP(ParamDescr[i])); 432 vector<RingElem> reducer; 433 vector<RingElem> polyx; 434 map<PPMonoidElem,int> FindReducerIndex; 435 map<PPMonoidElem,int> FindFactorIndex; 436 437 reducer.push_back(one(Kt)); 438 polyx.push_back(one(Kx)); 439 FindReducerIndex[LPP(reducer[0])] = 1; // deliberately add 1 440 FindFactorIndex[LPP(polyx[0])] = 0; 441 QBGenerator QBG(PPM(Kx)); 442 PPMonoidElem PP1 = QBG.myCorners().front(); 443 QBG.myCornerPPIntoQB(PP1); 444 while (true) 445 { 446 list<PPMonoidElem>::const_iterator BestIter = QBG.myCorners().begin(); 447 PPMonoidElem BestLPP = *BestIter; 448 for(list<PPMonoidElem>::const_iterator it=QBG.myCorners().begin(); it != QBG.myCorners().end(); ++it) 449 { 450 const PPMonoidElem tmp = *it; 451 if (tmp < BestLPP) { BestLPP = tmp; BestIter = it; } 452 } 453 const PPMonoidElem PPx = *BestIter; 454 // if (deg(PPx) != CurrentDeg) 455 // { 456 // CurrentDeg = deg(PPx); 457 // // std::cout << "deg(PPx) = " << deg(PPx) << std::endl; 458 // } 459 QBG.myCornerPPIntoQB(PPx); 460 const int var = ChooseIndet(PPx, ParamDescr, FindFactorIndex, polyx); 461 const int ReducerIndex = FindFactorIndex[PPx/indet(PPM(Kx),var)]; 462 RingElem NewPolyx = polyx[ReducerIndex] *indet(Kx,var); 463 RingElem NewReducer = reducer[ReducerIndex]*ParamDescr[var]; 464 reduce(NewReducer, NewPolyx, FindReducerIndex, reducer, polyx); 465 if (IsZero(NewReducer)) return NewPolyx; 466 reducer.push_back(NewReducer); 467 polyx.push_back(NewPolyx); 468 FindReducerIndex[LPP(NewReducer)] = len(reducer); 469 FindFactorIndex[PPx] = len(reducer)-1; 470 } 471 } 472 473 ImplicitDirectLPP2HomogBis(const std::vector<RingElem> & ParamDescrOrig)474 RingElem ImplicitDirectLPP2HomogBis(const std::vector<RingElem>& ParamDescrOrig) 475 { 476 using std::list; 477 VerboseLog VERBOSE("ImplicitDirectLPP2HomogBis"); 478 VERBOSE(90) << "start" << std::endl; 479 const ring& Porig = owner(ParamDescrOrig[0]); 480 const int n = len(ParamDescrOrig); 481 482 ring Kt = NewPolyRingForImplicit(CoeffRing(Porig), "t", NumIndets(Porig)); 483 RingHom phi = PolyAlgebraHom(Porig, Kt, indets(Kt)); 484 vector<RingElem> ParamDescr = apply(phi, ParamDescrOrig); 485 ring Kx = NewPolyRingForImplicit(CoeffRing(Porig), "x",n); 486 vector<PPMonoidElem> VecLPP; VecLPP.reserve(len(ParamDescr)); for(int i=0;i<len(ParamDescr);++i)VecLPP.push_back(LPP(ParamDescr[i])); 487 vector<RingElem> reducerPrev; 488 vector<RingElem> polyxPrev; 489 map<PPMonoidElem,int> FindReducerIndexPrev; 490 map<PPMonoidElem,int> FindFactorIndexPrev; 491 vector<RingElem> reducer; 492 vector<RingElem> polyx; 493 map<PPMonoidElem,int> FindReducerIndex; 494 map<PPMonoidElem,int> FindFactorIndex; 495 496 reducer.push_back(one(Kt)); 497 polyx.push_back(one(Kx)); 498 reducerPrev.push_back(one(Kt)); 499 polyxPrev.push_back(one(Kx)); 500 FindReducerIndex[LPP(reducer[0])] = 1; // deliberately add 1 501 FindFactorIndex[LPP(polyx[0])] = 0; 502 FindReducerIndexPrev[LPP(reducer[0])] = 1; // deliberately add 1 503 FindFactorIndexPrev[LPP(polyx[0])] = 0; 504 QBGenerator QBG(PPM(Kx)); 505 PPMonoidElem PP1 = QBG.myCorners().front(); 506 QBG.myCornerPPIntoQB(PP1); 507 long CurrentDeg = 0; 508 while (true) 509 { 510 list<PPMonoidElem>::const_iterator BestIter = QBG.myCorners().begin(); 511 PPMonoidElem BestLPP = *BestIter; 512 for(list<PPMonoidElem>::const_iterator it=QBG.myCorners().begin(); it != QBG.myCorners().end(); ++it) 513 { 514 const PPMonoidElem tmp = *it; 515 if (tmp < BestLPP) { BestLPP = tmp; BestIter = it; } 516 } 517 const PPMonoidElem PPx = *BestIter; 518 if (deg(PPx) != CurrentDeg) 519 { 520 CurrentDeg = deg(PPx); 521 reducerPrev.clear(); 522 polyxPrev.clear(); 523 FindReducerIndexPrev.clear(); 524 FindFactorIndexPrev.clear(); 525 reducerPrev.swap(reducer); 526 polyxPrev.swap(polyx); 527 FindReducerIndexPrev.swap(FindReducerIndex); 528 FindFactorIndexPrev.swap(FindFactorIndex); 529 reducer.push_back(one(Kt)); 530 polyx.push_back(one(Kx)); 531 FindReducerIndex[LPP(reducer[0])] = 1; // deliberately add 1 532 FindFactorIndex[LPP(polyx[0])] = 0; 533 } 534 QBG.myCornerPPIntoQB(PPx); 535 const int var = ChooseIndet(PPx, ParamDescr, FindFactorIndexPrev, polyxPrev); 536 const int ReducerIndex = FindFactorIndexPrev[PPx/indet(PPM(Kx),var)]; 537 RingElem NewPolyx = polyxPrev[ReducerIndex] *indet(Kx,var); 538 RingElem NewReducer = reducerPrev[ReducerIndex]*ParamDescr[var]; 539 reduce(NewReducer, NewPolyx, FindReducerIndex, reducer, polyx); 540 if (IsZero(NewReducer)) return NewPolyx; 541 reducer.push_back(NewReducer); 542 polyx.push_back(NewPolyx); 543 FindReducerIndex[LPP(NewReducer)] = len(reducer); 544 FindFactorIndex[PPx] = len(reducer)-1; 545 } 546 } 547 548 ImplicitDirectOrd2HomogBis(const std::vector<RingElem> & ParamDescrOrig)549 RingElem ImplicitDirectOrd2HomogBis(const std::vector<RingElem>& ParamDescrOrig) 550 { 551 using std::list; 552 VerboseLog VERBOSE("ImplicitDirectOrd2HomogBis"); 553 VERBOSE(90) << "start" << std::endl; 554 555 const ring& Porig = owner(ParamDescrOrig[0]); 556 const int n = len(ParamDescrOrig); 557 558 ring Kt = NewPolyRingForImplicit(CoeffRing(Porig), "t", NumIndets(Porig)); 559 RingHom phi = PolyAlgebraHom(Porig, Kt, indets(Kt)); 560 vector<RingElem> ParamDescr = apply(phi, ParamDescrOrig); 561 // ring Kx = NewPolyRingForImplicit(CoeffRing(Porig), "x",n); 562 vector<PPMonoidElem> VecLPP; VecLPP.reserve(len(ParamDescr)); for(int i=0;i<len(ParamDescr);++i)VecLPP.push_back(LPP(ParamDescr[i])); 563 // odering compatible with LPP in Kt 564 vector<long> V(n,0); 565 vector<vector<long> > VV(n,V); 566 for (long i=0; i<n; ++i) exponents(VV[i], LPP(ParamDescr[i])); 567 matrix M = MakeTermOrd(StdDegRevLexMat(n-1) * transpose(NewDenseMat(RingZZ(),VV))); 568 // std::cout << M << std::endl; 569 570 /// we should always check that there are no constants in ParamDescr 571 ring Kx = NewPolyRingForImplicit(CoeffRing(Porig), "x", NewMatrixOrdering(M, 1)); 572 573 vector<RingElem> reducerPrev; 574 vector<RingElem> polyxPrev; 575 map<PPMonoidElem,int> FindReducerIndexPrev; 576 map<PPMonoidElem,int> FindFactorIndexPrev; 577 vector<RingElem> reducer; 578 vector<RingElem> polyx; 579 map<PPMonoidElem,int> FindReducerIndex; 580 map<PPMonoidElem,int> FindFactorIndex; 581 582 reducer.push_back(one(Kt)); 583 polyx.push_back(one(Kx)); 584 reducerPrev.push_back(one(Kt)); 585 polyxPrev.push_back(one(Kx)); 586 FindReducerIndex[LPP(reducer[0])] = 1; // deliberately add 1 587 FindFactorIndex[LPP(polyx[0])] = 0; 588 FindReducerIndexPrev[LPP(reducer[0])] = 1; // deliberately add 1 589 FindFactorIndexPrev[LPP(polyx[0])] = 0; 590 QBGenerator QBG(PPM(Kx)); 591 PPMonoidElem PP1 = QBG.myCorners().front(); 592 QBG.myCornerPPIntoQB(PP1); 593 degree CurrentDeg = wdeg(PP1); 594 while (true) 595 { 596 // list<PPMonoidElem>::const_iterator BestIter = QBG.myCorners().begin(); 597 // PPMonoidElem BestLPP = *BestIter; 598 // for(list<PPMonoidElem>::const_iterator it=QBG.myCorners().begin(); it != QBG.myCorners().end(); ++it) 599 // { 600 // const PPMonoidElem tmp = *it; 601 // if (tmp < BestLPP) { BestLPP = tmp; BestIter = it; } 602 // } 603 // // const PPMonoidElem PPx = *BestIter; 604 // if (*BestIter != QBG.myCorners().front()) 605 // { 606 // std::cout << "ppx is " << QBG.myCorners().front() << 607 // << " *BestIter is " << *BestIter << std::endl; 608 // } 609 const PPMonoidElem PPx = QBG.myCorners().front(); // special ord 610 if (wdeg(PPx) != CurrentDeg) 611 { 612 CurrentDeg = wdeg(PPx); 613 reducerPrev.clear(); 614 polyxPrev.clear(); 615 FindReducerIndexPrev.clear(); 616 FindFactorIndexPrev.clear(); 617 reducerPrev.swap(reducer); 618 polyxPrev.swap(polyx); 619 FindReducerIndexPrev.swap(FindReducerIndex); 620 FindFactorIndexPrev.swap(FindFactorIndex); 621 reducer.push_back(one(Kt)); 622 polyx.push_back(one(Kx)); 623 FindReducerIndex[LPP(reducer[0])] = 1; // deliberately add 1 624 FindFactorIndex[LPP(polyx[0])] = 0; 625 } 626 QBG.myCornerPPIntoQB(PPx); 627 const int var = ChooseIndet(PPx, ParamDescr, FindFactorIndexPrev, polyxPrev); 628 const int ReducerIndex = FindFactorIndexPrev[PPx/indet(PPM(Kx),var)]; 629 RingElem NewPolyx = polyxPrev[ReducerIndex] *indet(Kx,var); 630 RingElem NewReducer = reducerPrev[ReducerIndex]*ParamDescr[var]; 631 reduce(NewReducer, NewPolyx, FindReducerIndex, reducer, polyx); 632 if (IsZero(NewReducer)) return NewPolyx; 633 reducer.push_back(NewReducer); 634 polyx.push_back(NewPolyx); 635 FindReducerIndex[LPP(NewReducer)] = len(reducer); 636 FindFactorIndex[PPx] = len(reducer)-1; 637 } 638 } 639 640 641 642 643 // Procedure for curves. but there is no stopping criterion :-( 644 // std::vector<RingElem> BM_param(const std::vector<RingElem>& ParamDescrOrig) 645 // { 646 // if (ParamDescrOrig.empty()) 647 // CoCoA_THROW_ERROR(ERR::BadArg, "BM_param"); 648 // const ring& Porig = owner(ParamDescrOrig[0]); 649 // const int n = len(ParamDescrOrig); 650 651 // ring Kt = NewPolyRingForImplicit(CoeffRing(Porig), "t", NumIndets(Porig)); 652 // RingHom phi = PolyAlgebraHom(Porig, Kt, indets(Kt)); 653 // vector<RingElem> ParamDescr = apply(phi, ParamDescrOrig); 654 // ring Kx = NewPolyRingForImplicit(CoeffRing(Porig), "x",n); 655 // vector<PPMonoidElem> VecLPP; VecLPP.reserve(len(ParamDescr)); 656 // for (int i=0;i<n;++i) VecLPP.push_back(LPP(ParamDescr[i])); 657 // vector<RingElem> reducer; 658 // vector<RingElem> polyx; 659 // map<PPMonoidElem,int> FindReducerIndex; 660 // map<PPMonoidElem,int> FindFactorIndex; 661 662 // reducer.push_back(one(Kt)); 663 // polyx.push_back(one(Kx)); 664 // FindReducerIndex[LPP(reducer[0])] = 1; // deliberately add 1 665 // FindFactorIndex[LPP(polyx[0])] = 0; 666 // QBGenerator QBG(PPM(Kx)); 667 // PPMonoidElem PP1 = QBG.myCorners().front(); 668 // QBG.myCornerPPIntoQB(PP1); 669 // vector<RingElem> GB(0); 670 // while (!QBG.myCorners().empty()) 671 // { 672 // std::cout << "QBG = " << QBG << std::endl; 673 674 // using std::list; 675 // const PPMonoidElem t = QBG.myCorners().front(); 676 // const int var = ChooseIndet(t, ParamDescr, FindFactorIndex, polyx); 677 // const int ReducerIndex = FindFactorIndex[t/indet(PPM(Kx),var)]; 678 // RingElem NewPolyx = indet(Kx,var)*polyx[ReducerIndex]; 679 // RingElem NewReducer = reducer[ReducerIndex]*ParamDescr[var]; 680 // reduce(NewReducer, NewPolyx, FindReducerIndex, reducer, polyx); 681 // // if (IsZero(NewReducer)) return NewPolyx; 682 // if (IsZero(NewReducer)) 683 // { 684 // QBG.myCornerPPIntoAvoidSet(t); 685 // GB.push_back(NewPolyx); 686 // std::cout << "--GB = " << GB << std::endl; 687 // } 688 // else 689 // { 690 // QBG.myCornerPPIntoQB(t); 691 // reducer.push_back(NewReducer); 692 // polyx.push_back(NewPolyx); 693 // FindReducerIndex[LPP(NewReducer)] = len(reducer); 694 // FindFactorIndex[t] = len(reducer)-1; 695 // } 696 // } 697 // return GB; 698 // } 699 700 ImplicitDirectWithCond(const std::vector<RingElem> & ParamDescrOrig,const std::vector<RingElem> & cond)701 RingElem ImplicitDirectWithCond(const std::vector<RingElem>& ParamDescrOrig, const std::vector<RingElem>& cond) 702 { 703 if (ParamDescrOrig.empty() || len(ParamDescrOrig)+len(cond) != 1+NumIndets(owner(ParamDescrOrig[0]))) 704 CoCoA_THROW_ERROR(ERR::BadArg, "ImplicitDirectWithCond"); 705 if (cond.empty()) return ImplicitDirect(ParamDescrOrig); 706 const ring& Porig = owner(ParamDescrOrig[0]); 707 const int n = len(ParamDescrOrig); 708 709 // ring Kt = NewPolyRingForImplicit(CoeffRing(Porig), "t",NumIndets(Porig)); 710 // RingHom phi = PolyAlgebraHom(Porig, Kt, indets(Kt)); 711 // vector<RingElem> ParamDescr = apply(phi, ParamDescrOrig); 712 // vector<RingElem> relations = apply(phi, cond); 713 ring Kt = Porig; 714 vector<RingElem> ParamDescr = ParamDescrOrig; 715 vector<RingElem> relations = cond; 716 ideal RelationIdeal(relations); 717 // Originally the code was like this (ie. compute in the originally given ring) 718 // if (ParamDescr.empty() || len(ParamDescr) != 1+NumIndets(owner(ParamDescr[0]))) 719 // CoCoA_THROW_ERROR(ERR::BadArg, "ImplicitDirect"); 720 // const ring& P = owner(ParamDescr[0]); 721 // const int n = len(ParamDescr); 722 ring Kx = NewPolyRingForImplicit(CoeffRing(Kt), "x",n); 723 vector<RingElem> reducer; 724 vector<RingElem> polyx; 725 map<PPMonoidElem,int> FindReducerIndex; 726 727 reducer.push_back(one(Kt)); 728 polyx.push_back(one(Kx)); 729 FindReducerIndex[LPP(reducer[0])] = 1; // deliberately add 1 730 QBGenerator QBG(PPM(Kx)); 731 // PPMonoidElem PP1 = QBG.myCorners().front(); 732 // QBG.myCornerPPIntoQB(PP1); 733 QBG.myCornerPPIntoQB(QBG.myCorners().front()); 734 while (true) 735 { 736 const PPMonoidElem t = QBG.myCorners().front(); 737 QBG.myCornerPPIntoQB(t); 738 RingElem NewReducer = NF(ComputeImage(t,ParamDescr), RelationIdeal); 739 RingElem NewPolyx = monomial(Kx, t); 740 reduce(NewReducer, NewPolyx, FindReducerIndex, reducer, polyx); 741 if (IsZero(NewReducer)) return NewPolyx; 742 reducer.push_back(NewReducer); 743 polyx.push_back(NewPolyx); 744 FindReducerIndex[LPP(NewReducer)] = len(reducer); 745 } 746 } 747 748 ImplicitDirectWithCondLPP(const std::vector<RingElem> & ParamDescrOrig,const std::vector<RingElem> & cond)749 RingElem ImplicitDirectWithCondLPP(const std::vector<RingElem>& ParamDescrOrig, const std::vector<RingElem>& cond) 750 { 751 if (ParamDescrOrig.empty()) 752 CoCoA_THROW_ERROR(ERR::Empty, "ImplicitDirectWithCondLPP"); 753 if (len(ParamDescrOrig)+len(cond) != 1+NumIndets(owner(ParamDescrOrig[0]))) 754 CoCoA_THROW_ERROR(ERR::BadArg, "ImplicitDirectWithCondLPP"); 755 if (cond.empty()) return ImplicitDirectLPP(ParamDescrOrig); 756 ring Porig = owner(ParamDescrOrig[0]); 757 const int n = len(ParamDescrOrig); 758 759 ring Kt = NewPolyRingForImplicit(CoeffRing(Porig), "t", NumIndets(Porig)); 760 RingHom phi = PolyAlgebraHom(Porig, Kt, indets(Kt)); 761 vector<RingElem> ParamDescr = apply(phi, ParamDescrOrig); 762 vector<RingElem> relations = apply(phi, cond); 763 ideal RelationIdeal(relations); 764 // Originally the code was like this (ie. compute in the originally given ring) 765 // if (ParamDescr.empty() || len(ParamDescr) != 1+NumIndets(owner(ParamDescr[0]))) 766 // CoCoA_THROW_ERROR(ERR::BadArg, "ImplicitDirect"); 767 // const ring& P = owner(ParamDescr[0]); 768 // const int n = len(ParamDescr); 769 ring Kx = NewPolyRingForImplicit(CoeffRing(Kt), "x",n); 770 vector<PPMonoidElem> VecLPP; VecLPP.reserve(len(ParamDescr)); for(int i=0;i<len(ParamDescr);++i)VecLPP.push_back(LPP(ParamDescr[i])); 771 vector<RingElem> reducer; 772 vector<RingElem> polyx; 773 map<PPMonoidElem,int> FindReducerIndex; 774 775 reducer.push_back(one(Kt)); 776 polyx.push_back(one(Kx)); 777 FindReducerIndex[LPP(reducer[0])] = 1; // deliberately add 1 778 QBGenerator QBG(PPM(Kx)); 779 PPMonoidElem PP1 = QBG.myCorners().front(); 780 QBG.myCornerPPIntoQB(PP1); 781 while (true) 782 { 783 using std::list; 784 // pick corner elem giving smallest LPP in image 785 list<PPMonoidElem>::const_iterator BestIter = QBG.myCorners().begin(); 786 PPMonoidElem BestLPP = ComputeLPP(*BestIter, VecLPP); 787 788 for(list<PPMonoidElem>::const_iterator it=QBG.myCorners().begin(); it != QBG.myCorners().end(); ++it) 789 // for(int i=1;i<len(QBG.myCorners());++i) 790 { 791 const PPMonoidElem tmp = ComputeLPP(*it, VecLPP); 792 if (tmp < BestLPP) { BestLPP = tmp; BestIter = it; } 793 } 794 const PPMonoidElem t = *BestIter; 795 QBG.myCornerPPIntoQB(t); 796 RingElem NewReducer = NF(ComputeImage(t,ParamDescr), RelationIdeal); 797 RingElem NewPolyx = monomial(Kx, t); 798 reduce(NewReducer, NewPolyx, FindReducerIndex, reducer, polyx); 799 if (IsZero(NewReducer)) return NewPolyx; 800 reducer.push_back(NewReducer); 801 polyx.push_back(NewPolyx); 802 FindReducerIndex[LPP(NewReducer)] = len(reducer); 803 } 804 } 805 806 ImplicitDirectWithCondOrd2(const std::vector<RingElem> & ParamDescrOrig,const std::vector<RingElem> & cond)807 RingElem ImplicitDirectWithCondOrd2(const std::vector<RingElem>& ParamDescrOrig, const std::vector<RingElem>& cond) 808 { 809 if (ParamDescrOrig.empty()) 810 CoCoA_THROW_ERROR(ERR::Empty, "ImplicitDirectWithCondOrd2"); 811 if (len(ParamDescrOrig)+len(cond) != 1+NumIndets(owner(ParamDescrOrig[0]))) 812 CoCoA_THROW_ERROR(ERR::BadArg, "ImplicitDirectWithCondOrd2"); 813 if (cond.empty()) return ImplicitDirectOrd2(ParamDescrOrig); 814 CoCoA_THROW_ERROR(ERR::NYI, "ImplicitDirectOrd2"); 815 return zero(RingZZ()); // just to keep the compiler quiet 816 } 817 818 ImplicitDirectWithCondLPP2(const std::vector<RingElem> & ParamDescrOrig,const std::vector<RingElem> & cond)819 RingElem ImplicitDirectWithCondLPP2(const std::vector<RingElem>& ParamDescrOrig, const std::vector<RingElem>& cond) 820 { 821 if (ParamDescrOrig.empty()) 822 CoCoA_THROW_ERROR(ERR::Empty, "ImplicitDirectWithCondLPP2"); 823 if (len(ParamDescrOrig)+len(cond) != 1+NumIndets(owner(ParamDescrOrig[0]))) 824 CoCoA_THROW_ERROR(ERR::BadArg, "ImplicitDirectWithCondLPP2"); 825 if (cond.empty()) return ImplicitDirectLPP2(ParamDescrOrig); 826 const ring& Porig = owner(ParamDescrOrig[0]); 827 const int n = len(ParamDescrOrig); 828 829 // faster if creating a new ring (for locality?!?) 830 ring Kt = NewPolyRingForImplicit(CoeffRing(Porig), "t", NumIndets(Porig)); 831 RingHom phi = PolyAlgebraHom(Porig, Kt, indets(Kt)); 832 vector<RingElem> ParamDescr = apply(phi, ParamDescrOrig); 833 vector<RingElem> relations = apply(phi, cond); 834 // ring Kt = Porig; 835 // const vector<RingElem>& ParamDescr = ParamDescrOrig; 836 // const vector<RingElem>& relations = cond; 837 ideal RelationIdeal(relations); 838 ring Kx = NewPolyRingForImplicit(CoeffRing(Kt), "x", n); 839 vector<PPMonoidElem> VecLPP; 840 VecLPP.reserve(len(ParamDescr)); 841 for(int i=0;i<len(ParamDescr);++i) VecLPP.push_back(LPP(ParamDescr[i])); 842 vector<RingElem> reducer; 843 vector<RingElem> polyx; 844 map<PPMonoidElem,int> FindReducerIndex; 845 map<PPMonoidElem,int> FindFactorIndex; // -- LPP2 846 847 polyx.push_back(one(Kx)); 848 reducer.push_back(one(Kt)); 849 FindReducerIndex[LPP(reducer[0])] = 1; // deliberately add 1 850 QBGenerator QBG(PPM(Kx)); 851 PPMonoidElem PP1 = QBG.myCorners().front(); 852 QBG.myCornerPPIntoQB(PP1); 853 while (true) 854 { 855 using std::list; 856 // pick corner elem giving smallest LPP in image 857 list<PPMonoidElem>::const_iterator BestIter = QBG.myCorners().begin(); 858 PPMonoidElem BestLPP = ComputeLPP(*BestIter, VecLPP); 859 860 for(list<PPMonoidElem>::const_iterator it=QBG.myCorners().begin(); it != QBG.myCorners().end(); ++it) 861 // for(int i=1;i<len(QBG.myCorners());++i) 862 { 863 const PPMonoidElem tmp = ComputeLPP(*it, VecLPP); 864 if (tmp < BestLPP) { BestLPP = tmp; BestIter = it; } 865 } 866 const PPMonoidElem PP = *BestIter; 867 QBG.myCornerPPIntoQB(PP); 868 const int var = ChooseIndet(PP, ParamDescr, FindFactorIndex, reducer); 869 const int ReducerIndex = FindFactorIndex[PP/indet(PPM(Kx),var)]; 870 RingElem NewPolyx = polyx[ReducerIndex]*indet(Kx,var); 871 RingElem NewReducer = NF(reducer[ReducerIndex]*ParamDescr[var], RelationIdeal); 872 reduce(NewReducer, NewPolyx, FindReducerIndex, reducer, polyx); 873 if (IsZero(NewReducer)) return NewPolyx; 874 reducer.push_back(NewReducer); 875 polyx.push_back(NewPolyx); 876 FindReducerIndex[LPP(NewReducer)] = len(reducer); 877 FindFactorIndex[PP] = len(reducer)-1; 878 } 879 } 880 881 882 883 //------------------------------------------------------- 884 // ImplicitByPoints 885 886 namespace // anonymous for file local fns 887 { 888 eval(const RingElem & f,const vector<RingElem> & pt)889 RingElem eval(const RingElem& f, const vector<RingElem>& pt) 890 { 891 const ring& P = owner(f); 892 RingElem ans = zero(CoeffRing(P)); 893 RingHom phi = CoeffEmbeddingHom(P); 894 for (SparsePolyIter it = BeginIter(f); !IsEnded(it); ++it) 895 { 896 RingElem ThisTerm = coeff(it); 897 for (int i=0; i < NumIndets(P); ++i) 898 ThisTerm *= power(pt[i], exponent(PP(it),i)); 899 ans += ThisTerm; 900 } 901 return ans; 902 } 903 904 // currently just use random points GetNewSamplePt(const vector<RingElem> & ParamDescr,const vector<vector<RingElem>> & SamplePts)905 vector<RingElem> GetNewSamplePt(const vector<RingElem>& ParamDescr, const vector< vector<RingElem> >& SamplePts) 906 { 907 const ring& R = owner(ParamDescr[0]); 908 const int NumParams = NumIndets(R); 909 vector<RingElem> params(NumParams); 910 vector<RingElem> NewSamplePt(len(ParamDescr)); 911 TryAgain: 912 for (int i=0; i < NumParams; ++i) 913 { 914 params[i] = RingElem(CoeffRing(R),RandomLong(-999,999)); 915 } 916 for (int i=0; i < len(ParamDescr); ++i) 917 NewSamplePt[i] = eval(ParamDescr[i], params); 918 if (find(SamplePts.begin(), SamplePts.end(), NewSamplePt) != SamplePts.end()) goto TryAgain; 919 return NewSamplePt; 920 } 921 922 } // end of anonymous namespace 923 ImplicitByPoints(const std::vector<RingElem> & ParamDescr)924 RingElem ImplicitByPoints(const std::vector<RingElem>& ParamDescr) 925 { 926 if (ParamDescr.empty() || len(ParamDescr) != 1+NumIndets(owner(ParamDescr[0]))) 927 CoCoA_THROW_ERROR(ERR::BadArg, "ImplicitByPoints"); 928 const ring& P = owner(ParamDescr[0]); 929 const long n = len(ParamDescr); 930 ring Kx = NewPolyRing(CoeffRing(P), SymbolRange("x", 1, n)); 931 RingHom phi = CoeffEmbeddingHom(Kx); 932 933 QBGenerator QBG(PPM(Kx)); 934 const PPMonoidElem PP1 = QBG.myCorners().front(); 935 QBG.myCornerPPIntoQB(PP1); 936 937 const RingElem one = monomial(Kx, PP1); 938 vector<RingElem> polyx; polyx.push_back(one); 939 vector< vector<RingElem> > RowReducers(1); // essentially the BM matrix 940 941 vector< vector<RingElem> > SamplePts; // should be a std::set!!! 942 int NumPts = 0; 943 while (true) 944 { 945 ++NumPts; 946 const vector<RingElem> NewSamplePt = GetNewSamplePt(ParamDescr, SamplePts); 947 SamplePts.push_back(NewSamplePt); 948 949 for (long i=0; i < NumPts; ++i) 950 { 951 RowReducers[i].push_back(eval(polyx[i], NewSamplePt)); // clear but inefficient 952 } 953 if (IsZero(RowReducers[NumPts-1][NumPts-1])) return polyx[NumPts-1]; 954 const PPMonoidElem t = QBG.myCorners().front(); 955 QBG.myCornerPPIntoQB(t); 956 RingElem NewPolyx = monomial(Kx, t); 957 vector<RingElem> NewRow; 958 for(long i=0; i < NumPts ; ++i) 959 NewRow.push_back(eval(NewPolyx, SamplePts[i])); 960 961 // just gaussian reduction -- we know that matrix is upper triangular 962 for (long i=0; i < NumPts; ++i) 963 { 964 if (IsZero(NewRow[i])) continue; 965 const RingElem c = NewRow[i]/RowReducers[i][i]; 966 NewPolyx -= phi(c)*polyx[i]; 967 for (long j=i; j < NumPts; ++j) 968 NewRow[j] -= c*RowReducers[i][j]; 969 } 970 polyx.push_back(NewPolyx); 971 RowReducers.push_back(NewRow); 972 } 973 } 974 975 976 ImplicitByPoints2(const std::vector<RingElem> & ParamDescr)977 RingElem ImplicitByPoints2(const std::vector<RingElem>& ParamDescr) 978 { 979 using std::list; 980 if (ParamDescr.empty() || len(ParamDescr) != 1+NumIndets(owner(ParamDescr[0]))) 981 CoCoA_THROW_ERROR(ERR::BadArg, "ImplicitByPoints"); 982 const ring& P = owner(ParamDescr[0]); 983 const long p = ConvertTo<long>(characteristic(P)); 984 CoCoA_ASSERT(IsPrime(p)); 985 SmallFpImpl ModP(p); 986 const long n = len(ParamDescr); 987 ring Kx = NewPolyRing(CoeffRing(P), SymbolRange("x", 1, n)); 988 RingHom phi = CoeffEmbeddingHom(Kx); 989 vector<PPMonoidElem> VecLPP; VecLPP.reserve(len(ParamDescr)); for(int i=0;i<len(ParamDescr);++i)VecLPP.push_back(LPP(ParamDescr[i])); 990 991 QBGenerator QBG(PPM(Kx)); 992 const PPMonoidElem PP1 = QBG.myCorners().front(); 993 QBG.myCornerPPIntoQB(PP1); 994 995 const RingElem one = monomial(Kx, PP1); 996 vector<RingElem> polyx; polyx.push_back(one); 997 ///JAA vector< vector<SmallFpImpl::value_t> > RowReducers(1); // essentially the BM matrix 998 vector< vector<SmallFpImpl::value> > RowReducers2(1); // essentially the BM matrix 999 1000 vector< vector<RingElem> > SamplePts; // should be a std::set!!! 1001 int NumPts = 0; 1002 1003 while (true) 1004 { 1005 ++NumPts; 1006 const vector<RingElem> NewSamplePt = GetNewSamplePt(ParamDescr, SamplePts); 1007 SamplePts.push_back(NewSamplePt); 1008 1009 for (long i=0; i < NumPts; ++i) 1010 { 1011 ///JAA RowReducers[i].push_back(ModP.myReduce(ConvertTo<long>(eval(polyx[i], NewSamplePt)))); // clear but inefficient 1012 RowReducers2[i].push_back(ModP.myReduce(ConvertTo<long>(eval(polyx[i], NewSamplePt)))); // clear but inefficient 1013 } 1014 ///JAA CoCoA_ASSERT(ModP.myExport(RowReducers[NumPts-1][NumPts-1]) == ModP.myExport(RowReducers2[NumPts-1][NumPts-1])); 1015 if (IsZero(RowReducers2[NumPts-1][NumPts-1])) return polyx[NumPts-1]; 1016 1017 // pick corner elem giving smallest LPP in image 1018 list<PPMonoidElem>::const_iterator BestIter = QBG.myCorners().begin(); 1019 PPMonoidElem BestLPP = ComputeLPP(*BestIter, VecLPP); 1020 1021 for(list<PPMonoidElem>::const_iterator it=QBG.myCorners().begin(); it != QBG.myCorners().end(); ++it) 1022 // for(int i=1;i<len(QBG.myCorners());++i) 1023 { 1024 const PPMonoidElem tmp = ComputeLPP(*it, VecLPP); 1025 if (tmp < BestLPP) { BestLPP = tmp; BestIter = it; } 1026 } 1027 const PPMonoidElem t = *BestIter; 1028 QBG.myCornerPPIntoQB(t); 1029 RingElem NewPolyx = monomial(Kx, t); 1030 ///JAA vector<SmallFpImpl::value_t> NewRow; 1031 vector<SmallFpImpl::NonRedValue> NewRow2; 1032 for(long i=0; i < NumPts ; ++i) 1033 { 1034 ///JAA NewRow.push_back(ModP.myReduce(ConvertTo<long>(eval(NewPolyx, SamplePts[i])))); 1035 NewRow2.push_back(ModP.myReduce(ConvertTo<long>(eval(NewPolyx, SamplePts[i])))); 1036 ///JAA CoCoA_ASSERT(ModP.myExport(NewRow[i]) == ModP.myExport(ModP.myNormalize(NewRow2[i]))); 1037 } 1038 1039 // just gaussian reduction -- we know that matrix is upper triangular 1040 long IterCount = 0; 1041 for (long i=0; i < NumPts; ++i) 1042 { 1043 ///JAA NewRow[i] = ModP.myNormalize(NewRow[i]); 1044 const SmallFpImpl::value coordi = ModP.myNormalize(NewRow2[i]); 1045 ///JAA NewRow2[i] = coordi; /// USELESS??? 1046 if (IsZero(coordi)) continue; 1047 ///JAA const SmallFpImpl::value_t c = ModP.myDiv(NewRow[i], RowReducers[i][i]); 1048 const SmallFpImpl::value c2 = ModP.myNegate(ModP.myDiv(coordi, RowReducers2[i][i])); 1049 ///JAA CoCoA_ASSERT(ModP.myExport(ModP.myNegate(c)) == ModP.myExport(c2)); 1050 ///JAA NewRow[i] = 0; 1051 NewRow2[i] = zero(SmallFp); // USELESS???? 1052 ///JAA NewPolyx -= RingElem(Kx,c)*polyx[i]; 1053 NewPolyx += ModP.myExportNonNeg(c2)*polyx[i]; 1054 for (long j=i+1; j < NumPts; ++j) 1055 { 1056 ///JAA NewRow[j] += (p-c)*RowReducers[i][j]; 1057 NewRow2[j] += c2 * RowReducers2[i][j]; 1058 // NewRow2[j] = add(NewRow2[j], mul(c2,RowReducers2[i][j])); 1059 } 1060 ++IterCount; 1061 if (IterCount == ModP.myMaxIters()) 1062 { 1063 IterCount = 0; 1064 for (long j=i+1; j < NumPts; ++j) 1065 { 1066 ///JAA NewRow[j] = ModP.myHalfNormalize(NewRow[j]); 1067 NewRow2[j] = ModP.myHalfNormalize(NewRow2[j]); 1068 ///JAA CoCoA_ASSERT(ModP.myExport(ModP.myNormalize(NewRow[j])) == ModP.myExport(ModP.myNormalize(NewRow2[j]))); 1069 } 1070 } 1071 } 1072 polyx.push_back(NewPolyx); 1073 ///JAA NewRow[NumPts-1] = ModP.myNormalize(NewRow[NumPts-1]); 1074 NewRow2[NumPts-1] = ModP.myNormalize(NewRow2[NumPts-1]); 1075 ///JAA CoCoA_ASSERT(ModP.myExport(NewRow[NumPts-1]) == ModP.myExport(ModP.myNormalize(NewRow2[NumPts-1]))); 1076 // CoCoA_ASSERT(NewRow[NumPts-1] == 0); 1077 ///JAA RowReducers.push_back(NewRow); 1078 RowReducers2.push_back(vector<SmallFpImpl::value>(NumPts)); 1079 } 1080 } 1081 1082 1083 1084 1085 // currently just use random points GetNewSamplePt3(const SmallFpImpl & ModP,const vector<RingElem> & ParamDescr,const vector<vector<SmallFpImpl::value>> & SamplePts)1086 vector<SmallFpImpl::value> GetNewSamplePt3(const SmallFpImpl& ModP, const vector<RingElem>& ParamDescr, const vector< vector<SmallFpImpl::value> >& SamplePts) 1087 { 1088 const ring& R = owner(ParamDescr[0]); 1089 const int NumParams = NumIndets(R); 1090 vector<RingElem> params(NumParams); 1091 vector<RingElem> ImagePt(len(ParamDescr)); 1092 vector<SmallFpImpl::value> CandidatePt(len(ParamDescr)); 1093 TryAgain: 1094 for (int i=0; i < NumParams; ++i) 1095 { 1096 params[i] = RingElem(CoeffRing(R),RandomLong(0,9999)); 1097 } 1098 for (int i=0; i < len(ParamDescr); ++i) 1099 ImagePt[i] = eval(ParamDescr[i], params); 1100 for (int i=0; i < len(ParamDescr); ++i) 1101 { 1102 CandidatePt[i] = ModP.myReduce(ConvertTo<long>(ImagePt[i])); 1103 } 1104 if (find(SamplePts.begin(), SamplePts.end(), CandidatePt) != SamplePts.end()) goto TryAgain; 1105 return CandidatePt; 1106 } 1107 ImplicitByPoints3(const std::vector<RingElem> & ParamDescr)1108 RingElem ImplicitByPoints3(const std::vector<RingElem>& ParamDescr) 1109 { 1110 using std::list; 1111 if (ParamDescr.empty() || len(ParamDescr) != 1+NumIndets(owner(ParamDescr[0]))) 1112 CoCoA_THROW_ERROR(ERR::BadArg, "ImplicitByPoints"); 1113 const ring& P = owner(ParamDescr[0]); 1114 const long p = ConvertTo<long>(characteristic(P)); 1115 CoCoA_ASSERT(IsPrime(p)); 1116 SmallFpImpl ModP(p); 1117 const long n = len(ParamDescr); 1118 ring Kx = NewPolyRing(CoeffRing(P), SymbolRange("x", 1, n)); 1119 RingHom phi = CoeffEmbeddingHom(Kx); 1120 vector<PPMonoidElem> VecLPP; VecLPP.reserve(len(ParamDescr)); for(int i=0;i<len(ParamDescr);++i)VecLPP.push_back(LPP(ParamDescr[i])); 1121 map<PPMonoidElem,int> FindFactorIndex; 1122 1123 QBGenerator QBG(PPM(Kx)); 1124 const PPMonoidElem PP1 = QBG.myCorners().front(); 1125 QBG.myCornerPPIntoQB(PP1); 1126 vector<long> QBFactorPP; QBFactorPP.push_back(-1); 1127 vector<long> QBFactorVar; QBFactorVar.push_back(-1); 1128 1129 //POLY const RingElem one = monomial(Kx, PP1); 1130 //POLY vector<RingElem> polyx; polyx.push_back(one); 1131 FindFactorIndex[one(PPM(Kx))] = 0; 1132 vector< vector<SmallFpImpl::value> > PolyCoeffs; PolyCoeffs.push_back(vector<SmallFpImpl::value>(1,ModP.myReduce(1))); 1133 1134 vector< vector<SmallFpImpl::value> > RowReducers(1); // essentially the BM matrix 1135 1136 vector< vector<SmallFpImpl::value> > SamplePts; // should be a std::set!!! 1137 int NumPts = 0; 1138 1139 while (true) 1140 { 1141 ++NumPts; 1142 // if (NumPts%10 == 0) clog<<"NumPts="<<NumPts<<endl; 1143 const vector<SmallFpImpl::value> NewSamplePt = GetNewSamplePt3(ModP, ParamDescr, SamplePts); 1144 SamplePts.push_back(NewSamplePt); 1145 1146 //POLY vector<RingElem> NewPt(n); for(int i=0;i<n;++i)NewPt[i]=RingElem(CoeffRing(P),ModP.myExport(NewSamplePt[i])); 1147 1148 vector<SmallFpImpl::value> QBvalue(len(QBG.myQB())); 1149 QBvalue[0] = ModP.myReduce(1); 1150 for (long i=1; i < len(QBG.myQB()); ++i) 1151 QBvalue[i] = ModP.myMul(QBvalue[QBFactorPP[i]], NewSamplePt[QBFactorVar[i]]); 1152 1153 for (long i=0; i < NumPts; ++i) 1154 { 1155 //POLY RowReducers[i].push_back(ModP.myReduce(ConvertTo<long>(eval(polyx[i], NewSamplePt)))); // clear but inefficient 1156 SmallFpImpl::NonRedValue EvalAtNewPt = PolyCoeffs[i][0]; 1157 long IterCount=0; 1158 const long ReduceNow = ModP.myMaxIters(); 1159 for (int j=1;j<len(PolyCoeffs[i]);++j) 1160 { 1161 /// EvalAtNewPt = add(EvalAtNewPt, mul(QBvalue[j], PolyCoeffs[i][j])); 1162 EvalAtNewPt += QBvalue[j] * PolyCoeffs[i][j]; 1163 ++IterCount; 1164 if (IterCount == ReduceNow) { EvalAtNewPt = ModP.myHalfNormalize(EvalAtNewPt); IterCount = 0; } 1165 // EvalAtNewPt = ModP.myAdd(EvalAtNewPt, ModP.myMul(QBvalue[j],PolyCoeffs[i][j])); 1166 } 1167 1168 RowReducers[i].push_back(ModP.myNormalize(EvalAtNewPt)); 1169 //// RowReducers[i].push_back(ModP.myReduce(ConvertTo<long>(eval(polyx[i], NewPt)))); // clear but inefficient 1170 //// if (ModP.myExport(EvalAtNewPt) != ModP.myExport(RowReducers[i].back())) CoCoA_THROW_ERROR("BUNGLED", "eval check"); 1171 } 1172 ///JAA CoCoA_ASSERT(ModP.myExport(RowReducers[NumPts-1][NumPts-1]) == ModP.myExport(RowReducers2[NumPts-1][NumPts-1])); 1173 // if (IsZero(RowReducers[NumPts-1][NumPts-1])) return polyx[NumPts-1]; 1174 if (IsZero(RowReducers[NumPts-1][NumPts-1])) 1175 { 1176 // convert answer to a polynomial 1177 RingElem ans(Kx); 1178 for (long i=0; i < len(PolyCoeffs[NumPts-1]); ++i) 1179 ans += monomial(Kx, ModP.myExportNonNeg(PolyCoeffs[NumPts-1][i]), QBG.myQB()[i]); 1180 return ans; 1181 } 1182 1183 // pick corner elem giving smallest LPP in image 1184 list<PPMonoidElem>::const_iterator BestIter = QBG.myCorners().begin(); 1185 PPMonoidElem BestLPP = ComputeLPP(*BestIter, VecLPP); 1186 1187 for(list<PPMonoidElem>::const_iterator it=QBG.myCorners().begin(); it != QBG.myCorners().end(); ++it) 1188 // for(int i=1;i<len(QBG.myCorners());++i) 1189 { 1190 const PPMonoidElem tmp = ComputeLPP(*it, VecLPP); 1191 if (tmp < BestLPP) { BestLPP = tmp; BestIter = it; } 1192 } 1193 const PPMonoidElem PP = *BestIter; 1194 QBG.myCornerPPIntoQB(PP); 1195 const vector<PPMonoidElem>& QB = QBG.myQB(); 1196 FindFactorIndex[PP] = len(QB)-1; 1197 1198 ///JAA vector<SmallFpImpl::value_t> NewRow; 1199 vector<long> exp; exponents(exp, PP); 1200 int BestVar = -1; // just to keep compiler quiet 1201 long BestFactorPP = -1; 1202 for (int var=0; var < n; ++var) 1203 { 1204 if (exp[var] > 0) 1205 { 1206 const int ReducerIndex = FindFactorIndex[PP/indet(PPM(Kx),var)]; 1207 if (ReducerIndex > BestFactorPP) { BestFactorPP = ReducerIndex; BestVar = var; } 1208 } 1209 } 1210 QBFactorPP.push_back(BestFactorPP); 1211 QBFactorVar.push_back(BestVar); 1212 const PPMonoidElem x = indet(PPM(Kx), BestVar); 1213 vector<SmallFpImpl::NonRedValue> NewRow; 1214 for(long j=0; j < NumPts ; ++j) 1215 { 1216 // NewRow.push_back(ModP.myReduce(ConvertTo<long>(eval(NewPolyx, SamplePts[i])))); 1217 NewRow.push_back(ModP.myMul(SamplePts[j][BestVar], RowReducers[BestFactorPP][j])); 1218 } 1219 //POLY RingElem NewPolyx = polyx[BestFactorPP]*indet(Kx,BestVar); 1220 vector<SmallFpImpl::NonRedValue> NewPolyCoeffs2(1+NumPts); 1221 ////NOT NEEDED NewPolyCoeffs[NumPts] = ModP.myReduce(1); 1222 const vector<SmallFpImpl::value>& OldPolyCoeffs = PolyCoeffs[BestFactorPP]; 1223 for(int it=0;it<len(OldPolyCoeffs);++it) 1224 { 1225 if (IsZero(OldPolyCoeffs[it])) continue; 1226 const int j = FindFactorIndex[x*QB[it]]; 1227 NewPolyCoeffs2[j] = OldPolyCoeffs[it]; 1228 } 1229 1230 // just gaussian reduction -- we know that matrix is upper triangular 1231 long IterCount = 0; 1232 for (long i=0; i < NumPts; ++i) 1233 { 1234 const SmallFpImpl::value coordi = ModP.myNormalize(NewRow[i]); 1235 if (IsZero(coordi)) continue; 1236 const SmallFpImpl::value q = ModP.myNegate(ModP.myDiv(coordi, RowReducers[i][i])); 1237 NewRow[i] = zero(SmallFp); // USELESS, but reassuring???? 1238 //POLY NewPolyx += ModP.myExport(q)*polyx[i]; 1239 for (long j=0; j<len(PolyCoeffs[i]);++j) 1240 { 1241 /// NewPolyCoeffs2[j] = add(NewPolyCoeffs2[j], mul(q,PolyCoeffs[i][j])); /// SLLOOOOOOWWWWWW 1242 NewPolyCoeffs2[j] += q * PolyCoeffs[i][j]; 1243 } 1244 for (long j=i+1; j < NumPts; ++j) 1245 { 1246 /// NewRow[j] = add(NewRow[j], mul(q,RowReducers[i][j])); 1247 NewRow[j] += q * RowReducers[i][j]; 1248 } 1249 ++IterCount; 1250 if (IterCount == ModP.myMaxIters()) 1251 { 1252 IterCount = 0; 1253 for (long j=i+1; j < NumPts; ++j) 1254 { 1255 NewRow[j] = ModP.myHalfNormalize(NewRow[j]); 1256 } 1257 for (long j=0; j<=NumPts;++j) 1258 { 1259 NewPolyCoeffs2[j] = ModP.myHalfNormalize(NewPolyCoeffs2[j]); 1260 1261 } 1262 } 1263 } 1264 //POLY polyx.push_back(NewPolyx); 1265 vector<SmallFpImpl::value> NewPolyCoeffs(NumPts+1); for(long j=0;j<=NumPts;++j) NewPolyCoeffs[j]=ModP.myNormalize(NewPolyCoeffs2[j]); 1266 PolyCoeffs.push_back(NewPolyCoeffs); 1267 NewRow[NumPts-1] = ModP.myNormalize(NewRow[NumPts-1]); 1268 // CoCoA_ASSERT(IsZero(ModP.myNormalize(NewRow[NumPts-1]) == 0); 1269 RowReducers.push_back(vector<SmallFpImpl::value>(NumPts)); 1270 } 1271 } 1272 1273 1274 //---------------------------------------------------------------------- 1275 //-- SLICING ALGORITHM 1276 //---------------------------------------------------------------------- 1277 1278 namespace // anonymous for file local fns 1279 { 1280 1281 template <typename T> first(const std::vector<T> & v,long n)1282 std::vector<T> first(const std::vector<T>& v, long n) 1283 { 1284 std::vector<T> w; 1285 if (n>len(v)) CoCoA_THROW_ERROR("vector too short", "first(v, n)"); 1286 for (long i=0; i<n; ++i) w.push_back(v[i]); 1287 return w; 1288 } 1289 1290 template <typename T> last(const std::vector<T> & v,long n)1291 std::vector<T> last(const std::vector<T>& v, long n) 1292 { 1293 std::vector<T> w; 1294 if (n>len(v)) CoCoA_THROW_ERROR("vector too short", "first(v, n)"); 1295 for (long i=len(v)-n; i<len(v); ++i) w.push_back(v[i]); 1296 return w; 1297 } 1298 1299 template <typename T> concat(const std::vector<T> & v1,const std::vector<T> & v2)1300 std::vector<T> concat(const std::vector<T>& v1, const std::vector<T>& v2) 1301 { 1302 std::vector<T> w; 1303 for (long i=0; i<len(v1); ++i) w.push_back(v1[i]); 1304 for (long i=0; i<len(v2); ++i) w.push_back(v2[i]); 1305 return w; 1306 } 1307 1308 } // end of anonymous namespace 1309 1310 //------------------------------ 1311 class ImplicitMill 1312 { 1313 public: 1314 enum FinalCall_t { elim, elim1, elimth, IDWC, IDWCLPP, IDWCLPP2, IDWCOrd2 }; 1315 1316 1317 public: 1318 ImplicitMill(const vector<RingElem>& ParamDescr, 1319 long NumXEndRec, 1320 const string& FinalCall); 1321 1322 // private: 1323 vector<RingElem> myParamDescr; 1324 long myNumXEndRec; 1325 FinalCall_t myFinalCall; 1326 ring myKt; 1327 ring myKx; 1328 ring myKxt; // for elim(1)(t) with/without weights, or homogenizing indet 1329 RingHom myHomT_XT; 1330 RingHom myHomXT_X; 1331 vector<RingElem> myXEndRec; // myKx: x[0],..,x[myNumXEndRec] 1332 mutable vector<long> myMinNumSlices; 1333 mutable vector<long> myXValues; 1334 }; 1335 1336 ImplicitMill(const vector<RingElem> & ParamDescr,long NumXEndRec,const string & FinalCall)1337 ImplicitMill::ImplicitMill(const vector<RingElem>& ParamDescr, 1338 long NumXEndRec, 1339 const string& FinalCall): 1340 myHomT_XT(IdentityHom(owner(ParamDescr[0]))), 1341 myHomXT_X(IdentityHom(owner(ParamDescr[0]))) 1342 { 1343 myNumXEndRec = NumXEndRec; 1344 const ring& Ktorig = owner(ParamDescr[0]); 1345 // myKt = owner(ParamDescr[0]); 1346 long NumT = NumIndets(Ktorig); 1347 long NumX = len(ParamDescr); 1348 if (NumT != NumX-1) CoCoA_THROW_ERROR("NumT != NumX-1", "ImplicitMill"); 1349 myKx = NewPolyRingForImplicit(CoeffRing(Ktorig), "x",NumX); 1350 myKt = NewPolyRingForImplicit(CoeffRing(Ktorig), "t",NumT); 1351 myParamDescr = apply(PolyAlgebraHom(Ktorig,myKt,indets(myKt)), ParamDescr); 1352 1353 myMinNumSlices = vector<long>(NumX, 0); 1354 myXValues = vector<long>(NumX, 0); // x[i] = n 1355 1356 if (FinalCall == "elim") myFinalCall = elim; 1357 else if (FinalCall == "elimth") myFinalCall = elimth; 1358 else if (FinalCall == "elim1") myFinalCall = elim1; 1359 else if (FinalCall == "IDWC") myFinalCall = IDWC; 1360 else if (FinalCall == "IDWCLPP") myFinalCall = IDWCLPP; 1361 else if (FinalCall == "IDWCLPP2") myFinalCall = IDWCLPP2; 1362 else if (FinalCall == "IDWCOrd2") myFinalCall = IDWCOrd2; 1363 else CoCoA_THROW_ERROR("either elim(1)(th), IDWC, IDWCLPP2, or IDWCLPP ","ImplicitMill"); 1364 1365 myXEndRec = first(indets(myKx), myNumXEndRec); 1366 1367 switch (myFinalCall) 1368 { 1369 case ImplicitMill::IDWC: 1370 case ImplicitMill::IDWCLPP: 1371 case ImplicitMill::IDWCLPP2: 1372 case ImplicitMill::IDWCOrd2: 1373 break; 1374 case ImplicitMill::elim: // ring for elim --> 1375 case ImplicitMill::elim1: // ring for elim --> 1376 { 1377 matrix W=NewDenseMat(RingQQ(), 1, NumT+NumXEndRec); 1378 if (myFinalCall==elim1) 1379 for (long i=0; i<NumXEndRec; ++i) SetEntry(W,0,i,1); 1380 else 1381 for (long i=0; i<NumXEndRec; ++i) SetEntry(W,0,i,deg(ParamDescr[i])); 1382 for (long i=0; i<NumT; ++i) SetEntry(W, 0, NumXEndRec+i, 1); 1383 myKxt = NewPolyRingForImplicit(CoeffRing(myKt), 1384 NumXEndRec+NumT, 1385 NewMatrixOrdering(MakeTermOrd(W),1)); 1386 1387 myHomT_XT = PolyAlgebraHom(myKt, myKxt, last(indets(myKxt),NumT)); 1388 myHomXT_X = PolyAlgebraHom(myKxt, myKx, 1389 concat(myXEndRec, 1390 vector<RingElem>(NumT,zero(myKx)))); 1391 } 1392 break; 1393 case ImplicitMill::elimth: // ring for elim --> 1394 { 1395 matrix W=NewDenseMat(RingQQ(), 1, NumT+NumXEndRec+1); 1396 for (long i=0; i<NumXEndRec; ++i) SetEntry(W,0,i,deg(ParamDescr[i])); 1397 for (long i=0; i<=NumT; ++i) SetEntry(W, 0, NumXEndRec+i, 1); 1398 myKxt = NewPolyRingForImplicit(CoeffRing(myKt), 1399 NumXEndRec+NumT+1, 1400 NewMatrixOrdering(MakeTermOrd(W),1)); 1401 1402 myHomT_XT = PolyAlgebraHom(myKt, myKxt, first(last(indets(myKxt),NumT+1),NumT)); 1403 myHomXT_X = PolyAlgebraHom(myKxt, myKx, 1404 concat(myXEndRec, 1405 vector<RingElem>(NumT+1,one(myKx)))); 1406 } 1407 break; 1408 default: 1409 CoCoA_THROW_ERROR(ERR::ShouldNeverGetHere, "ImplicitMill ctor"); 1410 } 1411 } 1412 1413 1414 //------------------------------ 1415 reconstruction(const vector<RingElem> & F,const vector<RingElem> & L)1416 RingElem reconstruction(const vector<RingElem>& F, const vector<RingElem>& L) 1417 { 1418 long d = len(L); 1419 if (d==1) return one(owner(F[0])); 1420 RingElem s(owner(F[0])); 1421 RingElem PrLWithout_i(owner(F[0])); 1422 for (long i=0; i<d; ++i) 1423 { 1424 PrLWithout_i = one(owner(F[0])); 1425 for (long j=0; j<d; ++j) if (i!=j) PrLWithout_i *= L[j]; 1426 s += (PrLWithout_i*F[i]) / NR(PrLWithout_i, vector<RingElem>(1,L[i])); 1427 } 1428 // std::cout << " s = " << s << std::endl; 1429 1430 return s; 1431 } 1432 1433 1434 namespace{ // anonymous 1435 CallElim(const ImplicitMill & IM,long)1436 RingElem CallElim(const ImplicitMill& IM, long /*NumX*/) 1437 { 1438 vector<RingElem> v; 1439 long j=0; 1440 for (; j< IM.myNumXEndRec; ++j) 1441 v.push_back(indet(IM.myKxt,j) - IM.myHomT_XT(IM.myParamDescr[j])); 1442 for (; j<NumIndets(IM.myKx); ++j) 1443 v.push_back(IM.myXValues[j] - IM.myHomT_XT(IM.myParamDescr[j])); 1444 ideal J = ideal(v); // non empty 1445 MakeUnique(J)->myElim(apply(IM.myHomT_XT, indets(IM.myKt))); 1446 // //---------- 1447 // std::cout << "\n---------\n" << v << "\n---------\n" << std::endl; 1448 // RingElem ff = ComputeHypersurface(v, 1449 // LPP(product(apply(IM.myHomT_XT, indets(IM.myKt))))); 1450 // if (len(gens(J)) != 1) 1451 // std::cout << "\n len(gens(J)) != 1" << std::endl; 1452 // if (IM.myHomXT_X(monic(gens(J)[0])) 1453 // != IM.myHomXT_X(monic(ff))) 1454 // std::cout << IM.myHomXT_X(monic(gens(J)[0])) << "\n != \n" 1455 // << IM.myHomXT_X(monic(ff)) << std::endl; 1456 // RingElem ff = CallImplicitDirectWithCondLPP2(IM, NumX); 1457 // if (monic(IM.myHomXT_X(gens(J)[0])) != monic(ff)) 1458 // std::cout << monic(IM.myHomXT_X(gens(J)[0])) << "\n != \n" 1459 // << monic(ff) << std::endl; 1460 //---------- 1461 // return monic(IM.myHomXT_X(gens(J)[0])); 1462 return IM.myHomXT_X(monic(gens(J)[0])); 1463 } 1464 1465 CallElimTH(const ImplicitMill & IM,long)1466 RingElem CallElimTH(const ImplicitMill& IM, long /*NumX*/) 1467 { 1468 vector<RingElem> v; 1469 long j=0; 1470 RingElem h = indets(IM.myKxt)[NumIndets(IM.myKxt)-1]; 1471 for (; j< IM.myNumXEndRec; ++j) 1472 v.push_back(homog(indet(IM.myKxt,j) - IM.myHomT_XT(IM.myParamDescr[j]), h)); 1473 for (; j<NumIndets(IM.myKx); ++j) 1474 v.push_back(homog(IM.myXValues[j] - IM.myHomT_XT(IM.myParamDescr[j]), h)); 1475 // std::cout <<" v =" << v << std::endl; 1476 RingElem ff = ComputeElimFirst(v, 1477 LPP(product(apply(IM.myHomT_XT, indets(IM.myKt))))); 1478 return IM.myHomXT_X(monic(ff)); 1479 } 1480 1481 CallImplicitDirectWithCond(const ImplicitMill & IM,long NumX)1482 RingElem CallImplicitDirectWithCond(const ImplicitMill& IM, long NumX) 1483 { 1484 vector<RingElem> cond; 1485 for (long i=NumX; i<NumIndets(IM.myKx); ++i) 1486 cond.push_back(IM.myParamDescr[i] - IM.myXValues[i]); 1487 RingElem f = monic(ImplicitDirectWithCond(first(IM.myParamDescr,NumX),cond)); 1488 return PolyAlgebraHom(owner(f), IM.myKx, IM.myXEndRec)(f); 1489 } 1490 1491 CallImplicitDirectWithCondLPP(const ImplicitMill & IM,long NumX)1492 RingElem CallImplicitDirectWithCondLPP(const ImplicitMill& IM, long NumX) 1493 { 1494 vector<RingElem> cond; 1495 for (long i=NumX; i<NumIndets(IM.myKx); ++i) 1496 cond.push_back(IM.myParamDescr[i] - IM.myXValues[i]); 1497 RingElem f = monic(ImplicitDirectWithCondLPP(first(IM.myParamDescr,NumX), cond)); 1498 return PolyAlgebraHom(owner(f), IM.myKx, IM.myXEndRec)(f); 1499 } 1500 1501 CallImplicitDirectWithCondLPP2(const ImplicitMill & IM,long NumX)1502 RingElem CallImplicitDirectWithCondLPP2(const ImplicitMill& IM, long NumX) 1503 { 1504 vector<RingElem> cond; 1505 for (long i=NumX; i<NumIndets(IM.myKx); ++i) 1506 cond.push_back(IM.myParamDescr[i] - IM.myXValues[i]); 1507 RingElem f = monic(ImplicitDirectWithCondLPP2(first(IM.myParamDescr,NumX), cond)); 1508 return PolyAlgebraHom(owner(f), IM.myKx, IM.myXEndRec)(f); 1509 } 1510 1511 CallImplicitDirectWithCondOrd2(const ImplicitMill & IM,long NumX)1512 RingElem CallImplicitDirectWithCondOrd2(const ImplicitMill& IM, long NumX) 1513 { 1514 vector<RingElem> cond; 1515 for (long i=NumX; i<NumIndets(IM.myKx); ++i) 1516 cond.push_back(IM.myParamDescr[i] - IM.myXValues[i]); 1517 RingElem f = monic(ImplicitDirectWithCondOrd2(first(IM.myParamDescr,NumX), cond)); 1518 return PolyAlgebraHom(owner(f), IM.myKx, IM.myXEndRec)(f); 1519 } 1520 1521 SliceCoreRec(const ImplicitMill & IM,long NumX)1522 RingElem SliceCoreRec(const ImplicitMill& IM, long NumX) 1523 { // X = [x[0], .., x[NumX-1], a[NumX], .., a[N-1]] 1524 VerboseLog VERBOSE("SliceCoreRec"); 1525 ring KX = IM.myKx; 1526 ring KT = IM.myKt; 1527 if (NumX <= IM.myNumXEndRec) 1528 switch (IM.myFinalCall) 1529 { 1530 case ImplicitMill::elim: return CallElim(IM, NumX); 1531 case ImplicitMill::elim1: return CallElim(IM, NumX); 1532 case ImplicitMill::elimth: return CallElimTH(IM, NumX); 1533 case ImplicitMill::IDWC: return CallImplicitDirectWithCond(IM, NumX); 1534 case ImplicitMill::IDWCLPP: return CallImplicitDirectWithCondLPP(IM,NumX); 1535 case ImplicitMill::IDWCLPP2: return CallImplicitDirectWithCondLPP2(IM, NumX); 1536 case ImplicitMill::IDWCOrd2: return CallImplicitDirectWithCondOrd2(IM, NumX); 1537 default: CoCoA_THROW_ERROR(ERR::ShouldNeverGetHere, "SliceCoreRec"); 1538 } 1539 RingElem x_n = indet(KX, NumX-1); 1540 vector<RingElem> Slices; 1541 vector<RingElem> ResultSlices; 1542 for (long i=1; i<=50; ++i) 1543 { 1544 // std::cout << " " << i << "(" << NumX << ")" << std::flush; 1545 Slices.push_back(x_n - i); 1546 IM.myXValues[NumX-1] = i; // [0,.., 0, *i*, a[NumX], .., a[N-1]] 1547 ResultSlices.push_back(SliceCoreRec(IM, NumX-1)); 1548 if (i < IM.myMinNumSlices[NumX-1]) continue; 1549 RingElem candidate = monic(reconstruction(ResultSlices, Slices)); 1550 if (IM.myMinNumSlices[NumX-1]==0 && 1551 deg(candidate,NumX-1) == len(Slices)-1) 1552 { 1553 // std::cout << ":LowD" << std::flush; 1554 continue; 1555 } 1556 vector<RingElem> v, img(NumIndets(KX), zero(KT)); 1557 for (long j=0; j<NumX; ++j) img[j] = IM.myParamDescr[j]; 1558 for (long j=NumX; j<len(IM.myXValues); ++j) img[j] = IM.myXValues[j]; 1559 for (long j=NumX; j<len(IM.myXValues); ++j) v.push_back(IM.myParamDescr[j]-IM.myXValues[j]); 1560 RingElem sbst = NF(PolyAlgebraHom(KX,KT,img)(candidate), ideal(KT,v)); 1561 if (IsZero(sbst)) 1562 { 1563 if (IM.myMinNumSlices[NumX-1]==0) IM.myMinNumSlices[NumX-1] = i-1; 1564 // std::cout << "\n--> MinNumSlices = " << MinNumSlices; 1565 // std::cout << " time: " << CpuTime()-T << std::endl; 1566 return candidate; 1567 } 1568 VERBOSE(50) << " !!!wrong candidate!!! one more slice..."; 1569 } 1570 CoCoA_THROW_ERROR("More than 50 slices: probably bad slicer!","SliceCoreRec"); 1571 return zero(KX); 1572 } 1573 } 1574 1575 SliceCore(const vector<RingElem> & ParamDescr,long RecDepth,const string & FinalCall)1576 RingElem SliceCore(const vector<RingElem>& ParamDescr, 1577 long RecDepth, 1578 const string& FinalCall) 1579 { 1580 ImplicitMill IM(ParamDescr, len(ParamDescr)-RecDepth, FinalCall); 1581 return SliceCoreRec(IM, len(ParamDescr)); 1582 } 1583 1584 1585 namespace // anonymous for file local defns 1586 { SliceCoreModp(SparsePolyRing QQx,const vector<RingElem> & ParamDescr,long RecDepth,const string & FinalCall,long p)1587 RingElem SliceCoreModp(SparsePolyRing QQx, 1588 const vector<RingElem>& ParamDescr, 1589 long RecDepth, 1590 const string& FinalCall, long p) 1591 { 1592 const SparsePolyRing& QQt(owner(ParamDescr[0])); 1593 const long s = NumIndets(QQt); 1594 const SparsePolyRing Fpt(NewPolyRing_DMPII(NewZZmod(p), SymbolRange("t",1,s))); 1595 const RingHom phi = PolyRingHom(QQt,Fpt, QQEmbeddingHom(Fpt), indets(Fpt)); 1596 const RingElem f = SliceCore(apply(phi,ParamDescr), RecDepth, FinalCall); 1597 const SparsePolyRing& Fpx(owner(f)); 1598 const PPMonoidHom psi(GeneralHom(PPM(Fpx), first(indets(PPM(QQx)),NumIndets(Fpx)))); 1599 RingElem ans(QQx); 1600 for (SparsePolyIter it=BeginIter(f) ; !IsEnded(it) ; ++it ) 1601 ans += monomial(QQx, ConvertTo<BigInt>(coeff(it)), psi(PP(it))); 1602 return ans; 1603 } 1604 1605 IsZeroEvalHorner(ConstRefRingElem f,const vector<RingElem> & ParamDescr)1606 bool IsZeroEvalHorner(ConstRefRingElem f, const vector<RingElem>& ParamDescr) 1607 { 1608 VerboseLog VERBOSE("IsZeroEvalHorner"); 1609 const SparsePolyRing& P1 = owner(f); 1610 const SparsePolyRing& P2 = owner(ParamDescr[0]); 1611 SparsePolyRing P = NewPolyRing(CoeffRing(P1), 1612 NewSymbols(NumIndets(P1)+NumIndets(P2))); 1613 RingHom phi1 = PolyAlgebraHom(P1, P, first(indets(P), NumIndets(P1))); 1614 RingHom phi2 = PolyAlgebraHom(P2, P, last(indets(P), NumIndets(P2))); 1615 const vector<RingElem>& ParamDescr_P = apply(phi2, ParamDescr); 1616 RingElem eval_f = phi1(f); 1617 for (long i=0; i<NumIndets(P1); ++i) 1618 { 1619 VERBOSE(50) << "horner-"<<i << std::endl; 1620 std::vector<RingElem> c = CoeffVecWRT(eval_f, indet(P,i)); 1621 eval_f = zero(P); 1622 for (long d=len(c)-1; d>=0; --d) 1623 eval_f = eval_f*ParamDescr_P[i] + c[d]; 1624 } 1625 return IsZero(eval_f); 1626 } 1627 1628 1629 1630 } // namespace // anonymous 1631 1632 1633 1634 SliceCoreQQ(const vector<RingElem> & ParamDescr,long RecDepth,const string & FinalCall)1635 RingElem SliceCoreQQ(const vector<RingElem>& ParamDescr, 1636 long RecDepth, 1637 const string& FinalCall) 1638 { 1639 VerboseLog VERBOSE("SliceCoreQQ"); 1640 SparsePolyRing QQx(RingQQt(len(ParamDescr))); 1641 const SparsePolyRing& QQt = owner(ParamDescr[0]); 1642 RingHom eval = PolyAlgebraHom(QQx, QQt, ParamDescr); 1643 long PrimeCount = 0; 1644 long p = 46349; 1645 RingElem fCRT(QQx); 1646 BigInt fModulus(p); 1647 // RingElem ParDDenom = CommonDenom(ParamDescr[0]); 1648 // for (long i=1; i<len(ParamDescr); ++i) 1649 // ParDDenom = lcm(ParDDenom, CommonDenom(ParamDescr[i])); 1650 RingElem ParDDenom = CommonDenom(ParamDescr); 1651 while (IsZero(fCRT)) 1652 { 1653 CheckForInterrupt("SliceCoreQQ"); 1654 p = PrevPrime(p); 1655 VERBOSE(20) << ++PrimeCount << ": prime is " << p << std::endl; 1656 if (IsDivisible(ParDDenom, p)) 1657 { 1658 VERBOSE(80) << " UGLY PRIME: going to another prime" << std::endl; 1659 continue; 1660 } 1661 fCRT = SliceCoreModp(QQx, ParamDescr, RecDepth, FinalCall, p); 1662 fModulus = p; 1663 try 1664 { 1665 VERBOSE(80) << "---- RatReconstruct" << std::endl; 1666 const RingElem f = RatReconstructPoly(fCRT, fModulus); 1667 VERBOSE(80) << "---- IsZero(eval(f))" << std::endl; 1668 if (IsZero(eval(f))) return f; 1669 } 1670 catch (const CoCoA::ErrorInfo& err) 1671 { 1672 if (err != ERR::CannotReconstruct) throw; 1673 } 1674 } 1675 while (true) 1676 { 1677 CheckForInterrupt("SliceCoreQQ"); 1678 p = PrevPrime(p); 1679 VERBOSE(20) << ++PrimeCount << ": prime is " << p << std::endl; 1680 if (IsDivisible(ParDDenom, p)) 1681 { 1682 VERBOSE(80) << " UGLY PRIME: going to another prime" << std::endl; 1683 continue; 1684 } 1685 const RingElem fp = SliceCoreModp(QQx,ParamDescr, RecDepth, FinalCall, p); 1686 CRTPoly(fCRT, fModulus, fCRT, fModulus, fp, BigInt(p)); 1687 try 1688 { 1689 VERBOSE(80) << "---- RatReconstruct" << std::endl; 1690 const RingElem f = RatReconstructPoly(fCRT, fModulus); 1691 double t=CpuTime(); 1692 bool b = IsZero(eval(f)); 1693 VERBOSE(90) << "---- IsZero time: " << CpuTime()-t << std::endl; 1694 // t = CpuTime(); 1695 // b = IsZeroEvalHorner(f, ParamDescr); 1696 // std::cout << std::boolalpha; // so that bools print out as true/false 1697 // std::cout << "---- IsZeroEvalHorner " << b 1698 // << " time: " << CpuTime()-t << std::endl; 1699 if (b) return f; 1700 } 1701 catch (const CoCoA::ErrorInfo& err) 1702 { 1703 if (err != ERR::CannotReconstruct) throw; 1704 } 1705 } 1706 } 1707 1708 1709 1710 1711 } // end of namespace CoCoA 1712 1713 1714 // RCS header/log in the next few lines 1715 // $Header: /Volumes/Home_1/cocoa/cvs-repository/CoCoALib-0.99/src/AlgebraicCore/SparsePolyOps-implicit.C,v 1.5 2020/06/17 15:49:28 abbott Exp $ 1716 // $Log: SparsePolyOps-implicit.C,v $ 1717 // Revision 1.5 2020/06/17 15:49:28 abbott 1718 // Summary: Changed CoCoA_ERROR into CoCoA_THROW_ERROR 1719 // 1720 // Revision 1.4 2019/12/21 16:43:43 abbott 1721 // Summary: Cleaned (most ring& instead of ring); removed some cruft 1722 // 1723 // Revision 1.3 2018/05/18 16:38:51 bigatti 1724 // -- added include SparsePolyOps-RingElem.H 1725 // 1726 // Revision 1.2 2018/05/17 15:50:13 bigatti 1727 // -- renamed MatrixOperations --> MatrixOps 1728 // -- sorted includes 1729 // -- renamed VectorOperations --> VectorOps 1730 // -- added include SparsePolyIter 1731 // 1732 // Revision 1.1 2018/04/06 15:16:11 bigatti 1733 // -- renamed TmpImplicit.C 1734 // 1735 // Revision 1.50 2018/03/29 09:36:40 bigatti 1736 // -- added member functions myTestIsRadical, myTestIsPrimary and flags 1737 // 1738 // Revision 1.49 2018/02/27 17:30:23 abbott 1739 // Summary: Renamed NumTheory_prime to NumTheory-prime; changed includes 1740 // 1741 // Revision 1.48 2018/02/27 10:55:37 abbott 1742 // Summary: Added include NumTheory_prime 1743 // 1744 // Revision 1.47 2017/09/06 11:56:30 abbott 1745 // Summary: Changed ERR::SERIOUS into ERR::ShouldNeverGetHere 1746 // 1747 // Revision 1.46 2017/07/07 09:52:59 bigatti 1748 // -- changed NextPrime --> PrevPrime 1749 // 1750 // Revision 1.45 2017/06/28 12:48:18 bigatti 1751 // -- all cout now are VERBOSE(..) 1752 // 1753 // Revision 1.44 2017/05/15 10:24:37 bigatti 1754 // -- detecting UGLY PRIME using divisibility (instead of try/catch) 1755 // 1756 // Revision 1.43 2017/05/11 10:58:10 bigatti 1757 // -- now catching ugly primes (as for MinPolyModular) 1758 // 1759 // Revision 1.42 2017/05/11 08:46:46 bigatti 1760 // -- added error ERR::CannotReconstruct 1761 // -- cleaned up code accordingly 1762 // 1763 // Revision 1.41 2017/04/18 09:52:51 bigatti 1764 // -- added verbosity 1765 // 1766 // Revision 1.40 2017/04/10 15:51:19 bigatti 1767 // -- debug printing set to false (to be substituted by verbose) 1768 // 1769 // Revision 1.39 2017/03/08 15:35:09 bigatti 1770 // -- added verbosity (90) 1771 // 1772 // Revision 1.38 2016/11/03 12:25:25 abbott 1773 // Summary: Changed IsRadical (for PPMonoidElem) into IsSqFree 1774 // 1775 // Revision 1.37 2016/07/19 11:35:57 bigatti 1776 // -- static bool for printing algorithm name 1777 // 1778 // Revision 1.36 2016/06/24 14:27:41 bigatti 1779 // -- renamed CRT_poly --> CRTPoly 1780 // 1781 // Revision 1.35 2016/04/26 14:33:00 abbott 1782 // Summary: commented out some debug print commands 1783 // 1784 // Revision 1.34 2016/04/22 16:06:59 bigatti 1785 // -- added first implementation of IsZeroEvalHorner 1786 // 1787 // Revision 1.33 2016/03/09 10:17:10 abbott 1788 // Summary: Commented out some logging print commands 1789 // 1790 // Revision 1.32 2016/02/17 10:30:37 bigatti 1791 // -- added a return to keep compiler quiet (NYI) 1792 // 1793 // Revision 1.31 2016/01/26 13:51:00 bigatti 1794 // -- added ImplicitDirectOrd2 1795 // 1796 // Revision 1.30 2015/12/08 13:56:09 abbott 1797 // Summary: Updated Mario's code! Very many changes! 1798 // 1799 // Revision 1.29 2015/11/04 12:14:58 abbott 1800 // Summary: Several consequential changes (after revising SmallFpImpl) 1801 // 1802 // Revision 1.28 2015/10/01 10:15:05 bigatti 1803 // -- moved CRT_poly, RatReconstructPoly --> SparsePolyRing.C 1804 // 1805 // Revision 1.27 2015/06/26 14:59:51 abbott 1806 // Summary: Revised after changing interface to ErrorInfo (inherited from CoCoA::exception) 1807 // Author: JAA 1808 // 1809 // Revision 1.26 2015/06/11 16:57:24 bigatti 1810 // -- using new functions monomial(ring, pp) and monomial(ring, expv) 1811 // 1812 // Revision 1.25 2015/05/04 13:17:16 bigatti 1813 // -- just new names: FindReducerIndex, NewPolyx, polyx,... 1814 // 1815 // Revision 1.24 2015/04/24 15:34:02 bigatti 1816 // -- added ImplicitDirectLPP2Homog, ImplicitDirectLPP2HomogBis 1817 // -- changed ring names into Kx and Kt 1818 // 1819 // Revision 1.23 2015/03/09 13:13:27 bigatti 1820 // -- added SliceCoreQQ (and all needed functions in anonymous namespace) 1821 // 1822 // Revision 1.22 2015/03/05 16:49:32 bigatti 1823 // -- changed back to creating a new ring (faster) in ImplicitDirectWithCondLPP2 1824 // 1825 // Revision 1.21 2015/03/04 11:18:19 bigatti 1826 // -- added elim1 (no weights), elimth (truncated + homogeneous) 1827 // 1828 // Revision 1.20 2015/02/12 15:59:03 bigatti 1829 // -- in ImplicitDirectWithCondLPP2: ChooseIndet(..poly) --> ChooseIndet(..reducer) 1830 // 1831 // Revision 1.19 2015/02/12 09:38:11 bigatti 1832 // -- changed all ring names into Kt and Kx 1833 // 1834 // Revision 1.18 2015/02/06 13:48:03 bigatti 1835 // -- allowing direct computation in QQ[...] using DMPI (instead of just DMPII) 1836 // 1837 // Revision 1.17 2014/12/18 15:57:45 bigatti 1838 // -- renamed myRingX and T into myKx myKt in ImplicitMill 1839 // -- myKt is now DMPII 1840 // -- fixed usage of myKt in ImplicitDirectWithCondLPP2 1841 // 1842 // Revision 1.16 2014/12/18 14:20:40 bigatti 1843 // -- renamed rings into Kx, Kt in ImplicitDirectWithCondLPP2 1844 // -- no new ring Kt (using the input one) 1845 // -- indexing x from 1 in SliceCore 1846 // 1847 // Revision 1.15 2014/12/11 15:36:34 bigatti 1848 // -- SliceCore output in "x" (instead of anonymous symbols) 1849 // 1850 // Revision 1.14 2014/12/10 11:55:21 bigatti 1851 // -- new base ring for SliceCore 1852 // 1853 // Revision 1.13 2014/12/04 08:55:43 bigatti 1854 // -- minor optimization to reduce (use myAddMul instead?) 1855 // 1856 // Revision 1.12 2014/12/03 15:55:31 bigatti 1857 // -- trivial change in ChooseIndet 1858 // 1859 // Revision 1.11 2014/12/03 13:55:02 bigatti 1860 // -- added ImplicitDirectWithCondLPP2 1861 // 1862 // Revision 1.10 2014/11/28 15:22:55 bigatti 1863 // -- added arg to SliceCore for algorithm of final calls 1864 // -- now SliceCore uses ImplicitMill 1865 // -- minor fix in ImplicitDirectWithCondLPP 1866 // (now calling ImplicitDirectLPP is cond is empty) 1867 // 1868 // Revision 1.9 2014/11/27 11:32:25 abbott 1869 // Summary: Added ImplicitDirectWithCondLPP 1870 // Author: JAA 1871 // 1872 // Revision 1.8 2014/11/27 09:45:37 bigatti 1873 // -- SliceCoreRec: choose CallElim or CallImplicitDirect (at compiletime) 1874 // 1875 // Revision 1.7 2014/11/17 10:34:53 abbott 1876 // Summary: Added ImplicitDirectLPP2 and ImplicitDirectWithCond 1877 // Author: JAA 1878 // 1879 // Revision 1.6 2014/11/13 16:39:52 bigatti 1880 // -- SliceCore: creating NewPolyRing_DMPII, removed printouts 1881 // 1882 // Revision 1.5 2014/11/13 16:30:40 abbott 1883 // Summary: Now ImplicitDirect uses DMPII 1884 // Author: JAA 1885 // 1886 // Revision 1.4 2014/11/13 15:45:05 abbott 1887 // Summary: update to ImplicitDirectLPP 1888 // Author: JAA 1889 // 1890 // Revision 1.3 2014/11/11 09:19:59 bigatti 1891 // -- added SliceCore 1892 // 1893 // Revision 1.2 2014/11/04 18:20:29 abbott 1894 // Summary: Added anon namespaces 1895 // Author: JAA 1896 // 1897 // Revision 1.1 2014/11/04 18:10:13 abbott 1898 // Summary: new tmp code for implicitization 1899 // Author: JAA 1900 // 1901 // 1902