1 // Copyright (c) 2005-2010,2014 John Abbott 2 3 // This file is part of the source of CoCoALib, the CoCoA Library. 4 5 // CoCoALib is free software: you can redistribute it and/or modify 6 // it under the terms of the GNU General Public License as published by 7 // the Free Software Foundation, either version 3 of the License, or 8 // (at your option) any later version. 9 10 // CoCoALib is distributed in the hope that it will be useful, 11 // but WITHOUT ANY WARRANTY; without even the implied warranty of 12 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 13 // GNU General Public License for more details. 14 15 // You should have received a copy of the GNU General Public License 16 // along with CoCoALib. If not, see <http://www.gnu.org/licenses/>. 17 18 19 // Source code for RingWeylImpl 20 21 #include "CoCoA/RingWeyl.H" 22 23 #include "CoCoA/BigIntOps.H" 24 #include "CoCoA/DistrMPolyClean.H" 25 #include "CoCoA/MatrixForOrdering.H" 26 #include "CoCoA/MatrixOps.H" 27 #include "CoCoA/OpenMath.H" 28 #include "CoCoA/PPMonoid.H" 29 #include "CoCoA/PPMonoidEv.H" 30 #include "CoCoA/PPOrdering.H" 31 #include "CoCoA/RingDistrMPolyClean.H" 32 #include "CoCoA/RingHom.H" 33 #include "CoCoA/TmpGReductor.H" 34 #include "CoCoA/TmpUniversalInvolutiveBasisContainer.H" 35 #include "CoCoA/assert.H" 36 #include "CoCoA/convert.H" 37 #include "CoCoA/ideal.H" 38 #include "CoCoA/matrix.H" 39 #include "CoCoA/symbol.H" 40 #include "CoCoA/utils.H" 41 42 #include <memory> 43 using std::unique_ptr; 44 #include <iostream> 45 using std::ostream; 46 using std::endl; // just for debugging 47 //#include <vector> 48 using std::vector; 49 50 51 52 namespace CoCoA 53 { 54 55 class RingWeylImpl: public SparsePolyRingBase //???? should be NCRingBase non-commutative!!! 56 { 57 58 public: 59 // RingWeylImpl(const ring& CoeffRing, long NumIndets, const std::vector<long>& ElimIndets); 60 RingWeylImpl(const ring& CoeffRing, const std::vector<symbol>& names, const std::vector<long>& ElimIndets); 61 ~RingWeylImpl(); 62 63 private: // data members of RingWeyl 64 const SparsePolyRing myReprRing; 65 long myNumTrueIndetsValue; 66 std::unique_ptr<RingElem> myZeroPtr; ///< Every ring stores its own zero. 67 std::unique_ptr<RingElem> myOnePtr; ///< Every ring stores its own one. 68 std::vector<RingElem> myIndetVector; 69 std::vector<RingElem> myDerivationVector; 70 71 public: 72 //---------------------------------------------------------------------- 73 // Functions which every ring must implement: 74 //---------------------------------------------------------------------- 75 virtual bool IamCommutative() const; 76 virtual bool3 IamIntegralDomain3(bool) const; 77 virtual bool IamTrueGCDDomain() const; 78 virtual ConstRefRingElem myZero() const; 79 virtual ConstRefRingElem myOne() const; 80 using RingBase::myNew; // disable warnings of overloading 81 using RingBase::myAssign; // disable warnings of overloading 82 using PolyRingBase::myIndets; // disable warnings of overloading 83 virtual RingElemRawPtr myNew() const; 84 virtual RingElemRawPtr myNew(const MachineInt& n) const; 85 virtual RingElemRawPtr myNew(const BigInt& N) const; 86 virtual RingElemRawPtr myNew(ConstRawPtr rawcopy) const; 87 virtual void myDelete(RawPtr rawx) const; // destroys x (incl all resources) 88 virtual void mySwap(RawPtr rawx, RawPtr rawy) const; // swap(x, y) 89 virtual void myAssign(RawPtr rawlhs, ConstRawPtr rawx) const; // lhs = x 90 virtual void myAssign(RawPtr rawlhs, const MachineInt& n) const; // lhs = n 91 virtual void myAssign(RawPtr rawlhs, const BigInt& N) const; // lhs = N 92 virtual void myAssignZero(RawPtr rawlhs) const; // lhs = 0 93 virtual void myRecvTwinFloat(RawPtr rawlhs, ConstRawPtr rawx) const; 94 virtual void myNegate(RawPtr rawlhs, ConstRawPtr rawx) const; // lhs = -x 95 virtual void myAdd(RawPtr rawlhs, ConstRawPtr rawx, ConstRawPtr rawy) const; // lhs = x+y 96 virtual void mySub(RawPtr rawlhs, ConstRawPtr rawx, ConstRawPtr rawy) const; // lhs = x-y 97 // virtual void myMul(RawPtr rawlhs, ConstRawPtr rawx, ConstRawPtr rawy) const; // lhs = x*y 98 // virtual void myDiv(RawPtr rawlhs, ConstRawPtr rawx, ConstRawPtr rawy) const; // lhs = x/y 99 virtual bool myIsDivisible(RawPtr rawlhs, ConstRawPtr rawx, ConstRawPtr rawy) const;// lhs = x/y, if divisible 100 virtual bool myIsInvertible(ConstRawPtr rawx) const; // true iff x is invertible 101 // virtual bool myIsPrintAtom(ConstRawPtr rawx) const; 102 virtual void myGcd(RawPtr rawlhs, ConstRawPtr rawx, ConstRawPtr rawy) const; // not a TrueGCDDomain 103 virtual void myPowerSmallExp(RawPtr rawlhs, ConstRawPtr rawx, long n) const;// lhs = x^n, n>1, x not -1,0,1 104 virtual void mySymbols(std::vector<symbol>& SymList) const; // appends ring's symbols to SymList 105 virtual void myOutput(std::ostream& out, ConstRawPtr rawx) const; // out << x 106 using PolyRingBase::myOutputSelf; // disable warnings of overloading 107 // virtual void myOutputSelf(OpenMathOutput& OMOut) const; // OMOut << R myImplDetails()108 virtual std::string myImplDetails() const {return "RingWeyl";} 109 virtual void myOutput(OpenMathOutput& OMOut, ConstRawPtr rawx) const; // OMOut << x 110 virtual bool myIsZero(ConstRawPtr rawx) const; // x == 0 111 virtual bool myIsOne(ConstRawPtr rawx) const; // x == 1 112 virtual bool myIsMinusOne(ConstRawPtr rawx) const; // x == -1 113 // virtual bool myIsZeroAddMul: use default definition 114 virtual bool myIsEqual(ConstRawPtr rawx, ConstRawPtr rawy) const; // x == y 115 116 virtual ideal myIdealCtor(const std::vector<RingElem>& gens) const; 117 118 virtual RingHom myCompose(const RingHom& phi, const RingHom& theta) const; // phi(theta(...)) 119 120 121 /*---------------------------------------------------------------------- 122 Member functions every PolyRing must have 123 in addition to those of RingBase. 124 ----------------------------------------------------------------------*/ 125 virtual long myNumIndets() const; 126 virtual const ring& myCoeffRing() const; 127 virtual const std::vector<RingElem>& myIndets() const; 128 virtual void myIndetPower(RawPtr rawf, long var, long exp) const; 129 130 ///@name Simple functions on polynomials 131 //@{ 132 virtual long myNumTerms(ConstRawPtr rawf) const; 133 virtual bool myIsConstant(ConstRawPtr rawf) const; 134 virtual bool myIsIndet(long& index, ConstRawPtr rawf) const; 135 virtual bool myIsMonomial(ConstRawPtr rawf) const; 136 virtual long myStdDeg(ConstRawPtr rawf) const; 137 virtual long myDeg(ConstRawPtr rawf, long var) const; 138 virtual RingElemAlias myLC(ConstRawPtr rawf) const; 139 virtual void myContent(RawPtr rawcontent, ConstRawPtr rawf) const; 140 virtual void myRemoveBigContent(RawPtr rawf) const; 141 virtual void myMulByCoeff(RawPtr rawf, ConstRawPtr rawc) const; ///< WEAK EXCEPTION GUARANTEE 142 virtual bool myDivByCoeff(RawPtr rawf, ConstRawPtr rawc) const; ///< WEAK EXCEPTION GUARANTEE 143 virtual void myDeriv(RawPtr rawlhs, ConstRawPtr rawf, ConstRawPtr rawx) const; ///< lhs = deriv(f, x) 144 //@} 145 146 // friend const RingElem& indet(const RingWeyl& RW, long var); 147 // friend const RingElem& derivation(const RingWeyl& RW, long var); 148 149 virtual RingHom myHomCtor(const ring& codomain, const RingHom& CoeffHom, const std::vector<RingElem>& IndetImages) const; 150 151 //---------------------------------------------------------------------- 152 // Functions which every SparsePolyRing must implement: 153 //---------------------------------------------------------------------- 154 155 virtual const PPMonoid& myPPM() const; 156 157 ///@name Functions for creating/building polynomials 158 //@{ 159 virtual RingElem myMonomial(ConstRawPtr rawc, PPMonoidElemConstRawPtr rawpp) const; 160 virtual SparsePolyIter myBeginIter(ConstRawPtr rawf) const; 161 virtual SparsePolyIter myEndIter(ConstRawPtr rawf) const; 162 163 virtual void myPushFront(RawPtr rawf, ConstRawPtr rawc, const std::vector<long>& expv) const; 164 virtual void myPushBack(RawPtr rawf, ConstRawPtr rawc, const std::vector<long>& expv) const; 165 virtual void myPushFront(RawPtr rawf, ConstRawPtr rawc, PPMonoidElemConstRawPtr rawpp) const; 166 virtual void myPushBack(RawPtr rawf, ConstRawPtr rawc, PPMonoidElemConstRawPtr rawpp) const; 167 //@} 168 169 ///@name Special functions on polynomials needed for implementing Buchberger's Algorithm 170 //@{ 171 // virtual void myWDeg(degree& d, ConstRawPtr rawf) const; 172 // virtual int myCmpWDeg(ConstRawPtr rawf, ConstRawPtr rawg) const; 173 virtual ConstRefPPMonoidElem myLPP(ConstRawPtr rawf) const; 174 virtual void myMulByPP(RawPtr rawf, PPMonoidElemConstRawPtr rawpp) const; 175 virtual bool myIsZeroAddLCs(RawPtr rawf, RawPtr rawg) const; ///< f+=LM(g); g-=LM(g); assumes LPP(f)==LPP(g); returns LC(f)+LC(g)==0 176 virtual void myMoveLMToFront(RawPtr rawf, RawPtr rawg) const; 177 virtual void myMoveLMToBack(RawPtr rawf, RawPtr rawg) const; 178 virtual void myDeleteLM(RawPtr rawf) const; // ????? right interface 179 virtual void myDivLM(RawPtr rawlhs, ConstRawPtr rawf, ConstRawPtr rawg) const; ///< lhs=div(LM(f),LM(g)); assumes f!=0,g!=0 180 virtual int myCmpLPP(ConstRawPtr rawf, ConstRawPtr rawg) const; ///< cmp(LPP(f),LPP(g)); assumes f!=0,g!=0 181 virtual void myAddClear(RawPtr rawf, RawPtr rawg) const; ///< f+=g; g=0; 182 virtual void myAppendClear(RawPtr rawf, RawPtr rawg) const; ///< f+=g; g=0; appends g to f with no checks 183 184 virtual void myAddMulLM(RawPtr rawf, ConstRawPtr rawh, ConstRawPtr rawg) const; ///< f += LM(h)*g 185 virtual void myAddMulLM(RawPtr rawf, ConstRawPtr rawh, ConstRawPtr rawg, SkipLMFlag) const; ///< f += LM(h)*g 186 virtual void myReductionStep(RawPtr rawf, ConstRawPtr rawg) const; 187 // ??? aggiungere coefficiente 188 virtual void myReductionStepGCD(RawPtr rawf, ConstRawPtr rawg, RingElem& FScale) const; 189 // should it all be in ReductionStep ??? ANNA 190 //@} 191 192 193 private: // HomImpl class: overwrite some implementations of SparsePolyRing functions 194 class HomImpl: public SparsePolyRingBase::HomImpl 195 { 196 public: 197 HomImpl(const SparsePolyRing& domain, const ring& codomain, const RingHom& CoeffHom, const std::vector<RingElem>& IndetImages); 198 }; 199 200 private: // Ideal class: overwrite some implementations of SparsePolyRing functions 201 class IdealImpl: public SparsePolyRingBase::IdealImpl 202 { 203 public: 204 IdealImpl(const SparsePolyRing& P, const std::vector<RingElem>& gens); 205 // default copy ctor is OK 206 virtual void myReduceMod(RingElemRawPtr rawr) const; // r elem of R, where I is ideal in R 207 virtual bool IhaveElem(RingElemConstRawPtr rawr) const; 208 virtual void myIntersect(const ideal&); 209 virtual void myColon(const ideal&); 210 virtual bool myDivMod(RingElemRawPtr rawlhs, RingElemConstRawPtr rawnum, RingElemConstRawPtr rawden) const; // result is true iff quot exists & is unique; if true set lhs = num/den modulo I 211 212 virtual void myTestIsMaximal() const; 213 virtual void myTestIsPrimary() const; 214 virtual void myTestIsPrime() const; 215 virtual void myTestIsRadical() const; 216 virtual const std::vector<RingElem>& myGBasis(const CpuTimeLimit&) const; 217 }; 218 }; 219 220 //---------------------------------------------------------------------- 221 WANAMES(long NumIndets)222 std::vector<symbol> WANAMES(long NumIndets) 223 { 224 vector<symbol> ans = SymbolRange("x", 0, NumIndets-1); 225 vector<symbol> d = SymbolRange("d", 0, NumIndets-1); 226 ans.insert(ans.end(), d.begin(), d.end()); 227 return ans; 228 } 229 WANAMES(const std::vector<symbol> & names)230 std::vector<symbol> WANAMES(const std::vector<symbol>& names) 231 { 232 vector<symbol> ans=names; 233 const long NumNames = len(names); 234 ans.reserve(2*NumNames); 235 for ( long i=0 ; i < NumNames ; ++i ) 236 { 237 if ( NumSubscripts(names[i])!=0 ) 238 CoCoA_THROW_ERROR("names must have no subscripts", 239 "WANAMES(const vector<symbol>& names)"); 240 ans.push_back(symbol("d"+head(names[i]))); 241 } 242 return ans; 243 } 244 245 RingWeylImpl(const ring & CoeffRing,const std::vector<symbol> & WANames,const std::vector<long> & ElimIndets)246 RingWeylImpl::RingWeylImpl(const ring& CoeffRing, const std::vector<symbol>& WANames, const std::vector<long>& ElimIndets): 247 myReprRing(NewPolyRing(CoeffRing, 248 WANames, 249 NewMatrixOrdering(ElimMat(ElimIndets, len(WANames)), 250 0))), 251 myNumTrueIndetsValue(len(WANames)/2) 252 { 253 const long NumTrueIndets=len(WANames)/2; 254 CoCoA_ASSERT(IsCommutative(CoeffRing)); 255 myRefCountInc(); // this is needed for exception cleanliness, in case one of the lines below throws 256 myZeroPtr.reset(new RingElem(ring(this))); 257 myOnePtr.reset(new RingElem(ring(this), 1)); 258 // myIndetVector.resize(myNumIndetsValue, *myZeroPtr); 259 // myDerivationVector.resize(myNumIndetsValue, *myZeroPtr); 260 // for (long i=0; i < myNumIndetsValue; ++i) 261 myIndetVector.resize(2*NumTrueIndets, *myZeroPtr); 262 vector<long> expv(2*NumTrueIndets); 263 264 for (long i=0; i < 2*NumTrueIndets; ++i) 265 { 266 expv[i] = 1; 267 myPushFront(raw(myIndetVector[i]), raw(one(CoeffRing)), expv); 268 expv[i] = 0; 269 } 270 myRefCountZero(); // otherwise it is 2 + NumIndets and won't be destroyed 271 } 272 273 ~RingWeylImpl()274 RingWeylImpl::~RingWeylImpl() 275 {} 276 277 278 // inline const RingWeylImpl* RingWeyl::operator->() const 279 // { 280 // return static_cast<const RingWeylImpl*>(myRingPtr()); 281 // } 282 283 284 // inline WeylPoly& ring_weyl::AsWeylPoly(RawPtr rawx) const 285 // { 286 // return *x.WeylPolyPtr; 287 // // return *static_cast<ring_weyl::WeylPoly*>(x.WeylPolyPtr); 288 // } 289 290 291 // inline const WeylPoly& ring_weyl::AsWeylPoly(ConstRawPtr rawx) const 292 // { 293 // return *x.WeylPolyPtr; 294 // // return *static_cast<ring_weyl::WeylPoly*>(x.WeylPolyPtr); 295 // } 296 297 298 // RingWeyl::RingWeyl(const RingWeylImpl* RingPtr): 299 // SparsePolyRing(RingPtr) 300 // {} 301 302 303 //---------------------------------------------------------------------- 304 // Functions which every ring must implement: 305 //---------------------------------------------------------------------- 306 IamCommutative()307 bool RingWeylImpl::IamCommutative() const 308 { 309 return myNumIndets() == 0; // we are assuming the coeff ring is commutative 310 } 311 312 IamIntegralDomain3(bool QuickMode)313 bool3 RingWeylImpl::IamIntegralDomain3(bool QuickMode) const 314 { 315 return myReprRing->IamIntegralDomain3(QuickMode); //??? I think this is right 316 } 317 318 IamTrueGCDDomain()319 bool RingWeylImpl::IamTrueGCDDomain() const 320 { 321 return false; // I have no clue how to compute GCDs even if they do exist 322 } 323 324 myZero()325 ConstRefRingElem RingWeylImpl::myZero() const 326 { 327 return *myZeroPtr; 328 } 329 330 myOne()331 ConstRefRingElem RingWeylImpl::myOne() const 332 { 333 return *myOnePtr; 334 } 335 336 myNew()337 RingElemRawPtr RingWeylImpl::myNew() const 338 { 339 return myReprRing->myNew(); 340 } 341 342 myNew(const MachineInt & n)343 RingElemRawPtr RingWeylImpl::myNew(const MachineInt& n) const 344 { 345 return myReprRing->myNew(n); 346 } 347 348 myNew(const BigInt & N)349 RingElemRawPtr RingWeylImpl::myNew(const BigInt& N) const 350 { 351 return myReprRing->myNew(N); 352 } 353 354 myNew(ConstRawPtr rawcopy)355 RingElemRawPtr RingWeylImpl::myNew(ConstRawPtr rawcopy) const 356 { 357 return myReprRing->myNew(rawcopy); 358 } 359 360 myDelete(RawPtr rawx)361 void RingWeylImpl::myDelete(RawPtr rawx) const 362 { 363 myReprRing->myDelete(rawx); 364 } 365 366 mySwap(RawPtr rawx,RawPtr rawy)367 void RingWeylImpl::mySwap(RawPtr rawx, RawPtr rawy) const 368 { 369 myReprRing->mySwap(rawx, rawy); 370 } 371 372 myAssign(RawPtr rawlhs,ConstRawPtr rawx)373 void RingWeylImpl::myAssign(RawPtr rawlhs, ConstRawPtr rawx) const 374 { 375 myReprRing->myAssign(rawlhs, rawx); 376 } 377 378 myAssign(RawPtr rawlhs,const MachineInt & n)379 void RingWeylImpl::myAssign(RawPtr rawlhs, const MachineInt& n) const 380 { 381 myReprRing->myAssign(rawlhs, n); 382 } 383 384 myAssign(RawPtr rawlhs,const BigInt & N)385 void RingWeylImpl::myAssign(RawPtr rawlhs, const BigInt& N) const 386 { 387 myReprRing->myAssign(rawlhs, N); 388 } 389 390 myAssignZero(RawPtr rawlhs)391 void RingWeylImpl::myAssignZero(RawPtr rawlhs) const 392 { 393 myReprRing->myAssignZero(rawlhs); 394 } 395 396 myRecvTwinFloat(RawPtr rawlhs,ConstRawPtr rawx)397 void RingWeylImpl::myRecvTwinFloat(RawPtr rawlhs, ConstRawPtr rawx) const 398 { 399 CoCoA_ASSERT(!IsExact(myReprRing)); 400 myReprRing->myRecvTwinFloat(rawlhs, rawx); 401 } 402 myNegate(RawPtr rawlhs,ConstRawPtr rawx)403 void RingWeylImpl::myNegate(RawPtr rawlhs, ConstRawPtr rawx) const 404 { 405 myReprRing->myNegate(rawlhs, rawx); 406 } 407 408 myAdd(RawPtr rawlhs,ConstRawPtr rawx,ConstRawPtr rawy)409 void RingWeylImpl::myAdd(RawPtr rawlhs, ConstRawPtr rawx, ConstRawPtr rawy) const 410 { 411 myReprRing->myAdd(rawlhs, rawx, rawy); 412 } 413 414 mySub(RawPtr rawlhs,ConstRawPtr rawx,ConstRawPtr rawy)415 void RingWeylImpl::mySub(RawPtr rawlhs, ConstRawPtr rawx, ConstRawPtr rawy) const 416 { 417 myReprRing->mySub(rawlhs, rawx, rawy); 418 } 419 420 421 //??? ANNA: I believe the general implementation in SparsePolyRing works for RingWeyl 422 // void RingWeylImpl::myMul(RawPtr rawlhs, ConstRawPtr rawx, ConstRawPtr rawy) const 423 // { 424 // if (myReprRing->myIsZero(x) || myReprRing->myIsZero(y)) { myReprRing->myAssignZero(rawlhs); return; } 425 // RingElem f(myReprRing); 426 // myReprRing->myAssign(raw(f), rawx); 427 // ConstRefRingElem g(myReprRing, rawy); 428 // // std::clog<<"f="<<f<<std::endl; 429 // // std::clog<<"g="<<g<<std::endl; 430 // RingElem ans(myReprRing); // compute answer in a temporary for exception safety 431 // while (!CoCoA::IsZero(f)) 432 // { 433 // RingElem LMf(myReprRing); 434 // myMoveLMToFront(raw(LMf), raw(f)); 435 // PPMonoidElem LPPf = LPP(LMf); 436 // RingElem LMfg(g); 437 // myMulByPP(raw(LMfg), raw(LPPf)); 438 // // std::clog<<"LMf="<<LMf<<std::endl; 439 // // std::clog<<"LMfg="<<LMfg<<std::endl; 440 // myReprRing->myMulByCoeff(raw(LMfg), raw(LC(LMf))); //????UGLY!!!! 441 // ans += LMfg; 442 // } 443 // mySwap(rawlhs, raw(ans));// really an assignment -- is this safe???? 444 // } 445 446 447 // void RingWeylImpl::myDiv(RawPtr /*rawlhs*/, ConstRawPtr /*rawx*/, ConstRawPtr /*rawy*/) const 448 // { 449 // CoCoA_THROW_ERROR(ERR::NYI, "RingWeylImpl::myDiv"); 450 // return;//??? 451 // // bool OK; // Two lines o/w compiler complains that 452 // // OK = CoCoA::WeylDiv(AsDMPI(rawlhs), AsDMPI(x), AsDMPI(y)); // OK is unused when debugging is off. 453 // // CoCoA_ASSERT(OK); 454 // } 455 456 myIsDivisible(RawPtr,ConstRawPtr,ConstRawPtr)457 bool RingWeylImpl::myIsDivisible(RawPtr /*rawlhs*/, ConstRawPtr /*rawx*/, ConstRawPtr /*rawy*/) const 458 { 459 CoCoA_THROW_ERROR(ERR::NYI, "RingWeylImpl::myIsDivisible"); 460 return false; 461 // return CoCoA::WeylDiv(AsDMPI(rawlhs), AsDMPI(x), AsDMPI(y)); 462 } 463 464 myIsInvertible(ConstRawPtr)465 bool RingWeylImpl::myIsInvertible(ConstRawPtr /*rawx*/) const 466 { 467 CoCoA_THROW_ERROR(ERR::NYI, "RingWeylImpl::myIsInvertible"); //??? d[1]*x[1] = 1 468 return false;//??? 469 } 470 471 472 // bool RingWeylImpl::myIsPrintAtom(ConstRawPtr rawx) const 473 // { 474 // return myReprRing->myIsPrintAtom(rawx); //??? default always false 475 // } 476 477 myGcd(RawPtr,ConstRawPtr,ConstRawPtr)478 void RingWeylImpl::myGcd(RawPtr /*rawlhs*/, ConstRawPtr /*rawx*/, ConstRawPtr /*rawy*/) const 479 { 480 CoCoA_THROW_ERROR(ERR::NYI, "RingWeylImpl::myGcd"); 481 } 482 483 myPowerSmallExp(RawPtr rawlhs,ConstRawPtr rawx,long n)484 void RingWeylImpl::myPowerSmallExp(RawPtr rawlhs, ConstRawPtr rawx, long n) const 485 { 486 // Assert that we have a genuinely non-trivial case. 487 CoCoA_ASSERT(n > 1); 488 CoCoA_ASSERT(!myIsZero(rawx) && !myIsOne(rawx) && !myIsMinusOne(rawx)); 489 490 mySequentialPower(rawlhs, rawx, n); // probably better than myBinaryPower, I guess. 491 } 492 493 mySymbols(std::vector<symbol> & SymList)494 void RingWeylImpl::mySymbols(std::vector<symbol>& SymList) const 495 { //??? ANNA: which symbols are in a RingWeylImpl?? 496 myReprRing->mySymbols(SymList); 497 } 498 499 myOutput(std::ostream & out,ConstRawPtr rawx)500 void RingWeylImpl::myOutput(std::ostream& out, ConstRawPtr rawx) const 501 { 502 if (!out) return; // short-cut for bad ostreams 503 myReprRing->myOutput(out, rawx); 504 } 505 506 507 // void RingWeylImpl::myOutputSelf(std::ostream& out) const 508 // { 509 // out << "RingWeylImpl(" << CoeffRing(myReprRing) << ", " << NumIndets(myReprRing) << ")"; 510 // //????????? WHAT ABOUT THE ORDERING AND GRADING???? 511 // } 512 513 514 // void RingWeylImpl::myOutputSelf(OpenMathOutput& OMOut) const 515 // { 516 // OMOut->mySendApplyStart(); 517 // OMOut << OpenMathSymbol("polyd", "poly_ring_d"); //??? weyl? 518 // OMOut << CoeffRing(myReprRing); 519 // OMOut << NumIndets(myReprRing); //???? losing the ordering and grading info here!!! 520 // OMOut->mySendApplyEnd(); 521 // } 522 523 myOutput(OpenMathOutput & OMOut,ConstRawPtr rawx)524 void RingWeylImpl::myOutput(OpenMathOutput& OMOut, ConstRawPtr rawx) const 525 { 526 myReprRing->myOutput(OMOut, rawx); 527 } 528 529 myIsZero(ConstRawPtr rawx)530 bool RingWeylImpl::myIsZero(ConstRawPtr rawx) const 531 { 532 return myReprRing->myIsZero(rawx); 533 } 534 535 myIsOne(ConstRawPtr rawx)536 bool RingWeylImpl::myIsOne(ConstRawPtr rawx) const 537 { 538 return myReprRing->myIsOne(rawx); 539 } 540 541 myIsMinusOne(ConstRawPtr rawx)542 bool RingWeylImpl::myIsMinusOne(ConstRawPtr rawx) const 543 { 544 return myReprRing->myIsMinusOne(rawx); 545 } 546 547 548 // bool RingWeylImpl::myIsZeroAddMul(RawPtr rawlhs, ConstRawPtr rawy, ConstRawPtr rawz) const; // use default definition 549 550 myIsEqual(ConstRawPtr rawx,ConstRawPtr rawy)551 bool RingWeylImpl::myIsEqual(ConstRawPtr rawx, ConstRawPtr rawy) const 552 { 553 return myReprRing->myIsEqual(rawx, rawy); 554 } 555 556 myCompose(const RingHom &,const RingHom & theta)557 RingHom RingWeylImpl::myCompose(const RingHom& /*phi*/, const RingHom& theta) const 558 { 559 CoCoA_THROW_ERROR(ERR::NYI, "RingWeylImpl::compose"); 560 return theta; // just to keep compiler quiet 561 } 562 563 //---------------------------------------------------------------------- 564 565 // long NumIndets(const RingWeyl& RW) 566 // { 567 // return RW->myNumIndets(); 568 // } 569 570 //---------------------------------------------------------------------- 571 // Functions which every PolyRing must implement 572 //---------------------------------------------------------------------- 573 myNumIndets()574 long RingWeylImpl::myNumIndets() const 575 { 576 return myReprRing->myNumIndets(); 577 } 578 579 myCoeffRing()580 const ring& RingWeylImpl::myCoeffRing() const 581 { 582 return myReprRing->myCoeffRing(); 583 } 584 585 myIndets()586 const std::vector<RingElem>& RingWeylImpl::myIndets() const 587 { 588 return myIndetVector; 589 } 590 591 myIndetPower(RawPtr rawf,long var,long exp)592 void RingWeylImpl::myIndetPower(RawPtr rawf, long var, long exp) const 593 { 594 myReprRing->myIndetPower(rawf, var, exp); 595 } 596 597 myNumTerms(ConstRawPtr rawx)598 long RingWeylImpl::myNumTerms(ConstRawPtr rawx) const 599 { 600 return myReprRing->myNumTerms(rawx); 601 } 602 603 myIsConstant(ConstRawPtr rawf)604 bool RingWeylImpl::myIsConstant(ConstRawPtr rawf) const 605 { 606 return myReprRing->myIsConstant(rawf); 607 } 608 609 myIsIndet(long & index,ConstRawPtr rawf)610 bool RingWeylImpl::myIsIndet(long& index, ConstRawPtr rawf) const 611 { 612 return myReprRing->myIsIndet(index, rawf); 613 } 614 615 myIsMonomial(ConstRawPtr rawf)616 bool RingWeylImpl::myIsMonomial(ConstRawPtr rawf) const 617 { 618 return myReprRing->myIsMonomial(rawf); 619 } 620 621 myStdDeg(ConstRawPtr rawf)622 long RingWeylImpl::myStdDeg(ConstRawPtr rawf) const 623 { 624 return myReprRing->myStdDeg(rawf); 625 } 626 627 myDeg(ConstRawPtr rawf,long var)628 long RingWeylImpl::myDeg(ConstRawPtr rawf, long var) const 629 { 630 CoCoA_ASSERT(!myIsZero(rawf)); 631 return myReprRing->myDeg(rawf, var);//???????? 632 } 633 634 myLC(ConstRawPtr rawf)635 RingElemAlias RingWeylImpl::myLC(ConstRawPtr rawf) const 636 { 637 CoCoA_ASSERT(!myIsZero(rawf)); 638 return myReprRing->myLC(rawf);//??????????? 639 } 640 641 myContent(RawPtr rawdest,ConstRawPtr rawf)642 void RingWeylImpl::myContent(RawPtr rawdest, ConstRawPtr rawf) const 643 { 644 myReprRing->myContent(rawdest, rawf); 645 } 646 647 myRemoveBigContent(RawPtr rawf)648 void RingWeylImpl::myRemoveBigContent(RawPtr rawf) const 649 { 650 myReprRing->myRemoveBigContent(rawf); 651 } 652 653 myMulByCoeff(RawPtr rawf,ConstRawPtr rawc)654 void RingWeylImpl::myMulByCoeff(RawPtr rawf, ConstRawPtr rawc) const // WEAK EXCEPTION GUARANTEE 655 { 656 myReprRing->myMulByCoeff(rawf, rawc); 657 } 658 659 myDivByCoeff(RawPtr rawf,ConstRawPtr rawc)660 bool RingWeylImpl::myDivByCoeff(RawPtr rawf, ConstRawPtr rawc) const // WEAK EXCEPTION GUARANTEE 661 { 662 return myReprRing->myDivByCoeff(rawf, rawc); 663 } 664 665 myDeriv(RawPtr,ConstRawPtr,ConstRawPtr)666 void RingWeylImpl::myDeriv(RawPtr /*rawlhs*/, ConstRawPtr /*rawf*/, ConstRawPtr /*rawx*/) const 667 { 668 CoCoA_THROW_ERROR(ERR::NYI, "RingWeylImpl::myDeriv"); 669 } 670 671 myHomCtor(const ring &,const RingHom &,const std::vector<RingElem> &)672 RingHom RingWeylImpl::myHomCtor(const ring& /*codomain*/, const RingHom& /*CoeffHom*/, const std::vector<RingElem>& /*IndetImages*/) const 673 { 674 CoCoA_THROW_ERROR("DOES NOT EXIST", "RingWeylImpl::myHomCtor"); 675 return IdentityHom(myReprRing); // just to keep compiler quiet 676 } 677 678 679 680 //---------------------------------------------------------------------- 681 // Functions which every SparsePolyRing must implement: 682 //---------------------------------------------------------------------- 683 myPPM()684 const PPMonoid& RingWeylImpl::myPPM() const 685 { 686 return myReprRing->myPPM(); 687 } 688 689 myMonomial(ConstRawPtr rawc,PPMonoidElemConstRawPtr rawpp)690 RingElem RingWeylImpl::myMonomial(ConstRawPtr rawc, PPMonoidElemConstRawPtr rawpp) const 691 { 692 RingElem ans(ring(this)); 693 RingElem m = myReprRing->myMonomial(rawc, rawpp); 694 mySwap(raw(ans), raw(m)); 695 return ans; 696 } 697 698 myBeginIter(ConstRawPtr rawf)699 SparsePolyIter RingWeylImpl::myBeginIter(ConstRawPtr rawf) const 700 { 701 return myReprRing->myBeginIter(rawf); 702 } 703 704 myEndIter(ConstRawPtr rawf)705 SparsePolyIter RingWeylImpl::myEndIter(ConstRawPtr rawf) const 706 { 707 return myReprRing->myEndIter(rawf); 708 } 709 710 myPushFront(RawPtr rawf,ConstRawPtr rawc,const std::vector<long> & expv)711 void RingWeylImpl::myPushFront(RawPtr rawf, ConstRawPtr rawc, const std::vector<long>& expv) const 712 { 713 myReprRing->myPushFront(rawf, rawc, expv);//??? how many elements does expv have??? 714 } 715 716 myPushBack(RawPtr rawf,ConstRawPtr rawc,const std::vector<long> & expv)717 void RingWeylImpl::myPushBack(RawPtr rawf, ConstRawPtr rawc, const std::vector<long>& expv) const 718 { 719 myReprRing->myPushBack(rawf, rawc, expv);//??? how many elements does expv have??? 720 } 721 722 myPushFront(RawPtr rawf,ConstRawPtr rawc,PPMonoidElemConstRawPtr rawpp)723 void RingWeylImpl::myPushFront(RawPtr rawf, ConstRawPtr rawc, PPMonoidElemConstRawPtr rawpp) const 724 { 725 myReprRing->myPushFront(rawf, rawc, rawpp); 726 } 727 728 myPushBack(RawPtr rawf,ConstRawPtr rawc,PPMonoidElemConstRawPtr rawpp)729 void RingWeylImpl::myPushBack(RawPtr rawf, ConstRawPtr rawc, PPMonoidElemConstRawPtr rawpp) const 730 { 731 myReprRing->myPushBack(rawf, rawc, rawpp); 732 } 733 734 myLPP(ConstRawPtr rawf)735 ConstRefPPMonoidElem RingWeylImpl::myLPP(ConstRawPtr rawf) const 736 { 737 CoCoA_ASSERT(!myIsZero(rawf)); 738 return myReprRing->myLPP(rawf); 739 } 740 741 742 // f = pp*f NOTE: pp is on the LEFT and f is on the RIGHT myMulByPP(RawPtr rawf,PPMonoidElemConstRawPtr rawpp)743 void RingWeylImpl::myMulByPP(RawPtr rawf, PPMonoidElemConstRawPtr rawpp) const 744 { 745 if (myReprRing->myIsZero(rawf)) return; 746 747 RingElem g(myReprRing); 748 myReprRing->myAssign(raw(g), rawf); 749 //??? std::clog<<"Alias(f)="<<g<<std::endl; 750 // const long nvars = myNumIndetsValue; 751 for (long idx=0; idx < myNumTrueIndetsValue; ++idx) 752 { 753 const long Didx = idx + myNumTrueIndetsValue; 754 const long d = PPM(myReprRing)->myExponent(rawpp, Didx); 755 //??? std::clog<<"deg of D["<<idx<<"]="<<d<<std::endl; 756 RingElem der(g); // copy of f in myReprRing 757 RingElem sum = g*IndetPower(myReprRing, Didx, d); 758 for (long i=1; i <= d; ++i) 759 { 760 der = deriv(der, idx); 761 if (IsZero(der)) break; 762 // std::clog<<"bin(" << d << ", " << i << ")="<<binomial(d, i)<<std::endl; 763 sum += binomial(d, i)*der*IndetPower(myReprRing, Didx, d-i); 764 // std::clog<<"summand="<<binomial(d, i)*der*IndetPower(myReprRing, Didx, d-i)<<std::endl; 765 } 766 g = sum; 767 // std::clog<<"g="<<g<<std::endl; 768 } 769 vector<long> expv(2*myNumTrueIndetsValue); 770 PPM(myReprRing)->myExponents(expv, rawpp); 771 for (long idx=0; idx < myNumTrueIndetsValue; ++idx) 772 expv[idx+myNumTrueIndetsValue] = 0; 773 PPMonoidElem pp2(PPM(myReprRing)); 774 PPM(myReprRing)->myAssign(raw(pp2), expv); 775 myReprRing->myMulByPP(raw(g), raw(pp2)); 776 //??? std::clog<<"mul gives "<<g<<std::endl; 777 mySwap(rawf, raw(g));// really an assignment -- is this safe???? 778 } 779 780 myIsZeroAddLCs(RawPtr rawf,RawPtr rawg)781 bool RingWeylImpl::myIsZeroAddLCs(RawPtr rawf, RawPtr rawg) const 782 { 783 return myReprRing->myIsZeroAddLCs(rawf, rawg); 784 } 785 786 myMoveLMToFront(RawPtr rawf,RawPtr rawg)787 void RingWeylImpl::myMoveLMToFront(RawPtr rawf, RawPtr rawg) const 788 { 789 CoCoA_ASSERT(!myIsZero(rawg)); 790 myReprRing->myMoveLMToFront(rawf, rawg); 791 } 792 793 myMoveLMToBack(RawPtr rawf,RawPtr rawg)794 void RingWeylImpl::myMoveLMToBack(RawPtr rawf, RawPtr rawg) const 795 { 796 CoCoA_ASSERT(!myIsZero(rawg)); 797 myReprRing->myMoveLMToBack(rawf, rawg); 798 } 799 800 myDeleteLM(RawPtr rawf)801 void RingWeylImpl::myDeleteLM(RawPtr rawf) const 802 { 803 CoCoA_ASSERT(!myIsZero(rawf)); 804 myReprRing->myDeleteLM(rawf); 805 } 806 807 myDivLM(RawPtr rawlhs,ConstRawPtr rawf,ConstRawPtr rawg)808 void RingWeylImpl::myDivLM(RawPtr rawlhs, ConstRawPtr rawf, ConstRawPtr rawg) const 809 { 810 CoCoA_ASSERT(!myIsZero(rawf) && !myIsZero(rawg)); 811 myReprRing->myDivLM(rawlhs, rawf, rawg); 812 } 813 814 myCmpLPP(ConstRawPtr rawf,ConstRawPtr rawg)815 int RingWeylImpl::myCmpLPP(ConstRawPtr rawf, ConstRawPtr rawg) const 816 { 817 CoCoA_ASSERT(!myIsZero(rawf) && !myIsZero(rawg)); 818 return myReprRing->myCmpLPP(rawf, rawg); 819 } 820 821 myAddClear(RawPtr rawf,RawPtr rawg)822 void RingWeylImpl::myAddClear(RawPtr rawf, RawPtr rawg) const 823 { 824 myReprRing->myAddClear(rawf, rawg); 825 } 826 827 myAppendClear(RawPtr rawf,RawPtr rawg)828 void RingWeylImpl::myAppendClear(RawPtr rawf, RawPtr rawg) const 829 { 830 myReprRing->myAppendClear(rawf, rawg); 831 } 832 833 myAddMulLM(RawPtr rawf,ConstRawPtr rawh,ConstRawPtr rawg)834 void RingWeylImpl::myAddMulLM(RawPtr rawf, ConstRawPtr rawh, ConstRawPtr rawg) const //???? delete me??? 835 { 836 myAddMulLM(rawf, rawh, rawg, DontSkipLMg); 837 } 838 839 myAddMulLM(RawPtr rawf,ConstRawPtr rawh,ConstRawPtr rawg,SkipLMFlag skip)840 void RingWeylImpl::myAddMulLM(RawPtr rawf, ConstRawPtr rawh, ConstRawPtr rawg, SkipLMFlag skip) const 841 { 842 CoCoA_ASSERT(myNumTerms(rawh)==1); 843 RingElem prod(ring(this), myNew(rawg)); 844 myMulByCoeff(raw(prod), raw(myLC(rawh))); 845 myMulByPP(raw(prod), raw(myLPP(rawh))); 846 myReprRing->myAddMulLM(rawf, raw(myOne()), raw(prod), skip); 847 } 848 849 myReductionStep(RawPtr rawf,ConstRawPtr rawg)850 void RingWeylImpl::myReductionStep(RawPtr rawf, ConstRawPtr rawg) const 851 { 852 PPMonoidElem q = myLPP(rawf)/myLPP(rawg); 853 RingElem c = -myLC(rawf)/myLC(rawg); // RingElem c = -LC(repf)/LC(repg); 854 RingElem qg(myReprRing); myReprRing->myAssign(raw(qg), rawg); // qg is copy of rawg 855 myMulByPP(raw(qg), raw(q)); // qg = q * g //??? should we use myAddMul 856 myMulByCoeff(raw(qg), raw(c)); 857 myAdd(rawf, rawf, raw(qg)); // repf += qg; 858 859 860 /* 861 //??? ConstRefRingElem aliasg(ring(this), rawg); 862 //??? std::clog<<"g="<<aliasg<<std::endl; 863 PPMonoidElem LPPf = myReprRing->myLPP(f); 864 PPMonoidElem LPPg = myReprRing->myLPP(g); 865 PPMonoidElem q = LPPf/LPPg; // quotient exists 866 RingElem LCf(CoeffRing(myReprRing)); 867 CoeffRing(myReprRing)->myAssign(raw(LCf), myReprRing->myRawLC(f)); 868 RingElem LCg(CoeffRing(myReprRing)); 869 CoeffRing(myReprRing)->myAssign(raw(LCg), myReprRing->myRawLC(g)); 870 RingElem c = -LCf/LCg; 871 /// RingElem c(myCoeffRing()); 872 /// myCoeffRing()->myDiv(raw(c), myReprRing->myRawLC(f), myReprRing->myRawLC(g)); 873 RingElem g1(myReprRing); 874 //??? std::clog<<"ReductionStep: q="<<q<<std::endl; 875 //??? std::clog<<"ReductionStep: c="<<c<<std::endl; 876 myReprRing->myAssign(raw(g1), rawg); 877 //??? std::clog<<"ReductionStep: g="<<g1<<std::endl; 878 myMul(raw(g1), raw(q)); // g1 = q*g; 879 myReprRing->myMulByCoeff(raw(g1), raw(c)); 880 //??? std::clog<<"ReductionStep: prod="<<g1<<std::endl; 881 { 882 //??? ConstRefRingElem aliasf(ring(this),f); 883 //??? std::clog<<"REDUCING f="<<aliasf<<std::endl; 884 //??? std::clog<<"REDN by g="<<aliasg<<std::endl; 885 myAdd(f, f, raw(g1)); 886 //??? std::clog<<"RESULT is "<<aliasf<<std::endl<<std::endl; 887 } 888 */ 889 } 890 891 myReductionStepGCD(RawPtr,ConstRawPtr,RingElem &)892 void RingWeylImpl::myReductionStepGCD(RawPtr /*rawf*/, ConstRawPtr /*rawg*/, RingElem& /*fscale*/) const 893 { 894 CoCoA_THROW_ERROR(ERR::NYI, "RingWeylImpl::ReductionStepGCD"); 895 } 896 897 898 //--------------------------------------------------------------------------- 899 // Functions to do with RingWeylImpl::HomImpl 900 901 HomImpl(const SparsePolyRing & domain,const ring & codomain,const RingHom & CoeffHom,const std::vector<RingElem> & IndetImages)902 RingWeylImpl::HomImpl::HomImpl(const SparsePolyRing& domain, const ring& codomain, const RingHom& CoeffHom, const std::vector<RingElem>& IndetImages): 903 SparsePolyRingBase::HomImpl(domain, codomain, CoeffHom, IndetImages) 904 { 905 CoCoA_THROW_ERROR(ERR::NYI, "RingWeylImpl::HomImpl::HomImpl"); 906 } 907 908 909 910 //--------------------------------------------------------------------------- 911 // Functions for the class RingWeylImpl::IdealImpl 912 913 // inheritance is delicate here: this function is **necessary** myIdealCtor(const std::vector<RingElem> & gens)914 ideal RingWeylImpl::myIdealCtor(const std::vector<RingElem>& gens) const 915 { 916 return ideal(new IdealImpl(SparsePolyRing(this), gens)); //??? ugly ??? 917 } 918 919 IdealImpl(const SparsePolyRing & P,const std::vector<RingElem> & gens)920 RingWeylImpl::IdealImpl::IdealImpl(const SparsePolyRing& P, const std::vector<RingElem>& gens): 921 SparsePolyRingBase::IdealImpl(P, gens) 922 {} 923 924 myTestIsMaximal()925 void RingWeylImpl::IdealImpl::myTestIsMaximal() const 926 { CoCoA_THROW_ERROR(ERR::NYI, "RingWeylImpl::IdealImpl::myTestIsMaximal()"); } 927 928 myTestIsPrimary()929 void RingWeylImpl::IdealImpl::myTestIsPrimary() const 930 { CoCoA_THROW_ERROR(ERR::NYI, "RingWeylImpl::IdealImpl::myTestIsPrimary()"); } 931 932 myTestIsPrime()933 void RingWeylImpl::IdealImpl::myTestIsPrime() const 934 { CoCoA_THROW_ERROR(ERR::NYI, "RingWeylImpl::IdealImpl::myTestIsPrime()"); } 935 936 myTestIsRadical()937 void RingWeylImpl::IdealImpl::myTestIsRadical() const 938 { CoCoA_THROW_ERROR(ERR::NYI, "RingWeylImpl::IdealImpl::myTestIsRadical()"); } 939 940 myReduceMod(RingElemRawPtr)941 void RingWeylImpl::IdealImpl::myReduceMod(RingElemRawPtr /*rawr*/) const 942 { CoCoA_THROW_ERROR(ERR::NYI, "RingWeylImpl::WeylIdeal::IdealImpl::myReduceMod"); } 943 944 IhaveElem(RingElemConstRawPtr)945 bool RingWeylImpl::IdealImpl::IhaveElem(RingElemConstRawPtr /*rawr*/) const 946 { 947 CoCoA_THROW_ERROR(ERR::NYI, "RingWeylImpl::WeylIdeal::IhaveElem"); 948 return false; /* just to keep compiler quiet */ 949 } 950 951 myIntersect(const ideal &)952 void RingWeylImpl::IdealImpl::myIntersect(const ideal& /*J*/) 953 { 954 CoCoA_THROW_ERROR(ERR::NYI, "RingWeylImpl::WeylIdeal::myIntersect"); //??? 955 } 956 957 myColon(const ideal &)958 void RingWeylImpl::IdealImpl::myColon(const ideal& /*J*/) 959 { 960 CoCoA_THROW_ERROR(ERR::NYI, "RingWeylImpl::WeylIdeal::myColon"); //??? 961 } 962 963 myDivMod(RingElemRawPtr,RingElemConstRawPtr,RingElemConstRawPtr)964 bool RingWeylImpl::IdealImpl::myDivMod(RingElemRawPtr /*rawlhs*/, RingElemConstRawPtr /*rawnum*/, RingElemConstRawPtr /*rawden*/) const 965 { 966 CoCoA_THROW_ERROR(ERR::NYI, "RingWeylImpl::WeylIdeal::myDivMod"); 967 return false; /* just to keep compiler quiet */ 968 } 969 970 myGBasis(const CpuTimeLimit & CheckForTimeout)971 const std::vector<RingElem>& RingWeylImpl::IdealImpl::myGBasis(const CpuTimeLimit& CheckForTimeout) const 972 { 973 if (IhaveGBasis()) return myGBasisValue; 974 CoCoA_ASSERT(myGBasisValue.empty()); 975 vector<RingElem> GensList, GBList; 976 GensList.insert(GensList.end(), myGensValue.begin(), myGensValue.end()); 977 bool IsHomog=false; 978 bool IsSatAlg=false; 979 GRingInfo GRI(myP, IsHomog, IsSatAlg, NewDivMaskEvenPowers(), CheckForTimeout); 980 GBCriteria criteria(GBCriteria::DontUseCoprime, GBCriteria::UseGM, 981 GBCriteria::UseBack, GBCriteria::DontUseDiv); 982 GReductor GR(GRI, GensList, 983 GReductor::AffineAlg, 984 Reductors::DontUseBorel, GReductor::DontUseDynamicAlg, 985 criteria); 986 GR.myDoGBasis(); 987 // GR.myDoAFFGBasis(); 988 GR.myGBasis(GBList); 989 myGBasisValue.insert(myGBasisValue.end(), GBList.begin(), GBList.end()); 990 IhaveGBasisValue = true; 991 return myGBasisValue; 992 } 993 994 995 //---------------------------------------------------------------------- 996 // Pseudo-ctors for (sparse) polynomial rings. 997 998 NewWeylAlgebra(const ring & CoeffRing,long NumIndets,const vector<long> & ElimIndets)999 SparsePolyRing NewWeylAlgebra(const ring& CoeffRing, long NumIndets, const vector<long>& ElimIndets) 1000 { 1001 return SparsePolyRing(new RingWeylImpl(CoeffRing, WANAMES(NumIndets), ElimIndets)); 1002 } 1003 1004 NewWeylAlgebra(const ring & CoeffRing,const std::vector<symbol> & names,const vector<long> & ElimIndets)1005 SparsePolyRing NewWeylAlgebra(const ring& CoeffRing, const std::vector<symbol>& names, const vector<long>& ElimIndets) 1006 { 1007 return SparsePolyRing(new RingWeylImpl(CoeffRing, WANAMES(names), ElimIndets)); 1008 } 1009 1010 1011 // const RingElem& indet(const RingWeyl& RW, long var) 1012 // { 1013 // if (var > CoCoA::NumIndets(RW)) CoCoA_THROW_ERROR("Indeterminate index too large", "indet(RW,var)"); 1014 // return RW->myIndetVector[var]; 1015 // } 1016 1017 1018 // const RingElem& derivation(const RingWeyl& RW, long var) 1019 // { 1020 // if (var > CoCoA::NumIndets(RW)) CoCoA_THROW_ERROR("Indeterminate index too large", "derivation(RW,var)"); 1021 // return RW->myDerivationVector[var]; 1022 // } 1023 1024 1025 } // end of namespace CoCoA 1026 1027 // #error "JUNK AFTER THIS LINE" 1028 // ????????????????????????????????????????????????????????????????????????????? 1029 // #include <algorithm> 1030 // #include <vector> 1031 // #include <iostream> 1032 1033 // using std::max; 1034 // using std::swap; 1035 // using std::vector; 1036 // using std::ostream; 1037 // using std::endl; // just for debugging 1038 1039 1040 // #include "CoCoA/deriv.H" 1041 // #include "CoCoA/ring_weyl.H" 1042 1043 // namespace // unnamed, for file local functions 1044 // { 1045 1046 // using namespace CoCoA; 1047 1048 1049 // //---------------------------------------------------------------------- 1050 1051 1052 // // // f = pp*g where the product is in the Weyl algebra 1053 // // void act(RingElem& f, const RingElem& pp, const RingElem& g) 1054 // // { 1055 // // ASSERT(IsPolyRing(owner(f))); 1056 // // const PolyRing& P = PolyRing::owner(f); 1057 // // ASSERT(&owner(g) == &P && &owner(pp) == &P); 1058 // // ASSERT(!IsZero(pp) && NumTerms(pp) == 1); 1059 // // // std::clog<<"entered ucha's action" << endl; 1060 // // // std::clog << "pp=" << pp << endl; 1061 // // f = g; 1062 // // if (IsZero(f)) return; 1063 // // const long nvars = NumIndets(P); 1064 1065 // // vector<long> expv(nvars); 1066 // // PPM(P).exponents(expv, raw(LPP(pp))); 1067 1068 // // for (long var = nvars/2; var < nvars; ++var) 1069 // // { 1070 // // // std::clog << "doing D var=" << var << endl; 1071 // // const long n = expv[var]; 1072 // // // std::clog << "order = " << n << endl; 1073 // // RingElem derf = f; 1074 // // f *= P.IndetPower(var, n); 1075 // // for (long i=1; i <= n; ++i) 1076 // // { 1077 // // derf = deriv(derf, var-nvars/2); 1078 // // f += binomial(n, i)*derf*P.IndetPower(var, n-i); // *IndetPower(h, 2*i); // for homog case 1079 // // std::clog<<"binomial("<<n<<","<<i<<")="<<binomial(n,i)<<std::endl; 1080 // // } 1081 // // } 1082 // // { // f *= var^deg(pp, var); for the early vars 1083 // // for (long var = nvars/2; var < nvars; ++var) 1084 // // expv[var] = 0; 1085 // // PPMonoid::elem qq(PPM(P), expv); 1086 // // f *= P.monomial(RingElem(CoeffRing(P), 1), qq); 1087 // // } 1088 // // } 1089 // //---------------------------------------------------------------------- 1090 // } // end of unnamed namespace 1091 1092 1093 // namespace CoCoA 1094 // { 1095 1096 1097 1098 1099 // RingWeylImpl::ring_weyl(const AbstractRing& R, const PPMonoid& PPM): 1100 // myCoeffRing(R), 1101 // myPPM(PPM), 1102 // myWeylPolyPool(sizeof(WeylPoly), "RingWeylImpl::myWeylPolyPool"), 1103 // myNumIndets(NumIndets(PPM)), 1104 // myOrdvWords(OrdvWords(ordering(PPM))), 1105 // mySummandSize(sizeof(WeylPoly::summand) + sizeof(int)*(myOrdvWords-1)), 1106 // mySummandPool(mySummandSize, "RingWeylImpl::mySummandPool") 1107 // { 1108 // myNumPolys = 0; 1109 // myZero = new RingElem(*this); 1110 // myIndetVector.resize(myNumIndets, RingElem(*this)); 1111 // vector<long> expv(myNumIndets); 1112 // RingElem one(myCoeffRing, 1); 1113 // for (long i=0; i < myNumIndets; ++i) 1114 // { 1115 // expv[i] = 1; 1116 // PushFront(raw(myIndetVector[i]), raw(one), expv); 1117 // expv[i] = 0; 1118 // } 1119 // } 1120 1121 1122 // RingWeylImpl::~ring_weyl() 1123 // { 1124 // myIndetVector.clear(); 1125 // delete myZero; 1126 // ASSERT(myNumPolys == 0); 1127 // } 1128 1129 1130 // void RingWeylImpl::MakeWritable(RawPtr rawx) const 1131 // { 1132 // if (x.WeylPolyPtr->myRefCount == 1) return; 1133 // --x.WeylPolyPtr->myRefCount; 1134 // WeylPoly* copy = static_cast<WeylPoly*>(myWeylPolyPool.alloc(sizeof(WeylPoly))); 1135 // const WeylPoly* orig = x.WeylPolyPtr; 1136 // x.WeylPolyPtr = copy; 1137 // copy->myRefCount = 1; 1138 // copy->mySummands = nullptr; 1139 // copy->myEnd = ©->mySummands; 1140 // for (const WeylPoly::summand* it = orig->mySummands; it != nullptr; it = it->myNext) 1141 // { 1142 // // more or less hand inlined body of PushBack -- NB using PushBack is "circular" 1143 // *copy->myEnd = CopySummand(it); 1144 // copy->myEnd = &((*copy->myEnd)->myNext); 1145 // } 1146 1147 // /// myCoeffRing.init(copy->myDenom); 1148 // // new (©->myLC) RingElem(myCoeffRing, AbstractRing::elem::ALIAS); 1149 // } 1150 1151 1152 // inline WeylPoly::summand* RingWeylImpl::AllocSummand() const 1153 // { 1154 // return static_cast<WeylPoly::summand*>(mySummandPool.alloc(mySummandSize)); 1155 // } 1156 1157 1158 // WeylPoly::summand* RingWeylImpl::InitSummand(WeylPoly::summand* ptr) const 1159 // { 1160 // myCoeffRing.init(ptr->myCoeff); 1161 // ptr->myNext = nullptr; 1162 // ptr->myModulePosn = 0; 1163 // return ptr; 1164 // } 1165 1166 1167 // WeylPoly::summand* RingWeylImpl::InitSummand(WeylPoly::summand* ptr, ConstRawPtr rawc, const vector<long>& expv) const 1168 // { 1169 // myCoeffRing.init(ptr->myCoeff, c); 1170 // ptr->myNext = nullptr; 1171 // ptr->myModulePosn = 0; 1172 // ordering(myPPM).ComputeOrdv(ptr->myOrdv, &expv[0]); // &expv[0] converts vector<T> to T* 1173 // return ptr; 1174 // } 1175 1176 1177 // WeylPoly::summand* RingWeylImpl::CopySummand(const WeylPoly::summand* original) const 1178 // { 1179 // WeylPoly::summand* copy = AllocSummand(); 1180 // copy->myNext = nullptr; 1181 // myCoeffRing.init(copy->myCoeff, original->myCoeff); 1182 // copy->myModulePosn = original->myModulePosn; 1183 // for (long i=0; i < myOrdvWords; ++i) 1184 // copy->myOrdv[i] = original->myOrdv[i]; 1185 // return copy; 1186 // } 1187 1188 1189 // void RingWeylImpl::SetSummandMOrdv(WeylPoly::summand* dest, const WeylPoly::summand* src) const 1190 // { 1191 // dest->myModulePosn = src->myModulePosn; 1192 // for (long i=0; i < myOrdvWords; ++i) 1193 // dest->myOrdv[i] = src->myOrdv[i]; 1194 // } 1195 1196 1197 // void RingWeylImpl::DeleteSummands(WeylPoly::summand* ptr) const 1198 // { 1199 // WeylPoly::summand* next; 1200 // while (ptr != nullptr) 1201 // { 1202 // next = ptr->myNext; 1203 // myCoeffRing.kill(ptr->myCoeff); 1204 // mySummandPool.free(ptr, mySummandSize); 1205 // ptr = next; 1206 // } 1207 // } 1208 1209 1210 // bool RingWeylImpl::EqualSummands(const WeylPoly::summand& lhs, const WeylPoly::summand& x) const 1211 // { 1212 // if (lhs.myModulePosn != x.myModulePosn) return false; 1213 // const PPOrdering::OrdvElem* const lordv = lhs.myOrdv; 1214 // const PPOrdering::OrdvElem* const rordv = x.myOrdv; 1215 // for (long i = 0; i < myOrdvWords; ++i) 1216 // if (lordv[i] != rordv[i]) return false; 1217 // return myCoeffRing.IsEqual(lhs.myCoeff, x.myCoeff); 1218 // } 1219 1220 1221 // inline void RingWeylImpl::MulOrdv(PPOrdering::OrdvElem* ov, const PPOrdering::OrdvElem* ov1, const PPOrdering::OrdvElem* ov2) const 1222 // { 1223 // for (long i=0; i < myOrdvWords; ++i) 1224 // ov[i] = ov1[i]+ov2[i]; 1225 // } 1226 1227 1228 // inline void RingWeylImpl::DivOrdv(PPOrdering::OrdvElem* ov, const PPOrdering::OrdvElem* ov1, const PPOrdering::OrdvElem* ov2) const 1229 // { 1230 // for (long i=0; i < myOrdvWords; ++i) 1231 // ov[i] = ov1[i]-ov2[i]; 1232 // } 1233 1234 1235 1236 1237 // long RingWeylImpl::NumIndets() const 1238 // { 1239 // return myNumIndets; 1240 // } 1241 1242 1243 // const AbstractRing& RingWeylImpl::CoeffRing() const 1244 // { 1245 // return myCoeffRing; 1246 // } 1247 1248 1249 // const PPMonoid& RingWeylImpl::PPM() const 1250 // { 1251 // return myPPM; 1252 // } 1253 1254 1255 // long RingWeylImpl::NumTerms(ConstRawPtr rawx) const 1256 // { 1257 // long nsummands = 0; 1258 // for (const WeylPoly::summand* it = AsWeylPoly(rawx).mySummands; it != nullptr; it = it->myNext) ++nsummands; 1259 // return nsummands; 1260 // } 1261 1262 1263 // PolyIter RingWeylImpl::BeginIter(AbstractRing::RawPtr rawx) const 1264 // { 1265 // return PolyIter(*this, &AsWeylPoly(rawx).mySummands); 1266 // } 1267 1268 1269 // PolyIter RingWeylImpl::EndIter(AbstractRing::RawPtr rawx) const 1270 // { 1271 // return PolyIter(*this, AsWeylPoly(rawx).myEnd); 1272 // } 1273 1274 1275 // const RingElem& RingWeylImpl::indet(long var) const 1276 // { 1277 // ASSERT(var < myNumIndets); 1278 // return myIndetVector[var]; 1279 // } 1280 1281 1282 // RingElem RingWeylImpl::IndetPower(long var, long n) const 1283 // { 1284 // ASSERT(0 <= var && var < myNumIndets); 1285 // return monomial(CoCoA::IndetPower(myPPM, var, n)); 1286 // } 1287 1288 1289 // RingElem RingWeylImpl::monomial(const RingElem& c, const PPMonoid::alias& pp) const 1290 // { 1291 // vector<long> expv(myNumIndets); 1292 // myPPM.exponents(expv, raw(pp)); 1293 // RingElem ans(*this); 1294 // PushFront(raw(ans), raw(c), expv); 1295 // return ans; 1296 // } 1297 1298 1299 // RingElem RingWeylImpl::monomial(const PPMonoid::alias& pp) const 1300 // { 1301 // RingElem c(myCoeffRing, 1); 1302 // vector<long> expv(myNumIndets); 1303 // myPPM.exponents(expv, raw(pp)); 1304 // RingElem ans(*this); 1305 // PushFront(raw(ans), raw(c), expv); 1306 // return ans; 1307 // } 1308 1309 1310 // RingElem RingWeylImpl::monomial(const RingElem& c) const 1311 // { 1312 // vector<long> expv(myNumIndets); 1313 // RingElem ans(*this); 1314 // PushFront(raw(ans), raw(c), expv); 1315 // return ans; 1316 // } 1317 1318 1319 // long RingWeylImpl::deg(ConstRawPtr rawf) const 1320 // { 1321 // if (IsZero(f)) CoCoA_THROW_ERROR("RingWeylImpl::deg: cannot compute degree of zero polynomial"); 1322 // if (GradingDim(myPPM) > 0) 1323 // return ordering(myPPM).deg(AsWeylPoly(f).mySummands->myOrdv); //// BUG???? "valid" only if grading dim == 1 1324 // // return ???; the vector in Z^0 -- ring not graded!!! 1325 // return -1;//?????????? BUG BUG INCOMPLETE 1326 // } 1327 1328 1329 // int RingWeylImpl::deg(ConstRawPtr rawf, long var) const 1330 // { 1331 // if (IsZero(f)) CoCoA_THROW_ERROR("RingWeylImpl::deg: cannot compute degree of zero polynomial"); 1332 // long d = 0; 1333 // for (WeylPoly::summand* it = AsWeylPoly(f).mySummands; it != nullptr; it = it->myNext) 1334 // d = max(d, ordering(myPPM).exponent(it->myOrdv, var)); 1335 // return d; 1336 // } 1337 1338 1339 // // int RingWeylImpl::deg(ConstRawPtr rawf, const PPGrading& G) const 1340 // // { 1341 // // if (IsZero(f)) CoCoA_THROW_ERROR("RingWeylImpl::deg: cannot compute degree of zero polynomial"); 1342 // // int d = 0; 1343 // // for (WeylPoly::summand* it = AsWeylPoly(f).mySummands; it != nullptr; it = it->myNext) 1344 // // d = max(d, deg(it, G)); // CANNOT JUST USE MAX HERE!! MUST USE LEX MAX FOR VECTORS 1345 // // return d; 1346 // // } 1347 1348 1349 1350 // const RingElem RingWeylImpl::LC(ConstRawPtr rawf) const 1351 // { 1352 // ASSERT(!IsZero(f)); 1353 // // if (IsZero(f)) return ::zero(myCoeffRing); 1354 // return RingElem(myCoeffRing, AsWeylPoly(f).mySummands->myCoeff); 1355 // // AsWeylPoly(f).myLC.reseat(AsWeylPoly(f).mySummands->myCoeff); 1356 // // return AsWeylPoly(f).myLC; 1357 // } 1358 1359 1360 // AbstractRing::RawPtr RingWeylImpl::RawLC(RawPtr rawf) const 1361 // { 1362 // MakeWritable(f);//????????? 1363 // ASSERT(!IsZero(f)); 1364 // return AsWeylPoly(f).mySummands->myCoeff; 1365 // } 1366 1367 1368 // const AbstractRing::RawPtr RingWeylImpl::RawLC(ConstRawPtr rawf) const 1369 // { 1370 // ASSERT(!IsZero(f)); 1371 // return AsWeylPoly(f).mySummands->myCoeff; 1372 // } 1373 1374 1375 // RingElem RingWeylImpl::content(ConstRawPtr rawf) const 1376 // { 1377 // ASSERT(::IsTrueGCDDomain(myCoeffRing)); 1378 // ASSERT(!IsZero(f)); 1379 // RingElem cont(myCoeffRing); 1380 // for (WeylPoly::summand* it = AsWeylPoly(f).mySummands; it != nullptr; it = it->myNext) 1381 // myCoeffRing.gcd(raw(cont), raw(cont), it->myCoeff); // be clever if cont == 1?? 1382 // return cont; 1383 // } 1384 1385 1386 1387 // void RingWeylImpl::MoveLMToFront(RawPtr rawf, RawPtr rawg) const 1388 // { 1389 // ASSERT(!IsZero(g)); 1390 // ASSERT(f.WeylPolyPtr->myRefCount == 1); 1391 // // ASSERT(g.WeylPolyPtr->myRefCount == 1); 1392 // // MakeWritable(f);//???????? 1393 // MakeWritable(g);//???????? 1394 // WeylPoly& G = AsWeylPoly(g); 1395 // WeylPoly::summand* ltg = G.mySummands; 1396 // G.mySummands = G.mySummands->myNext; 1397 // if (G.mySummands == nullptr) G.myEnd = &(G.mySummands); 1398 // PushFront(f, ltg); 1399 // } 1400 1401 1402 // void RingWeylImpl::DeleteLM(RawPtr rawf) const 1403 // { 1404 // ASSERT(!IsZero(f)); 1405 // // ASSERT(f.WeylPolyPtr->myRefCount == 1); 1406 // MakeWritable(f);//???????? 1407 // WeylPoly& F = AsWeylPoly(f); 1408 // WeylPoly::summand* old_ltf = F.mySummands; 1409 // F.mySummands = old_ltf->myNext; 1410 // if (F.mySummands == nullptr) F.myEnd = &F.mySummands; 1411 // old_ltf->myNext = nullptr; 1412 // DeleteSummands(old_ltf); 1413 // } 1414 1415 1416 // void RingWeylImpl::AddMul(RawPtr rawf, const WeylPoly::summand* s, ConstRawPtr rawg, bool SkipLM) const 1417 // { 1418 // //// std::clog << "AddMul: Doing funny product of the following two polys" << endl; 1419 // //// output(std::clog, rawg); 1420 // ASSERT(f.WeylPolyPtr->myRefCount == 1); 1421 // // MakeWritable(f); 1422 1423 // RingElem ppg(*this); 1424 // assign(raw(ppg), rawg); 1425 // vector<long> expv(myNumIndets); 1426 // ordering(myPPM).ComputeExpv(&expv[0], s->myOrdv); 1427 // //// std::clog << "expv: "; for (int i=0; i<myNumIndets;++i) std::clog << expv[i] << " "; std::clog << endl; 1428 // for (long var = myNumIndets/2; var < myNumIndets; ++var) 1429 // { 1430 // const long n = expv[var]; 1431 // if (n == 0) continue; 1432 // //// std::clog << "AddMul: doing D variable with index " << n << endl; 1433 // RingElem der(*this); 1434 // assign(raw(der), raw(ppg)); 1435 // // ppg *= IndetPower(var, n); CANNOT DO THIS --> INFINITE RECURSION 1436 // { 1437 // PlainMul(raw(ppg),raw(ppg),raw(IndetPower(var, n))); 1438 // } 1439 // // mul(raw(ppg), raw(ppg), raw(IndetPower(var, n))); 1440 // for (long i=1; i <= n; ++i) 1441 // { 1442 // deriv(raw(der), raw(der), var-myNumIndets/2); 1443 // //// std::clog << "der(" << i << ")="; output(std::clog, raw(der)); std::clog << endl; 1444 // // ppg += binomial(n, i)*der*IndetPower(var, n-i); // *IndetPower(h, 2*i); // for homog case 1445 // { 1446 // vector<long> expv2(myNumIndets); 1447 // expv2[var] = n-i; 1448 // RingElem shift(*this); 1449 // PushBack(raw(shift), raw(RingElem(myCoeffRing, binomial(n, i))), expv2); 1450 // RingElem tmp(*this); 1451 // PlainMul(raw(tmp), raw(shift), raw(der)); 1452 // //// std::clog << "AddMul: adding "; output(std::clog, raw(tmp)); std::clog << endl; 1453 // add(raw(ppg),raw(ppg),raw(tmp)); 1454 // } 1455 // } 1456 // } 1457 // { // f *= var^deg(pp, var); for the early vars 1458 // for (long var = myNumIndets/2; var < myNumIndets; ++var) 1459 // expv[var] = 0; 1460 // PPMonoid::elem qq(myPPM, expv); 1461 // RingElem c(myCoeffRing); 1462 // myCoeffRing.assign(raw(c), s->myCoeff); 1463 // PlainMul(raw(ppg), raw(ppg), raw(monomial(c, qq))); 1464 // /// f *= P.monomial(RingElem(CoeffRing(P), 1), qq); 1465 // } 1466 // //// std::clog << "AddMul: OUTPUT "; output(std::clog, raw(ppg)); std::clog << endl; 1467 // add(f, f, raw(ppg)); 1468 // } 1469 1470 1471 // void RingWeylImpl::AddMul2(RawPtr rawf, ConstRawPtr rawh, ConstRawPtr rawg, bool SkipLM) const //???? delete me??? 1472 // { //??? 1473 // AddMul(f, (AsWeylPoly(h).mySummands), rawg, SkipLM); //??? 1474 // } //??? 1475 1476 // void RingWeylImpl::PlainMul(RawPtr rawlhs, ConstRawPtr rawx, ConstRawPtr rawy) const 1477 // { 1478 // if (NumTerms(x) > NumTerms(y)) { PlainMul(lhs, rawy, rawx); return; } 1479 1480 // RingElem ans(*this); 1481 1482 // for (const WeylPoly::summand* xterm = AsWeylPoly(x).mySummands; xterm; xterm = xterm->myNext) 1483 // PlainAddMul(raw(ans), xterm, rawy, false);//????? 1484 1485 // swap(raw(ans), lhs); // really an assignment 1486 // } 1487 // void RingWeylImpl::PlainAddMul(RawPtr rawf, const WeylPoly::summand* s, ConstRawPtr rawg, bool SkipLM) const 1488 // { 1489 // // ASSERT(f.DMPIPtr->myRefCount == 1); 1490 // // MakeWritable(f); 1491 // // std::clog << "\nF := "; output(std::clog,f); 1492 // // std::clog<<";\nG := "; output(std::clog,g); 1493 // // std::clog<<endl; 1494 // // std::clog<<"s = "; 1495 // // myCoeffRing.output(std::clog, s->myCoeff); 1496 // // std::clog<<"x^("; 1497 // // for(long i=0;i<myOrdvWords;++i)std::clog<<s->myOrdv[i]<<" "; 1498 // // std::clog<<")"<<endl; 1499 // const AbstractRing& R = myCoeffRing; 1500 // typedef WeylPoly::summand summand; 1501 // const summand* g_smnd = AsWeylPoly(g).mySummands; 1502 // if (SkipLM) g_smnd = g_smnd->myNext; 1503 // summand** f_prev = &(AsWeylPoly(f).mySummands); 1504 // summand* f_smnd = *f_prev; 1505 1506 // summand* tmp_smnd = InitSummand(AllocSummand()); // just sets coeff = 0 1507 1508 // int CMP=0; 1509 1510 // // bool qIsOne = (myPPM.cmp(q, tmpPP)==0); 1511 // bool qIsOne = false; 1512 1513 // for (; f_smnd != nullptr && g_smnd != nullptr; g_smnd = g_smnd->myNext) 1514 // { 1515 // if (qIsOne) 1516 // while (f_smnd != nullptr && (CMP=ordering(myPPM).CmpOrdvs(f_smnd->myOrdv,g_smnd->myOrdv)) >0) 1517 // f_smnd = *(f_prev = &f_smnd->myNext); 1518 // else 1519 // { 1520 // MulOrdv(tmp_smnd->myOrdv, s->myOrdv, g_smnd->myOrdv); 1521 // while (f_smnd != nullptr && (CMP=ordering(myPPM).CmpOrdvs(f_smnd->myOrdv,tmp_smnd->myOrdv)) >0) 1522 // f_smnd = *(f_prev = &f_smnd->myNext); 1523 // } 1524 // if (f_smnd == nullptr) 1525 // { 1526 // R.mul(tmp_smnd->myCoeff, s->myCoeff, g_smnd->myCoeff); 1527 // PushBack(f, tmp_smnd); 1528 // tmp_smnd = InitSummand(AllocSummand()); 1529 // g_smnd = g_smnd->myNext; 1530 // break; 1531 // } 1532 // if (CMP == 0) 1533 // { 1534 // if (R.IsZeroAddMul(f_smnd->myCoeff, s->myCoeff, g_smnd->myCoeff)) 1535 // RemoveSmnd(f, f_prev); // f_prev = f_prev; 1536 // else 1537 // f_prev = &f_smnd->myNext; 1538 // f_smnd = *f_prev; 1539 // } 1540 // else // (CMP < 0) 1541 // { 1542 // R.mul(tmp_smnd->myCoeff, s->myCoeff, g_smnd->myCoeff); 1543 // InsertSmnd(f, tmp_smnd, f_prev); 1544 // tmp_smnd = InitSummand(AllocSummand()); 1545 // f_prev = &(*f_prev)->myNext; 1546 // // f_smnd = f_smnd; 1547 // } 1548 // } 1549 // for (;g_smnd != nullptr; g_smnd = g_smnd->myNext) 1550 // { 1551 // MulOrdv(tmp_smnd->myOrdv, s->myOrdv, g_smnd->myOrdv); 1552 // R.mul(tmp_smnd->myCoeff, s->myCoeff, g_smnd->myCoeff); 1553 // PushBack(f, tmp_smnd); 1554 // tmp_smnd = InitSummand(AllocSummand()); 1555 // } 1556 // DeleteSummands(tmp_smnd); // next ptr will be zero (set by InitSummand) 1557 // // std::clog << "AddMul: produced f=";output(std::clog,f);std::clog<<endl; 1558 // // std::clog << "------------------------------------------------------"<<endl; 1559 1560 // } 1561 1562 1563 1564 1565 1566 // void RingWeylImpl::ReductionStep(RawPtr rawf, ConstRawPtr rawg) const 1567 // { 1568 // // std::clog << "\nRingWeylImpl::reduce" << endl; 1569 // // std::clog << "\nF := "; output(std::clog,f); 1570 // // std::clog<<";\nG := "; output(std::clog,g); 1571 // // std::clog<<";"<<endl; 1572 // ASSERT(&g!=&f); 1573 // const AbstractRing& R = myCoeffRing; 1574 // const WeylPoly::summand* g_smnd = AsWeylPoly(g).mySummands; 1575 // const WeylPoly::summand* f_smnd = AsWeylPoly(f).mySummands; 1576 // WeylPoly::summand* tmp_smnd = InitSummand(AllocSummand()); // just sets coeff = 0 1577 1578 // DivOrdv(tmp_smnd->myOrdv, f_smnd->myOrdv, g_smnd->myOrdv); 1579 // R.div(tmp_smnd->myCoeff, f_smnd->myCoeff, g_smnd->myCoeff); 1580 // R.negate(tmp_smnd->myCoeff, tmp_smnd->myCoeff); 1581 1582 // // std::clog<<"-- S := "; 1583 // // myCoeffRing.output(std::clog, tmp_smnd->myCoeff); 1584 // // std::clog<<"x^("; 1585 // // for(long i=0;i<myOrdvWords;++i)std::clog<<tmp_smnd->myOrdv[i]<<" "; 1586 // // std::clog<<")"<<endl; 1587 1588 // DeleteLM(f); 1589 // AddMul(f, tmp_smnd, g, true /*SkipLM*/); 1590 1591 // DeleteSummands(tmp_smnd); // next ptr will be zero (set by InitSummand) 1592 // // std::clog << "H := "; output(std::clog,f); 1593 // // std::clog << ";\n--------------------------------------------------"<<endl; 1594 // } 1595 1596 1597 // void RingWeylImpl::AddClear(RawPtr rawf, RawPtr rawg) const 1598 // { 1599 // ASSERT(f.WeylPolyPtr->myRefCount == 1); 1600 // // MakeWritable(f); 1601 // MakeWritable(g); 1602 // // input polynomial are copied 1603 // //if (g.WeylPolyPtr->myRefCount != 1) 1604 // // std::clog << "AddClear: g.myRefCount == " << g.WeylPolyPtr->myRefCount << endl; 1605 // const AbstractRing& R = myCoeffRing; 1606 // typedef WeylPoly::summand summand; 1607 // WeylPoly& F = AsWeylPoly(f); 1608 // WeylPoly& G = AsWeylPoly(g); 1609 // summand* g_smnd = G.mySummands; 1610 // summand** f_prev = &(F.mySummands); 1611 // summand* f_smnd = *f_prev; 1612 // int CMP=0; 1613 // ASSERT(*(G.myEnd)==nullptr);//BUG HUNTING ??? 1614 1615 // // std::clog << "input f = "; output(std::clog, f) ;std::clog << endl; 1616 // while ( f_smnd!=nullptr && g_smnd!=nullptr ) 1617 // { 1618 // while (f_smnd!=nullptr && 1619 // (CMP=ordering(myPPM).CmpOrdvs(f_smnd->myOrdv,g_smnd->myOrdv)) >0) 1620 // f_smnd = *(f_prev = &f_smnd->myNext); 1621 // if (f_smnd == nullptr) break; 1622 // //std::clog << "(AddClear error: should never happen for Basic Reduction)" << endl; 1623 // G.mySummands = G.mySummands->myNext; 1624 // g_smnd->myNext = nullptr; 1625 // if (CMP == 0) 1626 // { 1627 // R.add(f_smnd->myCoeff, f_smnd->myCoeff, g_smnd->myCoeff); 1628 // if (R.IsZero(f_smnd->myCoeff)) 1629 // RemoveSmnd(f, f_prev); 1630 // DeleteSummands(g_smnd); 1631 // } 1632 // else // (CMP < 0) 1633 // { 1634 // InsertSmnd(f, g_smnd, f_prev); 1635 // f_prev = &(*f_prev)->myNext; 1636 // } 1637 // f_smnd = *f_prev; 1638 // g_smnd = G.mySummands; 1639 // } 1640 // if (G.mySummands!=nullptr) 1641 // { 1642 // *(F.myEnd) = G.mySummands; 1643 // F.myEnd = G.myEnd; 1644 // G.mySummands = nullptr; 1645 // } 1646 // G.myEnd = &G.mySummands; 1647 // // if (rare) {std::clog << "f2 = "; output(std::clog, f) ;std::clog << endl;} 1648 // } 1649 1650 1651 // void RingWeylImpl::AppendClear(RawPtr rawf, RawPtr rawg) const 1652 // { 1653 // ASSERT(f.WeylPolyPtr->myRefCount == 1); 1654 // // MakeWritable(f); 1655 // MakeWritable(g); 1656 // // input polynomial are copied 1657 // //if (g.WeylPolyPtr->myRefCount != 1) 1658 // // std::clog << "AppendClear: g.myRefCount == " << g.WeylPolyPtr->myRefCount << endl; 1659 // WeylPoly& F = AsWeylPoly(f); 1660 // WeylPoly& G = AsWeylPoly(g); 1661 // if (G.mySummands!=nullptr) 1662 // { 1663 // *(F.myEnd) = G.mySummands; 1664 // F.myEnd = G.myEnd; 1665 // G.mySummands = nullptr; 1666 // } 1667 // G.myEnd = &G.mySummands; 1668 // // if (rare) {std::clog << "f2 = "; output(std::clog, f) ;std::clog << endl;} 1669 // } 1670 1671 1672 // int RingWeylImpl::CmpLPP(ConstRawPtr rawf, ConstRawPtr rawg) const 1673 // { 1674 // ASSERT(!IsZero(f)); 1675 // ASSERT(!IsZero(g)); 1676 // return ordering(myPPM).CmpOrdvs(AsWeylPoly(f).mySummands->myOrdv,AsWeylPoly(g).mySummands->myOrdv); 1677 // } 1678 1679 1680 // void RingWeylImpl::DivLM(RawPtr rawlhs, ConstRawPtr rawf, ConstRawPtr rawg) const 1681 // { 1682 // // std::clog << "DivLM" << endl; 1683 // const AbstractRing& R = myCoeffRing; 1684 // typedef WeylPoly::summand summand; 1685 // const summand* f_smnd = AsWeylPoly(f).mySummands; 1686 // const summand* g_smnd = AsWeylPoly(g).mySummands; 1687 // MakeWritable(lhs); 1688 // assign(lhs,0); 1689 // summand* SpareSummand = InitSummand(AllocSummand()); 1690 // R.div(SpareSummand->myCoeff, f_smnd->myCoeff, g_smnd->myCoeff); 1691 // DivOrdv(SpareSummand->myOrdv, f_smnd->myOrdv, g_smnd->myOrdv); 1692 // PushBack(lhs, SpareSummand); 1693 // } 1694 1695 1696 // void RingWeylImpl::mul(RawPtr rawf, const PPMonoid::elem& pp) const 1697 // { 1698 // // bool qIsOne = (myPPM.cmp(q, tmpPP)==0); 1699 // MakeWritable(f); 1700 // std::vector<long> v(myNumIndets); 1701 // myPPM.exponents(v, raw(pp)); 1702 1703 // WeylPoly::summand* f_smnd = AsWeylPoly(f).mySummands; 1704 // WeylPoly::summand* s = InitSummand(AllocSummand()); 1705 1706 // ordering(myPPM).ComputeOrdv(s->myOrdv, &v[0]); // &v[0] converts vector<T> to T* 1707 1708 // for (; f_smnd != nullptr ; f_smnd = f_smnd->myNext) 1709 // MulOrdv(f_smnd->myOrdv, f_smnd->myOrdv, s->myOrdv); 1710 // } 1711 1712 1713 // void RingWeylImpl::PushFront(RawPtr rawf, ConstRawPtr rawc, const std::vector<long>& expv) const 1714 // { 1715 // if (CoeffRing().IsZero(c)) return; 1716 // MakeWritable(f); 1717 // WeylPoly::summand* t = AllocSummand(); 1718 // InitSummand(t, c, expv); 1719 // t->myNext = f.WeylPolyPtr->mySummands; 1720 // if (f.WeylPolyPtr->mySummands == nullptr) f.WeylPolyPtr->myEnd = &t->myNext; 1721 // f.WeylPolyPtr->mySummands = t; 1722 // } 1723 1724 1725 // void RingWeylImpl::PushBack(RawPtr rawf, ConstRawPtr rawc, const std::vector<long>& expv) const 1726 // { 1727 // if (CoeffRing().IsZero(c)) return; 1728 // MakeWritable(f); 1729 // WeylPoly::summand* t = AllocSummand(); 1730 // InitSummand(t, c, expv); 1731 // *(f.WeylPolyPtr->myEnd) = t; 1732 // f.WeylPolyPtr->myEnd = &t->myNext; 1733 // } 1734 1735 1736 // void RingWeylImpl::PushFront(RawPtr rawx, WeylPoly::summand* t) const 1737 // { 1738 // MakeWritable(x); 1739 // WeylPoly& f = AsWeylPoly(x); 1740 // t->myNext = f.mySummands; 1741 // f.mySummands = t; 1742 // if (f.myEnd == &f.mySummands) f.myEnd = &t->myNext; 1743 // } 1744 1745 1746 // void RingWeylImpl::PushBack(RawPtr rawx, WeylPoly::summand* t) const 1747 // { 1748 // MakeWritable(x); 1749 // WeylPoly& f = AsWeylPoly(x); 1750 // *f.myEnd = t; 1751 // f.myEnd = &t->myNext; 1752 // } 1753 1754 1755 // void RingWeylImpl::RemoveSmnd(RawPtr rawx, WeylPoly::summand** prev_link) const 1756 // { 1757 // ASSERT(x.WeylPolyPtr->myRefCount == 1); 1758 // // MakeWritable(x); 1759 // WeylPoly::summand* tmp = *prev_link; 1760 // ASSERT(tmp != nullptr); 1761 // if (tmp->myNext==nullptr) // f.myEnd == &(tmp->myNext) 1762 // { 1763 // WeylPoly& f = AsWeylPoly(x); 1764 // f.myEnd = prev_link; 1765 // } 1766 // *prev_link = tmp->myNext; 1767 // tmp->myNext = nullptr; 1768 // DeleteSummands(tmp); 1769 // } 1770 1771 1772 // void RingWeylImpl::InsertSmnd(RawPtr rawx, WeylPoly::summand* s, WeylPoly::summand** prev_link) const 1773 // { 1774 // ASSERT(x.WeylPolyPtr->myRefCount == 1); 1775 // // MakeWritable(x); 1776 // WeylPoly& f = AsWeylPoly(x); 1777 // s->myNext = (*prev_link); 1778 // (*prev_link) = s; 1779 // if (f.myEnd == prev_link) f.myEnd = &(s->myNext); 1780 // } 1781 1782 1783 // PPMonoid::elem RingWeylImpl::LPP(ConstRawPtr rawf) const 1784 // { 1785 // ASSERT(!IsZero(f)); 1786 // return PPMonoid::elem(PPM(), PPMonoid::elem::FromOrdv, AsWeylPoly(f).mySummands->myOrdv); 1787 // } 1788 1789 1790 // bool RingWeylImpl::IsZeroAddLCs(RawPtr rawf, RawPtr rawg) const 1791 // { 1792 // ASSERT(!IsZero(f) && !IsZero(g)); 1793 // ASSERT( CmpLPP(f,g) == 0); 1794 // WeylPoly& F = AsWeylPoly(f); 1795 // WeylPoly& G = AsWeylPoly(g); 1796 // ASSERT(F.myRefCount==1 && G.myRefCount==1); 1797 // myCoeffRing.add(F.mySummands->myCoeff, F.mySummands->myCoeff, G.mySummands->myCoeff); 1798 // DeleteLM(g); 1799 // if (!myCoeffRing.IsZero(F.mySummands->myCoeff)) return false; 1800 // DeleteLM(f); 1801 // return true; 1802 // } 1803 1804 1805 // void RingWeylImpl::deriv(RawPtr rawdest, ConstRawPtr rawf, long var) const 1806 // { 1807 // const WeylPoly& F = AsWeylPoly(f); 1808 // RingElem ans(*this); 1809 // vector<long> expv(myNumIndets); 1810 // for (WeylPoly::summand* i=F.mySummands; i; i = i->myNext) 1811 // { 1812 // ordering(myPPM).ComputeExpv(&expv[0], i->myOrdv); 1813 // if (expv[var] == 0) continue; 1814 // RingElem c(myCoeffRing, expv[var]); 1815 // if (CoCoA::IsZero(c)) continue; 1816 // myCoeffRing.mul(raw(c), raw(c), i->myCoeff); 1817 // if (CoCoA::IsZero(c)) continue; 1818 // --expv[var]; 1819 // PushBack(raw(ans), raw(c), expv); 1820 // } 1821 // swap(raw(ans), dest); 1822 // } 1823 1824 1825 // void RingWeylImpl::negate(RawPtr rawlhs, ConstRawPtr rawx) const 1826 // { 1827 // if (lhs.WeylPolyPtr == x.WeylPolyPtr) 1828 // { 1829 // MakeWritable(lhs); 1830 // typedef WeylPoly::summand summand; 1831 // // std::clog << "-- negate"; output(std::clog, lhs); std::clog << endl; 1832 // for (summand* smnd = AsWeylPoly(lhs).mySummands; smnd!=nullptr; smnd=smnd->myNext ) 1833 // myCoeffRing.negate(smnd->myCoeff, smnd->myCoeff); 1834 // // std::clog << "-> negate"; output(std::clog, lhs); std::clog << endl; 1835 // } 1836 // else 1837 // std::clog << "RingWeylImpl::negate: not yet implemented" << endl; 1838 // } 1839 1840 1841 // void RingWeylImpl::add(RawPtr rawlhs, ConstRawPtr rawx, ConstRawPtr rawy) const 1842 // { 1843 // const AbstractRing& R = myCoeffRing; 1844 // RawValue ans; 1845 // // ans.WeylPolyPtr = static_cast<WeylPoly*>(myWeylPolyPool.alloc(sizeof(WeylPoly))); 1846 // init(ans); 1847 // /// WeylPoly& sum = AsWeylPoly(ans); 1848 // typedef WeylPoly::summand summand; 1849 // const summand* gterm = AsWeylPoly(x).mySummands; 1850 // const summand* hterm = AsWeylPoly(y).mySummands; 1851 // summand* SpareSummand = InitSummand(AllocSummand()); // just sets coeff = 0 1852 // while (gterm != nullptr && hterm != nullptr) 1853 // { 1854 // const int cmp = ordering(myPPM).CmpOrdvs(gterm->myOrdv, hterm->myOrdv); 1855 1856 // if (cmp < 0) 1857 // { 1858 // summand* hcopy = CopySummand(hterm); 1859 // PushBack(ans, hcopy); 1860 // hterm = hterm->myNext; 1861 // continue; 1862 // } 1863 1864 // if (cmp > 0) 1865 // { 1866 // summand* gcopy = CopySummand(gterm); 1867 // PushBack(ans, gcopy); 1868 // gterm = gterm->myNext; 1869 // continue; 1870 // } 1871 1872 // // Must have cmp == 0 here. 1873 // // The leading PPs are the same, so we must sum the coeffs. 1874 // R.add(SpareSummand->myCoeff, gterm->myCoeff, hterm->myCoeff); 1875 // if (!R.IsZero(SpareSummand->myCoeff)) 1876 // { 1877 // SetSummandMOrdv(SpareSummand, gterm); // set module posn and PP 1878 // PushBack(ans, SpareSummand); 1879 // SpareSummand = InitSummand(AllocSummand());// just sets coeff = 0 1880 // } 1881 // gterm = gterm->myNext; 1882 // hterm = hterm->myNext; 1883 // } 1884 // while (gterm != nullptr) 1885 // { 1886 // summand* gcopy = CopySummand(gterm); 1887 // PushBack(ans, gcopy); 1888 // gterm = gterm->myNext; 1889 // } 1890 // while (hterm != nullptr) 1891 // { 1892 // summand* hcopy = CopySummand(hterm); 1893 // PushBack(ans, hcopy); 1894 // hterm = hterm->myNext; 1895 // } 1896 // DeleteSummands(SpareSummand); // next ptr will be zero (set by InitSummand) 1897 // swap(lhs, ans); // really an assignment 1898 // kill(ans); 1899 // } 1900 1901 1902 // void RingWeylImpl::sub(RawPtr rawlhs, ConstRawPtr rawx, ConstRawPtr rawy) const 1903 // { 1904 // // This code copied from RingWeylImpl::add... 1905 1906 // const AbstractRing& R = myCoeffRing; 1907 // RawValue ans; 1908 // // ans.WeylPolyPtr = static_cast<WeylPoly*>(myWeylPolyPool.alloc(sizeof(WeylPoly))); 1909 // init(ans); 1910 // /// WeylPoly& sum = AsWeylPoly(ans); 1911 // typedef WeylPoly::summand summand; 1912 // const summand* gterm = AsWeylPoly(x).mySummands; 1913 // const summand* hterm = AsWeylPoly(y).mySummands; 1914 // summand* SpareSummand = InitSummand(AllocSummand()); // just sets coeff = 0 1915 // while (gterm != nullptr && hterm != nullptr) 1916 // { 1917 // const int ord = ordering(myPPM).CmpOrdvs(gterm->myOrdv, hterm->myOrdv); 1918 1919 // if (ord < 0) 1920 // { 1921 // summand* hcopy = CopySummand(hterm); 1922 // R.negate(hcopy->myCoeff, hcopy->myCoeff); 1923 // // sum.push_front(hcopy); 1924 // PushBack(ans, hcopy); 1925 // hterm = hterm->myNext; 1926 // continue; 1927 // } 1928 1929 // if (ord > 0) 1930 // { 1931 // summand* gcopy = CopySummand(gterm); 1932 // // sum.push_front(gcopy); 1933 // PushBack(ans, gcopy); 1934 // gterm = gterm->myNext; 1935 // continue; 1936 // } 1937 1938 // // The leading PPs are the same, so we must sum the coeffs. 1939 // R.sub(SpareSummand->myCoeff, gterm->myCoeff, hterm->myCoeff); 1940 // if (!R.IsZero(SpareSummand->myCoeff)) 1941 // { 1942 // SetSummandMOrdv(SpareSummand, gterm); // set module posn and PP 1943 // // sum.push_front(SpareSummand); 1944 // PushBack(ans, SpareSummand); 1945 // SpareSummand = InitSummand(AllocSummand());// just sets coeff = 0 1946 // } 1947 // gterm = gterm->myNext; 1948 // hterm = hterm->myNext; 1949 // } 1950 // while (gterm != nullptr) 1951 // { 1952 // summand* gcopy = CopySummand(gterm); 1953 // // sum.push_front(gcopy); 1954 // PushBack(ans, gcopy); 1955 // gterm = gterm->myNext; 1956 // } 1957 // while (hterm != nullptr) 1958 // { 1959 // summand* hcopy = CopySummand(hterm); 1960 // // sum.push_front(hcopy); 1961 // R.negate(hcopy->myCoeff, hcopy->myCoeff); 1962 // PushBack(ans, hcopy); 1963 // hterm = hterm->myNext; 1964 // } 1965 // DeleteSummands(SpareSummand); // next ptr will be zero (set by InitSummand) 1966 // /// NO LONGER NEEDED sum.reverse(); 1967 // swap(lhs, ans); // really an assignment 1968 // kill(ans); 1969 // } 1970 1971 1972 // void RingWeylImpl::mul(RawPtr rawlhs, ConstRawPtr rawx, ConstRawPtr rawy) const 1973 // { 1974 // // NO!! NOT COMMUTATIVE!! if (NumTerms(x) > NumTerms(y)) { mul(lhs, y, x); return; } 1975 1976 // //// std::clog << "MUL on "; output(std::clog, x); std::clog << " and "; output(std::clog, y); std::clog << endl; 1977 // RingElem ans(*this); 1978 1979 // for (const WeylPoly::summand* xterm = AsWeylPoly(x).mySummands; xterm; xterm = xterm->myNext) 1980 // AddMul(raw(ans), xterm, y, false);//????? 1981 1982 // swap(raw(ans), lhs); // really an assignment 1983 // } 1984 1985 1986 1987 1988 1989 // RCS header/log in the next few lines 1990 // $Header: /Volumes/Home_1/cocoa/cvs-repository/CoCoALib-0.99/src/AlgebraicCore/RingWeyl.C,v 1.67 2020/06/17 15:49:27 abbott Exp $ 1991 // $Log: RingWeyl.C,v $ 1992 // Revision 1.67 2020/06/17 15:49:27 abbott 1993 // Summary: Changed CoCoA_ERROR into CoCoA_THROW_ERROR 1994 // 1995 // Revision 1.66 2020/02/12 09:01:47 bigatti 1996 // -- changed myTestIsMaximal etc to return void (and consequences) 1997 // 1998 // Revision 1.65 2020/02/11 16:56:42 abbott 1999 // Summary: Corrected last update (see redmine 969) 2000 // 2001 // Revision 1.64 2020/02/11 16:12:19 abbott 2002 // Summary: Added some checks for bad ostream (even to mem fns myOutput); see redmine 969 2003 // 2004 // Revision 1.63 2019/10/15 11:54:08 abbott 2005 // Summary: Changed 0 into nullptr (where appropriate) 2006 // 2007 // Revision 1.62 2019/03/04 11:11:08 abbott 2008 // Summary: Changed auto_ptr into unique_ptr 2009 // 2010 // Revision 1.61 2018/06/27 08:50:39 abbott 2011 // Summary: Revised to work with new CpuTimeLimit 2012 // 2013 // Revision 1.60 2018/05/25 09:24:46 abbott 2014 // Summary: Major redesign of CpuTimeLimit (many consequences) 2015 // 2016 // Revision 1.59 2018/05/18 12:22:30 bigatti 2017 // -- renamed IntOperations --> BigIntOps 2018 // 2019 // Revision 1.58 2018/05/17 15:59:37 bigatti 2020 // -- renamed MatrixOperations --> MatrixOps 2021 // -- sorted #includes 2022 // 2023 // Revision 1.57 2018/04/18 14:32:03 abbott 2024 // Summary: Added missing return stmts in some NYI fns 2025 // 2026 // Revision 1.56 2018/03/29 09:36:40 bigatti 2027 // -- added member functions myTestIsRadical, myTestIsPrimary and flags 2028 // 2029 // Revision 1.55 2018/03/20 11:38:08 bigatti 2030 // -- changed iAm***Test --> myTestIs***; and it returns bool 2031 // 2032 // Revision 1.54 2017/04/18 09:55:03 bigatti 2033 // -- removed StatLevel from GBasis 2034 // 2035 // Revision 1.53 2016/11/07 12:22:21 bigatti 2036 // -- changed myGBasisIsValid into IhaveGBasis 2037 // 2038 // Revision 1.52 2015/12/08 14:06:22 abbott 2039 // Summary: Added include for UIBC (otherwise does not compile) 2040 // 2041 // Revision 1.51 2015/12/01 13:34:44 abbott 2042 // Summary: Changed arg order in ElimMat and HomogElimMat; doc is out-of-date!! 2043 // 2044 // Revision 1.50 2015/11/30 21:53:55 abbott 2045 // Summary: Major update to matrices for orderings (not yet complete, some tests fail) 2046 // 2047 // Revision 1.49 2015/07/23 11:58:03 bigatti 2048 // -- removed obsolete doxygen header 2049 // 2050 // Revision 1.48 2015/04/24 15:40:58 bigatti 2051 // -- renamed: myAddMul --> myAddMulLM 2052 // -- renamed: myMoveLM --> myMoveLMToFront 2053 // -- new myMoveLMToBack (used in ReductionCog --> bug in test-TmpMorseGraph??) 2054 // 2055 // Revision 1.47 2014/07/28 15:51:31 abbott 2056 // Summary: Redesign: ringhoms no longer cached in rings (caused ref count trouble) 2057 // Author: JAA 2058 // 2059 // Revision 1.46 2014/07/11 16:00:39 bigatti 2060 // -- commented out myOutputSelf(OpenMathOutput& OMOut) 2061 // 2062 // Revision 1.45 2014/07/11 15:49:25 bigatti 2063 // -- commented out myOutputSelf (default impl) and added myImplDetails 2064 // (is thins correct for Weyl?) 2065 // 2066 // Revision 1.44 2014/07/09 13:03:25 abbott 2067 // Summary: Removed AsSparsePolyRing from commented out code 2068 // Author: JAA 2069 // 2070 // Revision 1.43 2014/05/14 15:57:15 bigatti 2071 // -- added "using" for clang with superpedantic flag 2072 // 2073 // Revision 1.42 2014/05/06 13:20:41 abbott 2074 // Summary: Changed names (my)MaxExponents into (my)Deg 2075 // Author: JAA 2076 // 2077 // Revision 1.41 2014/04/11 15:44:27 abbott 2078 // Summary: Renamed MatrixArith to MatrixOperations (in includes) 2079 // Author: JAA 2080 // 2081 // Revision 1.40 2014/04/02 10:57:46 abbott 2082 // Summary: Revised design of IamIntegralDomain3 2083 // Author: JAA 2084 // 2085 // Revision 1.39 2013/06/28 17:03:51 abbott 2086 // Modified semantics of IdealBase::myDivMod; 2087 // it now returns a boolean. 2088 // Several consequential changes. 2089 // 2090 // Revision 1.38 2013/02/14 17:34:35 bigatti 2091 // -- cleaned up code for elimination matrices 2092 // 2093 // Revision 1.37 2012/10/24 12:20:43 abbott 2094 // Changed return type of myLC. 2095 // 2096 // Revision 1.36 2012/10/17 09:40:16 abbott 2097 // Replaced RefRingElem by RingElem& 2098 // (plus a few consequential changes) 2099 // 2100 // Revision 1.35 2012/10/11 14:29:41 abbott 2101 // Rewrote myMulByPP and myReductionStep so that they no longer use RefRingElem. 2102 // 2103 // Revision 1.34 2012/05/28 09:18:20 abbott 2104 // Created IntOperations which gathers together all operations on 2105 // integers (both big and small). Many consequential changes. 2106 // 2107 // Revision 1.33 2012/05/24 14:49:23 bigatti 2108 // -- changed symbol "index" into "subscripts" 2109 // 2110 // Revision 1.32 2012/05/22 10:02:37 abbott 2111 // Removed IsGCDDomain; substituted by IsTrueGCDDomain. 2112 // Added IsFractionFieldOfGCDDomain. 2113 // 2114 // Revision 1.31 2012/01/26 16:47:00 bigatti 2115 // -- changed back_inserter into insert 2116 // 2117 // Revision 1.30 2011/11/09 14:29:37 bigatti 2118 // -- renamed MachineInteger --> MachineInt 2119 // 2120 // Revision 1.29 2011/08/14 15:52:16 abbott 2121 // Changed ZZ into BigInt (phase 1: just the library sources). 2122 // 2123 // Revision 1.28 2011/06/23 16:04:47 abbott 2124 // Added IamExact mem fn for rings. 2125 // Added myRecvTwinFloat mem fn for rings. 2126 // Added first imple of RingHom from RingTwinFloat to other rings. 2127 // 2128 // Revision 1.27 2011/03/10 16:39:33 abbott 2129 // Replaced (very many) size_t by long in function interfaces (for rings, 2130 // PPMonoids and modules). Also replaced most size_t inside fn defns. 2131 // 2132 // Revision 1.26 2011/03/08 17:58:29 bigatti 2133 // -- changed: ElimIndets is now of long instead of size_t 2134 // 2135 // Revision 1.25 2010/10/08 11:39:53 abbott 2136 // Renamed DistrMPoly to DistrMPolyClean. 2137 // 2138 // Revision 1.24 2010/10/06 14:10:24 abbott 2139 // Added increments to the ref count in ring and PPMonoid ctors to make 2140 // them exception safe. 2141 // 2142 // Revision 1.23 2010/07/27 07:37:13 bigatti 2143 // -- new class GBCriteria, simplified GReductor ctor 2144 // 2145 // Revision 1.22 2010/06/10 08:00:02 bigatti 2146 // -- fixed naming conventions 2147 // 2148 // Revision 1.21 2010/05/14 09:53:09 bigatti 2149 // -- removed empty ctor for SugarDegree 2150 // -- added marker for SugarDegree(uninitialized) 2151 // -- SugarDegree for GBasis input is initialized by myPrepareGBasis 2152 // 2153 // Revision 1.20 2009/10/02 13:47:07 bigatti 2154 // -- myDivByCoeff now returns bool 2155 // -- unique implementation of myDiv in PolyRing.C 2156 // 2157 // Revision 1.19 2009/09/28 08:39:01 bigatti 2158 // -- fixed (debugging) parameter by M.Caboara 2159 // 2160 // Revision 1.18 2009/09/25 13:02:43 bigatti 2161 // -- just a comment 2162 // 2163 // Revision 1.17 2009/09/22 13:35:55 bigatti 2164 // -- following coding conventions in function names Matrix --> Mat 2165 // -- forced all matrices to be over RingZ 2166 // 2167 // Revision 1.16 2009/01/30 13:41:50 bigatti 2168 // -- enum instead of bool arguments 2169 // 2170 // Revision 1.15 2008/12/17 12:11:52 abbott 2171 // Changed type from long to MachineInt in operations which use a machine integer 2172 // in place of a RingElem. The change is "superficial" but affects many files. 2173 // 2174 // Revision 1.14 2008/11/18 16:34:46 bigatti 2175 // -- removed debugging printout 2176 // 2177 // Revision 1.13 2008/11/18 15:20:10 bigatti 2178 // -- added const to myGBasis return value 2179 // -- added myIdealCtor to RingWeyl for proper inheritance 2180 // 2181 // Revision 1.12 2008/09/19 13:33:42 bigatti 2182 // -- added: Sat algorithm (M.Caboara) 2183 // 2184 // Revision 1.11 2008/04/21 11:23:11 abbott 2185 // Separated functions dealing with matrices and PPOrderings into a new file. 2186 // Added matrix norms, and completed adjoint. 2187 // 2188 // Revision 1.10 2008/04/10 15:13:21 bigatti 2189 // -- added myPushBack/Front(RawPtr, ConstRawPtr, PPMonoidElemConstRawPtr) 2190 // 2191 // Revision 1.9 2007/12/05 11:06:24 bigatti 2192 // -- changed "size_t StdDeg/myStdDeg(f)" into "long" (and related functions) 2193 // -- changed "log/myLog(f, i)" into "MaxExponent/myMaxExponent(f, i)" 2194 // -- fixed bug in "IsOne(ideal)" in SparsePolyRing.C 2195 // 2196 // Revision 1.8 2007/12/04 14:27:06 bigatti 2197 // -- changed "log(pp, i)" into "exponent(pp, i)" 2198 // 2199 // Revision 1.7 2007/10/30 17:14:07 abbott 2200 // Changed licence from GPL-2 only to GPL-3 or later. 2201 // New version for such an important change. 2202 // 2203 // Revision 1.6 2007/05/31 16:01:45 bigatti 2204 // -- default implementation for IamField, myCharacteristic in PolyRing 2205 // -- added !IsCommutative in test 2206 // 2207 // Revision 1.4 2007/05/22 22:45:14 abbott 2208 // Changed fn name IsUnit to IsInvertible. 2209 // 2210 // Revision 1.3 2007/05/21 12:45:13 abbott 2211 // Modified a pointless comment. 2212 // 2213 // Revision 1.2 2007/03/09 18:56:56 bigatti 2214 // -- added Tmp prefix to Groebner related files 2215 // 2216 // Revision 1.1.1.1 2007/03/09 15:16:11 abbott 2217 // Imported files 2218 // 2219 // Revision 1.26 2007/03/08 18:22:29 cocoa 2220 // Just whitespace cleaning. 2221 // 2222 // Revision 1.25 2007/03/08 17:43:10 cocoa 2223 // Swapped order of args to the NewPPMonoid pseudo ctors. 2224 // 2225 // Revision 1.24 2007/03/08 16:55:06 cocoa 2226 // Changed name of "range" function to "SymbolRange". 2227 // 2228 // Revision 1.23 2007/03/08 11:07:12 cocoa 2229 // Made pseudo ctors for polynomial rings more uniform. This allowed me to 2230 // remove an include of CoCoA/symbol.H from the RingDistrM*.H files, but then 2231 // I had to put the include in several .C files. 2232 // 2233 // Revision 1.22 2007/03/05 21:06:07 cocoa 2234 // New names for homomorphism pseudo-ctors: removed the "New" prefix. 2235 // 2236 // Revision 1.21 2007/02/10 18:44:03 cocoa 2237 // Added "const" twice to each test and example. 2238 // Eliminated dependency on io.H in several files. 2239 // Improved BuildInfo, and added an example about how to use it. 2240 // Some other minor cleaning. 2241 // 2242 // Revision 1.20 2007/01/20 14:07:25 bigatti 2243 // -- moved code for homomorphism into common implementation in SparsePolyRing 2244 // 2245 // Revision 1.19 2007/01/15 16:15:26 cocoa 2246 // -- added prefix "raw" to RawPtr arguments names 2247 // -- changed rhs into rawx, n, or N 2248 // 2249 // Revision 1.18 2007/01/13 14:14:34 cocoa 2250 // Overhaul of RingHom code: it nows uses SmartPtrIRC, and printing is more logical. 2251 // Have not yet updated the documentation. 2252 // 2253 // Revision 1.17 2006/12/21 13:48:32 cocoa 2254 // Made all increment/decrement calls prefix (except where the must be postfix). 2255 // 2256 // Revision 1.16 2006/12/07 17:23:46 cocoa 2257 // -- for compilation with _Wextra: commented out names of unused arguments 2258 // 2259 // Revision 1.15 2006/12/07 12:17:15 cocoa 2260 // -- style: RawPtr args are now called "raw.." 2261 // 2262 // Revision 1.14 2006/11/24 17:06:10 cocoa 2263 // -- reorganized includes of header files 2264 // 2265 // Revision 1.13 2006/11/21 18:09:23 cocoa 2266 // -- added myIsMonomial 2267 // -- implemented myIsOne, myIsMinusOne, myIsConstant, myIsIndet in SparsePolyRing 2268 // -- removed the 4 functions from DistrMPoly(..) and RingDistrMPoly(..) 2269 // -- changed all names of RawPtr arguments into "raw(..)" 2270 // 2271 // Revision 1.12 2006/11/17 12:04:29 cocoa 2272 // -- added (obvious) constructor for RingWeylImpl::IdealImpl 2273 // 2274 // Revision 1.11 2006/11/14 17:36:49 cocoa 2275 // -- fixed implementation for ideal in RingWeyl 2276 // 2277 // Revision 1.10 2006/11/09 17:46:58 cocoa 2278 // -- version 0.9712: 2279 // -- IdealImpl moved to SparsePolyRing from concrete rings 2280 // -- PolyList in GTypes is now vector<RingElem> (was list) 2281 // -- "my" coding convention applied to DistrMPoly 2282 // 2283 // Revision 1.9 2006/11/08 16:21:59 cocoa 2284 // Structural cleaning of RingHom; many consequential changes. 2285 // 2286 // Revision 1.8 2006/11/02 13:25:43 cocoa 2287 // Simplification of header files: the OpenMath classes have been renamed. 2288 // Many minor consequential changes. 2289 // 2290 // Revision 1.7 2006/10/16 23:18:59 cocoa 2291 // Corrected use of std::swap and various special swap functions. 2292 // Improved myApply memfn for homs of RingDistrMPolyInlPP. 2293 // 2294 // Revision 1.6 2006/10/06 14:04:14 cocoa 2295 // Corrected position of #ifndef in header files. 2296 // Separated CoCoA_ASSERT into assert.H from config.H; 2297 // many minor consequential changes (have to #include assert.H). 2298 // A little tidying of #include directives (esp. in Max's code). 2299 // 2300 // Revision 1.5 2006/08/17 09:41:33 cocoa 2301 // -- changed: I think it works again 2302 // 2303 // Revision 1.4 2006/08/07 21:23:25 cocoa 2304 // Removed almost all publicly visible references to SmallExponent_t; 2305 // changed to long in all PPMonoid functions and SparsePolyRing functions. 2306 // DivMask remains to sorted out. 2307 // 2308 // Revision 1.3 2006/07/20 14:11:27 cocoa 2309 // -- more stable version (more similar to RingDistrMPoly) 2310 // 2311 // Revision 1.1.1.1 2006/05/30 11:39:37 cocoa 2312 // Imported files 2313 // 2314 // Revision 1.6 2006/04/27 13:45:30 cocoa 2315 // Changed name of NewIdentityRingHom to NewIdentityHom. 2316 // Changed name of member functions which print out their own object 2317 // into myOutputSelf (to distinguish from "transitive" myOutput fns). 2318 // 2319 // Revision 1.5 2006/03/21 09:43:13 cocoa 2320 // Changed names of some member fns of ideals (dealing with setting and testing 2321 // the flags for primeness and maximality). Hope icc will complain less now. 2322 // 2323 // Revision 1.4 2006/03/14 15:01:49 cocoa 2324 // Improved the implementation of ring member fns for computing powers. 2325 // Should keep Intel C++ compiler quieter too. 2326 // 2327 // Revision 1.3 2006/03/12 21:28:33 cocoa 2328 // Major check in after many changes 2329 // 2330 // Revision 1.2 2006/02/20 22:41:19 cocoa 2331 // All forms of the log function for power products now return SmallExponent_t 2332 // (instead of int). exponents now resizes the vector rather than requiring 2333 // the user to pass in the correct size. 2334 // 2335 // Revision 1.1.1.1 2005/10/17 10:46:54 cocoa 2336 // Imported files 2337 // 2338 // Revision 1.1.1.1 2005/05/03 15:47:31 cocoa 2339 // Imported files 2340 // 2341 // Revision 1.3 2005/04/20 15:40:48 cocoa 2342 // Major change: modified the standard way errors are to be signalled 2343 // (now via a macro which records filename and line number). Updated 2344 // documentation in error.txt accordingly. 2345 // 2346 // Improved the documentation in matrix.txt (still more work to be done). 2347 // 2348 // Revision 1.2 2005/04/19 14:06:03 cocoa 2349 // Added GPL and GFDL licence stuff. 2350 // 2351 // Revision 1.1.1.1 2005/01/27 15:12:13 cocoa 2352 // Imported files 2353 // 2354 // Revision 1.16 2004/11/18 18:33:41 cocoa 2355 // Now every ring know its own "one" element (as well as "zero"). 2356 // Several consequential changes. 2357 // 2358 // Revision 1.15 2004/11/11 14:28:49 cocoa 2359 // -- minor changes for doxygen 2360 // -- change: cout --> std::clog 2361 // 2362 // Revision 1.14 2004/11/09 16:30:51 cocoa 2363 // Removed references to cout. 2364 // 2365 // Revision 1.13 2004/11/04 18:47:43 cocoa 2366 // (1) Ring member functions which previously expected mpz_t args 2367 // now expect ZZ args. Numerous minor consequential changes. 2368 // (2) Renamed function which gives access to the mpz_t value inside 2369 // a ZZ object: previously was raw(...), now is mpzref(...). 2370 // Plenty of calls had to be altered. 2371 // 2372 // Revision 1.12 2004/10/21 17:16:37 cocoa 2373 // Fairly major change: new OrdvArith namspace with various members, 2374 // new global typedef SmallExponent_t (defined in config.H). 2375 // 2376 // Revision 1.11 2004/07/27 16:03:39 cocoa 2377 // Added IsCommutative test and IamCommutative member function 2378 // to all rings. Tidied geobuckets a little. 2379 // 2380 // Revision 1.10 2004/07/20 15:37:08 cocoa 2381 // Minor fix for some errors which slipped through the net... 2382 // 2383 // Revision 1.9 2004/07/20 15:04:06 cocoa 2384 // The next step in the new "ring element" conversion process: 2385 // handling the case of creating a "const RefRingElem" object 2386 // (since C++ refuses to do this properly itself). 2387 // 2388 // Revision 1.8 2004/07/20 09:21:33 cocoa 2389 // -- fewer Stats are printed 2390 // 2391 // Revision 1.7 2004/06/17 12:44:58 cocoa 2392 // -- applied new coding conventions: compiles and seems to work 2393 // 2394 // Revision 1.6 2004/05/27 16:14:02 cocoa 2395 // Minor revision for new coding conventions. 2396 // 2397 // Revision 1.5 2004/05/24 15:52:13 cocoa 2398 // Major update: 2399 // new error mechanism 2400 // many fixes 2401 // RingHoms almost work now 2402 // RingFloat much improved 2403 // 2404 // Revision 1.4 2003/10/17 10:51:06 cocoa 2405 // Major cleaning, and new naming convention. 2406 // 2407 // Revision 1.3 2003/10/09 13:32:16 cocoa 2408 // A few glitches which slipped through the first major merge. 2409 // 2410 // Revision 1.2 2003/10/09 12:16:38 cocoa 2411 // New coding convention for rings. 2412 // 2413 // Revision 1.3 2003/06/23 16:58:43 abbott 2414 // Minor cleaning prior to public release. 2415 // Just consequential changes. 2416 // 2417 // Revision 1.2 2003/05/30 15:24:01 abbott 2418 // Changed beyond all recognition -- completely new version. 2419 // 2420 // Revision 1.1 2003/05/02 13:08:03 abbott 2421 // Initial revision 2422 // 2423 // 2424