001 /*
002 * $Id: HenselMultUtil.java 3732 2011-08-13 19:25:38Z kredel $
003 */
004
005 package edu.jas.ufd;
006
007
008 import java.util.ArrayList;
009 import java.util.List;
010 import java.util.Map;
011
012 import org.apache.log4j.Logger;
013
014 import edu.jas.arith.BigInteger;
015 import edu.jas.arith.ModIntegerRing;
016 import edu.jas.arith.ModLongRing;
017 import edu.jas.arith.Modular;
018 import edu.jas.arith.ModularRingFactory;
019 import edu.jas.poly.GenPolynomial;
020 import edu.jas.poly.GenPolynomialRing;
021 import edu.jas.poly.PolyUtil;
022 import edu.jas.poly.ExpVector;
023 import edu.jas.ps.PolynomialTaylorFunction;
024 import edu.jas.ps.TaylorFunction;
025 import edu.jas.ps.UnivPowerSeries;
026 import edu.jas.ps.UnivPowerSeriesRing;
027 import edu.jas.structure.GcdRingElem;
028 import edu.jas.structure.Power;
029 import edu.jas.structure.RingFactory;
030
031
032 /**
033 * Hensel multivariate lifting utilities.
034 * @author Heinz Kredel
035 */
036
037 public class HenselMultUtil {
038
039
040 private static final Logger logger = Logger.getLogger(HenselMultUtil.class);
041
042
043 private static final boolean debug = logger.isInfoEnabled();
044
045
046 /**
047 * Modular diophantine equation solution and lifting algorithm. Let p =
048 * A_i.ring.coFac.modul() and assume ggt(A,B) == 1 mod p.
049 * @param A modular GenPolynomial, mod p^k
050 * @param B modular GenPolynomial, mod p^k
051 * @param C modular GenPolynomial, mod p^k
052 * @param V list of substitution values, mod p^k
053 * @param d desired approximation exponent (x_i-v_i)^d.
054 * @param k desired approximation exponent p^k.
055 * @return [s, t] with s A' + t B' = C mod p^k, with A' = B, B' = A.
056 */
057 public static <MOD extends GcdRingElem<MOD> & Modular> List<GenPolynomial<MOD>> liftDiophant(
058 GenPolynomial<MOD> A, GenPolynomial<MOD> B, GenPolynomial<MOD> C, List<MOD> V, long d,
059 long k) throws NoLiftingException {
060 GenPolynomialRing<MOD> pkfac = C.ring;
061 if (pkfac.nvar == 1) { // V, d ignored
062 return HenselUtil.<MOD> liftDiophant(A, B, C, k);
063 }
064 if (!pkfac.equals(A.ring)) {
065 throw new IllegalArgumentException("A.ring != pkfac: " + A.ring + " != " + pkfac);
066 }
067
068 // evaluate at v_n:
069 List<MOD> Vp = new ArrayList<MOD>(V);
070 MOD v = Vp.remove(Vp.size() - 1);
071 GenPolynomial<MOD> zero = pkfac.getZERO();
072 // (x_n - v)
073 GenPolynomial<MOD> mon = pkfac.getONE();
074 GenPolynomial<MOD> xv = pkfac.univariate(0, 1);
075 xv = xv.subtract(pkfac.fromInteger(v.getSymmetricInteger().getVal()));
076 //System.out.println("xv = " + xv);
077 // A(v), B(v), C(v)
078 ModularRingFactory<MOD> cf = (ModularRingFactory<MOD>) pkfac.coFac;
079 MOD vp = cf.fromInteger(v.getSymmetricInteger().getVal());
080 //System.out.println("v = " + v + ", vp = " + vp);
081 GenPolynomialRing<MOD> ckfac = pkfac.contract(1);
082 GenPolynomial<MOD> Ap = PolyUtil.<MOD> evaluateMain(ckfac, A, vp);
083 GenPolynomial<MOD> Bp = PolyUtil.<MOD> evaluateMain(ckfac, B, vp);
084 GenPolynomial<MOD> Cp = PolyUtil.<MOD> evaluateMain(ckfac, C, vp);
085 //System.out.println("Ap = " + Ap);
086 //System.out.println("Bp = " + Bp);
087 //System.out.println("Cp = " + Cp);
088
089 // recursion:
090 List<GenPolynomial<MOD>> su = HenselMultUtil.<MOD> liftDiophant(Ap, Bp, Cp, Vp, d, k);
091 //System.out.println("su@p^" + k + " = " + su);
092 //System.out.println("coFac = " + su.get(0).ring.coFac.toScript());
093 if (pkfac.nvar == 2 && !HenselUtil.<MOD> isDiophantLift(Bp, Ap, su.get(0), su.get(1), Cp)) {
094 System.out.println("isDiophantLift: false");
095 }
096 if (!ckfac.equals(su.get(0).ring)) {
097 throw new IllegalArgumentException("qfac != ckfac: " + su.get(0).ring + " != " + ckfac);
098 }
099 GenPolynomialRing<BigInteger> ifac = new GenPolynomialRing<BigInteger>(new BigInteger(), pkfac);
100 GenPolynomialRing<BigInteger> cifac = new GenPolynomialRing<BigInteger>(new BigInteger(), ckfac);
101 //System.out.println("ifac = " + ifac.toScript());
102 String[] mn = new String[] { pkfac.getVars()[pkfac.nvar - 1] };
103 GenPolynomialRing<GenPolynomial<MOD>> qrfac = new GenPolynomialRing<GenPolynomial<MOD>>(ckfac, 1, mn);
104 //System.out.println("qrfac = " + qrfac);
105
106 List<GenPolynomial<MOD>> sup = new ArrayList<GenPolynomial<MOD>>(su.size());
107 List<GenPolynomial<BigInteger>> supi = new ArrayList<GenPolynomial<BigInteger>>(su.size());
108 for (GenPolynomial<MOD> s : su) {
109 GenPolynomial<MOD> sp = s.extend(pkfac, 0, 0L);
110 sup.add(sp);
111 GenPolynomial<BigInteger> spi = PolyUtil.integerFromModularCoefficients(ifac, sp);
112 supi.add(spi);
113 }
114 //System.out.println("sup = " + sup);
115 //System.out.println("supi = " + supi);
116 GenPolynomial<BigInteger> Ai = PolyUtil.integerFromModularCoefficients(ifac, A);
117 GenPolynomial<BigInteger> Bi = PolyUtil.integerFromModularCoefficients(ifac, B);
118 GenPolynomial<BigInteger> Ci = PolyUtil.integerFromModularCoefficients(ifac, C);
119 //System.out.println("Ai = " + Ai);
120 //System.out.println("Bi = " + Bi);
121 //System.out.println("Ci = " + Ci);
122 GenPolynomial<MOD> aq = PolyUtil.<MOD> fromIntegerCoefficients(pkfac, Ai);
123 GenPolynomial<MOD> bq = PolyUtil.<MOD> fromIntegerCoefficients(pkfac, Bi);
124 //System.out.println("aq = " + aq);
125 //System.out.println("bq = " + bq);
126
127 // compute error:
128 GenPolynomial<BigInteger> E = Ci; // - sum_i s_i b_i
129 E = E.subtract(Bi.multiply(supi.get(0)));
130 E = E.subtract(Ai.multiply(supi.get(1)));
131 //System.out.println("E = " + E);
132 if (E.isZERO()) {
133 logger.info("liftDiophant leaving on zero E");
134 return sup;
135 }
136 GenPolynomial<MOD> Ep = PolyUtil.<MOD> fromIntegerCoefficients(pkfac, E);
137 //System.out.println("Ep(0," + pkfac.nvar + ") = " + Ep);
138 logger.info("Ep(0," + pkfac.nvar + ") = " + Ep);
139 if (Ep.isZERO()) {
140 logger.info("liftDiophant leaving on zero Ep mod p^k");
141 return sup;
142 }
143 for (int e = 1; e <= d; e++) {
144 //System.out.println("\ne = " + e + " -------------------------------------- " + pkfac.nvar);
145 GenPolynomial<GenPolynomial<MOD>> Epr = PolyUtil.<MOD> recursive(qrfac, Ep);
146 //System.out.println("Epr = " + Epr);
147 UnivPowerSeriesRing<GenPolynomial<MOD>> psfac = new UnivPowerSeriesRing<GenPolynomial<MOD>>(qrfac);
148 //System.out.println("psfac = " + psfac);
149 TaylorFunction<GenPolynomial<MOD>> F = new PolynomialTaylorFunction<GenPolynomial<MOD>>(Epr);
150 //System.out.println("F = " + F);
151 List<GenPolynomial<MOD>> Vs = new ArrayList<GenPolynomial<MOD>>(1);
152 GenPolynomial<MOD> vq = ckfac.fromInteger(v.getSymmetricInteger().getVal());
153 Vs.add(vq);
154 //System.out.println("Vs = " + Vs);
155 UnivPowerSeries<GenPolynomial<MOD>> Epst = psfac.seriesOfTaylor(F, vq);
156 //System.out.println("Epst = " + Epst);
157 GenPolynomial<MOD> cm = Epst.coefficient(e);
158 //System.out.println("cm = " + cm + ", cm.ring = " + cm.ring.toScript());
159
160 // recursion:
161 List<GenPolynomial<MOD>> S = HenselMultUtil.<MOD> liftDiophant(Ap, Bp, cm, Vp, d, k);
162 //System.out.println("S = " + S);
163 if (!ckfac.coFac.equals(S.get(0).ring.coFac)) {
164 throw new IllegalArgumentException("ckfac != pkfac: " + ckfac.coFac + " != "
165 + S.get(0).ring.coFac);
166 }
167 if (pkfac.nvar == 2 && !HenselUtil.<MOD> isDiophantLift(Ap, Bp, S.get(1), S.get(0), cm)) {
168 System.out.println("isDiophantLift: false");
169 }
170 mon = mon.multiply(xv); // Power.<GenPolynomial<MOD>> power(pkfac,xv,e);
171 //System.out.println("mon = " + mon);
172 List<GenPolynomial<MOD>> Sp = new ArrayList<GenPolynomial<MOD>>(S.size());
173 int i = 0;
174 supi = new ArrayList<GenPolynomial<BigInteger>>(su.size());
175 for (GenPolynomial<MOD> dd : S) {
176 //System.out.println("dd = " + dd);
177 GenPolynomial<MOD> de = dd.extend(pkfac, 0, 0L);
178 GenPolynomial<MOD> dm = de.multiply(mon);
179 Sp.add(dm);
180 de = sup.get(i).sum(dm);
181 //System.out.println("dd = " + dd);
182 sup.set(i++, de);
183 GenPolynomial<BigInteger> spi = PolyUtil.integerFromModularCoefficients(ifac, dm);
184 supi.add(spi);
185 }
186 //System.out.println("Sp = " + Sp);
187 //System.out.println("sup = " + sup);
188 //System.out.println("supi = " + supi);
189 // compute new error
190 //E = E; // - sum_i s_i b_i
191 E = E.subtract(Bi.multiply(supi.get(0)));
192 E = E.subtract(Ai.multiply(supi.get(1)));
193 //System.out.println("E = " + E);
194 if (E.isZERO()) {
195 logger.info("liftDiophant leaving on zero E");
196 return sup;
197 }
198 Ep = PolyUtil.<MOD> fromIntegerCoefficients(pkfac, E);
199 //System.out.println("Ep(" + e + "," + pkfac.nvar + ") = " + Ep);
200 logger.info("Ep(" + e + "," + pkfac.nvar + ") = " + Ep);
201 if (Ep.isZERO()) {
202 logger.info("liftDiophant leaving on zero Ep mod p^k");
203 return sup;
204 }
205 }
206 //System.out.println("*** done: " + pkfac.nvar);
207 return sup;
208 }
209
210
211 /**
212 * Modular diophantine equation solution and lifting algorithm. Let p =
213 * A_i.ring.coFac.modul() and assume ggt(a,b) == 1 mod p, for a, b in A.
214 * @param A list of modular GenPolynomials, mod p^k
215 * @param C modular GenPolynomial, mod p^k
216 * @param V list of substitution values, mod p^k
217 * @param d desired approximation exponent (x_i-v_i)^d.
218 * @param k desired approximation exponent p^k.
219 * @return [s_1,..., s_n] with sum_i s_i A_i' = C mod p^k, with Ai' =
220 * prod_{j!=i} A_j.
221 */
222 public static <MOD extends GcdRingElem<MOD> & Modular> List<GenPolynomial<MOD>> liftDiophant(
223 List<GenPolynomial<MOD>> A, GenPolynomial<MOD> C, List<MOD> V, long d, long k)
224 throws NoLiftingException {
225 GenPolynomialRing<MOD> pkfac = C.ring;
226 if (pkfac.nvar == 1) { // V, d ignored
227 return HenselUtil.<MOD> liftDiophant(A, C, k);
228 }
229 if (!pkfac.equals(A.get(0).ring)) {
230 throw new IllegalArgumentException("A.ring != pkfac: " + A.get(0).ring + " != " + pkfac);
231 }
232 // co-products
233 GenPolynomial<MOD> As = pkfac.getONE();
234 for (GenPolynomial<MOD> a : A) {
235 As = As.multiply(a);
236 }
237 List<GenPolynomial<MOD>> Bp = new ArrayList<GenPolynomial<MOD>>(A.size());
238 for (GenPolynomial<MOD> a : A) {
239 GenPolynomial<MOD> b = PolyUtil.<MOD> basePseudoDivide(As, a);
240 Bp.add(b);
241 }
242
243 // evaluate at v_n:
244 List<MOD> Vp = new ArrayList<MOD>(V);
245 MOD v = Vp.remove(Vp.size() - 1);
246 // (x_n - v)
247 GenPolynomial<MOD> mon = pkfac.getONE();
248 GenPolynomial<MOD> xv = pkfac.univariate(0, 1);
249 xv = xv.subtract(pkfac.fromInteger(v.getSymmetricInteger().getVal()));
250 //System.out.println("xv = " + xv);
251 // A(v), B(v), C(v)
252 ModularRingFactory<MOD> cf = (ModularRingFactory<MOD>) pkfac.coFac;
253 MOD vp = cf.fromInteger(v.getSymmetricInteger().getVal());
254 //System.out.println("v = " + v + ", vp = " + vp);
255 GenPolynomialRing<MOD> ckfac = pkfac.contract(1);
256 List<GenPolynomial<MOD>> Ap = new ArrayList<GenPolynomial<MOD>>(A.size());
257 for (GenPolynomial<MOD> a : A) {
258 GenPolynomial<MOD> ap = PolyUtil.<MOD> evaluateMain(ckfac, a, vp);
259 Ap.add(ap);
260 }
261 GenPolynomial<MOD> Cp = PolyUtil.<MOD> evaluateMain(ckfac, C, vp);
262 //System.out.println("Ap = " + Ap);
263 //System.out.println("Cp = " + Cp);
264
265 // recursion:
266 List<GenPolynomial<MOD>> su = HenselMultUtil.<MOD> liftDiophant(Ap, Cp, Vp, d, k);
267 //System.out.println("su@p^" + k + " = " + su);
268 //System.out.println("coFac = " + su.get(0).ring.coFac.toScript());
269 if (pkfac.nvar == 2 && !HenselUtil.<MOD> isDiophantLift(Ap, su, Cp)) {
270 System.out.println("isDiophantLift: false");
271 }
272 if (!ckfac.equals(su.get(0).ring)) {
273 throw new IllegalArgumentException("qfac != ckfac: " + su.get(0).ring + " != " + ckfac);
274 }
275 GenPolynomialRing<BigInteger> ifac = new GenPolynomialRing<BigInteger>(new BigInteger(), pkfac);
276 GenPolynomialRing<BigInteger> cifac = new GenPolynomialRing<BigInteger>(new BigInteger(), ckfac);
277 //System.out.println("ifac = " + ifac.toScript());
278 String[] mn = new String[] { pkfac.getVars()[pkfac.nvar - 1] };
279 GenPolynomialRing<GenPolynomial<MOD>> qrfac = new GenPolynomialRing<GenPolynomial<MOD>>(ckfac, 1, mn);
280 //System.out.println("qrfac = " + qrfac);
281
282 List<GenPolynomial<MOD>> sup = new ArrayList<GenPolynomial<MOD>>(su.size());
283 List<GenPolynomial<BigInteger>> supi = new ArrayList<GenPolynomial<BigInteger>>(su.size());
284 for (GenPolynomial<MOD> s : su) {
285 GenPolynomial<MOD> sp = s.extend(pkfac, 0, 0L);
286 sup.add(sp);
287 GenPolynomial<BigInteger> spi = PolyUtil.integerFromModularCoefficients(ifac, sp);
288 supi.add(spi);
289 }
290 //System.out.println("sup = " + sup);
291 //System.out.println("supi = " + supi);
292 List<GenPolynomial<BigInteger>> Ai = new ArrayList<GenPolynomial<BigInteger>>(A.size());
293 for (GenPolynomial<MOD> a : A) {
294 GenPolynomial<BigInteger> ai = PolyUtil.integerFromModularCoefficients(ifac, a);
295 Ai.add(ai);
296 }
297 List<GenPolynomial<BigInteger>> Bi = new ArrayList<GenPolynomial<BigInteger>>(A.size());
298 for (GenPolynomial<MOD> b : Bp) {
299 GenPolynomial<BigInteger> bi = PolyUtil.integerFromModularCoefficients(ifac, b);
300 Bi.add(bi);
301 }
302 GenPolynomial<BigInteger> Ci = PolyUtil.integerFromModularCoefficients(ifac, C);
303 //System.out.println("Ai = " + Ai);
304 //System.out.println("Ci = " + Ci);
305
306 List<GenPolynomial<MOD>> Aq = new ArrayList<GenPolynomial<MOD>>(A.size());
307 for (GenPolynomial<BigInteger> ai : Ai) {
308 GenPolynomial<MOD> aq = PolyUtil.<MOD> fromIntegerCoefficients(pkfac, ai);
309 Aq.add(aq);
310 }
311 //System.out.println("Aq = " + Aq);
312
313 // compute error:
314 GenPolynomial<BigInteger> E = Ci; // - sum_i s_i b_i
315 int i = 0;
316 for (GenPolynomial<BigInteger> bi : Bi) {
317 E = E.subtract(bi.multiply(supi.get(i++)));
318 }
319 //System.out.println("E = " + E);
320 if (E.isZERO()) {
321 logger.info("liftDiophant leaving on zero E");
322 return sup;
323 }
324 GenPolynomial<MOD> Ep = PolyUtil.<MOD> fromIntegerCoefficients(pkfac, E);
325 //System.out.println("Ep(0," + pkfac.nvar + ") = " + Ep);
326 logger.info("Ep(0," + pkfac.nvar + ") = " + Ep);
327 if (Ep.isZERO()) {
328 logger.info("liftDiophant leaving on zero Ep mod p^k");
329 return sup;
330 }
331 for (int e = 1; e <= d; e++) {
332 //System.out.println("\ne = " + e + " -------------------------------------- " + pkfac.nvar);
333 GenPolynomial<GenPolynomial<MOD>> Epr = PolyUtil.<MOD> recursive(qrfac, Ep);
334 //System.out.println("Epr = " + Epr);
335 UnivPowerSeriesRing<GenPolynomial<MOD>> psfac = new UnivPowerSeriesRing<GenPolynomial<MOD>>(qrfac);
336 //System.out.println("psfac = " + psfac);
337 TaylorFunction<GenPolynomial<MOD>> F = new PolynomialTaylorFunction<GenPolynomial<MOD>>(Epr);
338 //System.out.println("F = " + F);
339 List<GenPolynomial<MOD>> Vs = new ArrayList<GenPolynomial<MOD>>(1);
340 GenPolynomial<MOD> vq = ckfac.fromInteger(v.getSymmetricInteger().getVal());
341 Vs.add(vq);
342 //System.out.println("Vs = " + Vs);
343 UnivPowerSeries<GenPolynomial<MOD>> Epst = psfac.seriesOfTaylor(F, vq);
344 //System.out.println("Epst = " + Epst);
345 GenPolynomial<MOD> cm = Epst.coefficient(e);
346 //System.out.println("cm = " + cm + ", cm.ring = " + cm.ring.toScript());
347 if (cm.isZERO()) {
348 continue;
349 }
350 // recursion:
351 List<GenPolynomial<MOD>> S = HenselMultUtil.<MOD> liftDiophant(Ap, cm, Vp, d, k);
352 //System.out.println("S = " + S);
353 if (!ckfac.coFac.equals(S.get(0).ring.coFac)) {
354 throw new IllegalArgumentException("ckfac != pkfac: " + ckfac.coFac + " != "
355 + S.get(0).ring.coFac);
356 }
357 if (pkfac.nvar == 2 && !HenselUtil.<MOD> isDiophantLift(Ap, S, cm)) {
358 System.out.println("isDiophantLift: false");
359 }
360 mon = mon.multiply(xv); // Power.<GenPolynomial<MOD>> power(pkfac,xv,e);
361 //System.out.println("mon = " + mon);
362 List<GenPolynomial<MOD>> Sp = new ArrayList<GenPolynomial<MOD>>(S.size());
363 i = 0;
364 supi = new ArrayList<GenPolynomial<BigInteger>>(su.size());
365 for (GenPolynomial<MOD> dd : S) {
366 //System.out.println("dd = " + dd);
367 GenPolynomial<MOD> de = dd.extend(pkfac, 0, 0L);
368 GenPolynomial<MOD> dm = de.multiply(mon);
369 Sp.add(dm);
370 de = sup.get(i).sum(dm);
371 //System.out.println("dd = " + dd);
372 sup.set(i++, de);
373 GenPolynomial<BigInteger> spi = PolyUtil.integerFromModularCoefficients(ifac, dm);
374 supi.add(spi);
375 }
376 //System.out.println("Sp = " + Sp);
377 //System.out.println("sup = " + sup);
378 //System.out.println("supi = " + supi);
379 // compute new error
380 //E = E; // - sum_i s_i b_i
381 i = 0;
382 for (GenPolynomial<BigInteger> bi : Bi) {
383 E = E.subtract(bi.multiply(supi.get(i++)));
384 }
385 //System.out.println("E = " + E);
386 if (E.isZERO()) {
387 logger.info("liftDiophant leaving on zero E");
388 return sup;
389 }
390 Ep = PolyUtil.<MOD> fromIntegerCoefficients(pkfac, E);
391 //System.out.println("Ep(" + e + "," + pkfac.nvar + ") = " + Ep);
392 logger.info("Ep(" + e + "," + pkfac.nvar + ") = " + Ep);
393 if (Ep.isZERO()) {
394 logger.info("liftDiophant leaving on zero Ep mod p^k");
395 return sup;
396 }
397 }
398 //System.out.println("*** done: " + pkfac.nvar);
399 return sup;
400 }
401
402
403 /**
404 * Modular Hensel lifting algorithm on coefficients test. Let p =
405 * f_i.ring.coFac.modul() and assume C == prod_{0,...,n-1} f_i mod p with
406 * gcd(f_i,f_j) == 1 mod p for i != j
407 * @param C integer polynomial
408 * @param Cp GenPolynomial mod p^k
409 * @param F = [f_0,...,f_{n-1}] list of monic modular polynomials.
410 * @param k approximation exponent.
411 * @param L = [g_0,...,g_{n-1}] list of lifted modular polynomials.
412 * @return true if C = prod_{0,...,n-1} g_i mod p^k, else false.
413 */
414 public static <MOD extends GcdRingElem<MOD> & Modular>
415 boolean isHenselLift(GenPolynomial<BigInteger> C, GenPolynomial<MOD> Cp,
416 List<GenPolynomial<MOD>> F, long k, List<GenPolynomial<MOD>> L) {
417 boolean t = true;
418 GenPolynomialRing<MOD> qfac = L.get(0).ring;
419 GenPolynomial<MOD> q = qfac.getONE();
420 for ( GenPolynomial<MOD> fi : L ) {
421 q = q.multiply(fi);
422 }
423 t = Cp.equals(q);
424 if ( ! t ) {
425 System.out.println("Cp = " + Cp);
426 System.out.println("q = " + q);
427 System.out.println("Cp != q: " + Cp.subtract(q));
428 return t;
429 }
430 GenPolynomialRing<BigInteger> dfac = C.ring;
431 GenPolynomial<BigInteger> Ci = PolyUtil.integerFromModularCoefficients(dfac, q);
432 t = C.equals(Ci);
433 if ( ! t ) {
434 System.out.println("C = " + C);
435 System.out.println("Ci = " + Ci);
436 System.out.println("C != Ci: " + C.subtract(Ci));
437 return t;
438 }
439 // test L mod id(V) == F
440 return t;
441 }
442
443
444 /**
445 * Modular Hensel lifting algorithm, monic case. Let p = A_i.ring.coFac.modul() and
446 * assume ggt(a,b) == 1 mod p, for a, b in A.
447 * @param C monic GenPolynomial with integer coefficients
448 * @param Cp GenPolynomial mod p^k
449 * @param F list of modular GenPolynomials, mod (I_v, p^k )
450 * @param V list of substitution values, mod p^k
451 * @param k desired approximation exponent p^k.
452 * @return [g'_1,..., g'_n] with prod_i g'_i = Cp mod p^k.
453 */
454 public static <MOD extends GcdRingElem<MOD> & Modular> List<GenPolynomial<MOD>> liftHenselMonic(
455 GenPolynomial<BigInteger> C, GenPolynomial<MOD> Cp, List<GenPolynomial<MOD>> F,
456 List<MOD> V, long k) throws NoLiftingException {
457 GenPolynomialRing<MOD> pkfac = Cp.ring;
458 //if (pkfac.nvar == 1) { // V ignored
459 // return HenselUtil.<MOD> liftHenselMonic(C,F,k);
460 //}
461 long d = C.degree();
462 //System.out.println("d = " + d);
463 // prepare stack of polynomial rings and polynomials
464 List<GenPolynomialRing<MOD>> Pfac = new ArrayList<GenPolynomialRing<MOD>>();
465 List<GenPolynomial<MOD>> Ap = new ArrayList<GenPolynomial<MOD>>();
466 List<MOD> Vb = new ArrayList<MOD>();
467 MOD v = V.get(0);
468 Pfac.add(pkfac);
469 Ap.add(Cp);
470 Vb.add(v);
471 GenPolynomialRing<MOD> pf = pkfac;
472 GenPolynomial<MOD> ap = Cp;
473 for (int j = pkfac.nvar; j > 2; j--) {
474 pf = pf.contract(1);
475 Pfac.add(0, pf);
476 MOD vp = pkfac.coFac.fromInteger(V.get(j - 2).getSymmetricInteger().getVal());
477 //System.out.println("vp = " + vp);
478 Vb.add(1, vp);
479 ap = PolyUtil.<MOD> evaluateMain(pf, ap, vp);
480 Ap.add(0, ap);
481 }
482 //System.out.println("Pfac = " + Pfac);
483 if (debug) {
484 logger.debug("Pfac = " + Pfac);
485 }
486 //System.out.println("Ap = " + Ap);
487 //System.out.println("V = " + V);
488 //System.out.println("Vb = " + Vb);
489 // setup bi-variate base case
490 GenPolynomialRing<MOD> pk1fac = F.get(0).ring;
491 if (!pkfac.coFac.equals(pk1fac.coFac)) {
492 throw new IllegalArgumentException("F.ring != pkfac: " + pk1fac + " != " + pkfac);
493 }
494 // TODO: adjust leading coefficients
495 pkfac = Pfac.get(0);
496 Cp = Ap.get(0);
497 //System.out.println("pkfac = " + pkfac.toScript());
498 //System.out.println("pk1fac = " + pk1fac.toScript());
499 GenPolynomialRing<BigInteger> i1fac = new GenPolynomialRing<BigInteger>(new BigInteger(), pk1fac);
500 //System.out.println("i1fac = " + i1fac.toScript());
501 List<GenPolynomial<BigInteger>> Bi = new ArrayList<GenPolynomial<BigInteger>>(F.size());
502 for (GenPolynomial<MOD> b : F) {
503 GenPolynomial<BigInteger> bi = PolyUtil.integerFromModularCoefficients(i1fac, b);
504 Bi.add(bi);
505 }
506 //System.out.println("Bi = " + Bi);
507 // evaluate Cp at v_n:
508 ModularRingFactory<MOD> cf = (ModularRingFactory<MOD>) pkfac.coFac;
509 MOD vp = cf.fromInteger(v.getSymmetricInteger().getVal());
510 //System.out.println("v = " + v + ", vp = " + vp);
511 GenPolynomialRing<MOD> ckfac = pkfac.contract(1);
512 GenPolynomial<MOD> Cs = PolyUtil.<MOD> evaluateMain(ckfac, Cp, vp);
513 //System.out.println("Cp = " + Cp);
514 //System.out.println("Cs = " + Cs);
515
516 List<GenPolynomial<MOD>> U = new ArrayList<GenPolynomial<MOD>>(F.size());
517 for (GenPolynomial<MOD> b : F) {
518 GenPolynomial<MOD> bi = b.extend(pkfac, 0, 0L);
519 U.add(bi);
520 }
521 //System.out.println("U = " + U);
522 List<GenPolynomial<MOD>> U1 = F;
523 //System.out.println("U1 = " + U1);
524
525 GenPolynomial<BigInteger> E = C.ring.getZERO();
526 List<MOD> Vh = new ArrayList<MOD>();
527
528 while (Pfac.size() > 0) { // loop through stack of polynomial rings
529 pkfac = Pfac.remove(0);
530 Cp = Ap.remove(0);
531 v = Vb.remove(0);
532 //Vh.add(0,v);
533 //System.out.println("\npkfac = " + pkfac.toScript() + " ================================== " + Vh);
534
535 // (x_n - v)
536 GenPolynomial<MOD> mon = pkfac.getONE();
537 GenPolynomial<MOD> xv = pkfac.univariate(0, 1);
538 xv = xv.subtract(pkfac.fromInteger(v.getSymmetricInteger().getVal()));
539 //System.out.println("xv = " + xv);
540
541 long deg = Cp.degree(pkfac.nvar - 1);
542 //System.out.println("deg = " + deg);
543
544 GenPolynomialRing<BigInteger> ifac = new GenPolynomialRing<BigInteger>(new BigInteger(), pkfac);
545 //System.out.println("ifac = " + ifac.toScript());
546 List<GenPolynomial<BigInteger>> Bip = new ArrayList<GenPolynomial<BigInteger>>(F.size());
547 for (GenPolynomial<BigInteger> b : Bi) {
548 GenPolynomial<BigInteger> bi = b.extend(ifac, 0, 0L);
549 Bip.add(bi);
550 }
551 Bi = Bip;
552 //System.out.println("Bi = " + Bi);
553 GenPolynomial<BigInteger> Ci = PolyUtil.integerFromModularCoefficients(ifac, Cp);
554 //System.out.println("Ci = " + Ci);
555
556 // compute error:
557 E = ifac.getONE();
558 for (GenPolynomial<BigInteger> bi : Bi) {
559 E = E.multiply(bi);
560 }
561 E = Ci.subtract(E);
562 //System.out.println("E = " + E);
563 GenPolynomial<MOD> Ep = PolyUtil.<MOD> fromIntegerCoefficients(pkfac, E);
564 //System.out.println("Ep(0," + pkfac.nvar + ") = " + Ep);
565 logger.info("Ep(0," + deg + "," + pkfac.nvar + ") = " + Ep);
566
567 String[] mn = new String[] { pkfac.getVars()[pkfac.nvar - 1] };
568 ckfac = pkfac.contract(1);
569 GenPolynomialRing<GenPolynomial<MOD>> pkrfac = new GenPolynomialRing<GenPolynomial<MOD>>(ckfac,
570 1, mn);
571 //System.out.println("pkrfac = " + pkrfac.toScript());
572
573 for (int e = 1; e <= deg && !Ep.isZERO(); e++) {
574 //System.out.println("\ne = " + e + " -------------------------------------- " + pkfac.nvar);
575 GenPolynomial<GenPolynomial<MOD>> Epr = PolyUtil.<MOD> recursive(pkrfac, Ep);
576 //System.out.println("Epr = " + Epr);
577 UnivPowerSeriesRing<GenPolynomial<MOD>> psfac = new UnivPowerSeriesRing<GenPolynomial<MOD>>(
578 pkrfac);
579 //System.out.println("psfac = " + psfac);
580 TaylorFunction<GenPolynomial<MOD>> T = new PolynomialTaylorFunction<GenPolynomial<MOD>>(Epr);
581 //System.out.println("T = " + T);
582 List<GenPolynomial<MOD>> Vs = new ArrayList<GenPolynomial<MOD>>(1);
583 GenPolynomial<MOD> vq = ckfac.fromInteger(v.getSymmetricInteger().getVal());
584 Vs.add(vq);
585 //System.out.println("Vs = " + Vs + ", Vh = " + Vh);
586 UnivPowerSeries<GenPolynomial<MOD>> Epst = psfac.seriesOfTaylor(T, vq);
587 //System.out.println("Epst = " + Epst);
588 logger.info("Epst(" + e + "," + deg + ", " + pkfac.nvar + ") = " + Epst);
589 GenPolynomial<MOD> cm = Epst.coefficient(e);
590 //System.out.println("cm = " + cm);
591 if (cm.isZERO()) {
592 continue;
593 }
594 List<GenPolynomial<MOD>> Ud = HenselMultUtil.<MOD> liftDiophant(U1, cm, Vh, d, k);
595 //System.out.println("Ud = " + Ud);
596
597 mon = mon.multiply(xv); // Power.<GenPolynomial<MOD>> power(pkfac,xv,e);
598 //System.out.println("mon = " + mon);
599 List<GenPolynomial<MOD>> Sd = new ArrayList<GenPolynomial<MOD>>(Ud.size());
600 int i = 0;
601 List<GenPolynomial<BigInteger>> Si = new ArrayList<GenPolynomial<BigInteger>>(Ud.size());
602 for (GenPolynomial<MOD> dd : Ud) {
603 //System.out.println("dd = " + dd);
604 GenPolynomial<MOD> de = dd.extend(pkfac, 0, 0L);
605 GenPolynomial<MOD> dm = de.multiply(mon);
606 Sd.add(dm);
607 de = U.get(i).sum(dm);
608 //System.out.println("de = " + de);
609 U.set(i++, de);
610 GenPolynomial<BigInteger> si = PolyUtil.integerFromModularCoefficients(ifac, de);
611 Si.add(si);
612 }
613 //System.out.println("Sd = " + Sd);
614 //System.out.println("U = " + U);
615 //System.out.println("Si = " + Si);
616
617 // compute new error:
618 E = ifac.getONE();
619 for (GenPolynomial<BigInteger> bi : Si) {
620 E = E.multiply(bi);
621 }
622 E = Ci.subtract(E);
623 //System.out.println("E = " + E);
624 Ep = PolyUtil.<MOD> fromIntegerCoefficients(pkfac, E);
625 //System.out.println("Ep(0," + pkfac.nvar + ") = " + Ep);
626 logger.info("Ep(" + e + "," + deg + "," + pkfac.nvar + ") = " + Ep);
627 }
628 Vh.add(v);
629 U1 = U;
630 if (Pfac.size() > 0) {
631 List<GenPolynomial<MOD>> U2 = new ArrayList<GenPolynomial<MOD>>(U.size());
632 pkfac = Pfac.get(0);
633 for (GenPolynomial<MOD> b : U) {
634 GenPolynomial<MOD> bi = b.extend(pkfac, 0, 0L);
635 U2.add(bi);
636 }
637 U = U2;
638 //System.out.println("U = " + U);
639 }
640 }
641 if (E.isZERO()) {
642 logger.info("liftHensel leaving with zero E");
643 }
644 return U;
645 }
646
647
648 /**
649 * Modular Hensel lifting algorithm. Let p = A_i.ring.coFac.modul() and
650 * assume ggt(a,b) == 1 mod p, for a, b in A.
651 * @param C GenPolynomial with integer coefficients
652 * @param Cp GenPolynomial C mod p^k
653 * @param F list of modular GenPolynomials, mod (I_v, p^k )
654 * @param V list of substitution values, mod p^k
655 * @param k desired approximation exponent p^k.
656 * @param G list of leading coefficients of the factors of C.
657 * @return [g'_1,..., g'_n] with prod_i g'_i = Cp mod p^k.
658 */
659 public static <MOD extends GcdRingElem<MOD> & Modular> List<GenPolynomial<MOD>> liftHensel(
660 GenPolynomial<BigInteger> C, GenPolynomial<MOD> Cp, List<GenPolynomial<MOD>> F,
661 List<MOD> V, long k, List<GenPolynomial<BigInteger>> G) throws NoLiftingException {
662 GenPolynomialRing<MOD> pkfac = Cp.ring;
663 long d = C.degree();
664 //System.out.println("C = " + C);
665 //System.out.println("Cp = " + Cp);
666 //System.out.println("G = " + G);
667
668 //GenPolynomial<BigInteger> cd = G.get(0); // 1
669 //System.out.println("cd = " + cd + ", ring = " + C.ring);
670 //if ( cd.equals(C.ring.univariate(0)) ) {
671 // System.out.println("cd == G[1]");
672 //}
673 // G mod p^k, in all variables
674 GenPolynomialRing<MOD> pkfac1 = new GenPolynomialRing<MOD>(pkfac.coFac, G.get(0).ring);
675 List<GenPolynomial<MOD>> Lp = new ArrayList<GenPolynomial<MOD>>(G.size());
676 for ( GenPolynomial<BigInteger> cd1 : G ) {
677 GenPolynomial<MOD> cdq = PolyUtil.<MOD> fromIntegerCoefficients(pkfac1, cd1);
678 cdq = cdq.extendLower(pkfac,0,0L); // reintroduce lower variable
679 Lp.add(cdq);
680 }
681 logger.info("G modulo p^k: " + Lp); // + ", ring = " + pkfac1);
682
683 // prepare stack of polynomial rings, polynomials and evaluated leading coefficients
684 List<GenPolynomialRing<MOD>> Pfac = new ArrayList<GenPolynomialRing<MOD>>();
685 List<GenPolynomial<MOD>> Ap = new ArrayList<GenPolynomial<MOD>>();
686 List<List<GenPolynomial<MOD>>> Gp = new ArrayList<List<GenPolynomial<MOD>>>();
687 List<MOD> Vb = new ArrayList<MOD>();
688 //MOD v = V.get(0);
689 Pfac.add(pkfac);
690 Ap.add(Cp);
691 Gp.add(Lp);
692 GenPolynomialRing<MOD> pf = pkfac;
693 GenPolynomialRing<MOD> pf1 = pkfac1;
694 GenPolynomial<MOD> ap = Cp;
695 List<GenPolynomial<MOD>> Lpp = Lp;
696 for (int j = pkfac.nvar; j > 2; j--) {
697 pf = pf.contract(1);
698 Pfac.add(0, pf);
699 MOD vp = pkfac.coFac.fromInteger(V.get(pkfac.nvar - j).getSymmetricInteger().getVal());
700 //System.out.println("vp = " + vp);
701 Vb.add(vp);
702 ap = PolyUtil.<MOD> evaluateMain(pf, ap, vp);
703 Ap.add(0, ap);
704 List<GenPolynomial<MOD>> Lps = new ArrayList<GenPolynomial<MOD>>(Lpp.size());
705 for ( GenPolynomial<MOD> qp : Lpp ) {
706 GenPolynomial<MOD> qpe = PolyUtil.<MOD> evaluateMain(pf, qp, vp);
707 Lps.add(qpe);
708 }
709 //System.out.println("Lps = " + Lps);
710 Lpp = Lps;
711 Gp.add(0, Lpp);
712 }
713 Vb.add( V.get(pkfac.nvar-2) );
714 //System.out.println("Pfac = " + Pfac);
715 if (debug) {
716 logger.debug("Pfac = " + Pfac);
717 }
718 //System.out.println("Ap = " + Ap);
719 //System.out.println("Gp = " + Gp);
720 //System.out.println("Gp[0] = " + Gp.get(0) + ", Gp[0].ring = " + Gp.get(0).get(0).ring);
721 //System.out.println("V = " + V);
722 //System.out.println("Vb = " + Vb + ", V == Vb: " + V.equals(Vb));
723
724 // check bi-variate base case
725 GenPolynomialRing<MOD> pk1fac = F.get(0).ring;
726 if (!pkfac.coFac.equals(pk1fac.coFac)) {
727 throw new IllegalArgumentException("F.ring != pkfac: " + pk1fac + " != " + pkfac);
728 }
729
730 // init recursion
731 List<GenPolynomial<MOD>> U = F;
732 //logger.info("to lift U = " + U); // + ", U1.ring = " + U1.get(0).ring);
733 GenPolynomial<BigInteger> E = C.ring.getZERO();
734 List<MOD> Vh = new ArrayList<MOD>();
735 List<GenPolynomial<BigInteger>> Si = new ArrayList<GenPolynomial<BigInteger>>(F.size());
736 MOD v = null;
737
738 while (Pfac.size() > 0) { // loop through stack of polynomial rings
739 pkfac = Pfac.remove(0);
740 Cp = Ap.remove(0);
741 Lpp = Gp.remove(0);
742 v = Vb.remove(Vb.size()-1); // last in stack
743 //System.out.println("\npkfac = " + pkfac.toScript() + " ================================== " + v);
744 logger.info("stack loop: pkfac = " + pkfac.toScript() + " v = " + v);
745
746 List<GenPolynomial<MOD>> U1 = U;
747 logger.info("to lift U1 = " + U1); // + ", U1.ring = " + U1.get(0).ring);
748 U = new ArrayList<GenPolynomial<MOD>>(U1.size());
749
750 // update U, replace leading coefficient if required
751 int j = 0;
752 for (GenPolynomial<MOD> b : U1) {
753 //System.out.println("b = " + b + ", b.ring = " + b.ring);
754 GenPolynomial<MOD> bi = b.extend(pkfac, 0, 0L);
755 GenPolynomial<MOD> li = Lpp.get(j);
756 if ( !li.isONE() ) {
757 //System.out.println("li = " + li + ", li.ring = " + li.ring);
758 //System.out.println("bi = " + bi);
759 GenPolynomialRing<GenPolynomial<MOD>> pkrfac = pkfac.recursive(pkfac.nvar-1);
760 //System.out.println("pkrfac = " + pkrfac);
761 GenPolynomial<GenPolynomial<MOD>> br = PolyUtil.<MOD>recursive(pkrfac,bi);
762 //System.out.println("br = " + br);
763 GenPolynomial<GenPolynomial<MOD>> bs = PolyUtil.<MOD>switchVariables(br);
764 //System.out.println("bs = " + bs + ", bs.ring = " + bs.ring);
765
766 GenPolynomial<GenPolynomial<MOD>> lr = PolyUtil.<MOD>recursive(pkrfac,li);
767 //System.out.println("lr = " + lr);
768 GenPolynomial<GenPolynomial<MOD>> ls = PolyUtil.<MOD>switchVariables(lr);
769 //System.out.println("ls = " + ls + ", ls.ring = " + ls.ring);
770 if ( !ls.isConstant() ) {
771 throw new RuntimeException("ls not constant " + ls);
772 }
773 bs.doPutToMap(bs.leadingExpVector(),ls.leadingBaseCoefficient());
774 //System.out.println("bs = " + bs + ", bs.ring = " + bs.ring);
775 br = PolyUtil.<MOD>switchVariables(bs);
776 //System.out.println("br = " + br);
777 bi = PolyUtil.<MOD>distribute(pkfac,br);
778 //System.out.println("bi = " + bi);
779 }
780 U.add(bi);
781 j++;
782 }
783 logger.info("U with leading coefficient replaced = " + U); // + ", U.ring = " + U.get(0).ring);
784
785 // (x_n - v)
786 GenPolynomial<MOD> mon = pkfac.getONE();
787 GenPolynomial<MOD> xv = pkfac.univariate(0, 1);
788 xv = xv.subtract(pkfac.fromInteger(v.getSymmetricInteger().getVal()));
789 //System.out.println("xv = " + xv);
790
791 long deg = Cp.degree(pkfac.nvar - 1);
792 //System.out.println("deg = " + deg + ", degv = " + Cp.degreeVector());
793
794 // convert to integer polynomials
795 GenPolynomialRing<BigInteger> ifac = new GenPolynomialRing<BigInteger>(new BigInteger(), pkfac);
796 //System.out.println("ifac = " + ifac.toScript());
797 List<GenPolynomial<BigInteger>> Bi = PolyUtil.integerFromModularCoefficients(ifac, U);
798 //System.out.println("Bi = " + Bi);
799 GenPolynomial<BigInteger> Ci = PolyUtil.integerFromModularCoefficients(ifac, Cp);
800 //System.out.println("Ci = " + Ci);
801
802 // compute error:
803 E = ifac.getONE();
804 for (GenPolynomial<BigInteger> bi : Bi) {
805 E = E.multiply(bi);
806 }
807 //System.out.println("E = " + E);
808 E = Ci.subtract(E);
809 //System.out.println("E = " + E);
810 GenPolynomial<MOD> Ep = PolyUtil.<MOD> fromIntegerCoefficients(pkfac, E);
811 logger.info("Ep(0," + deg + "," + pkfac.nvar + ") = " + Ep);
812
813 GenPolynomialRing<GenPolynomial<MOD>> pkrfac = pkfac.recursive(1);
814 GenPolynomialRing<MOD> ckfac = (GenPolynomialRing<MOD>) pkrfac.coFac;
815 //System.out.println("pkrfac = " + pkrfac.toScript());
816
817 for (int e = 1; e <= deg && !Ep.isZERO(); e++) {
818 //System.out.println("\ne = " + e + " -------------------------------------- " + deg);
819 logger.info("approximation loop: e = " + e + " of deg = " + deg);
820 GenPolynomial<GenPolynomial<MOD>> Epr = PolyUtil.<MOD> recursive(pkrfac, Ep);
821 //System.out.println("Epr = " + Epr);
822 UnivPowerSeriesRing<GenPolynomial<MOD>> psfac
823 = new UnivPowerSeriesRing<GenPolynomial<MOD>>(pkrfac);
824 //System.out.println("psfac = " + psfac);
825 TaylorFunction<GenPolynomial<MOD>> T = new PolynomialTaylorFunction<GenPolynomial<MOD>>(Epr);
826 //System.out.println("T = " + T);
827 GenPolynomial<MOD> vq = ckfac.fromInteger(v.getSymmetricInteger().getVal());
828 //System.out.println("vq = " + vq + ", Vh = " + Vh);
829 UnivPowerSeries<GenPolynomial<MOD>> Epst = psfac.seriesOfTaylor(T, vq);
830 //System.out.println("Epst = " + Epst);
831 logger.info("Epst(" + e + "," + deg + "," + pkfac.nvar + ") = " + Epst);
832 GenPolynomial<MOD> cm = Epst.coefficient(e);
833 if (cm.isZERO()) {
834 //System.out.println("cm = " + cm);
835 continue;
836 }
837 List<GenPolynomial<MOD>> Ud = HenselMultUtil.<MOD> liftDiophant(U1, cm, Vh, d, k);
838 //System.out.println("Ud = " + Ud);
839
840 mon = mon.multiply(xv); // Power.<GenPolynomial<MOD>> power(pkfac,xv,e);
841 //System.out.println("mon = " + mon);
842 List<GenPolynomial<MOD>> Sd = new ArrayList<GenPolynomial<MOD>>(Ud.size());
843 int i = 0;
844 Si = new ArrayList<GenPolynomial<BigInteger>>(Ud.size());
845 for (GenPolynomial<MOD> dd : Ud) {
846 //System.out.println("dd = " + dd);
847 GenPolynomial<MOD> de = dd.extend(pkfac, 0, 0L);
848 GenPolynomial<MOD> dm = de.multiply(mon);
849 Sd.add(dm);
850 de = U.get(i).sum(dm);
851 //System.out.println("de = " + de);
852 U.set(i++, de);
853 GenPolynomial<BigInteger> si = PolyUtil.integerFromModularCoefficients(ifac, de);
854 Si.add(si);
855 }
856 //System.out.println("Sd = " + Sd);
857 //System.out.println("U = " + U + ", U.ring = " + U.get(0).ring);
858 //System.out.println("Si = " + Si);
859
860 // compute new error:
861 E = ifac.getONE();
862 for (GenPolynomial<BigInteger> bi : Si) {
863 E = E.multiply(bi);
864 }
865 E = Ci.subtract(E);
866 //System.out.println("E = " + E);
867 Ep = PolyUtil.<MOD> fromIntegerCoefficients(pkfac, E);
868 //System.out.println("Ep(0," + pkfac.nvar + ") = " + Ep);
869 logger.info("Ep(" + e + "," + deg + "," + pkfac.nvar + ") = " + Ep);
870 }
871 Vh.add(v);
872 GenPolynomial<MOD> Uf = U.get(0).ring.getONE();
873 for ( GenPolynomial<MOD> Upp : U ) {
874 Uf = Uf.multiply(Upp);
875 }
876 if ( false && !Cp.equals(Uf) ) {
877 System.out.println("\nU = " + U);
878 System.out.println("Cp = " + Cp);
879 System.out.println("Uf = " + Uf);
880 System.out.println("isFactorization: " + Cp.equals(Uf) + ", Cp.ring = " + Cp.ring + ", Uf.ring = " + Uf.ring);
881 throw new RuntimeException("no factorization, something is wrong");
882 }
883 }
884 if (E.isZERO()) {
885 logger.info("liftHensel leaving with zero E, Ep");
886 }
887 if ( false && debug ) {
888 // remove normalization required ??
889 GreatestCommonDivisorAbstract<BigInteger> ufd = GCDFactory.getImplementation(new BigInteger());
890 List<GenPolynomial<BigInteger>> Fii = new ArrayList<GenPolynomial<BigInteger>>(U.size());
891 for ( GenPolynomial<BigInteger> bi : Si ) {
892 GenPolynomial<BigInteger> ci = ufd.content(bi); //ufd.primitivePart(bi); // ??
893 if ( !ci.isONE() ) {
894 System.out.println("bi = " + bi + ", cont(bi) = " + ci);
895 }
896 //Fii.add(ci);
897 }
898 //Si = Fii;
899 //System.out.println("Si = " + Si);
900 }
901 logger.info("multivariate lift = " + U + ", of " + F);
902 return U;
903 }
904
905
906
907 /**
908 * Modular Hensel full lifting algorithm. Let p = A_i.ring.coFac.modul() and
909 * assume ggt(a,b) == 1 mod p, for a, b in A.
910 * @param C GenPolynomial with integer coefficients
911 * @param F list of modular GenPolynomials, mod (I_v, p )
912 * @param V list of substitution values, mod p^k
913 * @param k desired approximation exponent p^k.
914 * @param G = [g_1,...,g_n] list of factors of leading coefficients.
915 * @return [c_1,..., c_n] with prod_i c_i = C mod p^k.
916 */
917 public static <MOD extends GcdRingElem<MOD> & Modular>
918 List<GenPolynomial<MOD>> liftHenselFull(
919 GenPolynomial<BigInteger> C, List<GenPolynomial<MOD>> F, List<MOD> V, long k,
920 List<GenPolynomial<BigInteger>> G) throws NoLiftingException {
921 if (F == null || F.size() == 0) {
922 return new ArrayList<GenPolynomial<MOD>>();
923 }
924 GenPolynomialRing<MOD> pkfac = F.get(0).ring;
925 long d = C.degree();
926 // setup q = p^k
927 RingFactory<MOD> cfac = pkfac.coFac;
928 ModularRingFactory<MOD> pcfac = (ModularRingFactory<MOD>) cfac;
929 //System.out.println("pcfac = " + pcfac);
930 BigInteger p = pcfac.getIntegerModul();
931 BigInteger q = Power.positivePower(p, k);
932 ModularRingFactory<MOD> mcfac;
933 if (ModLongRing.MAX_LONG.compareTo(q.getVal()) > 0) {
934 mcfac = (ModularRingFactory) new ModLongRing(q.getVal());
935 } else {
936 mcfac = (ModularRingFactory) new ModIntegerRing(q.getVal());
937 }
938 //System.out.println("mcfac = " + mcfac);
939
940 // convert C from Z[...] to Z_q[...]
941 GenPolynomialRing<MOD> qcfac = new GenPolynomialRing<MOD>(mcfac, C.ring);
942 GenPolynomial<MOD> Cq = PolyUtil.<MOD> fromIntegerCoefficients(qcfac, C);
943 //System.out.println("C = " + C);
944 //System.out.println("Cq = " + Cq);
945
946 // convert g_i from Z[...] to Z_q[...]
947 GenPolynomialRing<MOD> gcfac = new GenPolynomialRing<MOD>(mcfac, G.get(0).ring);
948 List<GenPolynomial<MOD>> GQ = new ArrayList<GenPolynomial<MOD>>();
949 boolean allOnes = true;
950 for ( GenPolynomial<BigInteger> g : G ) {
951 if ( !g.isONE() ) {
952 allOnes = false;
953 }
954 GenPolynomial<MOD> gq = PolyUtil.<MOD> fromIntegerCoefficients(gcfac, g);
955 GQ.add(gq);
956 }
957 //System.out.println("G = " + G);
958 //System.out.println("GQ = " + GQ);
959
960 // evaluate C to Z_q[x]
961 GenPolynomialRing<MOD> pf = qcfac;
962 GenPolynomial<MOD> ap = Cq;
963 for (int j = C.ring.nvar; j > 1; j--) {
964 pf = pf.contract(1);
965 MOD vp = mcfac.fromInteger(V.get(C.ring.nvar-j).getSymmetricInteger().getVal());
966 //System.out.println("vp = " + vp);
967 ap = PolyUtil.<MOD> evaluateMain(pf, ap, vp);
968 //System.out.println("ap = " + ap);
969 }
970 GenPolynomial<MOD> Cq1 = ap;
971 //System.out.println("Cq1 = " + Cq1);
972 if (Cq1.isZERO()) {
973 throw new NoLiftingException("C mod (I, p^k) == 0: " + C);
974 }
975 GenPolynomialRing<BigInteger> ifac = new GenPolynomialRing<BigInteger>(new BigInteger(), pf);
976 GenPolynomial<BigInteger> Ci = PolyUtil.integerFromModularCoefficients(ifac, Cq1);
977 //System.out.println("Ci = " + Ci);
978 GreatestCommonDivisorAbstract<BigInteger> ufd = GCDFactory.getImplementation(new BigInteger());
979 Ci = Ci.abs();
980 BigInteger cCi = ufd.baseContent(Ci);
981 Ci = Ci.divide(cCi);
982 //System.out.println("cCi = " + cCi);
983 //System.out.println("Ci = " + Ci);
984 ////System.out.println("F.fac = " + F.get(0).ring);
985
986 // evaluate G to Z_q
987 List<GenPolynomial<MOD>> GP = new ArrayList<GenPolynomial<MOD>>();
988 for ( GenPolynomial<MOD> gq : GQ ) {
989 GenPolynomialRing<MOD> gf = gcfac;
990 GenPolynomial<MOD> gp = gq;
991 for (int j = gcfac.nvar; j > 1; j--) {
992 gf = gf.contract(1);
993 MOD vp = mcfac.fromInteger(V.get(gcfac.nvar-j).getSymmetricInteger().getVal());
994 //System.out.println("vp = " + vp);
995 gp = PolyUtil.<MOD> evaluateMain(gf, gp, vp);
996 //System.out.println("gp = " + gp);
997 }
998 GP.add(gp);
999 }
1000 //System.out.println("GP = " + GP); // + ", GP.ring = " + GP.get(0).ring);
1001
1002 // leading coefficient for recursion base, for Cq1 and list GP
1003 BigInteger gi0 = Ci.leadingBaseCoefficient(); // gq0.getSymmetricInteger();
1004 //System.out.println("gi0 = " + gi0);
1005
1006 // lift F to Z_{p^k}[x]
1007 //System.out.println("Ci = " + Ci + ", F = " + F + ", k = " + k + ", p = " + F.get(0).ring + ", gi0 = " + gi0);
1008 List<GenPolynomial<MOD>> U1 = null;
1009 if ( gi0.isONE() ) {
1010 U1 = HenselUtil.<MOD> liftHenselMonic(Ci, F, k);
1011 } else {
1012 U1 = HenselUtil.<MOD> liftHensel(Ci, F, k, gi0); // GI0 TODO ??
1013 }
1014 logger.info("univariate lift: Ci = " + Ci + ", F = " + F + ", U1 = " + U1);
1015 //System.out.println("U1.fac = " + U1.get(0).ring);
1016
1017 // adjust leading coefficients of U1 with F
1018 List<GenPolynomial<BigInteger>> U1i = PolyUtil.<MOD> integerFromModularCoefficients(Ci.ring, U1);
1019 //System.out.println("U1i = " + U1i);
1020 boolean t = HenselUtil.isHenselLift(Ci, q, p, U1i);
1021 //System.out.println("isLift(U1) = " + t);
1022 if ( ! t ) {
1023 //System.out.println("NoLiftingException, Ci = " + Ci + ", U1i = " + U1i);
1024 throw new NoLiftingException("Ci = " + Ci + ", U1i = " + U1i);
1025 }
1026 MOD cC = mcfac.fromInteger(cCi.getVal());
1027 List<GenPolynomial<MOD>> U1f = PolyUtil.<MOD> fromIntegerCoefficients(F.get(0).ring, U1i);
1028 //System.out.println("U1f = " + U1f);
1029 List<GenPolynomial<MOD>> U1s = new ArrayList<GenPolynomial<MOD>>(U1.size());
1030 int j = 0;
1031 int s = 0;
1032 for ( GenPolynomial<MOD> u : U1 ) {
1033 GenPolynomial<MOD> uf = U1f.get(j);
1034 GenPolynomial<MOD> f = F.get(j);
1035 GenPolynomial<BigInteger> ui = U1i.get(j);
1036 GenPolynomial<BigInteger> gi = G.get(j);
1037 if ( ui.signum() != gi.signum() ) {
1038 //System.out.println("ui = " + ui + ", gi = " + gi);
1039 u = u.negate();
1040 uf = uf.negate();
1041 s++;
1042 }
1043 j++;
1044 if ( uf.isConstant() ) {
1045 //System.out.println("u = " + u);
1046 u = u.monic();
1047 //System.out.println("u = " + u);
1048 u = u.multiply(cC);
1049 cC = cC.divide(cC);
1050 //System.out.println("u = " + u);
1051 } else {
1052 MOD x = f.leadingBaseCoefficient().divide(uf.leadingBaseCoefficient());
1053 //System.out.println("x = " + x + ", xi = " + x.getSymmetricInteger());
1054 if ( !x.isONE() ) {
1055 MOD xq = mcfac.fromInteger(x.getSymmetricInteger().getVal());
1056 //System.out.println("xq = " + xq);
1057 u = u.multiply(xq);
1058 cC = cC.divide(xq);
1059 //System.out.println("cC = " + cC);
1060 }
1061 }
1062 U1s.add(u);
1063 }
1064 //if ( s % 2 != 0 || !cC.isONE()) {
1065 if ( !cC.isONE()) {
1066 throw new NoLiftingException("s = " + s + ", Ci = " + Ci + ", U1i = " + U1i + ", cC = " + cC);
1067 }
1068 U1 = U1s;
1069 U1i = PolyUtil.<MOD> integerFromModularCoefficients(Ci.ring, U1);
1070 //System.out.println("U1i = " + U1i);
1071 U1f = PolyUtil.<MOD> fromIntegerCoefficients(F.get(0).ring, U1i);
1072 if ( !F.equals(U1f) ) { // evtl loop until reached
1073 System.out.println("F = " + F);
1074 System.out.println("U1f = " + U1f);
1075 throw new NoLiftingException("F = " + F + ", U1f = " + U1f);
1076 }
1077 logger.info("multivariate lift: U1 = " + U1);
1078
1079 // lift U to Z_{p^k}[x,...]
1080 //System.out.println("C = " + C + ", U1 = " + U1 + ", V = " + V + ", k = " + k + ", q = " + U1.get(0).ring + ", G = " + G);
1081 List<GenPolynomial<MOD>> U = null;
1082 if ( allOnes ) {
1083 U = HenselMultUtil.<MOD> liftHenselMonic(C, Cq, U1, V, k);
1084 } else {
1085 U = HenselMultUtil.<MOD> liftHensel(C, Cq, U1, V, k, G);
1086 }
1087 logger.info("multivariate lift: C = " + C + ", U1 = " + U1 + ", U = " + U);
1088 //System.out.println("U = " + U);
1089 //System.out.println("U.fac = " + U.get(0).ring);
1090 return U;
1091 }
1092
1093 }