001/* 002 * $Id: HenselUtil.java 4115 2012-08-19 13:18:59Z kredel $ 003 */ 004 005package edu.jas.ufd; 006 007 008import java.util.ArrayList; 009import java.util.Arrays; 010import java.util.List; 011 012import org.apache.log4j.Logger; 013 014import edu.jas.arith.BigInteger; 015import edu.jas.arith.ModInteger; 016import edu.jas.arith.ModIntegerRing; 017import edu.jas.arith.ModLongRing; 018import edu.jas.arith.Modular; 019import edu.jas.arith.ModularRingFactory; 020import edu.jas.poly.ExpVector; 021import edu.jas.poly.GenPolynomial; 022import edu.jas.poly.GenPolynomialRing; 023import edu.jas.poly.Monomial; 024import edu.jas.poly.PolyUtil; 025import edu.jas.structure.GcdRingElem; 026import edu.jas.structure.RingFactory; 027import edu.jas.structure.Power; 028 029 030/** 031 * Hensel utilities for ufd. 032 * @author Heinz Kredel 033 */ 034 035public class HenselUtil { 036 037 038 private static final Logger logger = Logger.getLogger(HenselUtil.class); 039 040 041 private static final boolean debug = logger.isInfoEnabled(); 042 043 044 /** 045 * Modular Hensel lifting algorithm on coefficients. Let p = 046 * A.ring.coFac.modul() = B.ring.coFac.modul() and assume C == A*B mod p 047 * with gcd(A,B) == 1 mod p and S A + T B == 1 mod p. See Algorithm 6.1. in 048 * Geddes et.al.. Linear version, as it does not lift S A + T B == 1 mod 049 * p^{e+1}. 050 * @param C GenPolynomial 051 * @param A GenPolynomial 052 * @param B other GenPolynomial 053 * @param S GenPolynomial 054 * @param T GenPolynomial 055 * @param M bound on the coefficients of A1 and B1 as factors of C. 056 * @return [A1,B1,Am,Bm] = lift(C,A,B), with C = A1 * B1 mod p^e, Am = A1 057 * mod p^e, Bm = B1 mod p^e . 058 */ 059 @SuppressWarnings("unchecked") 060 public static <MOD extends GcdRingElem<MOD> & Modular> HenselApprox<MOD> liftHensel( 061 GenPolynomial<BigInteger> C, BigInteger M, GenPolynomial<MOD> A, GenPolynomial<MOD> B, 062 GenPolynomial<MOD> S, GenPolynomial<MOD> T) throws NoLiftingException { 063 if (C == null || C.isZERO()) { 064 return new HenselApprox<MOD>(C, C, A, B); 065 } 066 if (A == null || A.isZERO() || B == null || B.isZERO()) { 067 throw new IllegalArgumentException("A and B must be nonzero"); 068 } 069 GenPolynomialRing<BigInteger> fac = C.ring; 070 if (fac.nvar != 1) { // todo assert 071 throw new IllegalArgumentException("polynomial ring not univariate"); 072 } 073 // setup factories 074 GenPolynomialRing<MOD> pfac = A.ring; 075 RingFactory<MOD> p = pfac.coFac; 076 RingFactory<MOD> q = p; 077 ModularRingFactory<MOD> P = (ModularRingFactory<MOD>) p; 078 ModularRingFactory<MOD> Q = (ModularRingFactory<MOD>) q; 079 BigInteger Qi = Q.getIntegerModul(); 080 BigInteger M2 = M.multiply(M.fromInteger(2)); 081 BigInteger Mq = Qi; 082 083 // normalize c and a, b factors, assert p is prime 084 GenPolynomial<BigInteger> Ai; 085 GenPolynomial<BigInteger> Bi; 086 BigInteger c = C.leadingBaseCoefficient(); 087 C = C.multiply(c); // sic 088 MOD a = A.leadingBaseCoefficient(); 089 if (!a.isONE()) { // A = A.monic(); 090 A = A.divide(a); 091 S = S.multiply(a); 092 } 093 MOD b = B.leadingBaseCoefficient(); 094 if (!b.isONE()) { // B = B.monic(); 095 B = B.divide(b); 096 T = T.multiply(b); 097 } 098 MOD cm = P.fromInteger(c.getVal()); 099 A = A.multiply(cm); 100 B = B.multiply(cm); 101 T = T.divide(cm); 102 S = S.divide(cm); 103 Ai = PolyUtil.integerFromModularCoefficients(fac, A); 104 Bi = PolyUtil.integerFromModularCoefficients(fac, B); 105 // replace leading base coefficients 106 ExpVector ea = Ai.leadingExpVector(); 107 ExpVector eb = Bi.leadingExpVector(); 108 Ai.doPutToMap(ea, c); 109 Bi.doPutToMap(eb, c); 110 111 // polynomials mod p 112 GenPolynomial<MOD> Ap; 113 GenPolynomial<MOD> Bp; 114 GenPolynomial<MOD> A1p = A; 115 GenPolynomial<MOD> B1p = B; 116 GenPolynomial<MOD> Ep; 117 118 // polynomials over the integers 119 GenPolynomial<BigInteger> E; 120 GenPolynomial<BigInteger> Ea; 121 GenPolynomial<BigInteger> Eb; 122 GenPolynomial<BigInteger> Ea1; 123 GenPolynomial<BigInteger> Eb1; 124 125 while (Mq.compareTo(M2) < 0) { 126 // compute E=(C-AB)/q over the integers 127 E = C.subtract(Ai.multiply(Bi)); 128 if (E.isZERO()) { 129 logger.info("leaving on zero E"); 130 break; 131 } 132 try { 133 E = E.divide(Qi); 134 } catch (RuntimeException e) { 135 // useful in debuging 136 //System.out.println("C = " + C ); 137 //System.out.println("Ai = " + Ai ); 138 //System.out.println("Bi = " + Bi ); 139 //System.out.println("E = " + E ); 140 //System.out.println("Qi = " + Qi ); 141 throw e; 142 } 143 // E mod p 144 Ep = PolyUtil.<MOD> fromIntegerCoefficients(pfac, E); 145 //logger.info("Ep = " + Ep); 146 147 // construct approximation mod p 148 Ap = S.multiply(Ep); // S,T ++ T,S 149 Bp = T.multiply(Ep); 150 GenPolynomial<MOD>[] QR; 151 QR = Ap.quotientRemainder(B); 152 GenPolynomial<MOD> Qp; 153 GenPolynomial<MOD> Rp; 154 Qp = QR[0]; 155 Rp = QR[1]; 156 A1p = Rp; 157 B1p = Bp.sum(A.multiply(Qp)); 158 159 // construct q-adic approximation, convert to integer 160 Ea = PolyUtil.integerFromModularCoefficients(fac, A1p); 161 Eb = PolyUtil.integerFromModularCoefficients(fac, B1p); 162 Ea1 = Ea.multiply(Qi); 163 Eb1 = Eb.multiply(Qi); 164 165 Ea = Ai.sum(Eb1); // Eb1 and Ea1 are required 166 Eb = Bi.sum(Ea1); //-------------------------- 167 assert (Ea.degree(0) + Eb.degree(0) <= C.degree(0)); 168 //if ( Ea.degree(0)+Eb.degree(0) > C.degree(0) ) { // debug 169 // throw new RuntimeException("deg(A)+deg(B) > deg(C)"); 170 //} 171 172 // prepare for next iteration 173 Mq = Qi; 174 Qi = Q.getIntegerModul().multiply(P.getIntegerModul()); 175 // Q = new ModIntegerRing(Qi.getVal()); 176 if (ModLongRing.MAX_LONG.compareTo(Qi.getVal()) > 0) { 177 Q = (ModularRingFactory) new ModLongRing(Qi.getVal()); 178 } else { 179 Q = (ModularRingFactory) new ModIntegerRing(Qi.getVal()); 180 } 181 Ai = Ea; 182 Bi = Eb; 183 } 184 GreatestCommonDivisorAbstract<BigInteger> ufd = new GreatestCommonDivisorPrimitive<BigInteger>(); 185 186 // remove normalization 187 BigInteger ai = ufd.baseContent(Ai); 188 Ai = Ai.divide(ai); 189 BigInteger bi = null; 190 try { 191 bi = c.divide(ai); 192 Bi = Bi.divide(bi); // divide( c/a ) 193 } catch (RuntimeException e) { 194 //System.out.println("C = " + C ); 195 //System.out.println("Ai = " + Ai ); 196 //System.out.println("Bi = " + Bi ); 197 //System.out.println("c = " + c ); 198 //System.out.println("ai = " + ai ); 199 //System.out.println("bi = " + bi ); 200 //System.out.println("no exact lifting possible " + e); 201 throw new NoLiftingException("no exact lifting possible " + e); 202 } 203 return new HenselApprox<MOD>(Ai, Bi, A1p, B1p); 204 } 205 206 207 /** 208 * Modular Hensel lifting algorithm on coefficients. Let p = 209 * A.ring.coFac.modul() = B.ring.coFac.modul() and assume C == A*B mod p 210 * with gcd(A,B) == 1 mod p. See algorithm 6.1. in Geddes et.al. and 211 * algorithms 3.5.{5,6} in Cohen. 212 * @param C GenPolynomial 213 * @param A GenPolynomial 214 * @param B other GenPolynomial 215 * @param M bound on the coefficients of A1 and B1 as factors of C. 216 * @return [A1,B1] = lift(C,A,B), with C = A1 * B1. 217 */ 218 @SuppressWarnings("unchecked") 219 public static <MOD extends GcdRingElem<MOD> & Modular> HenselApprox<MOD> liftHensel( 220 GenPolynomial<BigInteger> C, BigInteger M, GenPolynomial<MOD> A, GenPolynomial<MOD> B) 221 throws NoLiftingException { 222 if (C == null || C.isZERO()) { 223 return new HenselApprox<MOD>(C, C, A, B); 224 } 225 if (A == null || A.isZERO() || B == null || B.isZERO()) { 226 throw new IllegalArgumentException("A and B must be nonzero"); 227 } 228 GenPolynomialRing<BigInteger> fac = C.ring; 229 if (fac.nvar != 1) { // todo assert 230 throw new IllegalArgumentException("polynomial ring not univariate"); 231 } 232 // one Hensel step on part polynomials 233 try { 234 GenPolynomial<MOD>[] gst = A.egcd(B); 235 if (!gst[0].isONE()) { 236 throw new NoLiftingException("A and B not coprime, gcd = " + gst[0] + ", A = " + A + ", B = " + B); 237 } 238 GenPolynomial<MOD> s = gst[1]; 239 GenPolynomial<MOD> t = gst[2]; 240 HenselApprox<MOD> ab = HenselUtil.<MOD> liftHensel(C, M, A, B, s, t); 241 return ab; 242 } catch (ArithmeticException e) { 243 throw new NoLiftingException("coefficient error " + e); 244 } 245 } 246 247 248 /** 249 * Modular quadratic Hensel lifting algorithm on coefficients. Let p = 250 * A.ring.coFac.modul() = B.ring.coFac.modul() and assume C == A*B mod p 251 * with gcd(A,B) == 1 mod p and S A + T B == 1 mod p. See algorithm 6.1. in 252 * Geddes et.al. and algorithms 3.5.{5,6} in Cohen. Quadratic version, as it 253 * also lifts S A + T B == 1 mod p^{e+1}. 254 * @param C GenPolynomial 255 * @param A GenPolynomial 256 * @param B other GenPolynomial 257 * @param S GenPolynomial 258 * @param T GenPolynomial 259 * @param M bound on the coefficients of A1 and B1 as factors of C. 260 * @return [A1,B1] = lift(C,A,B), with C = A1 * B1. 261 */ 262 @SuppressWarnings("unchecked") 263 public static <MOD extends GcdRingElem<MOD> & Modular> HenselApprox<MOD> liftHenselQuadratic( 264 GenPolynomial<BigInteger> C, BigInteger M, GenPolynomial<MOD> A, GenPolynomial<MOD> B, 265 GenPolynomial<MOD> S, GenPolynomial<MOD> T) throws NoLiftingException { 266 if (C == null || C.isZERO()) { 267 return new HenselApprox<MOD>(C, C, A, B); 268 } 269 if (A == null || A.isZERO() || B == null || B.isZERO()) { 270 throw new IllegalArgumentException("A and B must be nonzero"); 271 } 272 GenPolynomialRing<BigInteger> fac = C.ring; 273 if (fac.nvar != 1) { // todo assert 274 throw new IllegalArgumentException("polynomial ring not univariate"); 275 } 276 // setup factories 277 GenPolynomialRing<MOD> pfac = A.ring; 278 RingFactory<MOD> p = pfac.coFac; 279 RingFactory<MOD> q = p; 280 ModularRingFactory<MOD> P = (ModularRingFactory<MOD>) p; 281 ModularRingFactory<MOD> Q = (ModularRingFactory<MOD>) q; 282 BigInteger Qi = Q.getIntegerModul(); 283 BigInteger M2 = M.multiply(M.fromInteger(2)); 284 BigInteger Mq = Qi; 285 GenPolynomialRing<MOD> qfac; 286 qfac = new GenPolynomialRing<MOD>(Q, pfac); 287 288 // normalize c and a, b factors, assert p is prime 289 GenPolynomial<BigInteger> Ai; 290 GenPolynomial<BigInteger> Bi; 291 BigInteger c = C.leadingBaseCoefficient(); 292 C = C.multiply(c); // sic 293 MOD a = A.leadingBaseCoefficient(); 294 if (!a.isONE()) { // A = A.monic(); 295 A = A.divide(a); 296 S = S.multiply(a); 297 } 298 MOD b = B.leadingBaseCoefficient(); 299 if (!b.isONE()) { // B = B.monic(); 300 B = B.divide(b); 301 T = T.multiply(b); 302 } 303 MOD cm = P.fromInteger(c.getVal()); 304 A = A.multiply(cm); 305 B = B.multiply(cm); 306 T = T.divide(cm); 307 S = S.divide(cm); 308 Ai = PolyUtil.integerFromModularCoefficients(fac, A); 309 Bi = PolyUtil.integerFromModularCoefficients(fac, B); 310 // replace leading base coefficients 311 ExpVector ea = Ai.leadingExpVector(); 312 ExpVector eb = Bi.leadingExpVector(); 313 Ai.doPutToMap(ea, c); 314 Bi.doPutToMap(eb, c); 315 316 // polynomials mod p 317 GenPolynomial<MOD> Ap; 318 GenPolynomial<MOD> Bp; 319 GenPolynomial<MOD> A1p = A; 320 GenPolynomial<MOD> B1p = B; 321 GenPolynomial<MOD> Ep; 322 GenPolynomial<MOD> Sp = S; 323 GenPolynomial<MOD> Tp = T; 324 325 // polynomials mod q 326 GenPolynomial<MOD> Aq; 327 GenPolynomial<MOD> Bq; 328 GenPolynomial<MOD> Eq; 329 330 // polynomials over the integers 331 GenPolynomial<BigInteger> E; 332 GenPolynomial<BigInteger> Ea; 333 GenPolynomial<BigInteger> Eb; 334 GenPolynomial<BigInteger> Ea1; 335 GenPolynomial<BigInteger> Eb1; 336 GenPolynomial<BigInteger> Si; 337 GenPolynomial<BigInteger> Ti; 338 339 Si = PolyUtil.integerFromModularCoefficients(fac, S); 340 Ti = PolyUtil.integerFromModularCoefficients(fac, T); 341 342 Aq = PolyUtil.<MOD> fromIntegerCoefficients(qfac, Ai); 343 Bq = PolyUtil.<MOD> fromIntegerCoefficients(qfac, Bi); 344 345 while (Mq.compareTo(M2) < 0) { 346 // compute E=(C-AB)/q over the integers 347 E = C.subtract(Ai.multiply(Bi)); 348 if (E.isZERO()) { 349 logger.info("leaving on zero E"); 350 break; 351 } 352 E = E.divide(Qi); 353 // E mod p 354 Ep = PolyUtil.<MOD> fromIntegerCoefficients(qfac, E); 355 //logger.info("Ep = " + Ep + ", qfac = " + qfac); 356 //if (Ep.isZERO()) { 357 //System.out.println("leaving on zero error"); 358 //??logger.info("leaving on zero Ep"); 359 //??break; 360 //} 361 362 // construct approximation mod p 363 Ap = Sp.multiply(Ep); // S,T ++ T,S 364 Bp = Tp.multiply(Ep); 365 GenPolynomial<MOD>[] QR; 366 //logger.info("Ap = " + Ap + ", Bp = " + Bp + ", fac(Ap) = " + Ap.ring); 367 QR = Ap.quotientRemainder(Bq); 368 GenPolynomial<MOD> Qp; 369 GenPolynomial<MOD> Rp; 370 Qp = QR[0]; 371 Rp = QR[1]; 372 //logger.info("Qp = " + Qp + ", Rp = " + Rp); 373 A1p = Rp; 374 B1p = Bp.sum(Aq.multiply(Qp)); 375 376 // construct q-adic approximation, convert to integer 377 Ea = PolyUtil.integerFromModularCoefficients(fac, A1p); 378 Eb = PolyUtil.integerFromModularCoefficients(fac, B1p); 379 Ea1 = Ea.multiply(Qi); 380 Eb1 = Eb.multiply(Qi); 381 Ea = Ai.sum(Eb1); // Eb1 and Ea1 are required 382 Eb = Bi.sum(Ea1); //-------------------------- 383 assert (Ea.degree(0) + Eb.degree(0) <= C.degree(0)); 384 //if ( Ea.degree(0)+Eb.degree(0) > C.degree(0) ) { // debug 385 // throw new RuntimeException("deg(A)+deg(B) > deg(C)"); 386 //} 387 Ai = Ea; 388 Bi = Eb; 389 390 // gcd representation factors error -------------------------------- 391 // compute E=(1-SA-TB)/q over the integers 392 E = fac.getONE(); 393 E = E.subtract(Si.multiply(Ai)).subtract(Ti.multiply(Bi)); 394 E = E.divide(Qi); 395 // E mod q 396 Ep = PolyUtil.<MOD> fromIntegerCoefficients(qfac, E); 397 //logger.info("Ep2 = " + Ep); 398 399 // construct approximation mod q 400 Ap = Sp.multiply(Ep); // S,T ++ T,S 401 Bp = Tp.multiply(Ep); 402 QR = Bp.quotientRemainder(Aq); // Ai == A mod p ? 403 Qp = QR[0]; 404 Rp = QR[1]; 405 B1p = Rp; 406 A1p = Ap.sum(Bq.multiply(Qp)); 407 408 if (false && debug) { 409 Eq = A1p.multiply(Aq).sum(B1p.multiply(Bq)).subtract(Ep); 410 if (!Eq.isZERO()) { 411 System.out.println("A*A1p+B*B1p-Ep2 != 0 " + Eq); 412 throw new RuntimeException("A*A1p+B*B1p-Ep2 != 0 mod " + Q.getIntegerModul()); 413 } 414 } 415 416 // construct q-adic approximation, convert to integer 417 Ea = PolyUtil.integerFromModularCoefficients(fac, A1p); 418 Eb = PolyUtil.integerFromModularCoefficients(fac, B1p); 419 Ea1 = Ea.multiply(Qi); 420 Eb1 = Eb.multiply(Qi); 421 Ea = Si.sum(Ea1); // Eb1 and Ea1 are required 422 Eb = Ti.sum(Eb1); //-------------------------- 423 Si = Ea; 424 Ti = Eb; 425 426 // prepare for next iteration 427 Mq = Qi; 428 Qi = Q.getIntegerModul().multiply(Q.getIntegerModul()); 429 if (ModLongRing.MAX_LONG.compareTo(Qi.getVal()) > 0) { 430 Q = (ModularRingFactory) new ModLongRing(Qi.getVal()); 431 } else { 432 Q = (ModularRingFactory) new ModIntegerRing(Qi.getVal()); 433 } 434 //Q = new ModIntegerRing(Qi.getVal()); 435 //System.out.println("Q = " + Q + ", from Q = " + Mq); 436 437 qfac = new GenPolynomialRing<MOD>(Q, pfac); 438 439 Aq = PolyUtil.<MOD> fromIntegerCoefficients(qfac, Ai); 440 Bq = PolyUtil.<MOD> fromIntegerCoefficients(qfac, Bi); 441 Sp = PolyUtil.<MOD> fromIntegerCoefficients(qfac, Si); 442 Tp = PolyUtil.<MOD> fromIntegerCoefficients(qfac, Ti); 443 if (false && debug) { 444 E = Ai.multiply(Si).sum(Bi.multiply(Ti)); 445 Eq = PolyUtil.<MOD> fromIntegerCoefficients(qfac, E); 446 if (!Eq.isONE()) { 447 System.out.println("Ai*Si+Bi*Ti=1 " + Eq); 448 throw new RuntimeException("Ai*Si+Bi*Ti != 1 mod " + Q.getIntegerModul()); 449 } 450 } 451 } 452 GreatestCommonDivisorAbstract<BigInteger> ufd = new GreatestCommonDivisorPrimitive<BigInteger>(); 453 454 // remove normalization if possible 455 BigInteger ai = ufd.baseContent(Ai); 456 Ai = Ai.divide(ai); 457 BigInteger bi = null; 458 try { 459 bi = c.divide(ai); 460 Bi = Bi.divide(bi); // divide( c/a ) 461 } catch (RuntimeException e) { 462 //System.out.println("C = " + C ); 463 //System.out.println("Ai = " + Ai ); 464 //System.out.println("Bi = " + Bi ); 465 //System.out.println("c = " + c ); 466 //System.out.println("ai = " + ai ); 467 //System.out.println("bi = " + bi ); 468 //System.out.println("no exact lifting possible " + e); 469 throw new NoLiftingException("no exact lifting possible " + e); 470 } 471 return new HenselApprox<MOD>(Ai, Bi, A1p, B1p); 472 } 473 474 475 /** 476 * Modular quadratic Hensel lifting algorithm on coefficients. Let p = 477 * A.ring.coFac.modul() = B.ring.coFac.modul() and assume C == A*B mod p 478 * with gcd(A,B) == 1 mod p. See algorithm 6.1. in Geddes et.al. and 479 * algorithms 3.5.{5,6} in Cohen. Quadratic version. 480 * @param C GenPolynomial 481 * @param A GenPolynomial 482 * @param B other GenPolynomial 483 * @param M bound on the coefficients of A1 and B1 as factors of C. 484 * @return [A1,B1] = lift(C,A,B), with C = A1 * B1. 485 */ 486 @SuppressWarnings("unchecked") 487 public static <MOD extends GcdRingElem<MOD> & Modular> HenselApprox<MOD> liftHenselQuadratic( 488 GenPolynomial<BigInteger> C, BigInteger M, GenPolynomial<MOD> A, GenPolynomial<MOD> B) 489 throws NoLiftingException { 490 if (C == null || C.isZERO()) { 491 return new HenselApprox<MOD>(C, C, A, B); 492 } 493 if (A == null || A.isZERO() || B == null || B.isZERO()) { 494 throw new IllegalArgumentException("A and B must be nonzero"); 495 } 496 GenPolynomialRing<BigInteger> fac = C.ring; 497 if (fac.nvar != 1) { // todo assert 498 throw new IllegalArgumentException("polynomial ring not univariate"); 499 } 500 // one Hensel step on part polynomials 501 try { 502 GenPolynomial<MOD>[] gst = A.egcd(B); 503 if (!gst[0].isONE()) { 504 throw new NoLiftingException("A and B not coprime, gcd = " + gst[0] + ", A = " + A + ", B = " + B); 505 } 506 GenPolynomial<MOD> s = gst[1]; 507 GenPolynomial<MOD> t = gst[2]; 508 HenselApprox<MOD> ab = HenselUtil.<MOD> liftHenselQuadratic(C, M, A, B, s, t); 509 return ab; 510 } catch (ArithmeticException e) { 511 throw new NoLiftingException("coefficient error " + e); 512 } 513 } 514 515 516 /** 517 * Modular Hensel lifting algorithm on coefficients. Let p = 518 * A.ring.coFac.modul() = B.ring.coFac.modul() and assume C == A*B mod p 519 * with gcd(A,B) == 1 mod p. See algorithm 6.1. in Geddes et.al. and 520 * algorithms 3.5.{5,6} in Cohen. Quadratic version. 521 * @param C GenPolynomial 522 * @param A GenPolynomial 523 * @param B other GenPolynomial 524 * @param M bound on the coefficients of A1 and B1 as factors of C. 525 * @return [A1,B1] = lift(C,A,B), with C = A1 * B1. 526 */ 527 @SuppressWarnings("unchecked") 528 public static <MOD extends GcdRingElem<MOD> & Modular> HenselApprox<MOD> liftHenselQuadraticFac( 529 GenPolynomial<BigInteger> C, BigInteger M, GenPolynomial<MOD> A, GenPolynomial<MOD> B) 530 throws NoLiftingException { 531 if (C == null || C.isZERO()) { 532 throw new IllegalArgumentException("C must be nonzero"); 533 } 534 if (A == null || A.isZERO() || B == null || B.isZERO()) { 535 throw new IllegalArgumentException("A and B must be nonzero"); 536 } 537 GenPolynomialRing<BigInteger> fac = C.ring; 538 if (fac.nvar != 1) { // todo assert 539 throw new IllegalArgumentException("polynomial ring not univariate"); 540 } 541 // one Hensel step on part polynomials 542 try { 543 GenPolynomial<MOD>[] gst = A.egcd(B); 544 if (!gst[0].isONE()) { 545 throw new NoLiftingException("A and B not coprime, gcd = " + gst[0] + ", A = " + A + ", B = " + B); 546 } 547 GenPolynomial<MOD> s = gst[1]; 548 GenPolynomial<MOD> t = gst[2]; 549 HenselApprox<MOD> ab = HenselUtil.<MOD> liftHenselQuadraticFac(C, M, A, B, s, t); 550 return ab; 551 } catch (ArithmeticException e) { 552 throw new NoLiftingException("coefficient error " + e); 553 } 554 } 555 556 557 /** 558 * Modular Hensel lifting algorithm on coefficients. Let p = 559 * A.ring.coFac.modul() = B.ring.coFac.modul() and assume C == A*B mod p 560 * with gcd(A,B) == 1 mod p and S A + T B == 1 mod p. See algorithm 6.1. in 561 * Geddes et.al. and algorithms 3.5.{5,6} in Cohen. Quadratic version, as it 562 * also lifts S A + T B == 1 mod p^{e+1}. 563 * @param C primitive GenPolynomial 564 * @param A GenPolynomial 565 * @param B other GenPolynomial 566 * @param S GenPolynomial 567 * @param T GenPolynomial 568 * @param M bound on the coefficients of A1 and B1 as factors of C. 569 * @return [A1,B1] = lift(C,A,B), with C = A1 * B1. 570 */ 571 @SuppressWarnings("unchecked") 572 public static <MOD extends GcdRingElem<MOD> & Modular> HenselApprox<MOD> liftHenselQuadraticFac( 573 GenPolynomial<BigInteger> C, BigInteger M, GenPolynomial<MOD> A, GenPolynomial<MOD> B, 574 GenPolynomial<MOD> S, GenPolynomial<MOD> T) throws NoLiftingException { 575 //System.out.println("*** version for factorization *** "); 576 //GenPolynomial<BigInteger>[] AB = new GenPolynomial[2]; 577 if (C == null || C.isZERO()) { 578 throw new IllegalArgumentException("C must be nonzero"); 579 } 580 if (A == null || A.isZERO() || B == null || B.isZERO()) { 581 throw new IllegalArgumentException("A and B must be nonzero"); 582 } 583 GenPolynomialRing<BigInteger> fac = C.ring; 584 if (fac.nvar != 1) { // todo assert 585 throw new IllegalArgumentException("polynomial ring not univariate"); 586 } 587 // setup factories 588 GenPolynomialRing<MOD> pfac = A.ring; 589 RingFactory<MOD> p = pfac.coFac; 590 RingFactory<MOD> q = p; 591 ModularRingFactory<MOD> P = (ModularRingFactory<MOD>) p; 592 ModularRingFactory<MOD> Q = (ModularRingFactory<MOD>) q; 593 BigInteger PP = Q.getIntegerModul(); 594 BigInteger Qi = PP; 595 BigInteger M2 = M.multiply(M.fromInteger(2)); 596 if (debug) { 597 //System.out.println("M2 = " + M2); 598 logger.debug("M2 = " + M2); 599 } 600 BigInteger Mq = Qi; 601 GenPolynomialRing<MOD> qfac; 602 qfac = new GenPolynomialRing<MOD>(Q, pfac); // mod p 603 GenPolynomialRing<MOD> mfac; 604 BigInteger Mi = Q.getIntegerModul().multiply(Q.getIntegerModul()); 605 ModularRingFactory<MOD> Qmm; 606 // = new ModIntegerRing(Mi.getVal()); 607 if (ModLongRing.MAX_LONG.compareTo(Mi.getVal()) > 0) { 608 Qmm = (ModularRingFactory) new ModLongRing(Mi.getVal()); 609 } else { 610 Qmm = (ModularRingFactory) new ModIntegerRing(Mi.getVal()); 611 } 612 mfac = new GenPolynomialRing<MOD>(Qmm, qfac); // mod p^e 613 MOD Qm = Qmm.fromInteger(Qi.getVal()); 614 615 // partly normalize c and a, b factors, assert p is prime 616 GenPolynomial<BigInteger> Ai; 617 GenPolynomial<BigInteger> Bi; 618 BigInteger c = C.leadingBaseCoefficient(); 619 C = C.multiply(c); // sic 620 MOD a = A.leadingBaseCoefficient(); 621 if (!a.isONE()) { // A = A.monic(); 622 A = A.divide(a); 623 S = S.multiply(a); 624 } 625 MOD b = B.leadingBaseCoefficient(); 626 if (!b.isONE()) { // B = B.monic(); 627 B = B.divide(b); 628 T = T.multiply(b); 629 } 630 MOD cm = P.fromInteger(c.getVal()); 631 if (cm.isZERO()) { 632 System.out.println("c = " + c); 633 System.out.println("P = " + P); 634 throw new ArithmeticException("c mod p == 0 not meaningful"); 635 } 636 // mod p 637 A = A.multiply(cm); 638 S = S.divide(cm); 639 B = B.multiply(cm); 640 T = T.divide(cm); 641 Ai = PolyUtil.integerFromModularCoefficients(fac, A); 642 Bi = PolyUtil.integerFromModularCoefficients(fac, B); 643 // replace leading base coefficients 644 ExpVector ea = Ai.leadingExpVector(); 645 ExpVector eb = Bi.leadingExpVector(); 646 Ai.doPutToMap(ea, c); 647 Bi.doPutToMap(eb, c); 648 649 // polynomials mod p 650 GenPolynomial<MOD> Ap; 651 GenPolynomial<MOD> Bp; 652 GenPolynomial<MOD> A1p = A; 653 GenPolynomial<MOD> B1p = B; 654 GenPolynomial<MOD> Sp = S; 655 GenPolynomial<MOD> Tp = T; 656 657 // polynomials mod q 658 GenPolynomial<MOD> Aq; 659 GenPolynomial<MOD> Bq; 660 661 // polynomials mod p^e 662 GenPolynomial<MOD> Cm; 663 GenPolynomial<MOD> Am; 664 GenPolynomial<MOD> Bm; 665 GenPolynomial<MOD> Em; 666 GenPolynomial<MOD> Emp; 667 GenPolynomial<MOD> Sm; 668 GenPolynomial<MOD> Tm; 669 GenPolynomial<MOD> Ema; 670 GenPolynomial<MOD> Emb; 671 GenPolynomial<MOD> Ema1; 672 GenPolynomial<MOD> Emb1; 673 674 // polynomials over the integers 675 GenPolynomial<BigInteger> Ei; 676 GenPolynomial<BigInteger> Si; 677 GenPolynomial<BigInteger> Ti; 678 679 Si = PolyUtil.integerFromModularCoefficients(fac, S); 680 Ti = PolyUtil.integerFromModularCoefficients(fac, T); 681 682 Aq = PolyUtil.<MOD> fromIntegerCoefficients(qfac, Ai); 683 Bq = PolyUtil.<MOD> fromIntegerCoefficients(qfac, Bi); 684 685 // polynomials mod p^e 686 Cm = PolyUtil.<MOD> fromIntegerCoefficients(mfac, C); 687 Am = PolyUtil.<MOD> fromIntegerCoefficients(mfac, Ai); 688 Bm = PolyUtil.<MOD> fromIntegerCoefficients(mfac, Bi); 689 Sm = PolyUtil.<MOD> fromIntegerCoefficients(mfac, Si); 690 Tm = PolyUtil.<MOD> fromIntegerCoefficients(mfac, Ti); 691 //System.out.println("Cm = " + Cm); 692 //System.out.println("Am = " + Am); 693 //System.out.println("Bm = " + Bm); 694 //System.out.println("Ai = " + Ai); 695 //System.out.println("Bi = " + Bi); 696 //System.out.println("mfac = " + mfac); 697 698 while (Mq.compareTo(M2) < 0) { 699 // compute E=(C-AB)/p mod p^e 700 if (debug) { 701 //System.out.println("mfac = " + Cm.ring); 702 logger.debug("mfac = " + Cm.ring); 703 } 704 Em = Cm.subtract(Am.multiply(Bm)); 705 //System.out.println("Em = " + Em); 706 if (Em.isZERO()) { 707 if (C.subtract(Ai.multiply(Bi)).isZERO()) { 708 logger.info("leaving on zero E"); 709 break; 710 } 711 } 712 // Em = Em.divide( Qm ); 713 Ei = PolyUtil.integerFromModularCoefficients(fac, Em); 714 Ei = Ei.divide(Qi); 715 //System.out.println("Ei = " + Ei); 716 717 // Ei mod p 718 Emp = PolyUtil.<MOD> fromIntegerCoefficients(qfac, Ei); 719 // Emp = PolyUtil.<MOD>fromIntegerCoefficients(qfac, 720 // PolyUtil.integerFromModularCoefficients(fac,Em) ); 721 //System.out.println("Emp = " + Emp); 722 //logger.info("Emp = " + Emp); 723 if (Emp.isZERO()) { 724 if (C.subtract(Ai.multiply(Bi)).isZERO()) { 725 logger.info("leaving on zero Emp"); 726 break; 727 } 728 } 729 730 // construct approximation mod p 731 Ap = Sp.multiply(Emp); // S,T ++ T,S 732 Bp = Tp.multiply(Emp); 733 GenPolynomial<MOD>[] QR = null; 734 QR = Ap.quotientRemainder(Bq); // Bq ! 735 GenPolynomial<MOD> Qp = QR[0]; 736 GenPolynomial<MOD> Rp = QR[1]; 737 A1p = Rp; 738 B1p = Bp.sum(Aq.multiply(Qp)); // Aq ! 739 //System.out.println("A1p = " + A1p); 740 //System.out.println("B1p = " + B1p); 741 742 // construct q-adic approximation 743 Ema = PolyUtil.<MOD> fromIntegerCoefficients(mfac, PolyUtil.integerFromModularCoefficients(fac, 744 A1p)); 745 Emb = PolyUtil.<MOD> fromIntegerCoefficients(mfac, PolyUtil.integerFromModularCoefficients(fac, 746 B1p)); 747 //System.out.println("Ema = " + Ema); 748 //System.out.println("Emb = " + Emb); 749 Ema1 = Ema.multiply(Qm); 750 Emb1 = Emb.multiply(Qm); 751 Ema = Am.sum(Emb1); // Eb1 and Ea1 are required 752 Emb = Bm.sum(Ema1); //-------------------------- 753 assert (Ema.degree(0) + Emb.degree(0) <= Cm.degree(0)); 754 //if ( Ema.degree(0)+Emb.degree(0) > Cm.degree(0) ) { // debug 755 // throw new RuntimeException("deg(A)+deg(B) > deg(C)"); 756 //} 757 Am = Ema; 758 Bm = Emb; 759 Ai = PolyUtil.integerFromModularCoefficients(fac, Am); 760 Bi = PolyUtil.integerFromModularCoefficients(fac, Bm); 761 //System.out.println("Am = " + Am); 762 //System.out.println("Bm = " + Bm); 763 //System.out.println("Ai = " + Ai); 764 //System.out.println("Bi = " + Bi + "\n"); 765 766 // gcd representation factors error -------------------------------- 767 // compute E=(1-SA-TB)/p mod p^e 768 Em = mfac.getONE(); 769 Em = Em.subtract(Sm.multiply(Am)).subtract(Tm.multiply(Bm)); 770 //System.out.println("Em = " + Em); 771 // Em = Em.divide( Pm ); 772 773 Ei = PolyUtil.integerFromModularCoefficients(fac, Em); 774 Ei = Ei.divide(Qi); 775 //System.out.println("Ei = " + Ei); 776 // Ei mod p 777 Emp = PolyUtil.<MOD> fromIntegerCoefficients(qfac, Ei); 778 //Emp = PolyUtil.<MOD>fromIntegerCoefficients(qfac, 779 // PolyUtil.integerFromModularCoefficients(fac,Em) ); 780 //System.out.println("Emp = " + Emp); 781 782 // construct approximation mod p 783 Ap = Sp.multiply(Emp); // S,T ++ T,S // Ep Eqp 784 Bp = Tp.multiply(Emp); // Ep Eqp 785 QR = Bp.quotientRemainder(Aq); // Ap Aq ! // Ai == A mod p ? 786 Qp = QR[0]; 787 Rp = QR[1]; 788 B1p = Rp; 789 A1p = Ap.sum(Bq.multiply(Qp)); 790 //System.out.println("A1p = " + A1p); 791 //System.out.println("B1p = " + B1p); 792 793 // construct p^e-adic approximation 794 Ema = PolyUtil.<MOD> fromIntegerCoefficients(mfac, PolyUtil.integerFromModularCoefficients(fac, 795 A1p)); 796 Emb = PolyUtil.<MOD> fromIntegerCoefficients(mfac, PolyUtil.integerFromModularCoefficients(fac, 797 B1p)); 798 Ema1 = Ema.multiply(Qm); 799 Emb1 = Emb.multiply(Qm); 800 Ema = Sm.sum(Ema1); // Emb1 and Ema1 are required 801 Emb = Tm.sum(Emb1); //-------------------------- 802 Sm = Ema; 803 Tm = Emb; 804 Si = PolyUtil.integerFromModularCoefficients(fac, Sm); 805 Ti = PolyUtil.integerFromModularCoefficients(fac, Tm); 806 //System.out.println("Sm = " + Sm); 807 //System.out.println("Tm = " + Tm); 808 //System.out.println("Si = " + Si); 809 //System.out.println("Ti = " + Ti + "\n"); 810 811 // prepare for next iteration 812 Qi = Q.getIntegerModul().multiply(Q.getIntegerModul()); 813 Mq = Qi; 814 //lmfac = mfac; 815 // Q = new ModIntegerRing(Qi.getVal()); 816 if (ModLongRing.MAX_LONG.compareTo(Qi.getVal()) > 0) { 817 Q = (ModularRingFactory) new ModLongRing(Qi.getVal()); 818 } else { 819 Q = (ModularRingFactory) new ModIntegerRing(Qi.getVal()); 820 } 821 qfac = new GenPolynomialRing<MOD>(Q, pfac); 822 BigInteger Qmmi = Qmm.getIntegerModul().multiply(Qmm.getIntegerModul()); 823 //Qmm = new ModIntegerRing(Qmmi.getVal()); 824 if (ModLongRing.MAX_LONG.compareTo(Qmmi.getVal()) > 0) { 825 Qmm = (ModularRingFactory) new ModLongRing(Qmmi.getVal()); 826 } else { 827 Qmm = (ModularRingFactory) new ModIntegerRing(Qmmi.getVal()); 828 } 829 mfac = new GenPolynomialRing<MOD>(Qmm, qfac); 830 Qm = Qmm.fromInteger(Qi.getVal()); 831 832 Cm = PolyUtil.<MOD> fromIntegerCoefficients(mfac, C); 833 Am = PolyUtil.<MOD> fromIntegerCoefficients(mfac, Ai); 834 Bm = PolyUtil.<MOD> fromIntegerCoefficients(mfac, Bi); 835 Sm = PolyUtil.<MOD> fromIntegerCoefficients(mfac, Si); 836 Tm = PolyUtil.<MOD> fromIntegerCoefficients(mfac, Ti); 837 838 assert (isHenselLift(C, Mi, PP, Ai, Bi)); 839 Mi = Mi.fromInteger(Qmm.getIntegerModul().getVal()); 840 841 Aq = PolyUtil.<MOD> fromIntegerCoefficients(qfac, Ai); 842 Bq = PolyUtil.<MOD> fromIntegerCoefficients(qfac, Bi); 843 Sp = PolyUtil.<MOD> fromIntegerCoefficients(qfac, Si); 844 Tp = PolyUtil.<MOD> fromIntegerCoefficients(qfac, Ti); 845 846 //System.out.println("Am = " + Am); 847 //System.out.println("Bm = " + Bm); 848 //System.out.println("Sm = " + Sm); 849 //System.out.println("Tm = " + Tm); 850 //System.out.println("mfac = " + mfac); 851 //System.out.println("Qmm = " + Qmm + ", M2 = " + M2 + ", Mq = " + Mq + "\n"); 852 } 853 //System.out.println("*Ai = " + Ai); 854 //System.out.println("*Bi = " + Bi); 855 856 Em = Cm.subtract(Am.multiply(Bm)); 857 if (!Em.isZERO()) { 858 System.out.println("Em = " + Em); 859 //throw new NoLiftingException("no exact lifting possible"); 860 } 861 // remove normalization not possible when not exact factorization 862 GreatestCommonDivisorAbstract<BigInteger> ufd = new GreatestCommonDivisorPrimitive<BigInteger>(); 863 // remove normalization if possible 864 BigInteger ai = ufd.baseContent(Ai); 865 Ai = Ai.divide(ai); // Ai=pp(Ai) 866 BigInteger[] qr = c.quotientRemainder(ai); 867 BigInteger bi = null; 868 boolean exact = true; 869 if (qr[1].isZERO()) { 870 bi = qr[0]; 871 try { 872 Bi = Bi.divide(bi); // divide( c/a ) 873 } catch (RuntimeException e) { 874 System.out.println("*catch: no exact factorization: " + bi + ", e = " + e); 875 exact = false; 876 } 877 } else { 878 System.out.println("*remainder: no exact factorization: q = " + qr[0] + ", r = " + qr[1]); 879 exact = false; 880 } 881 if (!exact) { 882 System.out.println("*Ai = " + Ai); 883 System.out.println("*ai = " + ai); 884 System.out.println("*Bi = " + Bi); 885 System.out.println("*bi = " + bi); 886 System.out.println("*c = " + c); 887 throw new NoLiftingException("no exact lifting possible"); 888 } 889 return new HenselApprox<MOD>(Ai, Bi, Aq, Bq); 890 } 891 892 893 /** 894 * Modular Hensel lifting test. Let p be a prime number and assume C == 895 * prod_{0,...,n-1} g_i mod p with gcd(g_i,g_j) == 1 mod p for i != j. 896 * @param C GenPolynomial 897 * @param G = [g_0,...,g_{n-1}] list of GenPolynomial 898 * @param M bound on the coefficients of g_i as factors of C. 899 * @param p prime number. 900 * @return true if C = prod_{0,...,n-1} g_i mod p^e, else false. 901 */ 902 public static//<MOD extends GcdRingElem<MOD> & Modular> 903 boolean isHenselLift(GenPolynomial<BigInteger> C, BigInteger M, BigInteger p, 904 List<GenPolynomial<BigInteger>> G) { 905 if (C == null || C.isZERO()) { 906 return false; 907 } 908 GenPolynomialRing<BigInteger> pfac = C.ring; 909 ModIntegerRing pm = new ModIntegerRing(p.getVal(), true); 910 GenPolynomialRing<ModInteger> mfac = new GenPolynomialRing<ModInteger>(pm, pfac); 911 912 // check mod p 913 GenPolynomial<ModInteger> cl = mfac.getONE(); 914 GenPolynomial<ModInteger> hlp; 915 for (GenPolynomial<BigInteger> hl : G) { 916 //System.out.println("hl = " + hl); 917 hlp = PolyUtil.<ModInteger> fromIntegerCoefficients(mfac, hl); 918 //System.out.println("hl mod p = " + hlp); 919 cl = cl.multiply(hlp); 920 } 921 GenPolynomial<ModInteger> cp = PolyUtil.<ModInteger> fromIntegerCoefficients(mfac, C); 922 if (!cp.equals(cl)) { 923 System.out.println("Hensel precondition wrong!"); 924 System.out.println("cl = " + cl); 925 System.out.println("cp = " + cp); 926 System.out.println("mon(cl) = " + cl.monic()); 927 System.out.println("mon(cp) = " + cp.monic()); 928 System.out.println("cp-cl = " + cp.subtract(cl)); 929 System.out.println("M = " + M + ", p = " + p); 930 return false; 931 } 932 933 // check mod p^e 934 BigInteger mip = p; 935 while (mip.compareTo(M) < 0) { 936 mip = mip.multiply(mip); // p 937 } 938 // mip = mip.multiply(p); 939 pm = new ModIntegerRing(mip.getVal(), false); 940 mfac = new GenPolynomialRing<ModInteger>(pm, pfac); 941 cl = mfac.getONE(); 942 for (GenPolynomial<BigInteger> hl : G) { 943 //System.out.println("hl = " + hl); 944 hlp = PolyUtil.<ModInteger> fromIntegerCoefficients(mfac, hl); 945 //System.out.println("hl mod p^e = " + hlp); 946 cl = cl.multiply(hlp); 947 } 948 cp = PolyUtil.<ModInteger> fromIntegerCoefficients(mfac, C); 949 if (!cp.equals(cl)) { 950 System.out.println("Hensel post condition wrong!"); 951 System.out.println("cl = " + cl); 952 System.out.println("cp = " + cp); 953 System.out.println("cp-cl = " + cp.subtract(cl)); 954 System.out.println("M = " + M + ", p = " + p + ", p^e = " + mip); 955 return false; 956 } 957 return true; 958 } 959 960 961 /** 962 * Modular Hensel lifting test. Let p be a prime number and assume C == A * 963 * B mod p with gcd(A,B) == 1 mod p. 964 * @param C GenPolynomial 965 * @param A GenPolynomial 966 * @param B GenPolynomial 967 * @param M bound on the coefficients of A and B as factors of C. 968 * @param p prime number. 969 * @return true if C = A * B mod p**e, else false. 970 */ 971 public static//<MOD extends GcdRingElem<MOD> & Modular> 972 boolean isHenselLift(GenPolynomial<BigInteger> C, BigInteger M, BigInteger p, 973 GenPolynomial<BigInteger> A, GenPolynomial<BigInteger> B) { 974 List<GenPolynomial<BigInteger>> G = new ArrayList<GenPolynomial<BigInteger>>(2); 975 G.add(A); 976 G.add(B); 977 return isHenselLift(C, M, p, G); 978 } 979 980 981 /** 982 * Modular Hensel lifting test. Let p be a prime number and assume C == A * 983 * B mod p with gcd(A,B) == 1 mod p. 984 * @param C GenPolynomial 985 * @param Ha Hensel approximation. 986 * @param M bound on the coefficients of A and B as factors of C. 987 * @param p prime number. 988 * @return true if C = A * B mod p^e, else false. 989 */ 990 public static <MOD extends GcdRingElem<MOD> & Modular> boolean isHenselLift(GenPolynomial<BigInteger> C, 991 BigInteger M, BigInteger p, HenselApprox<MOD> Ha) { 992 List<GenPolynomial<BigInteger>> G = new ArrayList<GenPolynomial<BigInteger>>(2); 993 G.add(Ha.A); 994 G.add(Ha.B); 995 return isHenselLift(C, M, p, G); 996 } 997 998 999 /** 1000 * Constructing and lifting algorithm for extended Euclidean relation. Let p 1001 * = A.ring.coFac.modul() and assume gcd(A,B) == 1 mod p. 1002 * @param A modular GenPolynomial 1003 * @param B modular GenPolynomial 1004 * @param k desired approximation exponent p^k. 1005 * @return [s,t] with s A + t B = 1 mod p^k. 1006 */ 1007 public static <MOD extends GcdRingElem<MOD> & Modular> GenPolynomial<MOD>[] liftExtendedEuclidean( 1008 GenPolynomial<MOD> A, GenPolynomial<MOD> B, long k) throws NoLiftingException { 1009 if (A == null || A.isZERO() || B == null || B.isZERO()) { 1010 throw new IllegalArgumentException("A and B must be nonzero, A = " + A + ", B = " + B); 1011 } 1012 GenPolynomialRing<MOD> fac = A.ring; 1013 if (fac.nvar != 1) { // todo assert 1014 throw new IllegalArgumentException("polynomial ring not univariate"); 1015 } 1016 // start with extended Euclidean relation mod p 1017 GenPolynomial<MOD>[] gst = null; 1018 try { 1019 gst = A.egcd(B); 1020 if (!gst[0].isONE()) { 1021 throw new NoLiftingException("A and B not coprime, gcd = " + gst[0] + ", A = " + A + ", B = " + B); 1022 } 1023 } catch (ArithmeticException e) { 1024 throw new NoLiftingException("coefficient error " + e); 1025 } 1026 GenPolynomial<MOD> S = gst[1]; 1027 GenPolynomial<MOD> T = gst[2]; 1028 //System.out.println("eeS = " + S + ": " + S.ring.coFac); 1029 //System.out.println("eeT = " + T + ": " + T.ring.coFac); 1030 1031 // setup integer polynomial ring 1032 GenPolynomialRing<BigInteger> ifac = new GenPolynomialRing<BigInteger>(new BigInteger(), fac); 1033 GenPolynomial<BigInteger> one = ifac.getONE(); 1034 GenPolynomial<BigInteger> Ai = PolyUtil.integerFromModularCoefficients(ifac, A); 1035 GenPolynomial<BigInteger> Bi = PolyUtil.integerFromModularCoefficients(ifac, B); 1036 GenPolynomial<BigInteger> Si = PolyUtil.integerFromModularCoefficients(ifac, S); 1037 GenPolynomial<BigInteger> Ti = PolyUtil.integerFromModularCoefficients(ifac, T); 1038 //System.out.println("Ai = " + Ai); 1039 //System.out.println("Bi = " + Bi); 1040 //System.out.println("Si = " + Si); 1041 //System.out.println("Ti = " + Ti); 1042 1043 // approximate mod p^i 1044 ModularRingFactory<MOD> mcfac = (ModularRingFactory<MOD>) fac.coFac; 1045 BigInteger p = mcfac.getIntegerModul(); 1046 BigInteger modul = p; 1047 GenPolynomialRing<MOD> mfac; // = new GenPolynomialRing<MOD>(mcfac, fac); 1048 for (int i = 1; i < k; i++) { 1049 // e = 1 - s a - t b in Z[x] 1050 GenPolynomial<BigInteger> e = one.subtract(Si.multiply(Ai)).subtract(Ti.multiply(Bi)); 1051 //System.out.println("\ne = " + e); 1052 if (e.isZERO()) { 1053 logger.info("leaving on zero e in liftExtendedEuclidean"); 1054 break; 1055 } 1056 e = e.divide(modul); 1057 // move to Z_p[x] and compute next approximation 1058 GenPolynomial<MOD> c = PolyUtil.<MOD> fromIntegerCoefficients(fac, e); 1059 //System.out.println("c = " + c + ": " + c.ring.coFac); 1060 GenPolynomial<MOD> s = S.multiply(c); 1061 GenPolynomial<MOD> t = T.multiply(c); 1062 //System.out.println("s = " + s + ": " + s.ring.coFac); 1063 //System.out.println("t = " + t + ": " + t.ring.coFac); 1064 1065 GenPolynomial<MOD>[] QR = s.quotientRemainder(B); // watch for ordering 1066 GenPolynomial<MOD> q = QR[0]; 1067 s = QR[1]; 1068 t = t.sum(q.multiply(A)); 1069 //System.out.println("s = " + s + ": " + s.ring.coFac); 1070 //System.out.println("t = " + t + ": " + t.ring.coFac); 1071 1072 GenPolynomial<BigInteger> si = PolyUtil.integerFromModularCoefficients(ifac, s); 1073 GenPolynomial<BigInteger> ti = PolyUtil.integerFromModularCoefficients(ifac, t); 1074 //System.out.println("si = " + si); 1075 //System.out.println("ti = " + si); 1076 // add approximation to solution 1077 Si = Si.sum(si.multiply(modul)); 1078 Ti = Ti.sum(ti.multiply(modul)); 1079 //System.out.println("Si = " + Si); 1080 //System.out.println("Ti = " + Ti); 1081 modul = modul.multiply(p); 1082 //System.out.println("modul = " + modul + ", " + p + "^" + k + ", p^k = " + Power.power(new BigInteger(),p,k)); 1083 } 1084 //System.out.println("Si = " + Si + ", Ti = " + Ti); 1085 // setup ring mod p^i 1086 if (ModLongRing.MAX_LONG.compareTo(modul.getVal()) > 0) { 1087 mcfac = (ModularRingFactory) new ModLongRing(modul.getVal()); 1088 } else { 1089 mcfac = (ModularRingFactory) new ModIntegerRing(modul.getVal()); 1090 } 1091 //System.out.println("mcfac = " + mcfac); 1092 mfac = new GenPolynomialRing<MOD>(mcfac, fac); 1093 S = PolyUtil.<MOD> fromIntegerCoefficients(mfac, Si); 1094 T = PolyUtil.<MOD> fromIntegerCoefficients(mfac, Ti); 1095 //System.out.println("S = " + S + ": " + S.ring.coFac); 1096 //System.out.println("T = " + T + ": " + T.ring.coFac); 1097 if ( debug ) { 1098 List<GenPolynomial<MOD>> AP = new ArrayList<GenPolynomial<MOD>>(); 1099 AP.add(B); AP.add(A); 1100 List<GenPolynomial<MOD>> SP = new ArrayList<GenPolynomial<MOD>>(); 1101 SP.add(S); SP.add(T); 1102 if ( !HenselUtil.<MOD>isExtendedEuclideanLift(AP,SP) ) { 1103 System.out.println("isExtendedEuclideanLift: false"); 1104 } 1105 } 1106 GenPolynomial<MOD>[] rel = (GenPolynomial<MOD>[]) new GenPolynomial[2]; 1107 rel[0] = S; 1108 rel[1] = T; 1109 return rel; 1110 } 1111 1112 1113 /** 1114 * Constructing and lifting algorithm for extended Euclidean relation. Let p 1115 * = A_i.ring.coFac.modul() and assume gcd(A_i,A_j) == 1 mod p, i != j. 1116 * @param A list of modular GenPolynomials 1117 * @param k desired approximation exponent p^k. 1118 * @return [s_0,...,s_n-1] with sum_i s_i * B_i = 1 mod p^k, with B_i = 1119 * prod_{i!=j} A_j. 1120 */ 1121 public static <MOD extends GcdRingElem<MOD> & Modular> List<GenPolynomial<MOD>> liftExtendedEuclidean( 1122 List<GenPolynomial<MOD>> A, long k) throws NoLiftingException { 1123 if (A == null || A.size() == 0) { 1124 throw new IllegalArgumentException("A must be non null and non empty"); 1125 } 1126 GenPolynomialRing<MOD> fac = A.get(0).ring; 1127 if (fac.nvar != 1) { // todo assert 1128 throw new IllegalArgumentException("polynomial ring not univariate"); 1129 } 1130 GenPolynomial<MOD> zero = fac.getZERO(); 1131 int r = A.size(); 1132 List<GenPolynomial<MOD>> Q = new ArrayList<GenPolynomial<MOD>>(r); 1133 for (int i = 0; i < r; i++) { 1134 Q.add(zero); 1135 } 1136 //System.out.println("A = " + A); 1137 Q.set(r - 2, A.get(r - 1)); 1138 for (int j = r - 3; j >= 0; j--) { 1139 GenPolynomial<MOD> q = A.get(j + 1).multiply(Q.get(j + 1)); 1140 Q.set(j, q); 1141 } 1142 //System.out.println("Q = " + Q); 1143 List<GenPolynomial<MOD>> B = new ArrayList<GenPolynomial<MOD>>(r + 1); 1144 List<GenPolynomial<MOD>> lift = new ArrayList<GenPolynomial<MOD>>(r); 1145 for (int i = 0; i < r; i++) { 1146 B.add(zero); 1147 lift.add(zero); 1148 } 1149 GenPolynomial<MOD> one = fac.getONE(); 1150 GenPolynomialRing<BigInteger> ifac = new GenPolynomialRing<BigInteger>(new BigInteger(), fac); 1151 B.add(0, one); 1152 //System.out.println("B(0) = " + B.get(0)); 1153 GenPolynomial<MOD> b = one; 1154 for (int j = 0; j < r - 1; j++) { 1155 //System.out.println("Q("+(j)+") = " + Q.get(j)); 1156 //System.out.println("A("+(j)+") = " + A.get(j)); 1157 //System.out.println("B("+(j)+") = " + B.get(j)); 1158 List<GenPolynomial<MOD>> S = liftDiophant(Q.get(j), A.get(j), B.get(j), k); 1159 //System.out.println("\nSb = " + S); 1160 b = S.get(0); 1161 GenPolynomial<MOD> bb = PolyUtil.<MOD> fromIntegerCoefficients(fac, PolyUtil 1162 .integerFromModularCoefficients(ifac, b)); 1163 B.set(j + 1, bb); 1164 lift.set(j, S.get(1)); 1165 //System.out.println("B("+(j+1)+") = " + B.get(j+1)); 1166 if (debug) { 1167 logger.info("lift(" + j + ") = " + lift.get(j)); 1168 } 1169 } 1170 //System.out.println("liftb = " + lift); 1171 lift.set(r - 1, b); 1172 if (debug) { 1173 logger.info("lift(" + (r - 1) + ") = " + b); 1174 } 1175 //System.out.println("B("+(r-1)+") = " + B.get(r-1) + " : " + B.get(r-1).ring.coFac + ", b = " + b + " : " + b.ring.coFac); 1176 //System.out.println("B = " + B); 1177 //System.out.println("liftb = " + lift); 1178 return lift; 1179 } 1180 1181 1182 /** 1183 * Modular diophantine equation solution and lifting algorithm. Let p = 1184 * A_i.ring.coFac.modul() and assume gcd(A,B) == 1 mod p. 1185 * @param A modular GenPolynomial, mod p^k 1186 * @param B modular GenPolynomial, mod p^k 1187 * @param C modular GenPolynomial, mod p^k 1188 * @param k desired approximation exponent p^k. 1189 * @return [s, t] with s A' + t B' = C mod p^k, with A' = B, B' = A. 1190 */ 1191 public static <MOD extends GcdRingElem<MOD> & Modular> List<GenPolynomial<MOD>> liftDiophant( 1192 GenPolynomial<MOD> A, GenPolynomial<MOD> B, GenPolynomial<MOD> C, long k) 1193 throws NoLiftingException { 1194 if (A == null || A.isZERO() || B == null || B.isZERO()) { 1195 throw new IllegalArgumentException("A and B must be nonzero, A = " + A + ", B = " + B + ", C = " + C); 1196 } 1197 List<GenPolynomial<MOD>> sol = new ArrayList<GenPolynomial<MOD>>(); 1198 GenPolynomialRing<MOD> fac = C.ring; 1199 if (fac.nvar != 1) { // todo assert 1200 throw new IllegalArgumentException("polynomial ring not univariate"); 1201 } 1202 //System.out.println("C = " + C); 1203 GenPolynomial<MOD> zero = fac.getZERO(); 1204 for (int i = 0; i < 2; i++) { 1205 sol.add(zero); 1206 } 1207 GenPolynomialRing<BigInteger> ifac = new GenPolynomialRing<BigInteger>(new BigInteger(),fac); 1208 for (Monomial<MOD> m : C) { 1209 //System.out.println("monomial = " + m); 1210 long e = m.e.getVal(0); 1211 List<GenPolynomial<MOD>> S = liftDiophant(A, B, e, k); 1212 //System.out.println("Se = " + S); 1213 MOD a = m.c; 1214 //System.out.println("C.fac = " + fac.toScript()); 1215 a = fac.coFac.fromInteger(a.getSymmetricInteger().getVal()); 1216 int i = 0; 1217 for (GenPolynomial<MOD> d : S) { 1218 //System.out.println("d = " + d); 1219 d = PolyUtil.<MOD> fromIntegerCoefficients(fac, PolyUtil.integerFromModularCoefficients(ifac, d)); 1220 d = d.multiply(a); 1221 d = sol.get(i).sum(d); 1222 //System.out.println("d = " + d); 1223 sol.set(i++, d); 1224 } 1225 //System.out.println("sol = " + sol + ", for " + m); 1226 } 1227 if (true || debug) { 1228 //GenPolynomialRing<BigInteger> ifac = new GenPolynomialRing<BigInteger>(new BigInteger(), fac); 1229 A = PolyUtil.<MOD> fromIntegerCoefficients(fac, PolyUtil.integerFromModularCoefficients(ifac, A)); 1230 B = PolyUtil.<MOD> fromIntegerCoefficients(fac, PolyUtil.integerFromModularCoefficients(ifac, B)); 1231 C = PolyUtil.<MOD> fromIntegerCoefficients(fac, PolyUtil.integerFromModularCoefficients(ifac, C)); 1232 GenPolynomial<MOD> y = B.multiply(sol.get(0)).sum(A.multiply(sol.get(1))); 1233 if (!y.equals(C)) { 1234 System.out.println("A = " + A + ", B = " + B); 1235 System.out.println("s1 = " + sol.get(0) + ", s2 = " + sol.get(1)); 1236 System.out.println("Error: A*r1 + B*r2 = " + y + " : " + fac.coFac); 1237 } 1238 } 1239 return sol; 1240 } 1241 1242 1243 /** 1244 * Modular diophantine equation solution and lifting algorithm. Let p = 1245 * A_i.ring.coFac.modul() and assume gcd(a,b) == 1 mod p, for a, b in A. 1246 * @param A list of modular GenPolynomials, mod p^k 1247 * @param C modular GenPolynomial, mod p^k 1248 * @param k desired approximation exponent p^k. 1249 * @return [s_1,..., s_n] with sum_i s_i A_i' = C mod p^k, with Ai' = prod_{j!=i} A_j. 1250 */ 1251 public static <MOD extends GcdRingElem<MOD> & Modular> List<GenPolynomial<MOD>> liftDiophant( 1252 List<GenPolynomial<MOD>> A, GenPolynomial<MOD> C, long k) 1253 throws NoLiftingException { 1254 if ( false && A.size() <= 2 ) { 1255 return HenselUtil.<MOD> liftDiophant(A.get(0),A.get(1),C,k); 1256 } 1257 List<GenPolynomial<MOD>> sol = new ArrayList<GenPolynomial<MOD>>(); 1258 GenPolynomialRing<MOD> fac = C.ring; 1259 if (fac.nvar != 1) { // todo assert 1260 throw new IllegalArgumentException("polynomial ring not univariate"); 1261 } 1262 //System.out.println("C = " + C); 1263 GenPolynomial<MOD> zero = fac.getZERO(); 1264 for (int i = 0; i < A.size(); i++) { 1265 sol.add(zero); 1266 } 1267 GenPolynomialRing<BigInteger> ifac = new GenPolynomialRing<BigInteger>(new BigInteger(),fac); 1268 for (Monomial<MOD> m : C) { 1269 //System.out.println("monomial = " + m); 1270 long e = m.e.getVal(0); 1271 List<GenPolynomial<MOD>> S = liftDiophant(A, e, k); 1272 //System.out.println("Se = " + S); 1273 MOD a = m.c; 1274 //System.out.println("C.fac = " + fac.toScript()); 1275 a = fac.coFac.fromInteger(a.getSymmetricInteger().getVal()); 1276 int i = 0; 1277 for (GenPolynomial<MOD> d : S) { 1278 //System.out.println("d = " + d); 1279 d = PolyUtil.<MOD> fromIntegerCoefficients(fac, PolyUtil.integerFromModularCoefficients(ifac, d)); 1280 d = d.multiply(a); 1281 d = sol.get(i).sum(d); 1282 //System.out.println("d = " + d); 1283 sol.set(i++, d); 1284 } 1285 //System.out.println("sol = " + sol + ", for " + m); 1286 } 1287 /* 1288 if (true || debug) { 1289 //GenPolynomialRing<BigInteger> ifac = new GenPolynomialRing<BigInteger>(new BigInteger(), fac); 1290 A = PolyUtil.<MOD> fromIntegerCoefficients(fac, PolyUtil.integerFromModularCoefficients(ifac, A)); 1291 B = PolyUtil.<MOD> fromIntegerCoefficients(fac, PolyUtil.integerFromModularCoefficients(ifac, B)); 1292 C = PolyUtil.<MOD> fromIntegerCoefficients(fac, PolyUtil.integerFromModularCoefficients(ifac, C)); 1293 GenPolynomial<MOD> y = B.multiply(sol.get(0)).sum(A.multiply(sol.get(1))); 1294 if (!y.equals(C)) { 1295 System.out.println("A = " + A + ", B = " + B); 1296 System.out.println("s1 = " + sol.get(0) + ", s2 = " + sol.get(1)); 1297 System.out.println("Error: A*r1 + B*r2 = " + y + " : " + fac.coFac); 1298 } 1299 } 1300 */ 1301 return sol; 1302 } 1303 1304 1305 /** 1306 * Modular diophantine equation solution and lifting algorithm. Let p = 1307 * A_i.ring.coFac.modul() and assume gcd(A,B) == 1 mod p. 1308 * @param A modular GenPolynomial 1309 * @param B modular GenPolynomial 1310 * @param e exponent for x^e 1311 * @param k desired approximation exponent p^k. 1312 * @return [s, t] with s A' + t B' = x^e mod p^k, with A' = B, B' = A. 1313 */ 1314 public static <MOD extends GcdRingElem<MOD> & Modular> List<GenPolynomial<MOD>> liftDiophant( 1315 GenPolynomial<MOD> A, GenPolynomial<MOD> B, long e, long k) throws NoLiftingException { 1316 if (A == null || A.isZERO() || B == null || B.isZERO()) { 1317 throw new IllegalArgumentException("A and B must be nonzero, A = " + A + ", B = " + B); 1318 } 1319 List<GenPolynomial<MOD>> sol = new ArrayList<GenPolynomial<MOD>>(); 1320 GenPolynomialRing<MOD> fac = A.ring; 1321 if (fac.nvar != 1) { // todo assert 1322 throw new IllegalArgumentException("polynomial ring not univariate"); 1323 } 1324 // lift EE relation to p^k 1325 GenPolynomial<MOD>[] lee = liftExtendedEuclidean(B,A, k); 1326 GenPolynomial<MOD> s1 = lee[0]; 1327 GenPolynomial<MOD> s2 = lee[1]; 1328 if ( e == 0L ) { 1329 sol.add(s1); 1330 sol.add(s2); 1331 //System.out.println("sol@0 = " + sol); 1332 return sol; 1333 } 1334 fac = s1.ring; 1335 GenPolynomialRing<BigInteger> ifac = new GenPolynomialRing<BigInteger>(new BigInteger(), fac); 1336 A = PolyUtil.<MOD> fromIntegerCoefficients(fac, PolyUtil.integerFromModularCoefficients(ifac, A)); 1337 B = PolyUtil.<MOD> fromIntegerCoefficients(fac, PolyUtil.integerFromModularCoefficients(ifac, B)); 1338 1339 // this is the wrong sequence: 1340 // GenPolynomial<MOD> xe = fac.univariate(0,e); 1341 // GenPolynomial<MOD> q = s1.multiply(xe); 1342 // GenPolynomial<MOD>[] QR = q.quotientRemainder(B); 1343 // q = QR[0]; 1344 // GenPolynomial<MOD> r1 = QR[1]; 1345 // GenPolynomial<MOD> r2 = s2.multiply(xe).sum( q.multiply(A) ); 1346 1347 GenPolynomial<MOD> xe = fac.univariate(0, e); 1348 GenPolynomial<MOD> q = s1.multiply(xe); 1349 GenPolynomial<MOD>[] QR = q.quotientRemainder(A); 1350 q = QR[0]; 1351 //System.out.println("ee coeff qr = " + Arrays.toString(QR)); 1352 GenPolynomial<MOD> r1 = QR[1]; 1353 GenPolynomial<MOD> r2 = s2.multiply(xe).sum(q.multiply(B)); 1354 //System.out.println("r1 = " + r1 + ", r2 = " + r2); 1355 sol.add(r1); 1356 sol.add(r2); 1357 //System.out.println("sol@"+ e + " = " + sol); 1358 if (true || debug) { 1359 GenPolynomial<MOD> y = B.multiply(r1).sum(A.multiply(r2)); 1360 if (!y.equals(xe)) { 1361 System.out.println("A = " + A + ", B = " + B); 1362 System.out.println("r1 = " + r1 + ", r2 = " + r2); 1363 System.out.println("Error: A*r1 + B*r2 = " + y); 1364 } 1365 } 1366 return sol; 1367 } 1368 1369 1370 /** 1371 * Modular diophantine equation solution and lifting algorithm. Let p = 1372 * A_i.ring.coFac.modul() and assume gcd(a,b) == 1 mod p, for a, b in A. 1373 * @param A list of modular GenPolynomials 1374 * @param e exponent for x^e 1375 * @param k desired approximation exponent p^k. 1376 * @return [s_1,..., s_n] with sum_i s_i A_i' = x^e mod p^k, with Ai' = prod_{j!=i} A_j. 1377 */ 1378 public static <MOD extends GcdRingElem<MOD> & Modular> List<GenPolynomial<MOD>> liftDiophant( 1379 List<GenPolynomial<MOD>> A, long e, long k) throws NoLiftingException { 1380 if ( false && A.size() <= 2 ) { 1381 return HenselUtil.<MOD> liftDiophant(A.get(0),A.get(1),e,k); 1382 } 1383 List<GenPolynomial<MOD>> sol = new ArrayList<GenPolynomial<MOD>>(); 1384 GenPolynomialRing<MOD> fac = A.get(0).ring; 1385 if (fac.nvar != 1) { // todo assert 1386 throw new IllegalArgumentException("polynomial ring not univariate"); 1387 } 1388 // lift EE relation to p^k 1389 List<GenPolynomial<MOD>> lee = liftExtendedEuclidean(A, k); 1390 if ( e == 0L ) { 1391 //System.out.println("sol@0 = " + sol); 1392 return lee; 1393 } 1394 fac = lee.get(0).ring; 1395 GenPolynomialRing<BigInteger> ifac = new GenPolynomialRing<BigInteger>(new BigInteger(), fac); 1396 List<GenPolynomial<MOD>> S = new ArrayList<GenPolynomial<MOD>>(lee.size()); 1397 for ( GenPolynomial<MOD> a : lee ) { 1398 a = PolyUtil.<MOD> fromIntegerCoefficients(fac, PolyUtil.integerFromModularCoefficients(ifac, a)); 1399 S.add(a); 1400 } 1401 GenPolynomial<MOD> xe = fac.univariate(0, e); 1402 //List<GenPolynomial<MOD>> Sr = new ArrayList<GenPolynomial<MOD>>(lee.size()); 1403 int i = 0; 1404 for ( GenPolynomial<MOD> s : S ) { 1405 GenPolynomial<MOD> q = s.multiply(xe); 1406 GenPolynomial<MOD> r = q.remainder(A.get(i++)); 1407 //System.out.println("r = " + r); 1408 sol.add(r); 1409 } 1410 //System.out.println("sol@"+ e + " = " + sol); 1411 /* 1412 if (true || debug) { 1413 GenPolynomial<MOD> y = B.multiply(r1).sum(A.multiply(r2)); 1414 if (!y.equals(xe)) { 1415 System.out.println("A = " + A + ", B = " + B); 1416 System.out.println("r1 = " + r1 + ", r2 = " + r2); 1417 System.out.println("Error: A*r1 + B*r2 = " + y); 1418 } 1419 } 1420 */ 1421 return sol; 1422 } 1423 1424 1425 /** 1426 * Modular Diophant relation lifting test. 1427 * @param A modular GenPolynomial 1428 * @param B modular GenPolynomial 1429 * @param C modular GenPolynomial 1430 * @param S1 modular GenPolynomial 1431 * @param S2 modular GenPolynomial 1432 * @return true if A*S1 + B*S2 = C, else false. 1433 */ 1434 public static <MOD extends GcdRingElem<MOD> & Modular> boolean isDiophantLift( 1435 GenPolynomial<MOD> A, GenPolynomial<MOD> B, GenPolynomial<MOD> S1, GenPolynomial<MOD> S2, GenPolynomial<MOD> C) { 1436 GenPolynomialRing<MOD> fac = C.ring; 1437 GenPolynomialRing<BigInteger> ifac = new GenPolynomialRing<BigInteger>(new BigInteger(),fac); 1438 GenPolynomial<MOD> a = PolyUtil.<MOD> fromIntegerCoefficients(fac,PolyUtil.integerFromModularCoefficients(ifac,A)); 1439 GenPolynomial<MOD> b = PolyUtil.<MOD> fromIntegerCoefficients(fac,PolyUtil.integerFromModularCoefficients(ifac,B)); 1440 GenPolynomial<MOD> s1 = PolyUtil.<MOD> fromIntegerCoefficients(fac,PolyUtil.integerFromModularCoefficients(ifac,S1)); 1441 GenPolynomial<MOD> s2 = PolyUtil.<MOD> fromIntegerCoefficients(fac,PolyUtil.integerFromModularCoefficients(ifac,S2)); 1442 GenPolynomial<MOD> t = a.multiply(s1).sum( b.multiply(s2) ); 1443 if ( t.equals(C) ) { 1444 return true; 1445 } 1446 System.out.println("a = " + a); 1447 System.out.println("b = " + b); 1448 System.out.println("s1 = " + s1); 1449 System.out.println("s2 = " + s2); 1450 System.out.println("t = " + t); 1451 System.out.println("C = " + C); 1452 return false; 1453 } 1454 1455 1456 /** 1457 * Modular extended Euclidean relation lifting test. 1458 * @param A list of GenPolynomials 1459 * @param S = [s_0,...,s_{n-1}] list of GenPolynomial 1460 * @return true if prod_{0,...,n-1} s_i * B_i = 1 mod p^e, with B_i = 1461 * prod_{i!=j} A_j, else false. 1462 */ 1463 public static <MOD extends GcdRingElem<MOD> & Modular> boolean isExtendedEuclideanLift( 1464 List<GenPolynomial<MOD>> A, List<GenPolynomial<MOD>> S) { 1465 GenPolynomialRing<MOD> fac = A.get(0).ring; 1466 GenPolynomial<MOD> C = fac.getONE(); 1467 return isDiophantLift(A,S,C); 1468 } 1469 1470 1471 /** 1472 * Modular Diophant relation lifting test. 1473 * @param A list of GenPolynomials 1474 * @param S = [s_0,...,s_{n-1}] list of GenPolynomials 1475 * @param C = GenPolynomial 1476 * @return true if prod_{0,...,n-1} s_i * B_i = C mod p^k, with B_i = 1477 * prod_{i!=j} A_j, else false. 1478 */ 1479 public static <MOD extends GcdRingElem<MOD> & Modular> boolean isDiophantLift( 1480 List<GenPolynomial<MOD>> A, List<GenPolynomial<MOD>> S, GenPolynomial<MOD> C) { 1481 GenPolynomialRing<MOD> fac = A.get(0).ring; 1482 GenPolynomialRing<BigInteger> ifac = new GenPolynomialRing<BigInteger>(new BigInteger(),fac); 1483 List<GenPolynomial<MOD>> B = new ArrayList<GenPolynomial<MOD>>(A.size()); 1484 int i = 0; 1485 for (GenPolynomial<MOD> ai : A) { 1486 GenPolynomial<MOD> b = fac.getONE(); 1487 int j = 0; 1488 for (GenPolynomial<MOD> aj : A) { 1489 if (i != j /*!ai.equals(aj)*/) { 1490 b = b.multiply(aj); 1491 } 1492 j++; 1493 } 1494 //System.out.println("b = " + b); 1495 b = PolyUtil.<MOD> fromIntegerCoefficients(fac,PolyUtil.integerFromModularCoefficients(ifac,b)); 1496 B.add(b); 1497 i++; 1498 } 1499 //System.out.println("B = " + B); 1500 // check mod p^e 1501 GenPolynomial<MOD> t = fac.getZERO(); 1502 i = 0; 1503 for (GenPolynomial<MOD> a : B) { 1504 GenPolynomial<MOD> b = S.get(i++); 1505 b = PolyUtil.<MOD> fromIntegerCoefficients(fac,PolyUtil.integerFromModularCoefficients(ifac,b)); 1506 GenPolynomial<MOD> s = a.multiply(b); 1507 t = t.sum(s); 1508 } 1509 if (!t.equals(C)) { 1510 System.out.println("no diophant lift!"); 1511 System.out.println("A = " + A); 1512 System.out.println("B = " + B); 1513 System.out.println("S = " + S); 1514 System.out.println("C = " + C); 1515 System.out.println("t = " + t); 1516 return false; 1517 } 1518 return true; 1519 } 1520 1521 1522 /** 1523 * Modular Hensel lifting algorithm on coefficients. Let p = 1524 * f_i.ring.coFac.modul() and assume C == prod_{0,...,n-1} f_i mod p with 1525 * gcd(f_i,f_j) == 1 mod p for i != j 1526 * @param C monic integer polynomial 1527 * @param F = [f_0,...,f_{n-1}] list of monic modular polynomials. 1528 * @param k approximation exponent. 1529 * @return [g_0,...,g_{n-1}] with C = prod_{0,...,n-1} g_i mod p^k. 1530 */ 1531 public static <MOD extends GcdRingElem<MOD> & Modular> 1532 List<GenPolynomial<MOD>> liftHenselMonic(GenPolynomial<BigInteger> C, List<GenPolynomial<MOD>> F, long k) 1533 throws NoLiftingException { 1534 if (C == null || C.isZERO() || F == null || F.size() == 0) { 1535 throw new IllegalArgumentException("C must be nonzero and F must be nonempty"); 1536 } 1537 GenPolynomialRing<BigInteger> fac = C.ring; 1538 if (fac.nvar != 1) { // todo assert 1539 throw new IllegalArgumentException("polynomial ring not univariate"); 1540 } 1541 List<GenPolynomial<MOD>> lift = new ArrayList<GenPolynomial<MOD>>(F.size()); 1542 GenPolynomialRing<MOD> pfac = F.get(0).ring; 1543 RingFactory<MOD> pcfac = pfac.coFac; 1544 ModularRingFactory<MOD> PF = (ModularRingFactory<MOD>) pcfac; 1545 BigInteger P = PF.getIntegerModul(); 1546 int n = F.size(); 1547 if (n == 1) { // lift F_0, this case will probably never be used 1548 GenPolynomial<MOD> f = F.get(0); 1549 ModularRingFactory<MOD> mcfac; 1550 if (ModLongRing.MAX_LONG.compareTo(P.getVal()) > 0) { 1551 mcfac = (ModularRingFactory) new ModLongRing(P.getVal()); 1552 } else { 1553 mcfac = (ModularRingFactory) new ModIntegerRing(P.getVal()); 1554 } 1555 GenPolynomialRing<MOD> mfac = new GenPolynomialRing<MOD>(mcfac, fac); 1556 f = PolyUtil.fromIntegerCoefficients(mfac, PolyUtil.integerFromModularCoefficients(fac, f)); 1557 lift.add(f); 1558 return lift; 1559 } 1560 // if (n == 2) { // only one step 1561 // HenselApprox<MOD> ab = HenselUtil.<MOD> liftHenselQuadratic(C, M, F.get(0), F.get(1)); 1562 // lift.add(ab.Am); 1563 // lift.add(ab.Bm); 1564 // return lift; 1565 // } 1566 1567 // setup integer polynomial ring 1568 GenPolynomialRing<BigInteger> ifac = new GenPolynomialRing<BigInteger>(new BigInteger(), fac); 1569 List<GenPolynomial<BigInteger>> Fi = PolyUtil.integerFromModularCoefficients(ifac, F); 1570 //System.out.println("Fi = " + Fi); 1571 1572 List<GenPolynomial<MOD>> S = liftExtendedEuclidean(F, k + 1); // lift works for any k, TODO: use this 1573 //System.out.println("Sext = " + S); 1574 if (debug) { 1575 logger.info("EE lift = " + S); 1576 // adjust coefficients 1577 List<GenPolynomial<MOD>> Sx = PolyUtil.fromIntegerCoefficients(pfac, PolyUtil 1578 .integerFromModularCoefficients(ifac, S)); 1579 try { 1580 boolean il = HenselUtil.<MOD> isExtendedEuclideanLift(F, Sx); 1581 //System.out.println("islift = " + il); 1582 } catch (RuntimeException e) { 1583 e.printStackTrace(); 1584 } 1585 } 1586 List<GenPolynomial<BigInteger>> Si = PolyUtil.integerFromModularCoefficients(ifac, S); 1587 //System.out.println("Si = " + Si); 1588 //System.out.println("C = " + C); 1589 1590 // approximate mod p^i 1591 ModularRingFactory<MOD> mcfac = PF; 1592 BigInteger p = mcfac.getIntegerModul(); 1593 BigInteger modul = p; 1594 GenPolynomialRing<MOD> mfac = new GenPolynomialRing<MOD>(mcfac, fac); 1595 List<GenPolynomial<MOD>> Sp = PolyUtil.fromIntegerCoefficients(mfac, Si); 1596 //System.out.println("Sp = " + Sp); 1597 for (int i = 1; i < k; i++) { 1598 //System.out.println("i = " + i); 1599 GenPolynomial<BigInteger> e = fac.getONE(); 1600 for (GenPolynomial<BigInteger> fi : Fi) { 1601 e = e.multiply(fi); 1602 } 1603 e = C.subtract(e); 1604 //System.out.println("\ne = " + e); 1605 if (e.isZERO()) { 1606 logger.info("leaving on zero e"); 1607 break; 1608 } 1609 try { 1610 e = e.divide(modul); 1611 } catch (RuntimeException ex) { 1612 ex.printStackTrace(); 1613 throw ex; 1614 } 1615 //System.out.println("e = " + e); 1616 // move to in Z_p[x] 1617 GenPolynomial<MOD> c = PolyUtil.<MOD> fromIntegerCoefficients(mfac, e); 1618 //System.out.println("c = " + c + ": " + c.ring.coFac); 1619 1620 List<GenPolynomial<MOD>> s = new ArrayList<GenPolynomial<MOD>>(S.size()); 1621 int j = 0; 1622 for (GenPolynomial<MOD> f : Sp) { 1623 f = f.multiply(c); 1624 //System.out.println("f = " + f + " : " + f.ring.coFac); 1625 //System.out.println("F,i = " + F.get(j) + " : " + F.get(j).ring.coFac); 1626 f = f.remainder(F.get(j++)); 1627 //System.out.println("f = " + f + " : " + f.ring.coFac); 1628 s.add(f); 1629 } 1630 //System.out.println("s = " + s); 1631 List<GenPolynomial<BigInteger>> si = PolyUtil.integerFromModularCoefficients(ifac, s); 1632 //System.out.println("si = " + si); 1633 1634 List<GenPolynomial<BigInteger>> Fii = new ArrayList<GenPolynomial<BigInteger>>(F.size()); 1635 j = 0; 1636 for (GenPolynomial<BigInteger> f : Fi) { 1637 f = f.sum(si.get(j++).multiply(modul)); 1638 Fii.add(f); 1639 } 1640 //System.out.println("Fii = " + Fii); 1641 Fi = Fii; 1642 modul = modul.multiply(p); 1643 if (i >= k - 1) { 1644 logger.info("e != 0 for k = " + k); 1645 } 1646 } 1647 // setup ring mod p^k 1648 modul = Power.positivePower(p,k); 1649 if (ModLongRing.MAX_LONG.compareTo(modul.getVal()) > 0) { 1650 mcfac = (ModularRingFactory) new ModLongRing(modul.getVal()); 1651 } else { 1652 mcfac = (ModularRingFactory) new ModIntegerRing(modul.getVal()); 1653 } 1654 //System.out.println("mcfac = " + mcfac); 1655 mfac = new GenPolynomialRing<MOD>(mcfac, fac); 1656 lift = PolyUtil.<MOD> fromIntegerCoefficients(mfac, Fi); 1657 //System.out.println("lift = " + lift + ": " + lift.get(0).ring.coFac); 1658 return lift; 1659 } 1660 1661 1662 /** 1663 * Modular Hensel lifting algorithm on coefficients. Let p = 1664 * f_i.ring.coFac.modul() and assume C == prod_{0,...,n-1} f_i mod p with 1665 * gcd(f_i,f_j) == 1 mod p for i != j 1666 * @param C integer polynomial 1667 * @param F = [f_0,...,f_{n-1}] list of monic modular polynomials. 1668 * @param k approximation exponent. 1669 * @param g leading coefficient. 1670 * @return [g_0,...,g_{n-1}] with C = prod_{0,...,n-1} g_i mod p^k. 1671 */ 1672 public static <MOD extends GcdRingElem<MOD> & Modular> 1673 List<GenPolynomial<MOD>> liftHensel(GenPolynomial<BigInteger> C, List<GenPolynomial<MOD>> F, long k, BigInteger g) 1674 throws NoLiftingException { 1675 if (C == null || C.isZERO() || F == null || F.size() == 0) { 1676 throw new IllegalArgumentException("C must be nonzero and F must be nonempty"); 1677 } 1678 GenPolynomialRing<BigInteger> fac = C.ring; 1679 if (fac.nvar != 1) { // todo assert 1680 throw new IllegalArgumentException("polynomial ring not univariate"); 1681 } 1682 List<GenPolynomial<MOD>> lift = new ArrayList<GenPolynomial<MOD>>(F.size()); 1683 GenPolynomialRing<MOD> pfac = F.get(0).ring; 1684 RingFactory<MOD> pcfac = pfac.coFac; 1685 ModularRingFactory<MOD> PF = (ModularRingFactory<MOD>) pcfac; 1686 BigInteger P = PF.getIntegerModul(); 1687 int n = F.size(); 1688 if (n == 1) { // lift F_0, this case will probably never be used 1689 GenPolynomial<MOD> f = F.get(0); 1690 ModularRingFactory<MOD> mcfac; 1691 if (ModLongRing.MAX_LONG.compareTo(P.getVal()) > 0) { 1692 mcfac = (ModularRingFactory) new ModLongRing(P.getVal()); 1693 } else { 1694 mcfac = (ModularRingFactory) new ModIntegerRing(P.getVal()); 1695 } 1696 GenPolynomialRing<MOD> mfac = new GenPolynomialRing<MOD>(mcfac, fac); 1697 f = PolyUtil.fromIntegerCoefficients(mfac, PolyUtil.integerFromModularCoefficients(fac, f)); 1698 lift.add(f); 1699 return lift; 1700 } 1701 // if (n == 2) { // only one step 1702 // HenselApprox<MOD> ab = HenselUtil.<MOD> liftHenselQuadratic(C, M, F.get(0), F.get(1)); 1703 // lift.add(ab.Am); 1704 // lift.add(ab.Bm); 1705 // return lift; 1706 // } 1707 1708 // normalize C and F_i factors 1709 BigInteger cc = g; //C.leadingBaseCoefficient(); // == g ?? 1710 for ( int i = 1; i < F.size(); i++ ) { // #F-1 1711 C = C.multiply(cc); // sic 1712 } 1713 MOD cm = PF.fromInteger(cc.getVal()); 1714 List<GenPolynomial<MOD>> Fp = new ArrayList<GenPolynomial<MOD>>(F.size()); 1715 for ( GenPolynomial<MOD> fm : F ) { 1716 GenPolynomial<MOD> am = fm.monic(); 1717 am = am.multiply(cm); 1718 Fp.add(am); 1719 } 1720 F = Fp; 1721 1722 // setup integer polynomial ring 1723 GenPolynomialRing<BigInteger> ifac = new GenPolynomialRing<BigInteger>(new BigInteger(), fac); 1724 List<GenPolynomial<BigInteger>> Fi = PolyUtil.integerFromModularCoefficients(ifac, F); 1725 //System.out.println("Fi = " + Fi); 1726 1727 // inplace modify polynomials, replace leading coefficient 1728 for ( GenPolynomial<BigInteger> ai : Fi ) { 1729 ExpVector ea = ai.leadingExpVector(); 1730 ai.doPutToMap(ea, cc); 1731 } 1732 //System.out.println("Fi = " + Fi); 1733 1734 List<GenPolynomial<MOD>> S = liftExtendedEuclidean(F, k + 1); // lift works for any k, TODO: use this 1735 //System.out.println("Sext = " + S); 1736 if (debug) { 1737 logger.info("EE lift = " + S); 1738 // adjust coefficients 1739 List<GenPolynomial<MOD>> Sx = PolyUtil.fromIntegerCoefficients(pfac, PolyUtil 1740 .integerFromModularCoefficients(ifac, S)); 1741 try { 1742 boolean il = HenselUtil.<MOD> isExtendedEuclideanLift(F, Sx); 1743 //System.out.println("islift = " + il); 1744 } catch (RuntimeException e) { 1745 e.printStackTrace(); 1746 } 1747 } 1748 List<GenPolynomial<BigInteger>> Si = PolyUtil.integerFromModularCoefficients(ifac, S); 1749 //System.out.println("Si = " + Si); 1750 //System.out.println("C = " + C); 1751 1752 // approximate mod p^i 1753 ModularRingFactory<MOD> mcfac = PF; 1754 BigInteger p = mcfac.getIntegerModul(); 1755 BigInteger modul = p; 1756 GenPolynomialRing<MOD> mfac = new GenPolynomialRing<MOD>(mcfac, fac); 1757 List<GenPolynomial<MOD>> Sp = PolyUtil.fromIntegerCoefficients(mfac, Si); 1758 //System.out.println("Sp = " + Sp); 1759 for (int i = 1; i < k; i++) { 1760 //System.out.println("i = " + i); 1761 GenPolynomial<BigInteger> e = fac.getONE(); 1762 for (GenPolynomial<BigInteger> fi : Fi) { 1763 e = e.multiply(fi); 1764 } 1765 e = C.subtract(e); 1766 //System.out.println("\ne = " + e); 1767 if (e.isZERO()) { 1768 logger.info("leaving on zero e"); 1769 break; 1770 } 1771 try { 1772 e = e.divide(modul); 1773 } catch (RuntimeException ex) { 1774 ex.printStackTrace(); 1775 throw ex; 1776 } 1777 //System.out.println("e = " + e); 1778 // move to in Z_p[x] 1779 GenPolynomial<MOD> c = PolyUtil.<MOD> fromIntegerCoefficients(mfac, e); 1780 //System.out.println("c = " + c + ": " + c.ring.coFac); 1781 1782 List<GenPolynomial<MOD>> s = new ArrayList<GenPolynomial<MOD>>(S.size()); 1783 int j = 0; 1784 for (GenPolynomial<MOD> f : Sp) { 1785 f = f.multiply(c); 1786 //System.out.println("f = " + f + " : " + f.ring.coFac); 1787 //System.out.println("F,i = " + F.get(j) + " : " + F.get(j).ring.coFac); 1788 f = f.remainder(F.get(j++)); 1789 //System.out.println("f = " + f + " : " + f.ring.coFac); 1790 s.add(f); 1791 } 1792 //System.out.println("s = " + s); 1793 List<GenPolynomial<BigInteger>> si = PolyUtil.integerFromModularCoefficients(ifac, s); 1794 //System.out.println("si = " + si); 1795 1796 List<GenPolynomial<BigInteger>> Fii = new ArrayList<GenPolynomial<BigInteger>>(F.size()); 1797 j = 0; 1798 for (GenPolynomial<BigInteger> f : Fi) { 1799 f = f.sum(si.get(j++).multiply(modul)); 1800 Fii.add(f); 1801 } 1802 //System.out.println("Fii = " + Fii); 1803 Fi = Fii; 1804 modul = modul.multiply(p); 1805 if (i >= k - 1) { 1806 logger.info("e != 0 for k = " + k); 1807 } 1808 } 1809 //System.out.println("Fi = " + Fi); 1810 // remove normalization 1811 GreatestCommonDivisorAbstract<BigInteger> ufd = GCDFactory.getImplementation(cc); 1812 //BigInteger ai = ufd.baseContent(Fi.get(0)); 1813 //System.out.println("ai = " + ai + ", cc = " + cc); 1814 List<GenPolynomial<BigInteger>> Fii = new ArrayList<GenPolynomial<BigInteger>>(F.size()); 1815 //int j = 0; 1816 for ( GenPolynomial<BigInteger> bi : Fi ) { 1817 GenPolynomial<BigInteger> ci = null; 1818 //if ( j++ == 0 ) { 1819 // ci = bi.divide(ai); 1820 //} else { 1821 // BigInteger i = cc.divide(ai); 1822 // ci = bi.divide(i); 1823 //} 1824 ci = ufd.basePrimitivePart(bi); // ?? 1825 //System.out.println("bi = " + bi + ", ci = " + ci); 1826 Fii.add(ci); 1827 } 1828 Fi = Fii; 1829 1830 // setup ring mod p^k 1831 modul = Power.positivePower(p,k); 1832 if (ModLongRing.MAX_LONG.compareTo(modul.getVal()) > 0) { 1833 mcfac = (ModularRingFactory) new ModLongRing(modul.getVal()); 1834 } else { 1835 mcfac = (ModularRingFactory) new ModIntegerRing(modul.getVal()); 1836 } 1837 //System.out.println("mcfac = " + mcfac); 1838 mfac = new GenPolynomialRing<MOD>(mcfac, fac); 1839 lift = PolyUtil.<MOD> fromIntegerCoefficients(mfac, Fi); 1840 //System.out.println("lift = " + lift + ": " + lift.get(0).ring.coFac); 1841 return lift; 1842 } 1843 1844}