001 /*
002 * $Id: HenselUtil.java 3708 2011-08-01 11:49:59Z kredel $
003 */
004
005 package edu.jas.ufd;
006
007
008 import java.util.ArrayList;
009 import java.util.Arrays;
010 import java.util.List;
011
012 import org.apache.log4j.Logger;
013
014 import edu.jas.arith.BigInteger;
015 import edu.jas.arith.ModInteger;
016 import edu.jas.arith.ModIntegerRing;
017 import edu.jas.arith.ModLongRing;
018 import edu.jas.arith.Modular;
019 import edu.jas.arith.ModularRingFactory;
020 import edu.jas.poly.ExpVector;
021 import edu.jas.poly.GenPolynomial;
022 import edu.jas.poly.GenPolynomialRing;
023 import edu.jas.poly.Monomial;
024 import edu.jas.poly.PolyUtil;
025 import edu.jas.structure.GcdRingElem;
026 import edu.jas.structure.RingFactory;
027 import edu.jas.structure.Power;
028
029
030 /**
031 * Hensel utilities for ufd.
032 * @author Heinz Kredel
033 */
034
035 public 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 GenPolynomial<BigInteger> one = ifac.getONE();
1570 List<GenPolynomial<BigInteger>> Fi = PolyUtil.integerFromModularCoefficients(ifac, F);
1571 //System.out.println("Fi = " + Fi);
1572
1573 List<GenPolynomial<MOD>> S = liftExtendedEuclidean(F, k + 1); // lift works for any k, TODO: use this
1574 //System.out.println("Sext = " + S);
1575 if (debug) {
1576 logger.info("EE lift = " + S);
1577 // adjust coefficients
1578 List<GenPolynomial<MOD>> Sx = PolyUtil.fromIntegerCoefficients(pfac, PolyUtil
1579 .integerFromModularCoefficients(ifac, S));
1580 try {
1581 boolean il = HenselUtil.<MOD> isExtendedEuclideanLift(F, Sx);
1582 //System.out.println("islift = " + il);
1583 } catch (RuntimeException e) {
1584 e.printStackTrace();
1585 }
1586 }
1587 List<GenPolynomial<BigInteger>> Si = PolyUtil.integerFromModularCoefficients(ifac, S);
1588 //System.out.println("Si = " + Si);
1589 //System.out.println("C = " + C);
1590
1591 // approximate mod p^i
1592 ModularRingFactory<MOD> mcfac = PF;
1593 BigInteger p = mcfac.getIntegerModul();
1594 BigInteger modul = p;
1595 GenPolynomialRing<MOD> mfac = new GenPolynomialRing<MOD>(mcfac, fac);
1596 List<GenPolynomial<MOD>> Sp = PolyUtil.fromIntegerCoefficients(mfac, Si);
1597 //System.out.println("Sp = " + Sp);
1598 for (int i = 1; i < k; i++) {
1599 //System.out.println("i = " + i);
1600 GenPolynomial<BigInteger> e = fac.getONE();
1601 for (GenPolynomial<BigInteger> fi : Fi) {
1602 e = e.multiply(fi);
1603 }
1604 e = C.subtract(e);
1605 //System.out.println("\ne = " + e);
1606 if (e.isZERO()) {
1607 logger.info("leaving on zero e");
1608 break;
1609 }
1610 try {
1611 e = e.divide(modul);
1612 } catch (RuntimeException ex) {
1613 ex.printStackTrace();
1614 throw ex;
1615 }
1616 //System.out.println("e = " + e);
1617 // move to in Z_p[x]
1618 GenPolynomial<MOD> c = PolyUtil.<MOD> fromIntegerCoefficients(mfac, e);
1619 //System.out.println("c = " + c + ": " + c.ring.coFac);
1620
1621 List<GenPolynomial<MOD>> s = new ArrayList<GenPolynomial<MOD>>(S.size());
1622 int j = 0;
1623 for (GenPolynomial<MOD> f : Sp) {
1624 f = f.multiply(c);
1625 //System.out.println("f = " + f + " : " + f.ring.coFac);
1626 //System.out.println("F,i = " + F.get(j) + " : " + F.get(j).ring.coFac);
1627 f = f.remainder(F.get(j++));
1628 //System.out.println("f = " + f + " : " + f.ring.coFac);
1629 s.add(f);
1630 }
1631 //System.out.println("s = " + s);
1632 List<GenPolynomial<BigInteger>> si = PolyUtil.integerFromModularCoefficients(ifac, s);
1633 //System.out.println("si = " + si);
1634
1635 List<GenPolynomial<BigInteger>> Fii = new ArrayList<GenPolynomial<BigInteger>>(F.size());
1636 j = 0;
1637 for (GenPolynomial<BigInteger> f : Fi) {
1638 f = f.sum(si.get(j++).multiply(modul));
1639 Fii.add(f);
1640 }
1641 //System.out.println("Fii = " + Fii);
1642 Fi = Fii;
1643 modul = modul.multiply(p);
1644 if (i >= k - 1) {
1645 logger.info("e != 0 for k = " + k);
1646 }
1647 }
1648 // setup ring mod p^k
1649 modul = Power.positivePower(p,k);
1650 if (ModLongRing.MAX_LONG.compareTo(modul.getVal()) > 0) {
1651 mcfac = (ModularRingFactory) new ModLongRing(modul.getVal());
1652 } else {
1653 mcfac = (ModularRingFactory) new ModIntegerRing(modul.getVal());
1654 }
1655 //System.out.println("mcfac = " + mcfac);
1656 mfac = new GenPolynomialRing<MOD>(mcfac, fac);
1657 lift = PolyUtil.<MOD> fromIntegerCoefficients(mfac, Fi);
1658 //System.out.println("lift = " + lift + ": " + lift.get(0).ring.coFac);
1659 return lift;
1660 }
1661
1662
1663 /**
1664 * Modular Hensel lifting algorithm on coefficients. Let p =
1665 * f_i.ring.coFac.modul() and assume C == prod_{0,...,n-1} f_i mod p with
1666 * gcd(f_i,f_j) == 1 mod p for i != j
1667 * @param C integer polynomial
1668 * @param F = [f_0,...,f_{n-1}] list of monic modular polynomials.
1669 * @param k approximation exponent.
1670 * @param g leading coefficient.
1671 * @return [g_0,...,g_{n-1}] with C = prod_{0,...,n-1} g_i mod p^k.
1672 */
1673 public static <MOD extends GcdRingElem<MOD> & Modular>
1674 List<GenPolynomial<MOD>> liftHensel(GenPolynomial<BigInteger> C, List<GenPolynomial<MOD>> F, long k, BigInteger g)
1675 throws NoLiftingException {
1676 if (C == null || C.isZERO() || F == null || F.size() == 0) {
1677 throw new IllegalArgumentException("C must be nonzero and F must be nonempty");
1678 }
1679 GenPolynomialRing<BigInteger> fac = C.ring;
1680 if (fac.nvar != 1) { // todo assert
1681 throw new IllegalArgumentException("polynomial ring not univariate");
1682 }
1683 List<GenPolynomial<MOD>> lift = new ArrayList<GenPolynomial<MOD>>(F.size());
1684 GenPolynomialRing<MOD> pfac = F.get(0).ring;
1685 RingFactory<MOD> pcfac = pfac.coFac;
1686 ModularRingFactory<MOD> PF = (ModularRingFactory<MOD>) pcfac;
1687 BigInteger P = PF.getIntegerModul();
1688 int n = F.size();
1689 if (n == 1) { // lift F_0, this case will probably never be used
1690 GenPolynomial<MOD> f = F.get(0);
1691 ModularRingFactory<MOD> mcfac;
1692 if (ModLongRing.MAX_LONG.compareTo(P.getVal()) > 0) {
1693 mcfac = (ModularRingFactory) new ModLongRing(P.getVal());
1694 } else {
1695 mcfac = (ModularRingFactory) new ModIntegerRing(P.getVal());
1696 }
1697 GenPolynomialRing<MOD> mfac = new GenPolynomialRing<MOD>(mcfac, fac);
1698 f = PolyUtil.fromIntegerCoefficients(mfac, PolyUtil.integerFromModularCoefficients(fac, f));
1699 lift.add(f);
1700 return lift;
1701 }
1702 // if (n == 2) { // only one step
1703 // HenselApprox<MOD> ab = HenselUtil.<MOD> liftHenselQuadratic(C, M, F.get(0), F.get(1));
1704 // lift.add(ab.Am);
1705 // lift.add(ab.Bm);
1706 // return lift;
1707 // }
1708
1709 // normalize C and F_i factors
1710 BigInteger cc = g; //C.leadingBaseCoefficient(); // == g ??
1711 for ( int i = 1; i < F.size(); i++ ) { // #F-1
1712 C = C.multiply(cc); // sic
1713 }
1714 MOD cm = PF.fromInteger(cc.getVal());
1715 List<GenPolynomial<MOD>> Fp = new ArrayList<GenPolynomial<MOD>>(F.size());
1716 for ( GenPolynomial<MOD> fm : F ) {
1717 GenPolynomial<MOD> am = fm.monic();
1718 am = am.multiply(cm);
1719 Fp.add(am);
1720 }
1721 F = Fp;
1722
1723 // setup integer polynomial ring
1724 GenPolynomialRing<BigInteger> ifac = new GenPolynomialRing<BigInteger>(new BigInteger(), fac);
1725 GenPolynomial<BigInteger> one = ifac.getONE();
1726 List<GenPolynomial<BigInteger>> Fi = PolyUtil.integerFromModularCoefficients(ifac, F);
1727 //System.out.println("Fi = " + Fi);
1728
1729 // inplace modify polynomials, replace leading coefficient
1730 for ( GenPolynomial<BigInteger> ai : Fi ) {
1731 ExpVector ea = ai.leadingExpVector();
1732 ai.doPutToMap(ea, cc);
1733 }
1734 //System.out.println("Fi = " + Fi);
1735
1736 List<GenPolynomial<MOD>> S = liftExtendedEuclidean(F, k + 1); // lift works for any k, TODO: use this
1737 //System.out.println("Sext = " + S);
1738 if (debug) {
1739 logger.info("EE lift = " + S);
1740 // adjust coefficients
1741 List<GenPolynomial<MOD>> Sx = PolyUtil.fromIntegerCoefficients(pfac, PolyUtil
1742 .integerFromModularCoefficients(ifac, S));
1743 try {
1744 boolean il = HenselUtil.<MOD> isExtendedEuclideanLift(F, Sx);
1745 //System.out.println("islift = " + il);
1746 } catch (RuntimeException e) {
1747 e.printStackTrace();
1748 }
1749 }
1750 List<GenPolynomial<BigInteger>> Si = PolyUtil.integerFromModularCoefficients(ifac, S);
1751 //System.out.println("Si = " + Si);
1752 //System.out.println("C = " + C);
1753
1754 // approximate mod p^i
1755 ModularRingFactory<MOD> mcfac = PF;
1756 BigInteger p = mcfac.getIntegerModul();
1757 BigInteger modul = p;
1758 GenPolynomialRing<MOD> mfac = new GenPolynomialRing<MOD>(mcfac, fac);
1759 List<GenPolynomial<MOD>> Sp = PolyUtil.fromIntegerCoefficients(mfac, Si);
1760 //System.out.println("Sp = " + Sp);
1761 for (int i = 1; i < k; i++) {
1762 //System.out.println("i = " + i);
1763 GenPolynomial<BigInteger> e = fac.getONE();
1764 for (GenPolynomial<BigInteger> fi : Fi) {
1765 e = e.multiply(fi);
1766 }
1767 e = C.subtract(e);
1768 //System.out.println("\ne = " + e);
1769 if (e.isZERO()) {
1770 logger.info("leaving on zero e");
1771 break;
1772 }
1773 try {
1774 e = e.divide(modul);
1775 } catch (RuntimeException ex) {
1776 ex.printStackTrace();
1777 throw ex;
1778 }
1779 //System.out.println("e = " + e);
1780 // move to in Z_p[x]
1781 GenPolynomial<MOD> c = PolyUtil.<MOD> fromIntegerCoefficients(mfac, e);
1782 //System.out.println("c = " + c + ": " + c.ring.coFac);
1783
1784 List<GenPolynomial<MOD>> s = new ArrayList<GenPolynomial<MOD>>(S.size());
1785 int j = 0;
1786 for (GenPolynomial<MOD> f : Sp) {
1787 f = f.multiply(c);
1788 //System.out.println("f = " + f + " : " + f.ring.coFac);
1789 //System.out.println("F,i = " + F.get(j) + " : " + F.get(j).ring.coFac);
1790 f = f.remainder(F.get(j++));
1791 //System.out.println("f = " + f + " : " + f.ring.coFac);
1792 s.add(f);
1793 }
1794 //System.out.println("s = " + s);
1795 List<GenPolynomial<BigInteger>> si = PolyUtil.integerFromModularCoefficients(ifac, s);
1796 //System.out.println("si = " + si);
1797
1798 List<GenPolynomial<BigInteger>> Fii = new ArrayList<GenPolynomial<BigInteger>>(F.size());
1799 j = 0;
1800 for (GenPolynomial<BigInteger> f : Fi) {
1801 f = f.sum(si.get(j++).multiply(modul));
1802 Fii.add(f);
1803 }
1804 //System.out.println("Fii = " + Fii);
1805 Fi = Fii;
1806 modul = modul.multiply(p);
1807 if (i >= k - 1) {
1808 logger.info("e != 0 for k = " + k);
1809 }
1810 }
1811 //System.out.println("Fi = " + Fi);
1812 // remove normalization
1813 GreatestCommonDivisorAbstract<BigInteger> ufd = GCDFactory.getImplementation(cc);
1814 //BigInteger ai = ufd.baseContent(Fi.get(0));
1815 //System.out.println("ai = " + ai + ", cc = " + cc);
1816 List<GenPolynomial<BigInteger>> Fii = new ArrayList<GenPolynomial<BigInteger>>(F.size());
1817 //int j = 0;
1818 for ( GenPolynomial<BigInteger> bi : Fi ) {
1819 GenPolynomial<BigInteger> ci = null;
1820 //if ( j++ == 0 ) {
1821 // ci = bi.divide(ai);
1822 //} else {
1823 // BigInteger i = cc.divide(ai);
1824 // ci = bi.divide(i);
1825 //}
1826 ci = ufd.basePrimitivePart(bi); // ??
1827 //System.out.println("bi = " + bi + ", ci = " + ci);
1828 Fii.add(ci);
1829 }
1830 Fi = Fii;
1831
1832 // setup ring mod p^k
1833 modul = Power.positivePower(p,k);
1834 if (ModLongRing.MAX_LONG.compareTo(modul.getVal()) > 0) {
1835 mcfac = (ModularRingFactory) new ModLongRing(modul.getVal());
1836 } else {
1837 mcfac = (ModularRingFactory) new ModIntegerRing(modul.getVal());
1838 }
1839 //System.out.println("mcfac = " + mcfac);
1840 mfac = new GenPolynomialRing<MOD>(mcfac, fac);
1841 lift = PolyUtil.<MOD> fromIntegerCoefficients(mfac, Fi);
1842 //System.out.println("lift = " + lift + ": " + lift.get(0).ring.coFac);
1843 return lift;
1844 }
1845
1846 }