1 // Copyright (c) 2008 Anna Bigatti 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 #include "CoCoA/TmpPPVector.H" 19 #include "CoCoA/SparsePolyOps-RingElem.H" 20 //#include "CoCoA/SparsePolyRing.H" // from SparsePolyOps-RingElem.H 21 #include "CoCoA/VectorOps.H" 22 //#include "CoCoA/PPWithMask.H" -- included by CoCoA/TmpPPVector.H 23 //#include "CoCoA/DivMask.H" -- included by CoCoA/PPWithMask.H 24 //#include "CoCoA/PPMonoid.H" -- included by CoCoA/PPWithMask.H 25 //#include "CoCoA/ring.H" -- included by CoCoA/SparsePolyRing.H 26 27 #include <algorithm> 28 using std::sort; 29 //#include <vector> 30 using std::vector; 31 32 // ints MFI_Occurrences, MFI_Indets; 33 // int Ints120[121]; 34 // MixedPPVector MTL120[121]; 35 // PPVector GlobalSplitterTList; 36 37 /**********************************************************/ 38 namespace CoCoA 39 { 40 myInterreduce()41 void PPVector::myInterreduce() 42 { myInterreduceSort(); } // AMB: faster than without sorting 43 44 PushBack(PPVector & PPs,ConstRefPPMonoidElem pp)45 void PushBack(PPVector& PPs, ConstRefPPMonoidElem pp) 46 { 47 if (PPM(PPs) != owner(pp)) 48 { 49 vector<long> e(NumIndets(PPM(PPs))); 50 exponents(e, pp); 51 PPs.myPushBack(PPMonoidElem(PPM(PPs), e)); 52 } 53 else 54 PPs.myPushBack(pp); 55 } 56 57 PushBack(PPVector & PPs,const PPWithMask & pp)58 void PushBack(PPVector& PPs, const PPWithMask& pp) 59 { 60 CoCoA_ASSERT(PPM(PPs) == owner(PP(pp))); 61 PPs.myPushBack(pp); 62 } 63 64 PushBackPopBack(PPVector & ToPPs,PPVector & FromPPs)65 void PushBackPopBack(PPVector& ToPPs, PPVector& FromPPs) 66 { ToPPs.myPushBackPopBack(FromPPs); } 67 68 IsDivisible(const PPWithMask & pp,const PPVector & ByL)69 bool IsDivisible(const PPWithMask& pp, const PPVector& ByL) 70 { 71 // CoCoA_ASSERT(compatible) 72 return ByL.myDivides(pp); 73 } 74 75 IsDivisible(ConstRefPPMonoidElem pp,const PPVector & ByL)76 bool IsDivisible(ConstRefPPMonoidElem pp, const PPVector& ByL) 77 { 78 // CoCoA_ASSERT(compatible) 79 return ByL.myDivides(pp); 80 } 81 82 swap(PPVector & PPs1,PPVector & PPs2)83 void swap(PPVector& PPs1, PPVector& PPs2) 84 { 85 // CoCoA_ASSERT(compatible) 86 // PPs1.mySwap(PPs2); 87 swap(PPs1.myVec, PPs2.myVec); 88 } 89 90 interreduce(PPVector & PPs)91 void interreduce(PPVector& PPs) 92 { PPs.myInterreduce(); } 93 94 InterreduceSort(PPVector & PPs)95 void InterreduceSort(PPVector& PPs) 96 { PPs.myInterreduceSort(); } 97 98 lcms(PPVector & PPs,const PPVector & PPs1,const PPVector & PPs2)99 void lcms(PPVector& PPs, const PPVector& PPs1, const PPVector& PPs2) 100 { 101 // CoCoA_ASSERT(compatible) 102 PPs.myLcms(PPs1, PPs2); 103 } 104 105 convert(std::vector<RingElem> & v,ring P,const PPVector & PPs)106 void convert(std::vector<RingElem>& v, ring P, const PPVector& PPs) 107 { 108 // CoCoA_ASSERT(compatible) 109 PPs.myConvert(v, P); 110 } 111 112 convert(PPVector & PPs,const std::vector<RingElem> & v)113 void convert(PPVector& PPs, const std::vector<RingElem>& v) 114 { 115 PPVector tmpPPs(PPM(PPs), DMR(PPs)); 116 CoCoA_ASSERT(AreMonomials(v)); 117 for (const auto& f: v) PushBack(tmpPPs, LPP(f)); 118 swap(tmpPPs, PPs); 119 } 120 121 122 // const PPWithMask& PPVector::operator[](MachineInt i) const 123 // { 124 // // CoCoA_ASSERT(range) ??? anna 125 // // return myAt(AsUnsignedLong(i)); ??? 126 // return myVec[AsUnsignedLong(i)]; 127 // } 128 129 130 std::ostream& operator<<(std::ostream& out, PPVector PPs) 131 { 132 if (!out) return out; // short-cut for bad ostreams 133 out << PPs.myVec; 134 return out; 135 } 136 137 138 // class PPVector 139 // { 140 // public: 141 // PPVector(PPMonoid PPM, DivMaskRule DMR); 142 // virtual ~PPVector() {}; 143 // long myNumIndetsEnv() const { return NumIndets(myPPM); } 144 // long mySize() const { return len(myVec); } 145 // void myReserve(long n); 146 // bool IamEmpty() const; 147 // PPVectorElem myDefaultElem() const; 148 // const PPVectorElem myComp(unsigned long i) const; 149 // void myPushBack(ConstRefPPMonoidElem pp); 150 // void myPushBack(const PPWithMask& pp); 151 // void myPushBackPopBack(PPVector& FromPPs); 152 // void mySort(); 153 // bool myDivides(const PPWithMask& pp) const; 154 // bool myDivides(ConstRefPPMonoidElem pp) const; 155 // void myInterreduceSort(); 156 // void myLcms(const PPVector& PPs1, const PPVector& PPs2); 157 // void myConvert(std::vector<RingElem>& v, SparsePolyRing P) const; 158 // private: 159 // const PPWithMask& myNth(int n) const; 160 // private: // member fields 161 // PPMonoid myPPM; 162 // DivMaskRule myDMR; 163 // vector<PPWithMask> myVec; 164 // }; 165 166 167 //-------- auxiliary PPWithMaskLessThan(const PPWithMask & f,const PPWithMask & g)168 bool PPWithMaskLessThan(const PPWithMask& f, const PPWithMask& g) 169 { 170 return (owner(PP(f)))->myCmp(raw(PP(f)), raw(PP(g))) < 0; 171 } 172 173 //-------- auxiliary end 174 175 PPVector(PPMonoid PPM,DivMaskRule DMR)176 PPVector::PPVector(PPMonoid PPM, DivMaskRule DMR): 177 myPPM(PPM), myDMR(DMR) 178 { } 179 180 PPVector(PPMonoid PPM,DivMaskRule DMR,const std::vector<RingElem> & v)181 PPVector::PPVector(PPMonoid PPM, DivMaskRule DMR, const std::vector<RingElem>& v): 182 myPPM(PPM), myDMR(DMR) 183 { 184 convert(*this, v); 185 } 186 187 PPVector(const PPVector & PPs)188 PPVector::PPVector(const PPVector& PPs): 189 myPPM(PPs.myPPM), myDMR(PPs.myDMR) 190 { 191 myVec = PPs.myVec; 192 } 193 194 195 PPVector& PPVector::operator=(const PPVector& rhs) 196 { 197 if (this == &rhs) return *this; 198 PPVector copy(rhs); 199 swap(*this, copy); 200 return *this; 201 } 202 203 204 // PPVector* PPVector::myZeroClone() const 205 // { return new PPVector(myPPM, myDMR); } 206 207 myReserve(long n)208 void PPVector::myReserve(long n) // negative n check??? 209 { myVec.reserve(n); } 210 211 myClear()212 void PPVector::myClear() // negative n check??? 213 { myVec.clear(); } 214 215 IamEmpty()216 bool PPVector::IamEmpty() const 217 { return myVec.empty(); } 218 219 220 // const PPWithMask PPVector::myAt(unsigned long i) const 221 // { 222 // // anna: make all checks 223 // return myVec[i]; 224 // } 225 226 myPushBack(ConstRefPPMonoidElem pp)227 void PPVector::myPushBack(ConstRefPPMonoidElem pp) 228 { myVec.push_back(PPWithMask(pp, myDMR)); } 229 230 myPushBack(const PPWithMask & pp)231 void PPVector::myPushBack(const PPWithMask& pp) 232 { myVec.push_back(pp); } 233 234 myPushBackPopBack(PPVector & FromPPs)235 void PPVector::myPushBackPopBack(PPVector& FromPPs) 236 { 237 myPushBack((FromPPs.myVec).back()); 238 FromPPs.myVec.pop_back(); 239 } 240 241 // NOT exception clean mySort()242 void PPVector::mySort() 243 { 244 sort(myVec.begin(), myVec.end(), PPWithMaskLessThan); 245 } 246 247 myDivides(const PPWithMask & pp)248 bool PPVector::myDivides(const PPWithMask& pp) const 249 { 250 // const long n=myLen(); 251 // for (long i=0; i < n; ++i) 252 // if ( IsDivisibleFast(pp, myVec[i]) ) return true; 253 for (const auto& pp_i: myVec) 254 if ( IsDivisibleFast(pp, pp_i) ) return true; 255 return false; 256 } 257 258 myDivides(ConstRefPPMonoidElem pp)259 bool PPVector::myDivides(ConstRefPPMonoidElem pp) const 260 { 261 return myDivides(PPWithMask(pp, myDMR)); 262 } 263 264 265 // NOT exception clean myInterreduceSort()266 void PPVector::myInterreduceSort() 267 { 268 if (IamEmpty()) return; 269 mySort(); // ascending 270 vector<PPWithMask> v; 271 v.reserve(myLen()); 272 swap(myVec, v); //// ANNA: swap afterwards for exception safety? 273 for (const auto& pp: v) 274 if (!myDivides(pp)) myPushBack(pp); 275 } 276 277 278 // anna: careful with push_back and multiple copies of lcms myLcms(const PPVector & PPs1,const PPVector & PPs2)279 void PPVector::myLcms(const PPVector& PPs1, const PPVector& PPs2) 280 { 281 // anna: this calls the destructor too many times (2*mylcm calls) 282 // CoCoA_ASSERT(compatible) 283 if (IsEmpty(PPs1) || IsEmpty(PPs2)) 284 { 285 myVec.clear(); 286 return; 287 } 288 PPMonoidElem pp(myPPM); 289 vector<PPWithMask> res; 290 res.reserve(len(PPs1)*len(PPs2)); 291 for (const auto& pp1: PPs1.myVec) 292 for (const auto& pp2: PPs2.myVec) 293 { 294 myPPM->myLcm(raw(pp), raw(PP(pp1)), raw(PP(pp2))); 295 res.push_back(PPWithMask(pp, myDMR)); 296 } 297 swap(myVec, res); 298 } 299 300 myAlexanderDual(const PPVector & PPs)301 void PPVector::myAlexanderDual(const PPVector& PPs) 302 { 303 CoCoA_ASSERT(PPM(PPs) == myPPM); 304 CoCoA_ASSERT(DMR(PPs) == myDMR); 305 PPVector res(PPM(PPs), DMR(PPs)); 306 PPVector tmp(PPM(PPs), DMR(PPs)); 307 PushBack(res, one(PPM(PPs))); 308 // for (long i=0; i<len(PPs); ++i) 309 for (const auto& ppwm: PPs.myVec) 310 { 311 tmp.mySupport(PP(ppwm)); 312 res.myLcms(res, tmp); 313 res.myInterreduceSort(); 314 } 315 swap(myVec, res.myVec); 316 } 317 318 mySupport(ConstRefPPMonoidElem pp)319 void PPVector::mySupport(ConstRefPPMonoidElem pp) 320 { 321 CoCoA_ASSERT(owner(pp) == myPPM); 322 const PPMonoid& PPM = myPPM; 323 vector<PPWithMask> res; 324 for (long i=0 ; i < NumIndets(PPM) ; ++i ) 325 if ( exponent(pp, i) != 0 ) 326 res.push_back(PPWithMask(indet(PPM, i), myDMR)); 327 swap(myVec, res); 328 } 329 330 myConvert(std::vector<RingElem> & v,SparsePolyRing P)331 void PPVector::myConvert(std::vector<RingElem>& v, SparsePolyRing P) const 332 { 333 vector<RingElem> tmp; 334 if (myPPM == PPM(P)) 335 for (const auto& ppwm: myVec) 336 tmp.push_back(monomial(P, PP(ppwm))); 337 else 338 { 339 vector<long> e(NumIndets(P)); 340 for (const auto& ppwm: myVec) 341 { 342 exponents(e, PP(ppwm)); 343 tmp.push_back(monomial(P, e)); 344 } 345 } 346 swap(v, tmp); 347 } 348 349 350 /******************** list of terms ********************/ 351 352 /*************************************************************/ 353 354 //#define Random(len)( 1 ) 355 356 // void InsInTList(PPVector TL, eterm t, int *NewLen) 357 // { 358 // ints OccInd = Indets(t); 359 // if ( IntsGetLen(OccInd) == 1 ) 360 // { 361 // InsInSPList (OccInd[1], eterm_get_nth(t,OccInd[1]), SPList(TL) ); 362 // eterm_free(t); 363 // } 364 // else 365 // MTLPutLast (MTList(TL), *NewLen, t); 366 // } 367 368 369 370 /******** Simple Power List ********/ 371 372 // void InsInSPList(size_t Index, unsigned Exp, eterm SPL) 373 // { 374 // size_t CurrentExp =SPL[Index]; 375 376 // if (CurrentExp == 0) 377 // { 378 // eterm_put_non0_nth(SPL, Index, Exp); 379 // IntsPutLast(Indets(SPL), Index); 380 // } 381 // else if (CurrentExp > Exp) 382 // eterm_put_non0_nth(SPL, Index, Exp); 383 // } 384 385 // #define SPLDividesTerm(SPL,t) sp_SmallMult(t,SPL) 386 387 // bool TListSimplePowersDivide(PPVector PPs, eterm t) 388 // { 389 // eterm SPL =SPList(PPs); 390 // int i; 391 392 // for ( i=eterm_get_indetsNo(SPL) ; i>0 ; --i ) 393 // if ( SPL[i]!=0 && SPL[i]<=t[i] ) return TRUE; 394 // return FALSE; 395 // } 396 397 // /************************ SPLIT ************************/ 398 399 // void DelLinearSimplePowers(PPVector PPs, int *LPSNo) 400 // { 401 // eterm SPL = SPList(PPs); 402 // ints OccInd = Indets(SPL); 403 // size_t n = IntsGetLen(OccInd), sum = 0, index; 404 405 // for ( ; n>0 ; --n ) 406 // if ( SPL[(index=OccInd[n])] == 1 ) 407 // { 408 // ++sum; 409 // eterm_put0_nth(SPL, index); 410 // IntsMoveLastToNth(OccInd, n); 411 // } 412 // *LPSNo = sum; 413 // } 414 415 // eterm GetCoprimeSPList(PPVector PPs) 416 // { 417 // int MTLLen =MTListLen(PPs), *OccInd; 418 // MixedPPVector MTL =MTList(PPs); 419 // size_t i, n =MTLLen, index; 420 // eterm res, theSPL; 421 // shortbitset CoprimeSupp; 422 423 // if ( IntsGetLen(Indets(SPList(PPs))) == 0 ) return nullptr; 424 // CoprimeSupp = SqFr(SPList(PPs)); 425 // for ( ; n > 0 ; --n) 426 // { 427 // CoprimeSupp &=(~(SqFr(MTL[n]) ) ); 428 // if ( CoprimeSupp == 0 ) return nullptr; 429 // } 430 // theSPL = SPList(PPs); 431 // OccInd = Indets(theSPL); 432 // res = eterm_init(eterm_get_indetsNo(theSPL)); 433 // for ( i=IntsGetLen(OccInd) ; i>0 ; --i ) 434 // { 435 // if ( bs_get_n(CoprimeSupp, (size_t)OccInd[i]) ) 436 // { 437 // index = OccInd[i]; 438 // eterm_put_non0_nth(res, index, theSPL[index]); 439 // IntsPutLast(Indets(res), index); 440 // eterm_put0_nth(theSPL, index); 441 // IntsMoveLastToNth(OccInd, i); 442 // } 443 // } 444 // return res; 445 // } 446 447 448 // void MoveNotCoprimeSP(eterm FromSPList, eterm ToSPList, eterm theTerm) 449 // { 450 // ints OccInd = Indets(FromSPList); 451 // size_t i = IntsGetLen(OccInd), index; 452 453 // for ( ; i>0 ; --i ) 454 // if ( theTerm[(index = OccInd[i])] != 0 ) 455 // { 456 // eterm_put_non0_nth(ToSPList, index, FromSPList[index]); 457 // IntsPutLast(Indets(ToSPList), index); 458 // eterm_put0_nth(FromSPList, index); 459 // IntsMoveLastToNth(OccInd, i); 460 // } 461 // } 462 463 // void MoveNotCoprime(PPVector FromTList, PPVector ToTList, eterm theTerm) 464 // { 465 // int FromMTLLen =MTListLen(FromTList), ToMTLLen=0; 466 // MixedPPVector FromMTL =MTList(FromTList), ToMTL =MTList(ToTList); 467 // int n =FromMTLLen; 468 469 // for ( ; n > 0 ; --n) 470 // if ( !eterm_coprime(FromMTL[n], theTerm) ) 471 // { 472 // ToMTL.push_back(FromMTL[n]); 473 // FromMTL.myEraseNth(n); 474 // } 475 476 // MoveNotCoprimeSP(SPList(FromTList), SPList(ToTList), theTerm); 477 // } 478 479 // PPVector SplitIndets(PPVector PPs) 480 // { 481 // PPVector res =GlobalSplitterTList; 482 // int resLen = 0; 483 // MixedPPVector MTL = MTList(PPs), resMTL =MTList(res); 484 // int n = MTListLen(PPs), m; 485 // eterm tn, tm; 486 487 // for ( ; n > 0 ; --n) 488 // { 489 // tn = MTL[n]; 490 // for ( m=resLen ; m>0 ; --m) 491 // { 492 // tm = resMTL[m]; 493 // if ( !eterm_coprime(tn, tm) ) 494 // { 495 // eterm_union_and_assign(tm, tn); 496 // if ( tn != MTL[n] ) eterm_free(tn); 497 // tn = tm; 498 // resMTL.myEraseNth(m); 499 // } 500 // } 501 // if ( tn != MTL[n] ) MTLPutLast(resMTL, resLen, tn); 502 // else MTLPutLast(resMTL, resLen, eterm_dup(tn)); 503 // } 504 // if ( resLen == 1 ) 505 // { 506 // eterm_free(resMTL[1]); 507 // SetMTListLen(res, 0); 508 // return nullptr; 509 // } 510 // SetMTListLen(res, resLen); 511 512 // return res; 513 // } 514 515 // /************************ INTERREDUCE ************************/ 516 517 // void MTLOrderByDecrDegree(MixedPPVector theMTList, int Len, int MaxDeg) 518 // { 519 // eterm t; 520 // MixedPPVector *MTLDeg=(MixedPPVector*)mycalloc 521 // ((MaxDeg+1),sizeof(MixedPPVector),"*MixedPPVector"); 522 // int *MTLLenDeg=(int*)mycalloc((MaxDeg+1), sizeof(int),"*int"), auxLen =0; 523 // int i, d; 524 525 // for ( d=MaxDeg ; d>=2 ; --d) MTLDeg[d] = malloc_MTList(Len); 526 527 // for ( i=Len ; i>0 ; --i) 528 // { 529 // d = eterm_degree((t = theMTList[i]) ); 530 // MTLPutLast( MTLDeg[d], MTLLenDeg[d], t); 531 // } 532 // for ( d=MaxDeg ; d>1 ; --d) 533 // { 534 // if (MTLLenDeg[d] != 0) 535 // for ( i=MTLLenDeg[d]; i>0; --i) 536 // MTLPutLast(theMTList, auxLen, (MTLDeg[d])[i]); 537 // free_MTList(Len, MTLDeg[d]); 538 // } 539 // myfree((MaxDeg+1)*sizeof(MixedPPVector), MTLDeg, "*MixedPPVector"); 540 // myfree((MaxDeg+1)*sizeof(int), MTLLenDeg, "*int"); 541 542 // } 543 544 // bool MTLDividesTerm(MixedPPVector MTL, int MTLLen, eterm T) 545 // { 546 // int i = MTLLen; 547 548 // for ( ; i>0 ; --i ) if ( eterm_divides(MTL[i], T) ) return TRUE; 549 550 // return FALSE; 551 // } 552 553 // void InterreduceTList(PPVector PPs) 554 // { 555 // MixedPPVector MTL =MTList(PPs); 556 // int MTLLen =MTListLen(PPs); 557 // int n =MTLLen, j, MaxDeg =0; 558 // eterm auxT; 559 560 // for ( ; n> 0 ; --n) 561 // if ( TListSimplePowersDivide(PPs,(auxT=MTL[n]) ) ) 562 // MTL.myEraseNth(n); 563 // else 564 // MaxDeg = MAX(MaxDeg, eterm_degree(auxT)); 565 // if (MTLLen > 1) 566 // { 567 // MTLOrderByDecrDegree(MTL, MTLLen, MaxDeg); 568 // for ( n=MTLLen-1 ; n>0; --n) 569 // { 570 // auxT = MTL[n]; 571 // for ( j = MTLLen ; j>n ; --j) 572 // if ( eterm_divides(MTL[j], auxT) ) 573 // { 574 // MTL.myEraseNth(n); 575 // break; 576 // } 577 // } 578 // } 579 // SetMTListLen(PPs, MTLLen); 580 // } 581 582 // /************************ PIVOT ************************/ 583 584 // int MostFrequentIndet(PPVector PPs) 585 // { 586 // int IndNo = TListIndetsNo(PPs), exp, MFIndNo = 1, i, j; 587 // ints OccInd; 588 589 // for ( j=IndNo ; j>0 ; --j ) MFI_Occurrences[j] = 0; 590 // for ( i=len(PPs) ; i>0 ; --i) 591 // { 592 // OccInd = Indets(MTL[i]); 593 // for ( j=IntsGetLen(OccInd) ; j>0 ; --j ) MFI_Occurrences[OccInd[j]]++; 594 // } 595 // MFI_Indets[1] = (j=IndNo); 596 // i = MFI_Occurrences[j]; 597 // for ( --j ; j>0 ; --j ) 598 // if ((exp=MFI_Occurrences[j]) >= i) 599 // { 600 // if (exp > i) { MFIndNo = 0; i = exp; } 601 // MFI_Indets[++MFIndNo] = j; 602 // } 603 604 // if ( i == 1 ) return -1; 605 // if ( MFIndNo == 1 ) return MFI_Indets[1]; 606 // else return MFI_Indets[MFIndNo/2]; 607 // } 608 609 610 // void PPVector::myPivotSetIndet() 611 // { 612 // int MFIndet = MostFrequentIndet(PPs); 613 614 // AssignOne(myPivot()); 615 // if (MFIndet == -1) return; 616 617 // myPivot() = indet(myPivot().myPPM(), MFIndet); 618 // } 619 620 // void PPVector::myPivotSetGCD3() 621 // { 622 // int MFIndet = MostFrequentIndet(PPs); 623 624 // AssignOne(myPivot()); 625 // if (MFIndet == -1) return; 626 627 // // reverse order for compatibility with old C code 628 // for (int i=len(PPs), int count=0; --i>=0 && count<3; ) 629 // if ( PPs[i][MFIndet] != 0 ) 630 // { 631 // myPivot() = gcd(myPivot(), PPs[i]); 632 // count++; 633 // } 634 // } 635 636 // void PPVector::myPivotSetSimplePower() 637 // { 638 // int MFIndet = MostFrequentIndet(PPs); 639 640 // AssignOne(myPivot()); 641 // if (MFIndet == -1) return; 642 643 // int r = -1; 644 // while (myPPs[++r][MFIndet] == 0); 645 // int e1 = myPPs[r][MFIndet]; 646 // r = len(myPPs); 647 // while (myPPs[--r][MFIndet] == 0); 648 // int e2 = myPPs[r][MFIndet]; 649 650 // myPivot() = IndetPower(myPivot().myPPM(), MFIndet, min(e1, e2)); 651 // } 652 653 /* 654 #define PivotOf(PPs,MTLLen)\ 655 ((MTLLen<3)?GCD3PivotOf(PPs):BigPivotOf(PPs)) 656 */ 657 658 // #define PivotOf(PPs,MTLLen) BigPivotOf(PPs) 659 660 // eterm eterm_colon_SP(eterm theTerm, size_t index, unsigned exp) 661 // { 662 // unsigned OldExp; 663 664 // if ((OldExp =theTerm[index]) != 0 ) 665 // { 666 // if ( OldExp <= exp ) printf("OldExp <= exp"); 667 // eterm_put_non0_nth(theTerm, index, OldExp-exp ); 668 // } 669 // return theTerm; 670 // } 671 672 // void OccIndDel(eterm T, int index) 673 // { 674 // ints OccInd = Indets(T); 675 // int n = IntsGetLen(OccInd); 676 677 // while ( OccInd[n] != index ) --n; 678 // IntsMoveLastToNth(OccInd, n); 679 // } 680 681 // void ReduceAndDivideBySimplePower( PPVector PPs, PPVector *DivTList, 682 // size_t PIndex, unsigned PExp) 683 // { 684 // int MTLLen = MTListLen(PPs), auxLen = MAX(MTListLen(PPs),199); 685 // int DivMTLLen, MTLLenEe, i; 686 // unsigned TExp, e, index; 687 // ints DivMTLLenExp; 688 // eterm T; 689 // MixedPPVector *DivMTLExp, MTL =MTList(PPs), DivMTL, MTLEe; 690 // eterm DivSPL; 691 692 // if ( PExp>119 ) 693 // { 694 // DivMTLLenExp = (int*)mymalloc((PExp+2)*sizeof(int), "*int"); 695 // DivMTLExp = (MixedPPVector*)mymalloc((PExp+2)*sizeof(MixedPPVector), 696 // "*MixedPPVector"); 697 // } 698 // else 699 // { DivMTLLenExp = Ints120; DivMTLExp = MTL120; } 700 // for ( i=PExp+1 ; i>=0 ; --i ) 701 // { DivMTLExp[i] = malloc_MTList(auxLen); DivMTLLenExp[i] = 0; } 702 703 // /* DivTList */ 704 // DivSPL = eterm_colon_SP(eterm_dup(SPList(PPs)), PIndex, PExp); 705 // /* PPs */ 706 // InsInSPList(PIndex, PExp, SPList(PPs)); 707 708 // for ( i=MTLLen ; i > 0 ; --i) 709 // { 710 // TExp = (T=MTL[i])[PIndex]; 711 // if ( TExp > PExp ) 712 // { 713 // PPsSwapLastAndNth(MTL, MTLLen, i); 714 // eterm_put_non0_nth(T, PIndex, TExp-PExp); 715 // MTLPutLast(DivMTLExp[PExp+1], DivMTLLenExp[PExp+1], T); 716 // } 717 // else if ( TExp == PExp ) 718 // { 719 // PPsSwapLastAndNth(MTL, MTLLen, i); 720 // if ( eterm_get_OccIndNo(T) == 2 ) 721 // { 722 // if ((index = (Indets(T))[1]) == PIndex ) index = (Indets(T))[2]; 723 // InsInSPList(index, T[index], DivSPL); 724 // eterm_free(T); 725 // } 726 // else 727 // { 728 // eterm_put0_nth(T, PIndex); 729 // OccIndDel(T, PIndex); 730 // MTLPutLast(DivMTLExp[TExp], DivMTLLenExp[TExp], T); 731 // } 732 // } 733 // else if ( TExp != 0 ) 734 // { 735 // if ( eterm_get_OccIndNo(T) == 2 ) 736 // { 737 // if ((index = (Indets(T))[1]) == PIndex ) index = (Indets(T))[2]; 738 // InsInSPList(index, T[index], DivSPL); 739 // } 740 // else 741 // { 742 // eterm_put0_nth(T, PIndex); 743 // OccIndDel(T, PIndex); 744 // MTLPutLast(DivMTLExp[TExp], DivMTLLenExp[TExp], T); 745 // } 746 // } 747 // else /* TExp == 0 */ MTLPutLast(DivMTLExp[TExp], DivMTLLenExp[TExp], T); 748 // } 749 750 // /* DivTList(Interreduction) */ 751 // DivMTL = DivMTLExp[PExp]; 752 // DivMTLLen = DivMTLLenExp[PExp]; 753 // for ( e=PExp-1 ; e>0 ; --e ) 754 // { 755 // MTLEe = DivMTLExp[e]; 756 // MTLLenEe = DivMTLLenExp[e]; 757 // for ( i=MTLLenEe ; i>0 ; --i ) 758 // { 759 // T = MTLEe[i]; 760 // if ( SPLDividesTerm(DivSPL, T) ) PPsSwapLastAndNth(MTLEe, MTLLenEe, i); 761 // else if ( MTLDividesTerm(DivMTL, DivMTLLen, T) ) 762 // PPsSwapLastAndNth(MTLEe, MTLLenEe, i); 763 // else 764 // MTLPutNth(MTLEe, i, eterm_dup(T)); 765 // eterm_put_non0_nth(T, PIndex, e); 766 // IntsPutLast(Indets(T), PIndex); 767 // } 768 // if ( MTLLenEe!=0 ) MTLAppend(DivMTL, &DivMTLLen, MTLEe, MTLLenEe); 769 // free_MTList(auxLen, MTLEe); 770 // } 771 // MTLEe = DivMTLExp[0]; 772 // MTLLenEe = DivMTLLenExp[0]; 773 // for ( i=MTLLenEe ; i>0 ; --i ) 774 // { 775 // T = MTLEe[i]; 776 // if ( SPLDividesTerm(DivSPL, T) ) PPsSwapLastAndNth(MTLEe, MTLLenEe, i); 777 // else if ( MTLDividesTerm(DivMTL, DivMTLLen, T) ) 778 // PPsSwapLastAndNth(MTLEe, MTLLenEe, i); 779 // else 780 // MTLPutNth(MTLEe, i, eterm_dup(T)); 781 // } 782 // if ( MTLLenEe!=0 ) MTLAppend(DivMTL, &DivMTLLen, MTLEe, MTLLenEe); 783 // free_MTList(auxLen, MTLEe); 784 785 // MTLAppend(DivMTL, &DivMTLLen, DivMTLExp[PExp+1], DivMTLLenExp[PExp+1]); 786 // free_MTList(auxLen, DivMTLExp[PExp+1]); 787 // *DivTList = TListMake(DivSPL, DivMTL, auxLen, DivMTLLen); 788 789 // /* PPs */ 790 // TListReduceSize(PPs, MTLLen); 791 // if ( PExp>119 ) 792 // { 793 // myfree((PExp+2)*sizeof(int), DivMTLLenExp, "*int"); 794 // myfree((PExp+2)*sizeof(MixedPPVector), DivMTLExp, "*MixedPPVector"); 795 // } 796 // } 797 798 // void ReduceAndDivideByMixedTerm( PPVector PPs, PPVector *DivTList, 799 // eterm Pivot) 800 // { 801 // int MTLLen =MTListLen(PPs), OldMTLLen =MTLLen, 802 // BMLen =0, CMTLLen =0, DivMTLLen =0; 803 // int index, TDeg; 804 // eterm DivT, t, auxSPL; 805 // int i = MTLLen; 806 // MixedPPVector MTL =MTList(PPs), BigMultMTL, CoprimeMTL, DivMTL; 807 // eterm DivSPL; 808 809 // *DivTList = NewTList(MTLLen, TListIndetsNo(PPs)); 810 // DivMTL = MTList(*DivTList); 811 // DivSPL = eterm_colon(eterm_dup(SPList(PPs)), Pivot); 812 // auxSPL =(SPList(*DivTList)); /* eterm_init() */ 813 814 // BigMultMTL = malloc_MTList(MTLLen); 815 // CoprimeMTL = malloc_MTList(MTLLen); 816 817 // for ( ; i > 0 ; --i ) 818 // { 819 // t=MTL[i]; 820 // TDeg = eterm_degree(t); 821 822 // if ( eterm_divides(Pivot,t) ) 823 // { 824 // /* PPs */ 825 // PPsSwapLastAndNth(MTL, MTLLen, i); 826 // /* DivTList */ 827 // if (sp_BigMult(t,Pivot)) 828 // MTLPutLast(BigMultMTL, BMLen, eterm_colon(t,Pivot)); 829 // else 830 // { 831 // DivT = eterm_colon(t, Pivot); 832 // if ( IntsGetLen(Indets(DivT)) == 1) 833 // { 834 // index = (Indets(DivT))[ 1]; 835 // InsInSPList(index, DivT[index], DivSPL); 836 // InsInSPList(index, DivT[index] + Pivot[index], auxSPL); 837 // eterm_free(DivT); 838 // } 839 // else 840 // MTLPutLast(DivMTL, DivMTLLen, DivT); 841 // } 842 // } 843 // else if ( SPLDividesTerm(auxSPL, t) ) 844 // ; 845 // else if (TDeg==eterm_degree((DivT=eterm_colon(eterm_dup(t),Pivot) )) ) 846 // MTLPutLast(CoprimeMTL, CMTLLen, DivT); 847 // else if ( IntsGetLen(Indets(DivT)) == 1 ) 848 // { 849 // index = (Indets(DivT))[ 1]; 850 // InsInSPList(index, DivT[index], DivSPL); 851 // InsInSPList(index, DivT[index] + Pivot[index], auxSPL); 852 // eterm_free(DivT); 853 // } 854 // else 855 // MTLPutLast(DivMTL, DivMTLLen, DivT); 856 // } 857 // /* PPs */ 858 // MTLPutLast(MTL, MTLLen, Pivot); 859 // TListReduceSize(PPs, MTLLen); 860 // /* DivTList */ 861 // eterm_colon(auxSPL, Pivot); 862 // SetMTListLen(*DivTList, DivMTLLen); 863 // InterreduceTList(*DivTList); 864 // DivMTLLen = MTListLen(*DivTList); 865 866 // if ( DivMTLLen != 0 ) 867 868 // for ( i=CMTLLen ; i>0 ; --i ) 869 // { 870 // if ( SPLDividesTerm(DivSPL,(t=(CoprimeMTL)[i]) ) ) 871 // { 872 // eterm_free(t); 873 // PPsSwapLastAndNth(CoprimeMTL, CMTLLen, i); 874 // } 875 // else if ( MTLDividesTerm(DivMTL, DivMTLLen, t) ) 876 // { 877 // eterm_free(t); 878 // PPsSwapLastAndNth(CoprimeMTL, CMTLLen, i); 879 // } 880 // } 881 // else 882 // for ( i=CMTLLen ; i>0 ; --i ) 883 // if ( SPLDividesTerm(DivSPL,(t=(CoprimeMTL)[i]) ) ) 884 // PPs.myEraseNth(i); 885 // SPList(*DivTList) = DivSPL; 886 // if (CMTLLen!=0) MTLAppend(DivMTL, &DivMTLLen, CoprimeMTL, CMTLLen); 887 // if (BMLen!=0) MTLAppend(DivMTL, &DivMTLLen, BigMultMTL, BMLen); 888 // free_MTList(OldMTLLen, CoprimeMTL); 889 // free_MTList(OldMTLLen, BigMultMTL); 890 // eterm_free(auxSPL); 891 // TListReduceSize(*DivTList, DivMTLLen); 892 // } 893 894 // void ReduceAndDivideByPivot(PPVector PPs, PPVector *DivTList, eterm Pivot) 895 // { 896 // if ( eterm_get_OccIndNo(Pivot) == 1) 897 // { 898 // int index = (Indets(Pivot))[1], exp = eterm_degree(Pivot); 899 // eterm_free(Pivot); 900 // ReduceAndDivideBySimplePower(PPs, DivTList, index, exp); 901 // } 902 // else 903 // ReduceAndDivideByMixedTerm(PPs, DivTList, Pivot); 904 // } 905 906 907 //---------------------------------------------------------------------- 908 909 } // end of namespace CoCoA 910 911 //---------------------------------------------------------------------- 912 // RCS header/log in the next few lines 913 // $Header: /Volumes/Home_1/cocoa/cvs-repository/CoCoALib-0.99/src/AlgebraicCore/TmpPPVector.C,v 1.23 2019/10/21 16:31:45 bigatti Exp $ 914 // $Log: TmpPPVector.C,v $ 915 // Revision 1.23 2019/10/21 16:31:45 bigatti 916 // -- added ctor with vector<RingElem> 917 // -- using range-based for loop 918 // 919 // Revision 1.22 2019/03/19 11:07:08 abbott 920 // Summary: Replaced 0 by nullptr where appropriate 921 // 922 // Revision 1.21 2018/05/18 16:38:51 bigatti 923 // -- added include SparsePolyOps-RingElem.H 924 // 925 // Revision 1.20 2018/05/17 15:55:26 bigatti 926 // -- renamed VectorOperations --> VectorOps 927 // 928 // Revision 1.19 2016/11/11 14:15:34 abbott 929 // Summary: Added short-cut to operator<< when ostream is in bad state 930 // 931 // Revision 1.18 2016/05/23 12:48:40 bigatti 932 // -- just a comment (exception safety?) 933 // 934 // Revision 1.17 2015/06/11 16:57:24 bigatti 935 // -- using new functions monomial(ring, pp) and monomial(ring, expv) 936 // 937 // Revision 1.16 2014/07/31 14:45:19 abbott 938 // Summary: Merged io.H and UtilsTemplate.H into new header VectorOperations.H 939 // Author: JAA 940 // 941 // Revision 1.15 2014/07/07 13:15:57 abbott 942 // Summary: Removed AsSparsePolyRing 943 // Author: JAA 944 // 945 // Revision 1.14 2014/04/30 16:26:22 abbott 946 // Summary: Replaced X.size() by len(X); Replaced size_t by long; ONLY IN COMMENTED OUT CODE!! 947 // Author: JAA 948 // 949 // Revision 1.13 2011/12/05 16:56:22 bigatti 950 // -- changed: MachineInteger --> MachineInt (just in comment) 951 // 952 // Revision 1.12 2011/11/07 11:04:51 bigatti 953 // -- AreMonomials is now public 954 // 955 // Revision 1.11 2011/07/27 15:51:14 bigatti 956 // -- improved myDivides 957 // 958 // Revision 1.10 2011/06/27 12:50:56 bigatti 959 // -- added mySupport, myAlexanderDual 960 // 961 // Revision 1.9 2011/05/25 12:21:38 bigatti 962 // -- added myClear 963 // 964 // Revision 1.8 2011/03/11 11:05:44 bigatti 965 // -- changed size_t --> long 966 // -- changed size --> len 967 // 968 // Revision 1.7 2010/04/27 16:10:07 bigatti 969 // -- changed sorting (ascending) 970 // 971 // Revision 1.6 2010/02/03 11:56:38 bigatti 972 // -- more flexible conversion PPVector/vector<RingElem> 973 // 974 // Revision 1.5 2009/10/29 18:43:51 abbott 975 // Added necessary include directive for <algorithm>. 976 // 977 // Revision 1.4 2008/07/04 12:12:35 bigatti 978 // -- added operator<< 979 // 980 // Revision 1.3 2008/07/04 09:11:04 bigatti 981 // -- new PPVector class 982 // 983