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