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