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