001/* 002 * $Id: HenselMultUtil.java 6007 2020-03-29 13:34:49Z kredel $ 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 //System.out.println("Ep(0," + pkfac.nvar + ") = " + Ep); 136 logger.info("Ep(0," + pkfac.nvar + ") = " + Ep); 137 if (Ep.isZERO()) { 138 logger.info("liftDiophant leaving on zero Ep mod p^k"); 139 return sup; 140 } 141 for (int e = 1; e <= d; e++) { 142 //System.out.println("\ne = " + e + " -------------------------------------- " + pkfac.nvar); 143 GenPolynomial<GenPolynomial<MOD>> Epr = PolyUtil.<MOD> recursive(qrfac, Ep); 144 //System.out.println("Epr = " + Epr); 145 UnivPowerSeriesRing<GenPolynomial<MOD>> psfac = new UnivPowerSeriesRing<GenPolynomial<MOD>>( 146 qrfac); 147 //System.out.println("psfac = " + psfac); 148 TaylorFunction<GenPolynomial<MOD>> F = new PolynomialTaylorFunction<GenPolynomial<MOD>>(Epr); 149 //System.out.println("F = " + F); 150 //List<GenPolynomial<MOD>> Vs = new ArrayList<GenPolynomial<MOD>>(1); 151 GenPolynomial<MOD> vq = ckfac.fromInteger(v.getSymmetricInteger().getVal()); 152 //Vs.add(vq); 153 //System.out.println("Vs = " + Vs); 154 UnivPowerSeries<GenPolynomial<MOD>> Epst = psfac.seriesOfTaylor(F, vq); 155 //System.out.println("Epst = " + Epst); 156 GenPolynomial<MOD> cm = Epst.coefficient(e); 157 //System.out.println("cm = " + cm + ", cm.ring = " + cm.ring.toScript()); 158 159 // recursion: 160 List<GenPolynomial<MOD>> S = HenselMultUtil.<MOD> liftDiophant(Ap, Bp, cm, Vp, d, k); 161 //System.out.println("S = " + S); 162 if (!ckfac.coFac.equals(S.get(0).ring.coFac)) { 163 throw new IllegalArgumentException( 164 "ckfac != pkfac: " + ckfac.coFac + " != " + S.get(0).ring.coFac); 165 } 166 if (pkfac.nvar == 2 && !HenselUtil.<MOD> isDiophantLift(Ap, Bp, S.get(1), S.get(0), cm)) { 167 //System.out.println("isDiophantLift: false"); 168 throw new NoLiftingException("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 throw new NoLiftingException("isDiophantLift: false"); 272 } 273 if (!ckfac.equals(su.get(0).ring)) { 274 throw new IllegalArgumentException("qfac != ckfac: " + su.get(0).ring + " != " + ckfac); 275 } 276 GenPolynomialRing<BigInteger> ifac = new GenPolynomialRing<BigInteger>(new BigInteger(), pkfac); 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>>( 336 qrfac); 337 //System.out.println("psfac = " + psfac); 338 TaylorFunction<GenPolynomial<MOD>> F = new PolynomialTaylorFunction<GenPolynomial<MOD>>(Epr); 339 //System.out.println("F = " + F); 340 //List<GenPolynomial<MOD>> Vs = new ArrayList<GenPolynomial<MOD>>(1); 341 GenPolynomial<MOD> vq = ckfac.fromInteger(v.getSymmetricInteger().getVal()); 342 //Vs.add(vq); 343 //System.out.println("Vs = " + Vs); 344 UnivPowerSeries<GenPolynomial<MOD>> Epst = psfac.seriesOfTaylor(F, vq); 345 //System.out.println("Epst = " + Epst); 346 GenPolynomial<MOD> cm = Epst.coefficient(e); 347 //System.out.println("cm = " + cm + ", cm.ring = " + cm.ring.toScript()); 348 if (cm.isZERO()) { 349 continue; 350 } 351 // recursion: 352 List<GenPolynomial<MOD>> S = HenselMultUtil.<MOD> liftDiophant(Ap, cm, Vp, d, k); 353 //System.out.println("S = " + S); 354 if (!ckfac.coFac.equals(S.get(0).ring.coFac)) { 355 throw new IllegalArgumentException( 356 "ckfac != pkfac: " + ckfac.coFac + " != " + S.get(0).ring.coFac); 357 } 358 if (pkfac.nvar == 2 && !HenselUtil.<MOD> isDiophantLift(Ap, S, cm)) { 359 //System.out.println("isDiophantLift: false"); 360 throw new NoLiftingException("isDiophantLift: false"); 361 } 362 mon = mon.multiply(xv); // Power.<GenPolynomial<MOD>> power(pkfac,xv,e); 363 //System.out.println("mon = " + mon); 364 //List<GenPolynomial<MOD>> Sp = new ArrayList<GenPolynomial<MOD>>(S.size()); 365 i = 0; 366 supi = new ArrayList<GenPolynomial<BigInteger>>(su.size()); 367 for (GenPolynomial<MOD> dd : S) { 368 //System.out.println("dd = " + dd); 369 GenPolynomial<MOD> de = dd.extend(pkfac, 0, 0L); 370 GenPolynomial<MOD> dm = de.multiply(mon); 371 //Sp.add(dm); 372 de = sup.get(i).sum(dm); 373 //System.out.println("dd = " + dd); 374 sup.set(i++, de); 375 GenPolynomial<BigInteger> spi = PolyUtil.integerFromModularCoefficients(ifac, dm); 376 supi.add(spi); 377 } 378 //System.out.println("Sp = " + Sp); 379 //System.out.println("sup = " + sup); 380 //System.out.println("supi = " + supi); 381 // compute new error 382 //E = E; // - sum_i s_i b_i 383 i = 0; 384 for (GenPolynomial<BigInteger> bi : Bi) { 385 E = E.subtract(bi.multiply(supi.get(i++))); 386 } 387 //System.out.println("E = " + E); 388 if (E.isZERO()) { 389 logger.info("liftDiophant leaving on zero E"); 390 return sup; 391 } 392 Ep = PolyUtil.<MOD> fromIntegerCoefficients(pkfac, E); 393 //System.out.println("Ep(" + e + "," + pkfac.nvar + ") = " + Ep); 394 logger.info("Ep(" + e + "," + pkfac.nvar + ") = " + Ep); 395 if (Ep.isZERO()) { 396 logger.info("liftDiophant leaving on zero Ep mod p^k"); 397 return sup; 398 } 399 } 400 //System.out.println("*** done: " + pkfac.nvar); 401 return sup; 402 } 403 404 405 /** 406 * Modular Hensel lifting algorithm on coefficients test. Let p = 407 * f_i.ring.coFac.modul() and assume C == prod_{0,...,n-1} f_i mod p with 408 * gcd(f_i,f_j) == 1 mod p for i != j 409 * @param C integer polynomial 410 * @param Cp GenPolynomial mod p^k 411 * @param F = [f_0,...,f_{n-1}] list of monic modular polynomials. 412 * @param L = [g_0,...,g_{n-1}] list of lifted modular polynomials. 413 * @return true if C = prod_{0,...,n-1} g_i mod p^k, else false. 414 */ 415 @SuppressWarnings("unused") 416 public static <MOD extends GcdRingElem<MOD> & Modular> boolean isHenselLift(GenPolynomial<BigInteger> C, 417 GenPolynomial<MOD> Cp, List<GenPolynomial<MOD>> F, List<GenPolynomial<MOD>> L) { 418 boolean t = true; 419 GenPolynomialRing<MOD> qfac = L.get(0).ring; 420 GenPolynomial<MOD> q = qfac.getONE(); 421 for (GenPolynomial<MOD> fi : L) { 422 q = q.multiply(fi); 423 } 424 t = Cp.equals(q); 425 if (!t) { 426 System.out.println("Cp = " + Cp); 427 System.out.println("q = " + q); 428 System.out.println("Cp != q: " + Cp.subtract(q)); 429 return t; 430 } 431 GenPolynomialRing<BigInteger> dfac = C.ring; 432 GenPolynomial<BigInteger> Ci = PolyUtil.integerFromModularCoefficients(dfac, q); 433 t = C.equals(Ci); 434 if (!t) { 435 System.out.println("C = " + C); 436 System.out.println("Ci = " + Ci); 437 System.out.println("C != Ci: " + C.subtract(Ci)); 438 return t; 439 } 440 // test L mod id(V) == F 441 return t; 442 } 443 444 445 /** 446 * Modular Hensel lifting algorithm, monic case. Let p = 447 * A_i.ring.coFac.modul() and assume ggt(a,b) == 1 mod p, for a, b in A. 448 * @param C monic GenPolynomial with integer coefficients 449 * @param Cp GenPolynomial mod p^k 450 * @param F list of modular GenPolynomials, mod (I_v, p^k ) 451 * @param V list of integer substitution values 452 * @param k desired approximation exponent p^k. 453 * @return [g'_1,..., g'_n] with prod_i g'_i = Cp mod p^k. 454 */ 455 public static <MOD extends GcdRingElem<MOD> & Modular> List<GenPolynomial<MOD>> liftHenselMonic( 456 GenPolynomial<BigInteger> C, GenPolynomial<MOD> Cp, List<GenPolynomial<MOD>> F, 457 List<BigInteger> V, long k) throws NoLiftingException { 458 GenPolynomialRing<MOD> pkfac = Cp.ring; 459 //if (pkfac.nvar == 1) { // V ignored 460 // return HenselUtil.<MOD> liftHenselMonic(C,F,k); 461 //} 462 long d = C.degree(); 463 //System.out.println("d = " + d); 464 // prepare stack of polynomial rings and polynomials 465 List<GenPolynomialRing<MOD>> Pfac = new ArrayList<GenPolynomialRing<MOD>>(); 466 List<GenPolynomial<MOD>> Ap = new ArrayList<GenPolynomial<MOD>>(); 467 List<MOD> Vb = new ArrayList<MOD>(); 468 MOD v = pkfac.coFac.fromInteger(V.get(0).getVal()); 469 Pfac.add(pkfac); 470 Ap.add(Cp); 471 Vb.add(v); 472 GenPolynomialRing<MOD> pf = pkfac; 473 GenPolynomial<MOD> ap = Cp; 474 for (int j = pkfac.nvar; j > 2; j--) { 475 pf = pf.contract(1); 476 Pfac.add(0, pf); 477 //MOD vp = pkfac.coFac.fromInteger(V.get(j - 2).getSymmetricInteger().getVal()); 478 MOD vp = pkfac.coFac.fromInteger(V.get(j - 2).getVal()); 479 //System.out.println("vp = " + vp); 480 Vb.add(1, vp); 481 ap = PolyUtil.<MOD> evaluateMain(pf, ap, vp); 482 Ap.add(0, ap); 483 } 484 //System.out.println("Pfac = " + Pfac); 485 if (debug) { 486 logger.debug("Pfac = " + Pfac); 487 } 488 //System.out.println("Ap = " + Ap); 489 //System.out.println("V = " + V); 490 //System.out.println("Vb = " + Vb); 491 // setup bi-variate base case 492 GenPolynomialRing<MOD> pk1fac = F.get(0).ring; 493 if (!pkfac.coFac.equals(pk1fac.coFac)) { 494 throw new IllegalArgumentException("F.ring != pkfac: " + pk1fac + " != " + pkfac); 495 } 496 // TODO: adjust leading coefficients 497 pkfac = Pfac.get(0); 498 //Cp = Ap.get(0); 499 //System.out.println("pkfac = " + pkfac.toScript()); 500 //System.out.println("pk1fac = " + pk1fac.toScript()); 501 GenPolynomialRing<BigInteger> i1fac = new GenPolynomialRing<BigInteger>(new BigInteger(), pk1fac); 502 //System.out.println("i1fac = " + i1fac.toScript()); 503 List<GenPolynomial<BigInteger>> Bi = new ArrayList<GenPolynomial<BigInteger>>(F.size()); 504 for (GenPolynomial<MOD> b : F) { 505 GenPolynomial<BigInteger> bi = PolyUtil.integerFromModularCoefficients(i1fac, b); 506 Bi.add(bi); 507 } 508 //System.out.println("Bi = " + Bi); 509 // evaluate Cp at v_n: 510 //ModularRingFactory<MOD> cf = (ModularRingFactory<MOD>) pkfac.coFac; 511 //MOD vp = cf.fromInteger(v.getSymmetricInteger().getVal()); 512 //System.out.println("v = " + v + ", vp = " + vp); 513 GenPolynomialRing<MOD> ckfac; // = pkfac.contract(1); 514 //GenPolynomial<MOD> Cs = PolyUtil.<MOD> evaluateMain(ckfac, Cp, vp); 515 //System.out.println("Cp = " + Cp); 516 //System.out.println("Cs = " + Cs); 517 518 List<GenPolynomial<MOD>> U = new ArrayList<GenPolynomial<MOD>>(F.size()); 519 for (GenPolynomial<MOD> b : F) { 520 GenPolynomial<MOD> bi = b.extend(pkfac, 0, 0L); 521 U.add(bi); 522 } 523 //System.out.println("U = " + U); 524 List<GenPolynomial<MOD>> U1 = F; 525 //System.out.println("U1 = " + U1); 526 527 GenPolynomial<BigInteger> E = C.ring.getZERO(); 528 List<MOD> Vh = new ArrayList<MOD>(); 529 530 while (Pfac.size() > 0) { // loop through stack of polynomial rings 531 pkfac = Pfac.remove(0); 532 Cp = Ap.remove(0); 533 v = Vb.remove(0); 534 //Vh.add(0,v); 535 //System.out.println("\npkfac = " + pkfac.toScript() + " ================================== " + Vh); 536 537 // (x_n - v) 538 GenPolynomial<MOD> mon = pkfac.getONE(); 539 GenPolynomial<MOD> xv = pkfac.univariate(0, 1); 540 xv = xv.subtract(pkfac.fromInteger(v.getSymmetricInteger().getVal())); 541 //System.out.println("xv = " + xv); 542 543 long deg = Cp.degree(pkfac.nvar - 1); 544 //System.out.println("deg = " + deg); 545 546 GenPolynomialRing<BigInteger> ifac = new GenPolynomialRing<BigInteger>(new BigInteger(), pkfac); 547 //System.out.println("ifac = " + ifac.toScript()); 548 List<GenPolynomial<BigInteger>> Bip = new ArrayList<GenPolynomial<BigInteger>>(F.size()); 549 for (GenPolynomial<BigInteger> b : Bi) { 550 GenPolynomial<BigInteger> bi = b.extend(ifac, 0, 0L); 551 Bip.add(bi); 552 } 553 Bi = Bip; 554 //System.out.println("Bi = " + Bi); 555 GenPolynomial<BigInteger> Ci = PolyUtil.integerFromModularCoefficients(ifac, Cp); 556 //System.out.println("Ci = " + Ci); 557 558 // compute error: 559 E = ifac.getONE(); 560 for (GenPolynomial<BigInteger> bi : Bi) { 561 E = E.multiply(bi); 562 } 563 E = Ci.subtract(E); 564 //System.out.println("E = " + E); 565 GenPolynomial<MOD> Ep = PolyUtil.<MOD> fromIntegerCoefficients(pkfac, E); 566 //System.out.println("Ep(0," + pkfac.nvar + ") = " + Ep); 567 logger.info("Ep(0," + deg + "," + pkfac.nvar + ") = " + Ep); 568 569 String[] mn = new String[] { pkfac.getVars()[pkfac.nvar - 1] }; 570 ckfac = pkfac.contract(1); 571 GenPolynomialRing<GenPolynomial<MOD>> pkrfac = new GenPolynomialRing<GenPolynomial<MOD>>(ckfac, 1, 572 mn); 573 //System.out.println("pkrfac = " + pkrfac.toScript()); 574 575 for (int e = 1; e <= deg && !Ep.isZERO(); e++) { 576 //System.out.println("\ne = " + e + " -------------------------------------- " + pkfac.nvar); 577 GenPolynomial<GenPolynomial<MOD>> Epr = PolyUtil.<MOD> recursive(pkrfac, Ep); 578 //System.out.println("Epr = " + Epr); 579 UnivPowerSeriesRing<GenPolynomial<MOD>> psfac = new UnivPowerSeriesRing<GenPolynomial<MOD>>( 580 pkrfac); 581 //System.out.println("psfac = " + psfac); 582 TaylorFunction<GenPolynomial<MOD>> T = new PolynomialTaylorFunction<GenPolynomial<MOD>>(Epr); 583 //System.out.println("T = " + T); 584 //List<GenPolynomial<MOD>> Vs = new ArrayList<GenPolynomial<MOD>>(1); 585 GenPolynomial<MOD> vq = ckfac.fromInteger(v.getSymmetricInteger().getVal()); 586 //Vs.add(vq); 587 //System.out.println("Vs = " + Vs + ", Vh = " + Vh); 588 UnivPowerSeries<GenPolynomial<MOD>> Epst = psfac.seriesOfTaylor(T, vq); 589 //System.out.println("Epst = " + Epst); 590 logger.info("Epst(" + e + "," + deg + ", " + pkfac.nvar + ") = " + Epst); 591 GenPolynomial<MOD> cm = Epst.coefficient(e); 592 //System.out.println("cm = " + cm); 593 if (cm.isZERO()) { 594 continue; 595 } 596 List<GenPolynomial<MOD>> Ud = HenselMultUtil.<MOD> liftDiophant(U1, cm, Vh, d, k); 597 //System.out.println("Ud = " + Ud); 598 599 mon = mon.multiply(xv); 600 //System.out.println("mon = " + mon); 601 //List<GenPolynomial<MOD>> Sd = new ArrayList<GenPolynomial<MOD>>(Ud.size()); 602 int i = 0; 603 List<GenPolynomial<BigInteger>> Si = new ArrayList<GenPolynomial<BigInteger>>(Ud.size()); 604 for (GenPolynomial<MOD> dd : Ud) { 605 //System.out.println("dd = " + dd); 606 GenPolynomial<MOD> de = dd.extend(pkfac, 0, 0L); 607 GenPolynomial<MOD> dm = de.multiply(mon); 608 //Sd.add(dm); 609 de = U.get(i).sum(dm); 610 //System.out.println("de = " + de); 611 U.set(i++, de); 612 GenPolynomial<BigInteger> si = PolyUtil.integerFromModularCoefficients(ifac, de); 613 Si.add(si); 614 } 615 //System.out.println("Sd = " + Sd); 616 //System.out.println("U = " + U); 617 //System.out.println("Si = " + Si); 618 619 // compute new error: 620 E = ifac.getONE(); 621 for (GenPolynomial<BigInteger> bi : Si) { 622 E = E.multiply(bi); 623 } 624 E = Ci.subtract(E); 625 //System.out.println("E = " + E); 626 Ep = PolyUtil.<MOD> fromIntegerCoefficients(pkfac, E); 627 //System.out.println("Ep(0," + pkfac.nvar + ") = " + Ep); 628 logger.info("Ep(" + e + "," + deg + "," + pkfac.nvar + ") = " + Ep); 629 } 630 Vh.add(v); 631 U1 = U; 632 if (Pfac.size() > 0) { 633 List<GenPolynomial<MOD>> U2 = new ArrayList<GenPolynomial<MOD>>(U.size()); 634 pkfac = Pfac.get(0); 635 for (GenPolynomial<MOD> b : U) { 636 GenPolynomial<MOD> bi = b.extend(pkfac, 0, 0L); 637 U2.add(bi); 638 } 639 U = U2; 640 //System.out.println("U = " + U); 641 } 642 } 643 if (E.isZERO()) { 644 logger.info("liftHensel leaving with zero E"); 645 } 646 return U; 647 } 648 649 650 /** 651 * Modular Hensel lifting algorithm. Let p = A_i.ring.coFac.modul() and 652 * assume ggt(a,b) == 1 mod p, for a, b in A. 653 * @param C GenPolynomial with integer coefficients 654 * @param Cp GenPolynomial C mod p^k 655 * @param F list of modular GenPolynomials, mod (I_v, p^k ) 656 * @param V list of integral substitution values 657 * @param k desired approximation exponent p^k. 658 * @param G list of leading coefficients of the factors of C. 659 * @return [g'_1,..., g'_n] with prod_i g'_i = Cp mod p^k. 660 */ 661 public static <MOD extends GcdRingElem<MOD> & Modular> List<GenPolynomial<MOD>> liftHensel( 662 GenPolynomial<BigInteger> C, GenPolynomial<MOD> Cp, List<GenPolynomial<MOD>> F, 663 List<BigInteger> V, long k, List<GenPolynomial<BigInteger>> G) throws NoLiftingException { 664 GenPolynomialRing<MOD> pkfac = Cp.ring; 665 long d = C.degree(); 666 //System.out.println("C = " + C); 667 //System.out.println("Cp = " + Cp); 668 //System.out.println("G = " + G); 669 670 //GenPolynomial<BigInteger> cd = G.get(0); // 1 671 //System.out.println("cd = " + cd + ", ring = " + C.ring); 672 //if ( cd.equals(C.ring.univariate(0)) ) { 673 // System.out.println("cd == G[1]"); 674 //} 675 // G mod p^k, in all variables 676 GenPolynomialRing<MOD> pkfac1 = new GenPolynomialRing<MOD>(pkfac.coFac, G.get(0).ring); 677 List<GenPolynomial<MOD>> Lp = new ArrayList<GenPolynomial<MOD>>(G.size()); 678 for (GenPolynomial<BigInteger> cd1 : G) { 679 GenPolynomial<MOD> cdq = PolyUtil.<MOD> fromIntegerCoefficients(pkfac1, cd1); 680 cdq = cdq.extendLower(pkfac, 0, 0L); // reintroduce lower variable 681 Lp.add(cdq); 682 } 683 logger.info("G modulo p^k: " + Lp); // + ", ring = " + pkfac1); 684 685 // prepare stack of polynomial rings, polynomials and evaluated leading coefficients 686 List<GenPolynomialRing<MOD>> Pfac = new ArrayList<GenPolynomialRing<MOD>>(); 687 List<GenPolynomial<MOD>> Ap = new ArrayList<GenPolynomial<MOD>>(); 688 List<List<GenPolynomial<MOD>>> Gp = new ArrayList<List<GenPolynomial<MOD>>>(); 689 List<MOD> Vb = new ArrayList<MOD>(); 690 //MOD v = V.get(0); // fromInteger 691 Pfac.add(pkfac); 692 Ap.add(Cp); 693 Gp.add(Lp); 694 GenPolynomialRing<MOD> pf = pkfac; 695 //GenPolynomialRing<MOD> pf1 = pkfac1; 696 GenPolynomial<MOD> ap = Cp; 697 List<GenPolynomial<MOD>> Lpp = Lp; 698 for (int j = pkfac.nvar; j > 2; j--) { 699 pf = pf.contract(1); 700 Pfac.add(0, pf); 701 //MOD vp = pkfac.coFac.fromInteger(V.get(pkfac.nvar - j).getSymmetricInteger().getVal()); 702 MOD vp = pkfac.coFac.fromInteger(V.get(pkfac.nvar - j).getVal()); 703 //System.out.println("vp = " + vp); 704 Vb.add(vp); 705 ap = PolyUtil.<MOD> evaluateMain(pf, ap, vp); 706 Ap.add(0, ap); 707 List<GenPolynomial<MOD>> Lps = new ArrayList<GenPolynomial<MOD>>(Lpp.size()); 708 for (GenPolynomial<MOD> qp : Lpp) { 709 GenPolynomial<MOD> qpe = PolyUtil.<MOD> evaluateMain(pf, qp, vp); 710 Lps.add(qpe); 711 } 712 //System.out.println("Lps = " + Lps); 713 Lpp = Lps; 714 Gp.add(0, Lpp); 715 } 716 Vb.add(pkfac.coFac.fromInteger(V.get(pkfac.nvar - 2).getVal())); 717 //System.out.println("Pfac = " + Pfac); 718 if (debug) { 719 logger.debug("Pfac = " + Pfac); 720 } 721 //System.out.println("Ap = " + Ap); 722 //System.out.println("Gp = " + Gp); 723 //System.out.println("Gp[0] = " + Gp.get(0) + ", Gp[0].ring = " + Gp.get(0).get(0).ring); 724 //System.out.println("V = " + V); 725 //System.out.println("Vb = " + Vb + ", V == Vb: " + V.equals(Vb)); 726 727 // check bi-variate base case 728 GenPolynomialRing<MOD> pk1fac = F.get(0).ring; 729 if (!pkfac.coFac.equals(pk1fac.coFac)) { 730 throw new IllegalArgumentException("F.ring != pkfac: " + pk1fac + " != " + pkfac); 731 } 732 733 // init recursion 734 List<GenPolynomial<MOD>> U = F; 735 //logger.info("to lift U = " + U); // + ", U1.ring = " + U1.get(0).ring); 736 GenPolynomial<BigInteger> E = C.ring.getZERO(); 737 List<MOD> Vh = new ArrayList<MOD>(); 738 List<GenPolynomial<BigInteger>> Si; // = new ArrayList<GenPolynomial<BigInteger>>(F.size()); 739 MOD v = null; 740 741 while (Pfac.size() > 0) { // loop through stack of polynomial rings 742 pkfac = Pfac.remove(0); 743 Cp = Ap.remove(0); 744 Lpp = Gp.remove(0); 745 v = Vb.remove(Vb.size() - 1); // last in stack 746 //System.out.println("\npkfac = " + pkfac.toScript() + " ================================== " + v); 747 logger.info("stack loop: pkfac = " + pkfac.toScript() + " v = " + v); 748 749 List<GenPolynomial<MOD>> U1 = U; 750 logger.info("to lift U1 = " + U1); // + ", U1.ring = " + U1.get(0).ring); 751 U = new ArrayList<GenPolynomial<MOD>>(U1.size()); 752 753 // update U, replace leading coefficient if required 754 int j = 0; 755 for (GenPolynomial<MOD> b : U1) { 756 //System.out.println("b = " + b + ", b.ring = " + b.ring); 757 GenPolynomial<MOD> bi = b.extend(pkfac, 0, 0L); 758 GenPolynomial<MOD> li = Lpp.get(j); 759 if (!li.isONE()) { 760 //System.out.println("li = " + li + ", li.ring = " + li.ring); 761 //System.out.println("bi = " + bi); 762 GenPolynomialRing<GenPolynomial<MOD>> pkrfac = pkfac.recursive(pkfac.nvar - 1); 763 //System.out.println("pkrfac = " + pkrfac); 764 GenPolynomial<GenPolynomial<MOD>> br = PolyUtil.<MOD> recursive(pkrfac, bi); 765 //System.out.println("br = " + br); 766 GenPolynomial<GenPolynomial<MOD>> bs = PolyUtil.<MOD> switchVariables(br); 767 //System.out.println("bs = " + bs + ", bs.ring = " + bs.ring); 768 769 GenPolynomial<GenPolynomial<MOD>> lr = PolyUtil.<MOD> recursive(pkrfac, li); 770 //System.out.println("lr = " + lr); 771 GenPolynomial<GenPolynomial<MOD>> ls = PolyUtil.<MOD> switchVariables(lr); 772 //System.out.println("ls = " + ls + ", ls.ring = " + ls.ring); 773 if (!ls.isConstant() && !ls.isZERO()) { 774 throw new RuntimeException("ls not constant " + ls + ", li = " + li); 775 } 776 bs.doPutToMap(bs.leadingExpVector(), ls.leadingBaseCoefficient()); 777 //System.out.println("bs = " + bs + ", bs.ring = " + bs.ring); 778 br = PolyUtil.<MOD> switchVariables(bs); 779 //System.out.println("br = " + br); 780 bi = PolyUtil.<MOD> distribute(pkfac, br); 781 //System.out.println("bi = " + bi); 782 } 783 U.add(bi); 784 j++; 785 } 786 logger.info("U with leading coefficient replaced = " + U); // + ", U.ring = " + U.get(0).ring); 787 788 // (x_n - v) 789 GenPolynomial<MOD> mon = pkfac.getONE(); 790 GenPolynomial<MOD> xv = pkfac.univariate(0, 1); 791 xv = xv.subtract(pkfac.fromInteger(v.getSymmetricInteger().getVal())); 792 //System.out.println("xv = " + xv); 793 794 long deg = Cp.degree(pkfac.nvar - 1); 795 //System.out.println("deg = " + deg + ", degv = " + Cp.degreeVector()); 796 797 // convert to integer polynomials 798 GenPolynomialRing<BigInteger> ifac = new GenPolynomialRing<BigInteger>(new BigInteger(), pkfac); 799 //System.out.println("ifac = " + ifac.toScript()); 800 List<GenPolynomial<BigInteger>> Bi = PolyUtil.integerFromModularCoefficients(ifac, U); 801 //System.out.println("Bi = " + Bi); 802 GenPolynomial<BigInteger> Ci = PolyUtil.integerFromModularCoefficients(ifac, Cp); 803 //System.out.println("Ci = " + Ci); 804 805 // compute error: 806 E = ifac.getONE(); 807 for (GenPolynomial<BigInteger> bi : Bi) { 808 E = E.multiply(bi); 809 } 810 //System.out.println("E = " + E); 811 E = Ci.subtract(E); 812 //System.out.println("E = " + E); 813 GenPolynomial<MOD> Ep = PolyUtil.<MOD> fromIntegerCoefficients(pkfac, E); 814 logger.info("Ep(0," + deg + "," + pkfac.nvar + ") = " + Ep); 815 816 GenPolynomialRing<GenPolynomial<MOD>> pkrfac = pkfac.recursive(1); 817 GenPolynomialRing<MOD> ckfac = (GenPolynomialRing<MOD>) pkrfac.coFac; 818 //System.out.println("pkrfac = " + pkrfac.toScript()); 819 820 for (int e = 1; e <= deg && !Ep.isZERO(); e++) { 821 //System.out.println("\ne = " + e + " -------------------------------------- " + deg); 822 logger.info("approximation loop: e = " + e + " of deg = " + deg); 823 GenPolynomial<GenPolynomial<MOD>> Epr = PolyUtil.<MOD> recursive(pkrfac, Ep); 824 //System.out.println("Epr = " + Epr); 825 UnivPowerSeriesRing<GenPolynomial<MOD>> psfac = new UnivPowerSeriesRing<GenPolynomial<MOD>>( 826 pkrfac); 827 //System.out.println("psfac = " + psfac); 828 TaylorFunction<GenPolynomial<MOD>> T = new PolynomialTaylorFunction<GenPolynomial<MOD>>(Epr); 829 //System.out.println("T = " + T); 830 GenPolynomial<MOD> vq = ckfac.fromInteger(v.getSymmetricInteger().getVal()); 831 //System.out.println("vq = " + vq + ", Vh = " + Vh); 832 UnivPowerSeries<GenPolynomial<MOD>> Epst = psfac.seriesOfTaylor(T, vq); 833 //System.out.println("Epst = " + Epst); 834 logger.info("Epst(" + e + "," + deg + "," + pkfac.nvar + ") = " + Epst); 835 GenPolynomial<MOD> cm = Epst.coefficient(e); 836 if (cm.isZERO()) { 837 //System.out.println("cm = " + cm); 838 continue; 839 } 840 List<GenPolynomial<MOD>> Ud = HenselMultUtil.<MOD> liftDiophant(U1, cm, Vh, d, k); 841 //System.out.println("Ud = " + Ud); 842 843 mon = mon.multiply(xv); 844 //System.out.println("mon = " + mon); 845 //List<GenPolynomial<MOD>> Sd = new ArrayList<GenPolynomial<MOD>>(Ud.size()); 846 int i = 0; 847 Si = new ArrayList<GenPolynomial<BigInteger>>(Ud.size()); 848 for (GenPolynomial<MOD> dd : Ud) { 849 //System.out.println("dd = " + dd); 850 GenPolynomial<MOD> de = dd.extend(pkfac, 0, 0L); 851 GenPolynomial<MOD> dm = de.multiply(mon); 852 //Sd.add(dm); 853 de = U.get(i).sum(dm); 854 //System.out.println("de = " + de); 855 U.set(i++, de); 856 GenPolynomial<BigInteger> si = PolyUtil.integerFromModularCoefficients(ifac, de); 857 Si.add(si); 858 } 859 //System.out.println("Sd = " + Sd); 860 //System.out.println("U = " + U + ", U.ring = " + U.get(0).ring); 861 //System.out.println("Si = " + Si); 862 863 // compute new error: 864 E = ifac.getONE(); 865 for (GenPolynomial<BigInteger> bi : Si) { 866 E = E.multiply(bi); 867 } 868 E = Ci.subtract(E); 869 //System.out.println("E = " + E); 870 Ep = PolyUtil.<MOD> fromIntegerCoefficients(pkfac, E); 871 //System.out.println("Ep(0," + pkfac.nvar + ") = " + Ep); 872 logger.info("Ep(" + e + "," + deg + "," + pkfac.nvar + ") = " + Ep); 873 } 874 Vh.add(v); 875 GenPolynomial<MOD> Uf = U.get(0).ring.getONE(); 876 for (GenPolynomial<MOD> Upp : U) { 877 Uf = Uf.multiply(Upp); 878 } 879 if (false && !Cp.leadingExpVector().equals(Uf.leadingExpVector())) { // not meanigfull test 880 System.out.println("\nU = " + U); 881 System.out.println("Cp = " + Cp); 882 System.out.println("Uf = " + Uf); 883 //System.out.println("Cp.ring = " + Cp.ring.toScript() + ", Uf.ring = " + Uf.ring.toScript() + "\n"); 884 System.out.println(""); 885 //throw new NoLiftingException("no factorization, Cp != Uf"); 886 } 887 } 888 if (E.isZERO()) { 889 logger.info("liftHensel leaving with zero E, Ep"); 890 } 891 if (false && debug) { 892 // remove normalization required ?? 893 GreatestCommonDivisorAbstract<BigInteger> ufd = GCDFactory.getImplementation(new BigInteger()); 894 List<GenPolynomial<BigInteger>> Fii = new ArrayList<GenPolynomial<BigInteger>>(U.size()); 895 for (GenPolynomial<BigInteger> bi : Si) { 896 GenPolynomial<BigInteger> ci = ufd.content(bi); //ufd.primitivePart(bi); // ?? 897 if (!ci.isONE()) { 898 System.out.println("bi = " + bi + ", cont(bi) = " + ci); 899 } 900 //Fii.add(ci); 901 } 902 //Si = Fii; 903 //System.out.println("Si = " + Si); 904 } 905 logger.info("multivariate lift: U = " + U + ", of " + F); 906 return U; 907 } 908 909 910 /** 911 * Modular Hensel full lifting algorithm. Let p = A_i.ring.coFac.modul() and 912 * assume ggt(a,b) == 1 mod p, for a, b in A. 913 * @param C GenPolynomial with integer coefficients 914 * @param F list of modular GenPolynomials, mod (I_v, p ) 915 * @param V list of integer substitution values 916 * @param k desired approximation exponent p^k. 917 * @param G = [g_1,...,g_n] list of factors of leading coefficients. 918 * @return [c_1,..., c_n] with prod_i c_i = C mod p^k. 919 */ 920 @SuppressWarnings("unchecked") 921 public static <MOD extends GcdRingElem<MOD> & Modular> List<GenPolynomial<MOD>> liftHenselFull( 922 GenPolynomial<BigInteger> C, List<GenPolynomial<MOD>> F, List<BigInteger> V, long k, 923 List<GenPolynomial<BigInteger>> G) throws NoLiftingException { 924 if (F == null || F.size() == 0) { 925 return new ArrayList<GenPolynomial<MOD>>(); 926 } 927 GenPolynomialRing<MOD> pkfac = F.get(0).ring; 928 //long d = C.degree(); 929 // setup q = p^k 930 RingFactory<MOD> cfac = pkfac.coFac; 931 ModularRingFactory<MOD> pcfac = (ModularRingFactory<MOD>) cfac; 932 //System.out.println("pcfac = " + pcfac); 933 BigInteger p = pcfac.getIntegerModul(); 934 BigInteger q = p.power(k); 935 ModularRingFactory<MOD> mcfac; 936 if (ModLongRing.MAX_LONG.compareTo(q.getVal()) > 0) { 937 mcfac = (ModularRingFactory) new ModLongRing(q.getVal()); 938 } else { 939 mcfac = (ModularRingFactory) new ModIntegerRing(q.getVal()); 940 } 941 //System.out.println("mcfac = " + mcfac); 942 943 // convert C from Z[...] to Z_q[...] 944 GenPolynomialRing<MOD> qcfac = new GenPolynomialRing<MOD>(mcfac, C.ring); 945 GenPolynomial<MOD> Cq = PolyUtil.<MOD> fromIntegerCoefficients(qcfac, C); 946 //System.out.println("C = " + C); 947 //System.out.println("Cq = " + Cq); 948 949 // convert g_i from Z[...] to Z_q[...] 950 GenPolynomialRing<MOD> gcfac = new GenPolynomialRing<MOD>(mcfac, G.get(0).ring); 951 List<GenPolynomial<MOD>> GQ = new ArrayList<GenPolynomial<MOD>>(); 952 boolean allOnes = true; 953 for (GenPolynomial<BigInteger> g : G) { 954 if (!g.isONE()) { 955 allOnes = false; 956 } 957 GenPolynomial<MOD> gq = PolyUtil.<MOD> fromIntegerCoefficients(gcfac, g); 958 GQ.add(gq); 959 } 960 //System.out.println("G = " + G); 961 //System.out.println("GQ = " + GQ); 962 963 // evaluate C to Z_q[x] 964 GenPolynomialRing<MOD> pf = qcfac; 965 GenPolynomial<MOD> ap = Cq; 966 for (int j = C.ring.nvar; j > 1; j--) { 967 pf = pf.contract(1); 968 //MOD vp = mcfac.fromInteger(V.get(C.ring.nvar - j).getSymmetricInteger().getVal()); 969 MOD vp = mcfac.fromInteger(V.get(C.ring.nvar - j).getVal()); 970 //System.out.println("vp = " + vp); 971 ap = PolyUtil.<MOD> evaluateMain(pf, ap, vp); 972 //System.out.println("ap = " + ap); 973 } 974 GenPolynomial<MOD> Cq1 = ap; 975 //System.out.println("Cq1 = " + Cq1); 976 if (Cq1.isZERO()) { 977 throw new NoLiftingException("C mod (I, p^k) == 0: " + C); 978 } 979 GenPolynomialRing<BigInteger> ifac = new GenPolynomialRing<BigInteger>(new BigInteger(), pf); 980 GenPolynomial<BigInteger> Ci = PolyUtil.integerFromModularCoefficients(ifac, Cq1); 981 //System.out.println("Ci = " + Ci); 982 GreatestCommonDivisorAbstract<BigInteger> ufd = GCDFactory.getImplementation(new BigInteger()); 983 Ci = Ci.abs(); 984 BigInteger cCi = ufd.baseContent(Ci); 985 Ci = Ci.divide(cCi); 986 //System.out.println("cCi = " + cCi); 987 //System.out.println("Ci = " + Ci); 988 ////System.out.println("F.fac = " + F.get(0).ring); 989 990 // evaluate G to Z_q 991 //List<GenPolynomial<MOD>> GP = new ArrayList<GenPolynomial<MOD>>(); 992 for (GenPolynomial<MOD> gq : GQ) { 993 GenPolynomialRing<MOD> gf = gcfac; 994 GenPolynomial<MOD> gp = gq; 995 for (int j = gcfac.nvar; j > 1; j--) { 996 gf = gf.contract(1); 997 //MOD vp = mcfac.fromInteger(V.get(gcfac.nvar - j).getSymmetricInteger().getVal()); 998 MOD vp = mcfac.fromInteger(V.get(gcfac.nvar - j).getVal()); 999 //System.out.println("vp = " + vp); 1000 gp = PolyUtil.<MOD> evaluateMain(gf, gp, vp); 1001 //System.out.println("gp = " + gp); 1002 } 1003 //GP.add(gp); 1004 } 1005 //System.out.println("GP = " + GP); // + ", GP.ring = " + GP.get(0).ring); 1006 1007 // leading coefficient for recursion base, for Cq1 and list GP 1008 BigInteger gi0 = Ci.leadingBaseCoefficient(); // gq0.getSymmetricInteger(); 1009 //System.out.println("gi0 = " + gi0); 1010 1011 // lift F to Z_{p^k}[x] 1012 //System.out.println("Ci = " + Ci + ", F = " + F + ", k = " + k + ", p = " + F.get(0).ring + ", gi0 = " + gi0); 1013 List<GenPolynomial<MOD>> U1 = null; 1014 if (gi0.isONE()) { 1015 U1 = HenselUtil.<MOD> liftHenselMonic(Ci, F, k); 1016 } else { 1017 U1 = HenselUtil.<MOD> liftHensel(Ci, F, k, gi0); // gi0 TODO ?? 1018 } 1019 logger.info("univariate lift: Ci = " + Ci + ", F = " + F + ", U1 = " + U1); 1020 //System.out.println("U1.fac = " + U1.get(0).ring); 1021 1022 // adjust leading coefficients of U1 with F 1023 List<GenPolynomial<BigInteger>> U1i = PolyUtil.<MOD> integerFromModularCoefficients(Ci.ring, U1); 1024 //System.out.println("U1i = " + U1i); 1025 boolean t = HenselUtil.isHenselLift(Ci, q, p, U1i); 1026 //System.out.println("isLift(U1) = " + t); 1027 if (!t) { 1028 //System.out.println("NoLiftingException, Ci = " + Ci + ", U1i = " + U1i); 1029 throw new NoLiftingException("Ci = " + Ci + ", U1i = " + U1i); 1030 } 1031 MOD cC = mcfac.fromInteger(cCi.getVal()); 1032 List<GenPolynomial<MOD>> U1f = PolyUtil.<MOD> fromIntegerCoefficients(F.get(0).ring, U1i); 1033 //System.out.println("U1f = " + U1f); 1034 List<GenPolynomial<MOD>> U1s = new ArrayList<GenPolynomial<MOD>>(U1.size()); 1035 int j = 0; 1036 int s = 0; 1037 for (GenPolynomial<MOD> u : U1) { 1038 GenPolynomial<MOD> uf = U1f.get(j); 1039 GenPolynomial<MOD> f = F.get(j); 1040 GenPolynomial<BigInteger> ui = U1i.get(j); 1041 GenPolynomial<BigInteger> gi = G.get(j); 1042 if (ui.signum() != gi.signum()) { 1043 //System.out.println("ui = " + ui + ", gi = " + gi); 1044 u = u.negate(); 1045 uf = uf.negate(); 1046 s++; 1047 } 1048 j++; 1049 if (uf.isConstant()) { 1050 //System.out.println("u = " + u); 1051 u = u.monic(); 1052 //System.out.println("u = " + u); 1053 u = u.multiply(cC); 1054 cC = cC.divide(cC); 1055 //System.out.println("u = " + u); 1056 } else { 1057 MOD x = f.leadingBaseCoefficient().divide(uf.leadingBaseCoefficient()); 1058 //System.out.println("x = " + x + ", xi = " + x.getSymmetricInteger()); 1059 if (!x.isONE()) { 1060 MOD xq = mcfac.fromInteger(x.getSymmetricInteger().getVal()); 1061 //System.out.println("xq = " + xq); 1062 u = u.multiply(xq); 1063 cC = cC.divide(xq); 1064 //System.out.println("cC = " + cC); 1065 } 1066 } 1067 U1s.add(u); 1068 } 1069 //if ( s % 2 != 0 || !cC.isONE()) { 1070 if (!cC.isONE()) { 1071 throw new NoLiftingException("s = " + s + ", Ci = " + Ci + ", U1i = " + U1i + ", cC = " + cC); 1072 } 1073 U1 = U1s; 1074 U1i = PolyUtil.<MOD> integerFromModularCoefficients(Ci.ring, U1); 1075 //System.out.println("U1i = " + U1i); 1076 U1f = PolyUtil.<MOD> fromIntegerCoefficients(F.get(0).ring, U1i); 1077 if (!F.equals(U1f)) { // evtl loop until reached 1078 System.out.println("F = " + F); 1079 System.out.println("U1f = " + U1f); 1080 throw new NoLiftingException("F = " + F + ", U1f = " + U1f); 1081 } 1082 logger.info("multivariate lift: U1 = " + U1); 1083 1084 // lift U to Z_{p^k}[x,...] 1085 //System.out.println("C = " + C + ", U1 = " + U1 + ", V = " + V + ", k = " + k + ", q = " + U1.get(0).ring + ", G = " + G); 1086 List<GenPolynomial<MOD>> U = null; 1087 if (allOnes) { 1088 U = HenselMultUtil.<MOD> liftHenselMonic(C, Cq, U1, V, k); 1089 } else { 1090 U = HenselMultUtil.<MOD> liftHensel(C, Cq, U1, V, k, G); 1091 } 1092 logger.info("multivariate lift: C = " + C + ", U1 = " + U1 + ", U = " + U); 1093 //System.out.println("U = " + U); 1094 //System.out.println("U.fac = " + U.get(0).ring); 1095 return U; 1096 } 1097 1098}