001/* 002 * $Id: HenselMultUtil.java 5313 2015-10-26 22:46:38Z kredel $ 003 */ 004 005package edu.jas.ufd; 006 007 008import java.util.ArrayList; 009import java.util.List; 010 011import org.apache.log4j.Logger; 012 013import edu.jas.arith.BigInteger; 014import edu.jas.arith.ModIntegerRing; 015import edu.jas.arith.ModLongRing; 016import edu.jas.arith.Modular; 017import edu.jas.arith.ModularRingFactory; 018import edu.jas.poly.GenPolynomial; 019import edu.jas.poly.GenPolynomialRing; 020import edu.jas.poly.PolyUtil; 021import edu.jas.ps.PolynomialTaylorFunction; 022import edu.jas.ps.TaylorFunction; 023import edu.jas.ps.UnivPowerSeries; 024import edu.jas.ps.UnivPowerSeriesRing; 025import edu.jas.structure.GcdRingElem; 026import edu.jas.structure.Power; 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 = Logger.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 //GenPolynomialRing<BigInteger> cifac = new GenPolynomialRing<BigInteger>(new BigInteger(), ckfac); 100 //System.out.println("ifac = " + ifac.toScript()); 101 String[] mn = new String[] { pkfac.getVars()[pkfac.nvar - 1] }; 102 GenPolynomialRing<GenPolynomial<MOD>> qrfac = new GenPolynomialRing<GenPolynomial<MOD>>(ckfac, 1, mn); 103 //System.out.println("qrfac = " + qrfac); 104 105 List<GenPolynomial<MOD>> sup = new ArrayList<GenPolynomial<MOD>>(su.size()); 106 List<GenPolynomial<BigInteger>> supi = new ArrayList<GenPolynomial<BigInteger>>(su.size()); 107 for (GenPolynomial<MOD> s : su) { 108 GenPolynomial<MOD> sp = s.extend(pkfac, 0, 0L); 109 sup.add(sp); 110 GenPolynomial<BigInteger> spi = PolyUtil.integerFromModularCoefficients(ifac, sp); 111 supi.add(spi); 112 } 113 //System.out.println("sup = " + sup); 114 //System.out.println("supi = " + supi); 115 GenPolynomial<BigInteger> Ai = PolyUtil.integerFromModularCoefficients(ifac, A); 116 GenPolynomial<BigInteger> Bi = PolyUtil.integerFromModularCoefficients(ifac, B); 117 GenPolynomial<BigInteger> Ci = PolyUtil.integerFromModularCoefficients(ifac, C); 118 //System.out.println("Ai = " + Ai); 119 //System.out.println("Bi = " + Bi); 120 //System.out.println("Ci = " + Ci); 121 //GenPolynomial<MOD> aq = PolyUtil.<MOD> fromIntegerCoefficients(pkfac, Ai); 122 //GenPolynomial<MOD> bq = PolyUtil.<MOD> fromIntegerCoefficients(pkfac, Bi); 123 //System.out.println("aq = " + aq); 124 //System.out.println("bq = " + bq); 125 126 // compute error: 127 GenPolynomial<BigInteger> E = Ci; // - sum_i s_i b_i 128 E = E.subtract(Bi.multiply(supi.get(0))); 129 E = E.subtract(Ai.multiply(supi.get(1))); 130 //System.out.println("E = " + E); 131 if (E.isZERO()) { 132 logger.info("liftDiophant leaving on zero E"); 133 return sup; 134 } 135 GenPolynomial<MOD> Ep = PolyUtil.<MOD> fromIntegerCoefficients(pkfac, E); 136 //System.out.println("Ep(0," + pkfac.nvar + ") = " + Ep); 137 logger.info("Ep(0," + pkfac.nvar + ") = " + Ep); 138 if (Ep.isZERO()) { 139 logger.info("liftDiophant leaving on zero Ep mod p^k"); 140 return sup; 141 } 142 for (int e = 1; e <= d; e++) { 143 //System.out.println("\ne = " + e + " -------------------------------------- " + pkfac.nvar); 144 GenPolynomial<GenPolynomial<MOD>> Epr = PolyUtil.<MOD> recursive(qrfac, Ep); 145 //System.out.println("Epr = " + Epr); 146 UnivPowerSeriesRing<GenPolynomial<MOD>> psfac = new UnivPowerSeriesRing<GenPolynomial<MOD>>(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("ckfac != pkfac: " + ckfac.coFac + " != " 164 + 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 //GenPolynomialRing<BigInteger> cifac = new GenPolynomialRing<BigInteger>(new BigInteger(), ckfac); 278 //System.out.println("ifac = " + ifac.toScript()); 279 String[] mn = new String[] { pkfac.getVars()[pkfac.nvar - 1] }; 280 GenPolynomialRing<GenPolynomial<MOD>> qrfac = new GenPolynomialRing<GenPolynomial<MOD>>(ckfac, 1, mn); 281 //System.out.println("qrfac = " + qrfac); 282 283 List<GenPolynomial<MOD>> sup = new ArrayList<GenPolynomial<MOD>>(su.size()); 284 List<GenPolynomial<BigInteger>> supi = new ArrayList<GenPolynomial<BigInteger>>(su.size()); 285 for (GenPolynomial<MOD> s : su) { 286 GenPolynomial<MOD> sp = s.extend(pkfac, 0, 0L); 287 sup.add(sp); 288 GenPolynomial<BigInteger> spi = PolyUtil.integerFromModularCoefficients(ifac, sp); 289 supi.add(spi); 290 } 291 //System.out.println("sup = " + sup); 292 //System.out.println("supi = " + supi); 293 //List<GenPolynomial<BigInteger>> Ai = new ArrayList<GenPolynomial<BigInteger>>(A.size()); 294 //for (GenPolynomial<MOD> a : A) { 295 // GenPolynomial<BigInteger> ai = PolyUtil.integerFromModularCoefficients(ifac, a); 296 // Ai.add(ai); 297 //} 298 List<GenPolynomial<BigInteger>> Bi = new ArrayList<GenPolynomial<BigInteger>>(A.size()); 299 for (GenPolynomial<MOD> b : Bp) { 300 GenPolynomial<BigInteger> bi = PolyUtil.integerFromModularCoefficients(ifac, b); 301 Bi.add(bi); 302 } 303 GenPolynomial<BigInteger> Ci = PolyUtil.integerFromModularCoefficients(ifac, C); 304 //System.out.println("Ai = " + Ai); 305 //System.out.println("Ci = " + Ci); 306 307 //List<GenPolynomial<MOD>> Aq = new ArrayList<GenPolynomial<MOD>>(A.size()); 308 //for (GenPolynomial<BigInteger> ai : Ai) { 309 // GenPolynomial<MOD> aq = PolyUtil.<MOD> fromIntegerCoefficients(pkfac, ai); 310 // Aq.add(aq); 311 //} 312 //System.out.println("Aq = " + Aq); 313 314 // compute error: 315 GenPolynomial<BigInteger> E = Ci; // - sum_i s_i b_i 316 int i = 0; 317 for (GenPolynomial<BigInteger> bi : Bi) { 318 E = E.subtract(bi.multiply(supi.get(i++))); 319 } 320 //System.out.println("E = " + E); 321 if (E.isZERO()) { 322 logger.info("liftDiophant leaving on zero E"); 323 return sup; 324 } 325 GenPolynomial<MOD> Ep = PolyUtil.<MOD> fromIntegerCoefficients(pkfac, E); 326 //System.out.println("Ep(0," + pkfac.nvar + ") = " + Ep); 327 logger.info("Ep(0," + pkfac.nvar + ") = " + Ep); 328 if (Ep.isZERO()) { 329 logger.info("liftDiophant leaving on zero Ep mod p^k"); 330 return sup; 331 } 332 for (int e = 1; e <= d; e++) { 333 //System.out.println("\ne = " + e + " -------------------------------------- " + pkfac.nvar); 334 GenPolynomial<GenPolynomial<MOD>> Epr = PolyUtil.<MOD> recursive(qrfac, Ep); 335 //System.out.println("Epr = " + Epr); 336 UnivPowerSeriesRing<GenPolynomial<MOD>> psfac = new UnivPowerSeriesRing<GenPolynomial<MOD>>(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("ckfac != pkfac: " + ckfac.coFac + " != " 356 + 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 k approximation exponent. 413 * @param L = [g_0,...,g_{n-1}] list of lifted modular polynomials. 414 * @return true if C = prod_{0,...,n-1} g_i mod p^k, else false. 415 * @deprecated use isHenselLift() without parameter k 416 */ 417 @Deprecated 418 public static <MOD extends GcdRingElem<MOD> & Modular> boolean isHenselLift(GenPolynomial<BigInteger> C, 419 GenPolynomial<MOD> Cp, List<GenPolynomial<MOD>> F, long k, List<GenPolynomial<MOD>> L) { 420 return isHenselLift(C,Cp,F,L); 421 } 422 423 424 /** 425 * Modular Hensel lifting algorithm on coefficients test. Let p = 426 * f_i.ring.coFac.modul() and assume C == prod_{0,...,n-1} f_i mod p with 427 * gcd(f_i,f_j) == 1 mod p for i != j 428 * @param C integer polynomial 429 * @param Cp GenPolynomial mod p^k 430 * @param F = [f_0,...,f_{n-1}] list of monic modular polynomials. 431 * @param L = [g_0,...,g_{n-1}] list of lifted modular polynomials. 432 * @return true if C = prod_{0,...,n-1} g_i mod p^k, else false. 433 */ 434 public static <MOD extends GcdRingElem<MOD> & Modular> boolean isHenselLift(GenPolynomial<BigInteger> C, 435 GenPolynomial<MOD> Cp, List<GenPolynomial<MOD>> F, List<GenPolynomial<MOD>> L) { 436 boolean t = true; 437 GenPolynomialRing<MOD> qfac = L.get(0).ring; 438 GenPolynomial<MOD> q = qfac.getONE(); 439 for (GenPolynomial<MOD> fi : L) { 440 q = q.multiply(fi); 441 } 442 t = Cp.equals(q); 443 if (!t) { 444 System.out.println("Cp = " + Cp); 445 System.out.println("q = " + q); 446 System.out.println("Cp != q: " + Cp.subtract(q)); 447 return t; 448 } 449 GenPolynomialRing<BigInteger> dfac = C.ring; 450 GenPolynomial<BigInteger> Ci = PolyUtil.integerFromModularCoefficients(dfac, q); 451 t = C.equals(Ci); 452 if (!t) { 453 System.out.println("C = " + C); 454 System.out.println("Ci = " + Ci); 455 System.out.println("C != Ci: " + C.subtract(Ci)); 456 return t; 457 } 458 // test L mod id(V) == F 459 return t; 460 } 461 462 463 /** 464 * Modular Hensel lifting algorithm, monic case. Let p = 465 * A_i.ring.coFac.modul() and assume ggt(a,b) == 1 mod p, for a, b in A. 466 * @param C monic GenPolynomial with integer coefficients 467 * @param Cp GenPolynomial mod p^k 468 * @param F list of modular GenPolynomials, mod (I_v, p^k ) 469 * @param V list of integer substitution values 470 * @param k desired approximation exponent p^k. 471 * @return [g'_1,..., g'_n] with prod_i g'_i = Cp mod p^k. 472 */ 473 public static <MOD extends GcdRingElem<MOD> & Modular> List<GenPolynomial<MOD>> liftHenselMonic( 474 GenPolynomial<BigInteger> C, GenPolynomial<MOD> Cp, List<GenPolynomial<MOD>> F, 475 List<BigInteger> V, long k) throws NoLiftingException { 476 GenPolynomialRing<MOD> pkfac = Cp.ring; 477 //if (pkfac.nvar == 1) { // V ignored 478 // return HenselUtil.<MOD> liftHenselMonic(C,F,k); 479 //} 480 long d = C.degree(); 481 //System.out.println("d = " + d); 482 // prepare stack of polynomial rings and polynomials 483 List<GenPolynomialRing<MOD>> Pfac = new ArrayList<GenPolynomialRing<MOD>>(); 484 List<GenPolynomial<MOD>> Ap = new ArrayList<GenPolynomial<MOD>>(); 485 List<MOD> Vb = new ArrayList<MOD>(); 486 MOD v = pkfac.coFac.fromInteger( V.get(0).getVal() ); 487 Pfac.add(pkfac); 488 Ap.add(Cp); 489 Vb.add(v); 490 GenPolynomialRing<MOD> pf = pkfac; 491 GenPolynomial<MOD> ap = Cp; 492 for (int j = pkfac.nvar; j > 2; j--) { 493 pf = pf.contract(1); 494 Pfac.add(0, pf); 495 //MOD vp = pkfac.coFac.fromInteger(V.get(j - 2).getSymmetricInteger().getVal()); 496 MOD vp = pkfac.coFac.fromInteger(V.get(j - 2).getVal()); 497 //System.out.println("vp = " + vp); 498 Vb.add(1, vp); 499 ap = PolyUtil.<MOD> evaluateMain(pf, ap, vp); 500 Ap.add(0, ap); 501 } 502 //System.out.println("Pfac = " + Pfac); 503 if (debug) { 504 logger.debug("Pfac = " + Pfac); 505 } 506 //System.out.println("Ap = " + Ap); 507 //System.out.println("V = " + V); 508 //System.out.println("Vb = " + Vb); 509 // setup bi-variate base case 510 GenPolynomialRing<MOD> pk1fac = F.get(0).ring; 511 if (!pkfac.coFac.equals(pk1fac.coFac)) { 512 throw new IllegalArgumentException("F.ring != pkfac: " + pk1fac + " != " + pkfac); 513 } 514 // TODO: adjust leading coefficients 515 pkfac = Pfac.get(0); 516 //Cp = Ap.get(0); 517 //System.out.println("pkfac = " + pkfac.toScript()); 518 //System.out.println("pk1fac = " + pk1fac.toScript()); 519 GenPolynomialRing<BigInteger> i1fac = new GenPolynomialRing<BigInteger>(new BigInteger(), pk1fac); 520 //System.out.println("i1fac = " + i1fac.toScript()); 521 List<GenPolynomial<BigInteger>> Bi = new ArrayList<GenPolynomial<BigInteger>>(F.size()); 522 for (GenPolynomial<MOD> b : F) { 523 GenPolynomial<BigInteger> bi = PolyUtil.integerFromModularCoefficients(i1fac, b); 524 Bi.add(bi); 525 } 526 //System.out.println("Bi = " + Bi); 527 // evaluate Cp at v_n: 528 //ModularRingFactory<MOD> cf = (ModularRingFactory<MOD>) pkfac.coFac; 529 //MOD vp = cf.fromInteger(v.getSymmetricInteger().getVal()); 530 //System.out.println("v = " + v + ", vp = " + vp); 531 GenPolynomialRing<MOD> ckfac; // = pkfac.contract(1); 532 //GenPolynomial<MOD> Cs = PolyUtil.<MOD> evaluateMain(ckfac, Cp, vp); 533 //System.out.println("Cp = " + Cp); 534 //System.out.println("Cs = " + Cs); 535 536 List<GenPolynomial<MOD>> U = new ArrayList<GenPolynomial<MOD>>(F.size()); 537 for (GenPolynomial<MOD> b : F) { 538 GenPolynomial<MOD> bi = b.extend(pkfac, 0, 0L); 539 U.add(bi); 540 } 541 //System.out.println("U = " + U); 542 List<GenPolynomial<MOD>> U1 = F; 543 //System.out.println("U1 = " + U1); 544 545 GenPolynomial<BigInteger> E = C.ring.getZERO(); 546 List<MOD> Vh = new ArrayList<MOD>(); 547 548 while (Pfac.size() > 0) { // loop through stack of polynomial rings 549 pkfac = Pfac.remove(0); 550 Cp = Ap.remove(0); 551 v = Vb.remove(0); 552 //Vh.add(0,v); 553 //System.out.println("\npkfac = " + pkfac.toScript() + " ================================== " + Vh); 554 555 // (x_n - v) 556 GenPolynomial<MOD> mon = pkfac.getONE(); 557 GenPolynomial<MOD> xv = pkfac.univariate(0, 1); 558 xv = xv.subtract(pkfac.fromInteger(v.getSymmetricInteger().getVal())); 559 //System.out.println("xv = " + xv); 560 561 long deg = Cp.degree(pkfac.nvar - 1); 562 //System.out.println("deg = " + deg); 563 564 GenPolynomialRing<BigInteger> ifac = new GenPolynomialRing<BigInteger>(new BigInteger(), pkfac); 565 //System.out.println("ifac = " + ifac.toScript()); 566 List<GenPolynomial<BigInteger>> Bip = new ArrayList<GenPolynomial<BigInteger>>(F.size()); 567 for (GenPolynomial<BigInteger> b : Bi) { 568 GenPolynomial<BigInteger> bi = b.extend(ifac, 0, 0L); 569 Bip.add(bi); 570 } 571 Bi = Bip; 572 //System.out.println("Bi = " + Bi); 573 GenPolynomial<BigInteger> Ci = PolyUtil.integerFromModularCoefficients(ifac, Cp); 574 //System.out.println("Ci = " + Ci); 575 576 // compute error: 577 E = ifac.getONE(); 578 for (GenPolynomial<BigInteger> bi : Bi) { 579 E = E.multiply(bi); 580 } 581 E = Ci.subtract(E); 582 //System.out.println("E = " + E); 583 GenPolynomial<MOD> Ep = PolyUtil.<MOD> fromIntegerCoefficients(pkfac, E); 584 //System.out.println("Ep(0," + pkfac.nvar + ") = " + Ep); 585 logger.info("Ep(0," + deg + "," + pkfac.nvar + ") = " + Ep); 586 587 String[] mn = new String[] { pkfac.getVars()[pkfac.nvar - 1] }; 588 ckfac = pkfac.contract(1); 589 GenPolynomialRing<GenPolynomial<MOD>> pkrfac = new GenPolynomialRing<GenPolynomial<MOD>>(ckfac, 590 1, mn); 591 //System.out.println("pkrfac = " + pkrfac.toScript()); 592 593 for (int e = 1; e <= deg && !Ep.isZERO(); e++) { 594 //System.out.println("\ne = " + e + " -------------------------------------- " + pkfac.nvar); 595 GenPolynomial<GenPolynomial<MOD>> Epr = PolyUtil.<MOD> recursive(pkrfac, Ep); 596 //System.out.println("Epr = " + Epr); 597 UnivPowerSeriesRing<GenPolynomial<MOD>> psfac = new UnivPowerSeriesRing<GenPolynomial<MOD>>( 598 pkrfac); 599 //System.out.println("psfac = " + psfac); 600 TaylorFunction<GenPolynomial<MOD>> T = new PolynomialTaylorFunction<GenPolynomial<MOD>>(Epr); 601 //System.out.println("T = " + T); 602 //List<GenPolynomial<MOD>> Vs = new ArrayList<GenPolynomial<MOD>>(1); 603 GenPolynomial<MOD> vq = ckfac.fromInteger(v.getSymmetricInteger().getVal()); 604 //Vs.add(vq); 605 //System.out.println("Vs = " + Vs + ", Vh = " + Vh); 606 UnivPowerSeries<GenPolynomial<MOD>> Epst = psfac.seriesOfTaylor(T, vq); 607 //System.out.println("Epst = " + Epst); 608 logger.info("Epst(" + e + "," + deg + ", " + pkfac.nvar + ") = " + Epst); 609 GenPolynomial<MOD> cm = Epst.coefficient(e); 610 //System.out.println("cm = " + cm); 611 if (cm.isZERO()) { 612 continue; 613 } 614 List<GenPolynomial<MOD>> Ud = HenselMultUtil.<MOD> liftDiophant(U1, cm, Vh, d, k); 615 //System.out.println("Ud = " + Ud); 616 617 mon = mon.multiply(xv); // Power.<GenPolynomial<MOD>> power(pkfac,xv,e); 618 //System.out.println("mon = " + mon); 619 //List<GenPolynomial<MOD>> Sd = new ArrayList<GenPolynomial<MOD>>(Ud.size()); 620 int i = 0; 621 List<GenPolynomial<BigInteger>> Si = new ArrayList<GenPolynomial<BigInteger>>(Ud.size()); 622 for (GenPolynomial<MOD> dd : Ud) { 623 //System.out.println("dd = " + dd); 624 GenPolynomial<MOD> de = dd.extend(pkfac, 0, 0L); 625 GenPolynomial<MOD> dm = de.multiply(mon); 626 //Sd.add(dm); 627 de = U.get(i).sum(dm); 628 //System.out.println("de = " + de); 629 U.set(i++, de); 630 GenPolynomial<BigInteger> si = PolyUtil.integerFromModularCoefficients(ifac, de); 631 Si.add(si); 632 } 633 //System.out.println("Sd = " + Sd); 634 //System.out.println("U = " + U); 635 //System.out.println("Si = " + Si); 636 637 // compute new error: 638 E = ifac.getONE(); 639 for (GenPolynomial<BigInteger> bi : Si) { 640 E = E.multiply(bi); 641 } 642 E = Ci.subtract(E); 643 //System.out.println("E = " + E); 644 Ep = PolyUtil.<MOD> fromIntegerCoefficients(pkfac, E); 645 //System.out.println("Ep(0," + pkfac.nvar + ") = " + Ep); 646 logger.info("Ep(" + e + "," + deg + "," + pkfac.nvar + ") = " + Ep); 647 } 648 Vh.add(v); 649 U1 = U; 650 if (Pfac.size() > 0) { 651 List<GenPolynomial<MOD>> U2 = new ArrayList<GenPolynomial<MOD>>(U.size()); 652 pkfac = Pfac.get(0); 653 for (GenPolynomial<MOD> b : U) { 654 GenPolynomial<MOD> bi = b.extend(pkfac, 0, 0L); 655 U2.add(bi); 656 } 657 U = U2; 658 //System.out.println("U = " + U); 659 } 660 } 661 if (E.isZERO()) { 662 logger.info("liftHensel leaving with zero E"); 663 } 664 return U; 665 } 666 667 668 /** 669 * Modular Hensel lifting algorithm. Let p = A_i.ring.coFac.modul() and 670 * assume ggt(a,b) == 1 mod p, for a, b in A. 671 * @param C GenPolynomial with integer coefficients 672 * @param Cp GenPolynomial C mod p^k 673 * @param F list of modular GenPolynomials, mod (I_v, p^k ) 674 * @param V list of integral substitution values 675 * @param k desired approximation exponent p^k. 676 * @param G list of leading coefficients of the factors of C. 677 * @return [g'_1,..., g'_n] with prod_i g'_i = Cp mod p^k. 678 */ 679 public static <MOD extends GcdRingElem<MOD> & Modular> List<GenPolynomial<MOD>> liftHensel( 680 GenPolynomial<BigInteger> C, GenPolynomial<MOD> Cp, List<GenPolynomial<MOD>> F, 681 List<BigInteger> V, long k, List<GenPolynomial<BigInteger>> G) 682 throws NoLiftingException { 683 GenPolynomialRing<MOD> pkfac = Cp.ring; 684 long d = C.degree(); 685 //System.out.println("C = " + C); 686 //System.out.println("Cp = " + Cp); 687 //System.out.println("G = " + G); 688 689 //GenPolynomial<BigInteger> cd = G.get(0); // 1 690 //System.out.println("cd = " + cd + ", ring = " + C.ring); 691 //if ( cd.equals(C.ring.univariate(0)) ) { 692 // System.out.println("cd == G[1]"); 693 //} 694 // G mod p^k, in all variables 695 GenPolynomialRing<MOD> pkfac1 = new GenPolynomialRing<MOD>(pkfac.coFac, G.get(0).ring); 696 List<GenPolynomial<MOD>> Lp = new ArrayList<GenPolynomial<MOD>>(G.size()); 697 for (GenPolynomial<BigInteger> cd1 : G) { 698 GenPolynomial<MOD> cdq = PolyUtil.<MOD> fromIntegerCoefficients(pkfac1, cd1); 699 cdq = cdq.extendLower(pkfac, 0, 0L); // reintroduce lower variable 700 Lp.add(cdq); 701 } 702 logger.info("G modulo p^k: " + Lp); // + ", ring = " + pkfac1); 703 704 // prepare stack of polynomial rings, polynomials and evaluated leading coefficients 705 List<GenPolynomialRing<MOD>> Pfac = new ArrayList<GenPolynomialRing<MOD>>(); 706 List<GenPolynomial<MOD>> Ap = new ArrayList<GenPolynomial<MOD>>(); 707 List<List<GenPolynomial<MOD>>> Gp = new ArrayList<List<GenPolynomial<MOD>>>(); 708 List<MOD> Vb = new ArrayList<MOD>(); 709 //MOD v = V.get(0); // fromInteger 710 Pfac.add(pkfac); 711 Ap.add(Cp); 712 Gp.add(Lp); 713 GenPolynomialRing<MOD> pf = pkfac; 714 //GenPolynomialRing<MOD> pf1 = pkfac1; 715 GenPolynomial<MOD> ap = Cp; 716 List<GenPolynomial<MOD>> Lpp = Lp; 717 for (int j = pkfac.nvar; j > 2; j--) { 718 pf = pf.contract(1); 719 Pfac.add(0, pf); 720 //MOD vp = pkfac.coFac.fromInteger(V.get(pkfac.nvar - j).getSymmetricInteger().getVal()); 721 MOD vp = pkfac.coFac.fromInteger(V.get(pkfac.nvar - j).getVal()); 722 //System.out.println("vp = " + vp); 723 Vb.add(vp); 724 ap = PolyUtil.<MOD> evaluateMain(pf, ap, vp); 725 Ap.add(0, ap); 726 List<GenPolynomial<MOD>> Lps = new ArrayList<GenPolynomial<MOD>>(Lpp.size()); 727 for (GenPolynomial<MOD> qp : Lpp) { 728 GenPolynomial<MOD> qpe = PolyUtil.<MOD> evaluateMain(pf, qp, vp); 729 Lps.add(qpe); 730 } 731 //System.out.println("Lps = " + Lps); 732 Lpp = Lps; 733 Gp.add(0, Lpp); 734 } 735 Vb.add( pkfac.coFac.fromInteger( V.get(pkfac.nvar - 2).getVal() ) ); 736 //System.out.println("Pfac = " + Pfac); 737 if (debug) { 738 logger.debug("Pfac = " + Pfac); 739 } 740 //System.out.println("Ap = " + Ap); 741 //System.out.println("Gp = " + Gp); 742 //System.out.println("Gp[0] = " + Gp.get(0) + ", Gp[0].ring = " + Gp.get(0).get(0).ring); 743 //System.out.println("V = " + V); 744 //System.out.println("Vb = " + Vb + ", V == Vb: " + V.equals(Vb)); 745 746 // check bi-variate base case 747 GenPolynomialRing<MOD> pk1fac = F.get(0).ring; 748 if (!pkfac.coFac.equals(pk1fac.coFac)) { 749 throw new IllegalArgumentException("F.ring != pkfac: " + pk1fac + " != " + pkfac); 750 } 751 752 // init recursion 753 List<GenPolynomial<MOD>> U = F; 754 //logger.info("to lift U = " + U); // + ", U1.ring = " + U1.get(0).ring); 755 GenPolynomial<BigInteger> E = C.ring.getZERO(); 756 List<MOD> Vh = new ArrayList<MOD>(); 757 List<GenPolynomial<BigInteger>> Si; // = new ArrayList<GenPolynomial<BigInteger>>(F.size()); 758 MOD v = null; 759 760 while (Pfac.size() > 0) { // loop through stack of polynomial rings 761 pkfac = Pfac.remove(0); 762 Cp = Ap.remove(0); 763 Lpp = Gp.remove(0); 764 v = Vb.remove(Vb.size() - 1); // last in stack 765 //System.out.println("\npkfac = " + pkfac.toScript() + " ================================== " + v); 766 logger.info("stack loop: pkfac = " + pkfac.toScript() + " v = " + v); 767 768 List<GenPolynomial<MOD>> U1 = U; 769 logger.info("to lift U1 = " + U1); // + ", U1.ring = " + U1.get(0).ring); 770 U = new ArrayList<GenPolynomial<MOD>>(U1.size()); 771 772 // update U, replace leading coefficient if required 773 int j = 0; 774 for (GenPolynomial<MOD> b : U1) { 775 //System.out.println("b = " + b + ", b.ring = " + b.ring); 776 GenPolynomial<MOD> bi = b.extend(pkfac, 0, 0L); 777 GenPolynomial<MOD> li = Lpp.get(j); 778 if (!li.isONE()) { 779 //System.out.println("li = " + li + ", li.ring = " + li.ring); 780 //System.out.println("bi = " + bi); 781 GenPolynomialRing<GenPolynomial<MOD>> pkrfac = pkfac.recursive(pkfac.nvar - 1); 782 //System.out.println("pkrfac = " + pkrfac); 783 GenPolynomial<GenPolynomial<MOD>> br = PolyUtil.<MOD> recursive(pkrfac, bi); 784 //System.out.println("br = " + br); 785 GenPolynomial<GenPolynomial<MOD>> bs = PolyUtil.<MOD> switchVariables(br); 786 //System.out.println("bs = " + bs + ", bs.ring = " + bs.ring); 787 788 GenPolynomial<GenPolynomial<MOD>> lr = PolyUtil.<MOD> recursive(pkrfac, li); 789 //System.out.println("lr = " + lr); 790 GenPolynomial<GenPolynomial<MOD>> ls = PolyUtil.<MOD> switchVariables(lr); 791 //System.out.println("ls = " + ls + ", ls.ring = " + ls.ring); 792 if (!ls.isConstant() && !ls.isZERO()) { 793 throw new RuntimeException("ls not constant " + ls + ", li = " + li); 794 } 795 bs.doPutToMap(bs.leadingExpVector(), ls.leadingBaseCoefficient()); 796 //System.out.println("bs = " + bs + ", bs.ring = " + bs.ring); 797 br = PolyUtil.<MOD> switchVariables(bs); 798 //System.out.println("br = " + br); 799 bi = PolyUtil.<MOD> distribute(pkfac, br); 800 //System.out.println("bi = " + bi); 801 } 802 U.add(bi); 803 j++; 804 } 805 logger.info("U with leading coefficient replaced = " + U); // + ", U.ring = " + U.get(0).ring); 806 807 // (x_n - v) 808 GenPolynomial<MOD> mon = pkfac.getONE(); 809 GenPolynomial<MOD> xv = pkfac.univariate(0, 1); 810 xv = xv.subtract(pkfac.fromInteger(v.getSymmetricInteger().getVal())); 811 //System.out.println("xv = " + xv); 812 813 long deg = Cp.degree(pkfac.nvar - 1); 814 //System.out.println("deg = " + deg + ", degv = " + Cp.degreeVector()); 815 816 // convert to integer polynomials 817 GenPolynomialRing<BigInteger> ifac = new GenPolynomialRing<BigInteger>(new BigInteger(), pkfac); 818 //System.out.println("ifac = " + ifac.toScript()); 819 List<GenPolynomial<BigInteger>> Bi = PolyUtil.integerFromModularCoefficients(ifac, U); 820 //System.out.println("Bi = " + Bi); 821 GenPolynomial<BigInteger> Ci = PolyUtil.integerFromModularCoefficients(ifac, Cp); 822 //System.out.println("Ci = " + Ci); 823 824 // compute error: 825 E = ifac.getONE(); 826 for (GenPolynomial<BigInteger> bi : Bi) { 827 E = E.multiply(bi); 828 } 829 //System.out.println("E = " + E); 830 E = Ci.subtract(E); 831 //System.out.println("E = " + E); 832 GenPolynomial<MOD> Ep = PolyUtil.<MOD> fromIntegerCoefficients(pkfac, E); 833 logger.info("Ep(0," + deg + "," + pkfac.nvar + ") = " + Ep); 834 835 GenPolynomialRing<GenPolynomial<MOD>> pkrfac = pkfac.recursive(1); 836 GenPolynomialRing<MOD> ckfac = (GenPolynomialRing<MOD>) pkrfac.coFac; 837 //System.out.println("pkrfac = " + pkrfac.toScript()); 838 839 for (int e = 1; e <= deg && !Ep.isZERO(); e++) { 840 //System.out.println("\ne = " + e + " -------------------------------------- " + deg); 841 logger.info("approximation loop: e = " + e + " of deg = " + deg); 842 GenPolynomial<GenPolynomial<MOD>> Epr = PolyUtil.<MOD> recursive(pkrfac, Ep); 843 //System.out.println("Epr = " + Epr); 844 UnivPowerSeriesRing<GenPolynomial<MOD>> psfac = new UnivPowerSeriesRing<GenPolynomial<MOD>>( 845 pkrfac); 846 //System.out.println("psfac = " + psfac); 847 TaylorFunction<GenPolynomial<MOD>> T = new PolynomialTaylorFunction<GenPolynomial<MOD>>(Epr); 848 //System.out.println("T = " + T); 849 GenPolynomial<MOD> vq = ckfac.fromInteger(v.getSymmetricInteger().getVal()); 850 //System.out.println("vq = " + vq + ", Vh = " + Vh); 851 UnivPowerSeries<GenPolynomial<MOD>> Epst = psfac.seriesOfTaylor(T, vq); 852 //System.out.println("Epst = " + Epst); 853 logger.info("Epst(" + e + "," + deg + "," + pkfac.nvar + ") = " + Epst); 854 GenPolynomial<MOD> cm = Epst.coefficient(e); 855 if (cm.isZERO()) { 856 //System.out.println("cm = " + cm); 857 continue; 858 } 859 List<GenPolynomial<MOD>> Ud = HenselMultUtil.<MOD> liftDiophant(U1, cm, Vh, d, k); 860 //System.out.println("Ud = " + Ud); 861 862 mon = mon.multiply(xv); // Power.<GenPolynomial<MOD>> power(pkfac,xv,e); 863 //System.out.println("mon = " + mon); 864 //List<GenPolynomial<MOD>> Sd = new ArrayList<GenPolynomial<MOD>>(Ud.size()); 865 int i = 0; 866 Si = new ArrayList<GenPolynomial<BigInteger>>(Ud.size()); 867 for (GenPolynomial<MOD> dd : Ud) { 868 //System.out.println("dd = " + dd); 869 GenPolynomial<MOD> de = dd.extend(pkfac, 0, 0L); 870 GenPolynomial<MOD> dm = de.multiply(mon); 871 //Sd.add(dm); 872 de = U.get(i).sum(dm); 873 //System.out.println("de = " + de); 874 U.set(i++, de); 875 GenPolynomial<BigInteger> si = PolyUtil.integerFromModularCoefficients(ifac, de); 876 Si.add(si); 877 } 878 //System.out.println("Sd = " + Sd); 879 //System.out.println("U = " + U + ", U.ring = " + U.get(0).ring); 880 //System.out.println("Si = " + Si); 881 882 // compute new error: 883 E = ifac.getONE(); 884 for (GenPolynomial<BigInteger> bi : Si) { 885 E = E.multiply(bi); 886 } 887 E = Ci.subtract(E); 888 //System.out.println("E = " + E); 889 Ep = PolyUtil.<MOD> fromIntegerCoefficients(pkfac, E); 890 //System.out.println("Ep(0," + pkfac.nvar + ") = " + Ep); 891 logger.info("Ep(" + e + "," + deg + "," + pkfac.nvar + ") = " + Ep); 892 } 893 Vh.add(v); 894 GenPolynomial<MOD> Uf = U.get(0).ring.getONE(); 895 for (GenPolynomial<MOD> Upp : U) { 896 Uf = Uf.multiply(Upp); 897 } 898 if (false && !Cp.leadingExpVector().equals(Uf.leadingExpVector())) { // not meanigfull test 899 System.out.println("\nU = " + U); 900 System.out.println("Cp = " + Cp); 901 System.out.println("Uf = " + Uf); 902 //System.out.println("Cp.ring = " + Cp.ring.toScript() + ", Uf.ring = " + Uf.ring.toScript() + "\n"); 903 System.out.println(""); 904 //throw new NoLiftingException("no factorization, Cp != Uf"); 905 } 906 } 907 if (E.isZERO()) { 908 logger.info("liftHensel leaving with zero E, Ep"); 909 } 910 if (false && debug) { 911 // remove normalization required ?? 912 GreatestCommonDivisorAbstract<BigInteger> ufd = GCDFactory.getImplementation(new BigInteger()); 913 List<GenPolynomial<BigInteger>> Fii = new ArrayList<GenPolynomial<BigInteger>>(U.size()); 914 for (GenPolynomial<BigInteger> bi : Si) { 915 GenPolynomial<BigInteger> ci = ufd.content(bi); //ufd.primitivePart(bi); // ?? 916 if (!ci.isONE()) { 917 System.out.println("bi = " + bi + ", cont(bi) = " + ci); 918 } 919 //Fii.add(ci); 920 } 921 //Si = Fii; 922 //System.out.println("Si = " + Si); 923 } 924 logger.info("multivariate lift: U = " + U + ", of " + F); 925 return U; 926 } 927 928 929 /** 930 * Modular Hensel full lifting algorithm. Let p = A_i.ring.coFac.modul() and 931 * assume ggt(a,b) == 1 mod p, for a, b in A. 932 * @param C GenPolynomial with integer coefficients 933 * @param F list of modular GenPolynomials, mod (I_v, p ) 934 * @param V list of integer substitution values 935 * @param k desired approximation exponent p^k. 936 * @param G = [g_1,...,g_n] list of factors of leading coefficients. 937 * @return [c_1,..., c_n] with prod_i c_i = C mod p^k. 938 */ 939 @SuppressWarnings("unchecked") 940 public static <MOD extends GcdRingElem<MOD> & Modular> List<GenPolynomial<MOD>> liftHenselFull( 941 GenPolynomial<BigInteger> C, List<GenPolynomial<MOD>> F, List<BigInteger> V, long k, 942 List<GenPolynomial<BigInteger>> G) throws NoLiftingException { 943 if (F == null || F.size() == 0) { 944 return new ArrayList<GenPolynomial<MOD>>(); 945 } 946 GenPolynomialRing<MOD> pkfac = F.get(0).ring; 947 //long d = C.degree(); 948 // setup q = p^k 949 RingFactory<MOD> cfac = pkfac.coFac; 950 ModularRingFactory<MOD> pcfac = (ModularRingFactory<MOD>) cfac; 951 //System.out.println("pcfac = " + pcfac); 952 BigInteger p = pcfac.getIntegerModul(); 953 BigInteger q = Power.positivePower(p, k); 954 ModularRingFactory<MOD> mcfac; 955 if (ModLongRing.MAX_LONG.compareTo(q.getVal()) > 0) { 956 mcfac = (ModularRingFactory) new ModLongRing(q.getVal()); 957 } else { 958 mcfac = (ModularRingFactory) new ModIntegerRing(q.getVal()); 959 } 960 //System.out.println("mcfac = " + mcfac); 961 962 // convert C from Z[...] to Z_q[...] 963 GenPolynomialRing<MOD> qcfac = new GenPolynomialRing<MOD>(mcfac, C.ring); 964 GenPolynomial<MOD> Cq = PolyUtil.<MOD> fromIntegerCoefficients(qcfac, C); 965 //System.out.println("C = " + C); 966 //System.out.println("Cq = " + Cq); 967 968 // convert g_i from Z[...] to Z_q[...] 969 GenPolynomialRing<MOD> gcfac = new GenPolynomialRing<MOD>(mcfac, G.get(0).ring); 970 List<GenPolynomial<MOD>> GQ = new ArrayList<GenPolynomial<MOD>>(); 971 boolean allOnes = true; 972 for (GenPolynomial<BigInteger> g : G) { 973 if (!g.isONE()) { 974 allOnes = false; 975 } 976 GenPolynomial<MOD> gq = PolyUtil.<MOD> fromIntegerCoefficients(gcfac, g); 977 GQ.add(gq); 978 } 979 //System.out.println("G = " + G); 980 //System.out.println("GQ = " + GQ); 981 982 // evaluate C to Z_q[x] 983 GenPolynomialRing<MOD> pf = qcfac; 984 GenPolynomial<MOD> ap = Cq; 985 for (int j = C.ring.nvar; j > 1; j--) { 986 pf = pf.contract(1); 987 //MOD vp = mcfac.fromInteger(V.get(C.ring.nvar - j).getSymmetricInteger().getVal()); 988 MOD vp = mcfac.fromInteger(V.get(C.ring.nvar - j).getVal()); 989 //System.out.println("vp = " + vp); 990 ap = PolyUtil.<MOD> evaluateMain(pf, ap, vp); 991 //System.out.println("ap = " + ap); 992 } 993 GenPolynomial<MOD> Cq1 = ap; 994 //System.out.println("Cq1 = " + Cq1); 995 if (Cq1.isZERO()) { 996 throw new NoLiftingException("C mod (I, p^k) == 0: " + C); 997 } 998 GenPolynomialRing<BigInteger> ifac = new GenPolynomialRing<BigInteger>(new BigInteger(), pf); 999 GenPolynomial<BigInteger> Ci = PolyUtil.integerFromModularCoefficients(ifac, Cq1); 1000 //System.out.println("Ci = " + Ci); 1001 GreatestCommonDivisorAbstract<BigInteger> ufd = GCDFactory.getImplementation(new BigInteger()); 1002 Ci = Ci.abs(); 1003 BigInteger cCi = ufd.baseContent(Ci); 1004 Ci = Ci.divide(cCi); 1005 //System.out.println("cCi = " + cCi); 1006 //System.out.println("Ci = " + Ci); 1007 ////System.out.println("F.fac = " + F.get(0).ring); 1008 1009 // evaluate G to Z_q 1010 //List<GenPolynomial<MOD>> GP = new ArrayList<GenPolynomial<MOD>>(); 1011 for (GenPolynomial<MOD> gq : GQ) { 1012 GenPolynomialRing<MOD> gf = gcfac; 1013 GenPolynomial<MOD> gp = gq; 1014 for (int j = gcfac.nvar; j > 1; j--) { 1015 gf = gf.contract(1); 1016 //MOD vp = mcfac.fromInteger(V.get(gcfac.nvar - j).getSymmetricInteger().getVal()); 1017 MOD vp = mcfac.fromInteger(V.get(gcfac.nvar - j).getVal()); 1018 //System.out.println("vp = " + vp); 1019 gp = PolyUtil.<MOD> evaluateMain(gf, gp, vp); 1020 //System.out.println("gp = " + gp); 1021 } 1022 //GP.add(gp); 1023 } 1024 //System.out.println("GP = " + GP); // + ", GP.ring = " + GP.get(0).ring); 1025 1026 // leading coefficient for recursion base, for Cq1 and list GP 1027 BigInteger gi0 = Ci.leadingBaseCoefficient(); // gq0.getSymmetricInteger(); 1028 //System.out.println("gi0 = " + gi0); 1029 1030 // lift F to Z_{p^k}[x] 1031 //System.out.println("Ci = " + Ci + ", F = " + F + ", k = " + k + ", p = " + F.get(0).ring + ", gi0 = " + gi0); 1032 List<GenPolynomial<MOD>> U1 = null; 1033 if (gi0.isONE()) { 1034 U1 = HenselUtil.<MOD> liftHenselMonic(Ci, F, k); 1035 } else { 1036 U1 = HenselUtil.<MOD> liftHensel(Ci, F, k, gi0); // GI0 TODO ?? 1037 } 1038 logger.info("univariate lift: Ci = " + Ci + ", F = " + F + ", U1 = " + U1); 1039 //System.out.println("U1.fac = " + U1.get(0).ring); 1040 1041 // adjust leading coefficients of U1 with F 1042 List<GenPolynomial<BigInteger>> U1i = PolyUtil.<MOD> integerFromModularCoefficients(Ci.ring, U1); 1043 //System.out.println("U1i = " + U1i); 1044 boolean t = HenselUtil.isHenselLift(Ci, q, p, U1i); 1045 //System.out.println("isLift(U1) = " + t); 1046 if (!t) { 1047 //System.out.println("NoLiftingException, Ci = " + Ci + ", U1i = " + U1i); 1048 throw new NoLiftingException("Ci = " + Ci + ", U1i = " + U1i); 1049 } 1050 MOD cC = mcfac.fromInteger(cCi.getVal()); 1051 List<GenPolynomial<MOD>> U1f = PolyUtil.<MOD> fromIntegerCoefficients(F.get(0).ring, U1i); 1052 //System.out.println("U1f = " + U1f); 1053 List<GenPolynomial<MOD>> U1s = new ArrayList<GenPolynomial<MOD>>(U1.size()); 1054 int j = 0; 1055 int s = 0; 1056 for (GenPolynomial<MOD> u : U1) { 1057 GenPolynomial<MOD> uf = U1f.get(j); 1058 GenPolynomial<MOD> f = F.get(j); 1059 GenPolynomial<BigInteger> ui = U1i.get(j); 1060 GenPolynomial<BigInteger> gi = G.get(j); 1061 if (ui.signum() != gi.signum()) { 1062 //System.out.println("ui = " + ui + ", gi = " + gi); 1063 u = u.negate(); 1064 uf = uf.negate(); 1065 s++; 1066 } 1067 j++; 1068 if (uf.isConstant()) { 1069 //System.out.println("u = " + u); 1070 u = u.monic(); 1071 //System.out.println("u = " + u); 1072 u = u.multiply(cC); 1073 cC = cC.divide(cC); 1074 //System.out.println("u = " + u); 1075 } else { 1076 MOD x = f.leadingBaseCoefficient().divide(uf.leadingBaseCoefficient()); 1077 //System.out.println("x = " + x + ", xi = " + x.getSymmetricInteger()); 1078 if (!x.isONE()) { 1079 MOD xq = mcfac.fromInteger(x.getSymmetricInteger().getVal()); 1080 //System.out.println("xq = " + xq); 1081 u = u.multiply(xq); 1082 cC = cC.divide(xq); 1083 //System.out.println("cC = " + cC); 1084 } 1085 } 1086 U1s.add(u); 1087 } 1088 //if ( s % 2 != 0 || !cC.isONE()) { 1089 if (!cC.isONE()) { 1090 throw new NoLiftingException("s = " + s + ", Ci = " + Ci + ", U1i = " + U1i + ", cC = " + cC); 1091 } 1092 U1 = U1s; 1093 U1i = PolyUtil.<MOD> integerFromModularCoefficients(Ci.ring, U1); 1094 //System.out.println("U1i = " + U1i); 1095 U1f = PolyUtil.<MOD> fromIntegerCoefficients(F.get(0).ring, U1i); 1096 if (!F.equals(U1f)) { // evtl loop until reached 1097 System.out.println("F = " + F); 1098 System.out.println("U1f = " + U1f); 1099 throw new NoLiftingException("F = " + F + ", U1f = " + U1f); 1100 } 1101 logger.info("multivariate lift: U1 = " + U1); 1102 1103 // lift U to Z_{p^k}[x,...] 1104 //System.out.println("C = " + C + ", U1 = " + U1 + ", V = " + V + ", k = " + k + ", q = " + U1.get(0).ring + ", G = " + G); 1105 List<GenPolynomial<MOD>> U = null; 1106 if (allOnes) { 1107 U = HenselMultUtil.<MOD> liftHenselMonic(C, Cq, U1, V, k); 1108 } else { 1109 U = HenselMultUtil.<MOD> liftHensel(C, Cq, U1, V, k, G); 1110 } 1111 logger.info("multivariate lift: C = " + C + ", U1 = " + U1 + ", U = " + U); 1112 //System.out.println("U = " + U); 1113 //System.out.println("U.fac = " + U.get(0).ring); 1114 return U; 1115 } 1116 1117}