001 /*
002 * $Id: PolyUtilApp.java 3753 2011-08-27 20:34:30Z kredel $
003 */
004
005 package edu.jas.application;
006
007
008 import java.util.ArrayList;
009 import java.util.Arrays;
010 import java.util.List;
011 import java.util.Map;
012 import java.util.Set;
013 import java.util.SortedMap;
014 import java.util.TreeMap;
015 import java.util.TreeSet;
016
017 import org.apache.log4j.Logger;
018
019 import edu.jas.arith.BigDecimal;
020 import edu.jas.arith.Product;
021 import edu.jas.arith.ProductRing;
022 import edu.jas.arith.Rational;
023 import edu.jas.poly.AlgebraicNumber;
024 import edu.jas.poly.AlgebraicNumberRing;
025 import edu.jas.poly.Complex;
026 import edu.jas.poly.ComplexRing;
027 import edu.jas.poly.ExpVector;
028 import edu.jas.poly.GenPolynomial;
029 import edu.jas.poly.GenPolynomialRing;
030 import edu.jas.poly.PolyUtil;
031 import edu.jas.poly.PolynomialList;
032 import edu.jas.poly.TermOrder;
033 import edu.jas.root.ComplexRootsAbstract;
034 import edu.jas.root.ComplexRootsSturm;
035 import edu.jas.root.Interval;
036 import edu.jas.root.RealAlgebraicNumber;
037 import edu.jas.root.RealAlgebraicRing;
038 import edu.jas.root.RealRootAbstract;
039 import edu.jas.root.RealRootTuple;
040 import edu.jas.root.RealRootsSturm;
041 import edu.jas.root.RootFactory;
042 import edu.jas.structure.GcdRingElem;
043 import edu.jas.structure.RingElem;
044 import edu.jas.structure.RingFactory;
045 import edu.jas.structure.UnaryFunctor;
046 import edu.jas.util.ListUtil;
047
048
049 /**
050 * Polynomial utilities for applications, for example conversion ExpVector to
051 * Product or zero dimensional ideal root computation.
052 * @param <C> coefficient type
053 * @author Heinz Kredel
054 */
055 public class PolyUtilApp<C extends RingElem<C>> {
056
057
058 private static final Logger logger = Logger.getLogger(PolyUtilApp.class);
059
060
061 private static boolean debug = logger.isDebugEnabled();
062
063
064 /**
065 * Product representation.
066 * @param <C> coefficient type.
067 * @param pfac polynomial ring factory.
068 * @param L list of polynomials to be represented.
069 * @return Product represenation of L in the polynomial ring pfac.
070 */
071 public static <C extends GcdRingElem<C>> List<GenPolynomial<Product<Residue<C>>>> toProductRes(
072 GenPolynomialRing<Product<Residue<C>>> pfac, List<GenPolynomial<GenPolynomial<C>>> L) {
073
074 List<GenPolynomial<Product<Residue<C>>>> list = new ArrayList<GenPolynomial<Product<Residue<C>>>>();
075 if (L == null || L.size() == 0) {
076 return list;
077 }
078 GenPolynomial<Product<Residue<C>>> b;
079 for (GenPolynomial<GenPolynomial<C>> a : L) {
080 b = toProductRes(pfac, a);
081 list.add(b);
082 }
083 return list;
084 }
085
086
087 /**
088 * Product representation.
089 * @param <C> coefficient type.
090 * @param pfac polynomial ring factory.
091 * @param A polynomial to be represented.
092 * @return Product represenation of A in the polynomial ring pfac.
093 */
094 public static <C extends GcdRingElem<C>> GenPolynomial<Product<Residue<C>>> toProductRes(
095 GenPolynomialRing<Product<Residue<C>>> pfac, GenPolynomial<GenPolynomial<C>> A) {
096
097 GenPolynomial<Product<Residue<C>>> P = pfac.getZERO().clone();
098 if (A == null || A.isZERO()) {
099 return P;
100 }
101 RingFactory<Product<Residue<C>>> rpfac = pfac.coFac;
102 ProductRing<Residue<C>> fac = (ProductRing<Residue<C>>) rpfac;
103 Product<Residue<C>> p;
104 for (Map.Entry<ExpVector, GenPolynomial<C>> y : A.getMap().entrySet()) {
105 ExpVector e = y.getKey();
106 GenPolynomial<C> a = y.getValue();
107 p = toProductRes(fac, a);
108 if (p != null && !p.isZERO()) {
109 P.doPutToMap(e, p);
110 }
111 }
112 return P;
113 }
114
115
116 /**
117 * Product representation.
118 * @param <C> coefficient type.
119 * @param pfac product ring factory.
120 * @param c coefficient to be represented.
121 * @return Product represenation of c in the ring pfac.
122 */
123 public static <C extends GcdRingElem<C>> Product<Residue<C>> toProductRes(ProductRing<Residue<C>> pfac,
124 GenPolynomial<C> c) {
125
126 SortedMap<Integer, Residue<C>> elem = new TreeMap<Integer, Residue<C>>();
127 for (int i = 0; i < pfac.length(); i++) {
128 RingFactory<Residue<C>> rfac = pfac.getFactory(i);
129 ResidueRing<C> fac = (ResidueRing<C>) rfac;
130 Residue<C> u = new Residue<C>(fac, c);
131 //fac.fromInteger( c.getVal() );
132 if (u != null && !u.isZERO()) {
133 elem.put(i, u);
134 }
135 }
136 return new Product<Residue<C>>(pfac, elem);
137 }
138
139
140 /**
141 * Product residue representation.
142 * @param <C> coefficient type.
143 * @param CS list of ColoredSystems from comprehensive GB system.
144 * @return Product residue represenation of CS.
145 */
146 public static <C extends GcdRingElem<C>> List<GenPolynomial<Product<Residue<C>>>> toProductRes(
147 List<ColoredSystem<C>> CS) {
148
149 List<GenPolynomial<Product<Residue<C>>>> list = new ArrayList<GenPolynomial<Product<Residue<C>>>>();
150 if (CS == null || CS.size() == 0) {
151 return list;
152 }
153 GenPolynomialRing<GenPolynomial<C>> pr = null;
154 List<RingFactory<Residue<C>>> rrl = new ArrayList<RingFactory<Residue<C>>>(CS.size());
155 for (ColoredSystem<C> cs : CS) {
156 Ideal<C> id = cs.condition.zero;
157 ResidueRing<C> r = new ResidueRing<C>(id);
158 if (!rrl.contains(r)) {
159 rrl.add(r);
160 }
161 if (pr == null) {
162 if (cs.list.size() > 0) {
163 pr = cs.list.get(0).green.ring;
164 }
165 }
166 }
167 ProductRing<Residue<C>> pfac;
168 pfac = new ProductRing<Residue<C>>(rrl);
169 //System.out.println("pfac = " + pfac);
170 GenPolynomialRing<Product<Residue<C>>> rf = new GenPolynomialRing<Product<Residue<C>>>(pfac, pr.nvar,
171 pr.tord, pr.getVars());
172 GroebnerSystem<C> gs = new GroebnerSystem<C>(CS);
173 List<GenPolynomial<GenPolynomial<C>>> F = gs.getCGB();
174 list = PolyUtilApp.<C> toProductRes(rf, F);
175 return list;
176 }
177
178
179 /**
180 * Residue coefficient representation.
181 * @param pfac polynomial ring factory.
182 * @param L list of polynomials to be represented.
183 * @return Represenation of L in the polynomial ring pfac.
184 */
185 public static <C extends GcdRingElem<C>> List<GenPolynomial<Residue<C>>> toResidue(
186 GenPolynomialRing<Residue<C>> pfac, List<GenPolynomial<GenPolynomial<C>>> L) {
187 List<GenPolynomial<Residue<C>>> list = new ArrayList<GenPolynomial<Residue<C>>>();
188 if (L == null || L.size() == 0) {
189 return list;
190 }
191 GenPolynomial<Residue<C>> b;
192 for (GenPolynomial<GenPolynomial<C>> a : L) {
193 b = toResidue(pfac, a);
194 if (b != null && !b.isZERO()) {
195 list.add(b);
196 }
197 }
198 return list;
199 }
200
201
202 /**
203 * Residue coefficient representation.
204 * @param pfac polynomial ring factory.
205 * @param A polynomial to be represented.
206 * @return Represenation of A in the polynomial ring pfac.
207 */
208 public static <C extends GcdRingElem<C>> GenPolynomial<Residue<C>> toResidue(
209 GenPolynomialRing<Residue<C>> pfac, GenPolynomial<GenPolynomial<C>> A) {
210 GenPolynomial<Residue<C>> P = pfac.getZERO().clone();
211 if (A == null || A.isZERO()) {
212 return P;
213 }
214 RingFactory<Residue<C>> rpfac = pfac.coFac;
215 ResidueRing<C> fac = (ResidueRing<C>) rpfac;
216 Residue<C> p;
217 for (Map.Entry<ExpVector, GenPolynomial<C>> y : A.getMap().entrySet()) {
218 ExpVector e = y.getKey();
219 GenPolynomial<C> a = y.getValue();
220 p = new Residue<C>(fac, a);
221 if (p != null && !p.isZERO()) {
222 P.doPutToMap(e, p);
223 }
224 }
225 return P;
226 }
227
228
229 /**
230 * Product slice.
231 * @param <C> coefficient type.
232 * @param L list of polynomials with product coefficients.
233 * @return Slices represenation of L.
234 */
235 public static <C extends GcdRingElem<C>> Map<Ideal<C>, PolynomialList<GenPolynomial<C>>> productSlice(
236 PolynomialList<Product<Residue<C>>> L) {
237
238 Map<Ideal<C>, PolynomialList<GenPolynomial<C>>> map;
239 RingFactory<Product<Residue<C>>> fpr = L.ring.coFac;
240 ProductRing<Residue<C>> pr = (ProductRing<Residue<C>>) fpr;
241 int s = pr.length();
242 map = new TreeMap<Ideal<C>, PolynomialList<GenPolynomial<C>>>();
243 List<GenPolynomial<GenPolynomial<C>>> slist;
244
245 List<GenPolynomial<Product<Residue<C>>>> plist = L.list;
246 PolynomialList<GenPolynomial<C>> spl;
247
248 for (int i = 0; i < s; i++) {
249 RingFactory<Residue<C>> r = pr.getFactory(i);
250 ResidueRing<C> rr = (ResidueRing<C>) r;
251 Ideal<C> id = rr.ideal;
252 GenPolynomialRing<C> cof = rr.ring;
253 GenPolynomialRing<GenPolynomial<C>> pfc;
254 pfc = new GenPolynomialRing<GenPolynomial<C>>(cof, L.ring);
255 slist = fromProduct(pfc, plist, i);
256 spl = new PolynomialList<GenPolynomial<C>>(pfc, slist);
257 PolynomialList<GenPolynomial<C>> d = map.get(id);
258 if (d != null) {
259 throw new RuntimeException("ideal exists twice " + id);
260 }
261 map.put(id, spl);
262 }
263 return map;
264 }
265
266
267 /**
268 * Product slice at i.
269 * @param <C> coefficient type.
270 * @param L list of polynomials with product coeffients.
271 * @param i index of slice.
272 * @return Slice of of L at i.
273 */
274 public static <C extends GcdRingElem<C>> PolynomialList<GenPolynomial<C>> productSlice(
275 PolynomialList<Product<Residue<C>>> L, int i) {
276
277 RingFactory<Product<Residue<C>>> fpr = L.ring.coFac;
278 ProductRing<Residue<C>> pr = (ProductRing<Residue<C>>) fpr;
279 List<GenPolynomial<GenPolynomial<C>>> slist;
280
281 List<GenPolynomial<Product<Residue<C>>>> plist = L.list;
282 PolynomialList<GenPolynomial<C>> spl;
283
284 RingFactory<Residue<C>> r = pr.getFactory(i);
285 ResidueRing<C> rr = (ResidueRing<C>) r;
286 GenPolynomialRing<C> cof = rr.ring;
287 GenPolynomialRing<GenPolynomial<C>> pfc;
288 pfc = new GenPolynomialRing<GenPolynomial<C>>(cof, L.ring);
289 slist = fromProduct(pfc, plist, i);
290 spl = new PolynomialList<GenPolynomial<C>>(pfc, slist);
291 return spl;
292 }
293
294
295 /**
296 * From product representation.
297 * @param <C> coefficient type.
298 * @param pfac polynomial ring factory.
299 * @param L list of polynomials to be converted from product representation.
300 * @param i index of product representation to be taken.
301 * @return Represenation of i-slice of L in the polynomial ring pfac.
302 */
303 public static <C extends GcdRingElem<C>> List<GenPolynomial<GenPolynomial<C>>> fromProduct(
304 GenPolynomialRing<GenPolynomial<C>> pfac, List<GenPolynomial<Product<Residue<C>>>> L,
305 int i) {
306
307 List<GenPolynomial<GenPolynomial<C>>> list = new ArrayList<GenPolynomial<GenPolynomial<C>>>();
308
309 if (L == null || L.size() == 0) {
310 return list;
311 }
312 GenPolynomial<GenPolynomial<C>> b;
313 for (GenPolynomial<Product<Residue<C>>> a : L) {
314 b = fromProduct(pfac, a, i);
315 if (b != null && !b.isZERO()) {
316 b = b.abs();
317 if (!list.contains(b)) {
318 list.add(b);
319 }
320 }
321 }
322 return list;
323 }
324
325
326 /**
327 * From product representation.
328 * @param <C> coefficient type.
329 * @param pfac polynomial ring factory.
330 * @param P polynomial to be converted from product representation.
331 * @param i index of product representation to be taken.
332 * @return Represenation of i-slice of P in the polynomial ring pfac.
333 */
334 public static <C extends GcdRingElem<C>> GenPolynomial<GenPolynomial<C>> fromProduct(
335 GenPolynomialRing<GenPolynomial<C>> pfac, GenPolynomial<Product<Residue<C>>> P, int i) {
336
337 GenPolynomial<GenPolynomial<C>> b = pfac.getZERO().clone();
338 if (P == null || P.isZERO()) {
339 return b;
340 }
341
342 for (Map.Entry<ExpVector, Product<Residue<C>>> y : P.getMap().entrySet()) {
343 ExpVector e = y.getKey();
344 Product<Residue<C>> a = y.getValue();
345 Residue<C> r = a.get(i);
346 if (r != null && !r.isZERO()) {
347 GenPolynomial<C> p = r.val;
348 if (p != null && !p.isZERO()) {
349 b.doPutToMap(e, p);
350 }
351 }
352 }
353 return b;
354 }
355
356
357 /**
358 * Product slice to String.
359 * @param <C> coefficient type.
360 * @param L list of polynomials with to be represented.
361 * @return Product represenation of L in the polynomial ring pfac.
362 */
363 public static <C extends GcdRingElem<C>> String productSliceToString(
364 Map<Ideal<C>, PolynomialList<GenPolynomial<C>>> L) {
365 Set<GenPolynomial<GenPolynomial<C>>> sl = new TreeSet<GenPolynomial<GenPolynomial<C>>>();
366 PolynomialList<GenPolynomial<C>> pl = null;
367 StringBuffer sb = new StringBuffer(); //"\nproductSlice ----------------- begin");
368 for (Ideal<C> id : L.keySet()) {
369 sb.append("\n\ncondition == 0:\n");
370 sb.append(id.list.toScript());
371 pl = L.get(id);
372 sl.addAll(pl.list);
373 sb.append("\ncorresponding ideal:\n");
374 sb.append(pl.toScript());
375 }
376 //List<GenPolynomial<GenPolynomial<C>>> sll
377 // = new ArrayList<GenPolynomial<GenPolynomial<C>>>( sl );
378 //pl = new PolynomialList<GenPolynomial<C>>(pl.ring,sll);
379 // sb.append("\nunion = " + pl.toString());
380 //sb.append("\nproductSlice ------------------------- end\n");
381 return sb.toString();
382 }
383
384
385 /**
386 * Product slice to String.
387 * @param <C> coefficient type.
388 * @param L list of polynomials with product coefficients.
389 * @return string represenation of slices of L.
390 */
391 public static <C extends GcdRingElem<C>> String productToString(PolynomialList<Product<Residue<C>>> L) {
392 Map<Ideal<C>, PolynomialList<GenPolynomial<C>>> M;
393 M = productSlice(L);
394 String s = productSliceToString(M);
395 return s;
396 }
397
398
399 /**
400 * Construct superset of complex roots for zero dimensional ideal(G).
401 * @param I zero dimensional ideal.
402 * @param eps desired precision.
403 * @return list of coordinates of complex roots for ideal(G)
404 */
405 public static <C extends RingElem<C> & Rational, D extends GcdRingElem<D>> List<List<Complex<BigDecimal>>> complexRootTuples(
406 Ideal<D> I, C eps) {
407 List<GenPolynomial<D>> univs = I.constructUnivariate();
408 if (logger.isInfoEnabled()) {
409 logger.info("univs = " + univs);
410 }
411 return complexRoots(I, univs, eps);
412 }
413
414
415 /**
416 * Construct superset of complex roots for zero dimensional ideal(G).
417 * @param I zero dimensional ideal.
418 * @param univs list of univariate polynomials.
419 * @param eps desired precision.
420 * @return list of coordinates of complex roots for ideal(G)
421 */
422 public static <C extends RingElem<C> & Rational, D extends GcdRingElem<D>> List<List<Complex<BigDecimal>>> complexRoots(
423 Ideal<D> I, List<GenPolynomial<D>> univs, C eps) {
424 List<List<Complex<BigDecimal>>> croots = new ArrayList<List<Complex<BigDecimal>>>();
425 RingFactory<C> cf = (RingFactory<C>) I.list.ring.coFac;
426 ComplexRing<C> cr = new ComplexRing<C>(cf);
427 ComplexRootsAbstract<C> cra = new ComplexRootsSturm<C>(cr);
428 List<GenPolynomial<Complex<C>>> cunivs = new ArrayList<GenPolynomial<Complex<C>>>();
429 for (GenPolynomial<D> p : univs) {
430 GenPolynomialRing<Complex<C>> pfac = new GenPolynomialRing<Complex<C>>(cr, p.ring);
431 //System.out.println("pfac = " + pfac.toScript());
432 GenPolynomial<Complex<C>> cp = PolyUtil.<C> toComplex(pfac, (GenPolynomial<C>) p);
433 cunivs.add(cp);
434 //System.out.println("cp = " + cp);
435 }
436 for (int i = 0; i < I.list.ring.nvar; i++) {
437 List<Complex<BigDecimal>> cri = cra.approximateRoots(cunivs.get(i), eps);
438 //System.out.println("cri = " + cri);
439 croots.add(cri);
440 }
441 croots = ListUtil.<Complex<BigDecimal>> tupleFromList(croots);
442 return croots;
443 }
444
445
446 /**
447 * Construct superset of complex roots for zero dimensional ideal(G).
448 * @param Il list of zero dimensional ideals with univariate polynomials.
449 * @param eps desired precision.
450 * @return list of coordinates of complex roots for ideal(cap_i(G_i))
451 */
452 public static <C extends RingElem<C> & Rational, D extends GcdRingElem<D>> List<List<Complex<BigDecimal>>> complexRootTuples(
453 List<IdealWithUniv<D>> Il, C eps) {
454 List<List<Complex<BigDecimal>>> croots = new ArrayList<List<Complex<BigDecimal>>>();
455 for (IdealWithUniv<D> I : Il) {
456 List<List<Complex<BigDecimal>>> cr = complexRoots(I.ideal, I.upolys, eps);
457 croots.addAll(cr);
458 }
459 return croots;
460 }
461
462
463 /**
464 * Construct superset of complex roots for zero dimensional ideal(G).
465 * @param Il list of zero dimensional ideals with univariate polynomials.
466 * @param eps desired precision.
467 * @return list of ideals with coordinates of complex roots for
468 * ideal(cap_i(G_i))
469 */
470 public static <C extends RingElem<C> & Rational, D extends GcdRingElem<D>> List<IdealWithComplexRoots<D>> complexRoots(
471 List<IdealWithUniv<D>> Il, C eps) {
472 List<IdealWithComplexRoots<D>> Ic = new ArrayList<IdealWithComplexRoots<D>>(Il.size());
473 for (IdealWithUniv<D> I : Il) {
474 List<List<Complex<BigDecimal>>> cr = complexRoots(I.ideal, I.upolys, eps);
475 IdealWithComplexRoots<D> ic = new IdealWithComplexRoots<D>(I, cr);
476 Ic.add(ic);
477 }
478 return Ic;
479 }
480
481
482 /**
483 * Construct superset of complex roots for zero dimensional ideal(G).
484 * @param G list of polynomials of a of zero dimensional ideal.
485 * @param eps desired precision.
486 * @return list of ideals with coordinates of complex roots for ideal(G)
487 */
488 public static <C extends RingElem<C> & Rational, D extends GcdRingElem<D>> List<IdealWithComplexRoots<D>> complexRoots(
489 Ideal<D> G, C eps) {
490 List<IdealWithUniv<D>> Il = G.zeroDimDecomposition();
491 return complexRoots(Il, eps);
492 }
493
494
495 /**
496 * Construct superset of real roots for zero dimensional ideal(G).
497 * @param I zero dimensional ideal.
498 * @param eps desired precision.
499 * @return list of coordinates of real roots for ideal(G)
500 */
501 public static <C extends RingElem<C> & Rational, D extends GcdRingElem<D>> List<List<BigDecimal>> realRootTuples(
502 Ideal<D> I, C eps) {
503 List<GenPolynomial<D>> univs = I.constructUnivariate();
504 if (logger.isInfoEnabled()) {
505 logger.info("univs = " + univs);
506 }
507 return realRoots(I, univs, eps);
508 }
509
510
511 /**
512 * Construct superset of real roots for zero dimensional ideal(G).
513 * @param I zero dimensional ideal.
514 * @param univs list of univariate polynomials.
515 * @param eps desired precision.
516 * @return list of coordinates of real roots for ideal(G)
517 */
518 public static <C extends RingElem<C> & Rational, D extends GcdRingElem<D>> List<List<BigDecimal>> realRoots(
519 Ideal<D> I, List<GenPolynomial<D>> univs, C eps) {
520 List<List<BigDecimal>> roots = new ArrayList<List<BigDecimal>>();
521 //RingFactory<C> cf = (RingFactory<C>) I.list.ring.coFac;
522 RealRootAbstract<C> rra = new RealRootsSturm<C>();
523 for (int i = 0; i < I.list.ring.nvar; i++) {
524 List<BigDecimal> rri = rra.approximateRoots((GenPolynomial<C>) univs.get(i), eps);
525 //System.out.println("rri = " + rri);
526 roots.add(rri);
527 }
528 //System.out.println("roots-1 = " + roots);
529 roots = ListUtil.<BigDecimal> tupleFromList(roots);
530 //System.out.println("roots-2 = " + roots);
531 return roots;
532 }
533
534
535 /**
536 * Construct superset of real roots for zero dimensional ideal(G).
537 * @param Il list of zero dimensional ideals with univariate polynomials.
538 * @param eps desired precision.
539 * @return list of coordinates of real roots for ideal(cap_i(G_i))
540 */
541 public static <C extends RingElem<C> & Rational, D extends GcdRingElem<D>> List<List<BigDecimal>> realRootTuples(
542 List<IdealWithUniv<D>> Il, C eps) {
543 List<List<BigDecimal>> rroots = new ArrayList<List<BigDecimal>>();
544 for (IdealWithUniv<D> I : Il) {
545 List<List<BigDecimal>> rr = realRoots(I.ideal, I.upolys, eps);
546 rroots.addAll(rr);
547 }
548 return rroots;
549 }
550
551
552 /**
553 * Construct superset of real roots for zero dimensional ideal(G).
554 * @param Il list of zero dimensional ideals with univariate polynomials.
555 * @param eps desired precision.
556 * @return list of ideals with coordinates of real roots for
557 * ideal(cap_i(G_i))
558 */
559 public static <C extends RingElem<C> & Rational, D extends GcdRingElem<D>> List<IdealWithRealRoots<D>> realRoots(
560 List<IdealWithUniv<D>> Il, C eps) {
561 List<IdealWithRealRoots<D>> Ir = new ArrayList<IdealWithRealRoots<D>>(Il.size());
562 for (IdealWithUniv<D> I : Il) {
563 List<List<BigDecimal>> rr = realRoots(I.ideal, I.upolys, eps);
564 IdealWithRealRoots<D> ir = new IdealWithRealRoots<D>(I, rr);
565 Ir.add(ir);
566 }
567 return Ir;
568 }
569
570
571 /**
572 * Construct superset of real roots for zero dimensional ideal(G).
573 * @param G list of polynomials of a of zero dimensional ideal.
574 * @param eps desired precision.
575 * @return list of ideals with coordinates of real roots for ideal(G)
576 */
577 public static <C extends RingElem<C> & Rational, D extends GcdRingElem<D>> List<IdealWithRealRoots<D>> realRoots(
578 Ideal<D> G, C eps) {
579 List<IdealWithUniv<D>> Il = G.zeroDimDecomposition();
580 return realRoots(Il, eps);
581 }
582
583
584 /**
585 * Test for real roots of zero dimensional ideal(L).
586 * @param L list of polynomials.
587 * @param roots list of real roots for ideal(G).
588 * @param eps desired precision.
589 * @return true if root is a list of coordinates of real roots for ideal(L)
590 */
591 public static boolean isRealRoots(List<GenPolynomial<BigDecimal>> L, List<List<BigDecimal>> roots,
592 BigDecimal eps) {
593 if (L == null || L.size() == 0) {
594 return true;
595 }
596 // polynomials with decimal coefficients
597 BigDecimal dc = BigDecimal.ONE;
598 GenPolynomialRing<BigDecimal> dfac = L.get(0).ring;
599 //System.out.println("dfac = " + dfac);
600 for (GenPolynomial<BigDecimal> dp : L) {
601 //System.out.println("dp = " + dp);
602 for (List<BigDecimal> r : roots) {
603 //System.out.println("r = " + r);
604 BigDecimal ev = PolyUtil.<BigDecimal> evaluateAll(dc, dfac, dp, r);
605 if (ev.abs().compareTo(eps) > 0) {
606 System.out.println("ev = " + ev);
607 return false;
608 }
609 }
610 }
611 return true;
612 }
613
614
615 /**
616 * Test for complex roots of zero dimensional ideal(L).
617 * @param L list of polynomials.
618 * @param roots list of real roots for ideal(G).
619 * @param eps desired precision.
620 * @return true if root is a list of coordinates of complex roots for
621 * ideal(L)
622 */
623 public static boolean isComplexRoots(List<GenPolynomial<Complex<BigDecimal>>> L,
624 List<List<Complex<BigDecimal>>> roots, BigDecimal eps) {
625 if (L == null || L.size() == 0) {
626 return true;
627 }
628 // polynomials with decimal coefficients
629 BigDecimal dc = BigDecimal.ONE;
630 ComplexRing<BigDecimal> dcc = new ComplexRing<BigDecimal>(dc);
631 GenPolynomialRing<Complex<BigDecimal>> dfac = L.get(0).ring;
632 //System.out.println("dfac = " + dfac);
633 for (GenPolynomial<Complex<BigDecimal>> dp : L) {
634 //System.out.println("dp = " + dp);
635 for (List<Complex<BigDecimal>> r : roots) {
636 //System.out.println("r = " + r);
637 Complex<BigDecimal> ev = PolyUtil.<Complex<BigDecimal>> evaluateAll(dcc, dfac, dp, r);
638 if (ev.norm().getRe().compareTo(eps) > 0) {
639 System.out.println("ev = " + ev);
640 return false;
641 }
642 }
643 }
644 return true;
645 }
646
647
648 /**
649 * Construct real roots for zero dimensional ideal(G).
650 * @param I zero dimensional ideal with univariate irreducible polynomials
651 * and bi-variate polynomials.
652 * @return real algebraic roots for ideal(G)
653 */
654 public static <C extends RingElem<C> & Rational, D extends GcdRingElem<D> & Rational> IdealWithRealAlgebraicRoots<C, D> realAlgebraicRoots(
655 IdealWithUniv<D> I) {
656 List<List<RealAlgebraicNumber<D>>> ran = new ArrayList<List<RealAlgebraicNumber<D>>>();
657 if (I == null || I.ideal == null || I.ideal.isZERO() || I.upolys == null || I.upolys.size() == 0) {
658 return new IdealWithRealAlgebraicRoots<C, D>(I, ran);
659 }
660 GenPolynomialRing<D> fac = I.ideal.list.ring;
661 // case i == 0:
662 GenPolynomial<D> p0 = I.upolys.get(0);
663 GenPolynomial<D> p0p = PolyUtil.<D> selectWithVariable(I.ideal.list.list, fac.nvar - 1);
664 if (p0p == null) {
665 throw new RuntimeException("no polynomial found in " + (fac.nvar - 1) + " of " + I.ideal);
666 }
667 //System.out.println("p0 = " + p0);
668 if (logger.isInfoEnabled()) {
669 logger.info("p0p = " + p0p);
670 }
671 int[] dep0 = p0p.degreeVector().dependencyOnVariables();
672 //System.out.println("dep0 = " + Arrays.toString(dep0));
673 if (dep0.length != 1) {
674 throw new RuntimeException("wrong number of variables " + Arrays.toString(dep0));
675 }
676 List<RealAlgebraicNumber<D>> rra = RootFactory.<D> realAlgebraicNumbersIrred(p0);
677 for (RealAlgebraicNumber<D> rr : rra) {
678 List<RealAlgebraicNumber<D>> rl = new ArrayList<RealAlgebraicNumber<D>>();
679 rl.add(rr);
680 ran.add(rl);
681 }
682 // case i > 0:
683 for (int i = 1; i < I.upolys.size(); i++) {
684 List<List<RealAlgebraicNumber<D>>> rn = new ArrayList<List<RealAlgebraicNumber<D>>>();
685 GenPolynomial<D> pi = I.upolys.get(i);
686 GenPolynomial<D> pip = PolyUtil.selectWithVariable(I.ideal.list.list, fac.nvar - 1 - i);
687 if (pip == null) {
688 throw new RuntimeException("no polynomial found in " + (fac.nvar - 1 - i) + " of " + I.ideal);
689 }
690 //System.out.println("i = " + i);
691 //System.out.println("pi = " + pi);
692 if (logger.isInfoEnabled()) {
693 logger.info("pi = " + pi);
694 logger.info("pip = " + pip);
695 }
696 int[] depi = pip.degreeVector().dependencyOnVariables();
697 //System.out.println("depi = " + Arrays.toString(depi));
698 if (depi.length < 1 || depi.length > 2) {
699 throw new RuntimeException("wrong number of variables " + Arrays.toString(depi));
700 }
701 rra = RootFactory.<D> realAlgebraicNumbersIrred(pi);
702 if (depi.length == 1) {
703 // all combinations are roots of the ideal I
704 for (RealAlgebraicNumber<D> rr : rra) {
705 //System.out.println("rr.ring = " + rr.ring);
706 for (List<RealAlgebraicNumber<D>> rx : ran) {
707 //System.out.println("rx = " + rx);
708 List<RealAlgebraicNumber<D>> ry = new ArrayList<RealAlgebraicNumber<D>>();
709 ry.addAll(rx);
710 ry.add(rr);
711 rn.add(ry);
712 }
713 }
714 } else { // depi.length == 2
715 // select roots of the ideal I
716 GenPolynomial<D> pip2 = PolyUtil.removeUnusedUpperVariables(pip);
717 //System.out.println("pip2 = " + pip2.ring);
718 GenPolynomialRing<D> ufac = pip2.ring.contract(1);
719 TermOrder to = new TermOrder(TermOrder.INVLEX);
720 GenPolynomialRing<GenPolynomial<D>> rfac = new GenPolynomialRing<GenPolynomial<D>>(ufac, 1,
721 to);
722 GenPolynomial<GenPolynomial<D>> pip2r = PolyUtil.<D> recursive(rfac, pip2);
723 int ix = fac.nvar - 1 - depi[depi.length - 1];
724 //System.out.println("ix = " + ix);
725 for (RealAlgebraicNumber<D> rr : rra) {
726 //System.out.println("rr.ring = " + rr.ring);
727 Interval<D> rroot = rr.ring.getRoot();
728 GenPolynomial<D> pip2el = PolyUtil.<D> evaluateMainRecursive(ufac, pip2r, rroot.left);
729 GenPolynomial<D> pip2er = PolyUtil.<D> evaluateMainRecursive(ufac, pip2r, rroot.right);
730 GenPolynomialRing<D> upfac = I.upolys.get(ix).ring;
731 GenPolynomial<D> pip2elc = convert(upfac, pip2el);
732 GenPolynomial<D> pip2erc = convert(upfac, pip2er);
733 //System.out.println("pip2elc = " + pip2elc);
734 //System.out.println("pip2erc = " + pip2erc);
735 for (List<RealAlgebraicNumber<D>> rx : ran) {
736 //System.out.println("rx = " + rx);
737 RealAlgebraicRing<D> rar = rx.get(ix).ring;
738 //System.out.println("rar = " + rar.toScript());
739 RealAlgebraicNumber<D> rel = new RealAlgebraicNumber<D>(rar, pip2elc);
740 RealAlgebraicNumber<D> rer = new RealAlgebraicNumber<D>(rar, pip2erc);
741 int sl = rel.signum();
742 int sr = rer.signum();
743 //System.out.println("sl = " + sl + ", sr = " + sr + ", sl*sr = " + (sl*sr));
744 if (sl * sr <= 0) {
745 List<RealAlgebraicNumber<D>> ry = new ArrayList<RealAlgebraicNumber<D>>();
746 ry.addAll(rx);
747 ry.add(rr);
748 rn.add(ry);
749 }
750 }
751 }
752 }
753 ran = rn;
754 }
755 IdealWithRealAlgebraicRoots<C, D> Ir = new IdealWithRealAlgebraicRoots<C, D>(I, ran);
756 return Ir;
757 }
758
759
760 /**
761 * Construct real roots for zero dimensional ideal(G).
762 * @param I list of zero dimensional ideal with univariate irreducible
763 * polynomials and bi-variate polynomials.
764 * @return list of real algebraic roots for all ideal(I_i)
765 */
766 public static <C extends RingElem<C> & Rational, D extends GcdRingElem<D> & Rational> List<IdealWithRealAlgebraicRoots<C, D>> realAlgebraicRoots(
767 List<IdealWithUniv<D>> I) {
768 List<IdealWithRealAlgebraicRoots<C, D>> lir = new ArrayList<IdealWithRealAlgebraicRoots<C, D>>(
769 I.size());
770 for (IdealWithUniv<D> iu : I) {
771 IdealWithRealAlgebraicRoots<C, D> iur = PolyUtilApp.<C, D> realAlgebraicRoots(iu);
772 //System.out.println("iur = " + iur);
773 lir.add(iur);
774 }
775 return lir;
776 }
777
778
779 /**
780 * Construct complex roots for zero dimensional ideal(G).
781 * @param I zero dimensional ideal with univariate irreducible polynomials
782 * and bi-variate polynomials.
783 * @return complex algebraic roots for ideal(G)
784 */
785 public static <C extends RingElem<C> & Rational, D extends GcdRingElem<D> & Rational> IdealWithComplexAlgebraicRoots<C, D> complexAlgebraicRoots(
786 IdealWithUniv<D> I) {
787 List<List<Complex<edu.jas.application.RealAlgebraicNumber<D>>>> can = new ArrayList<List<Complex<edu.jas.application.RealAlgebraicNumber<D>>>>();
788 if (I == null || I.ideal == null || I.ideal.isZERO() || I.upolys == null || I.upolys.size() == 0) {
789 return new IdealWithComplexAlgebraicRoots<C, D>(I, can);
790 }
791 GenPolynomialRing<D> fac = I.ideal.list.ring;
792 if (fac.nvar == 0) {
793 return new IdealWithComplexAlgebraicRoots<C, D>(I, can);
794 }
795 // case i == 0:
796 GenPolynomial<D> p0 = I.upolys.get(0);
797 GenPolynomial<D> p0p = PolyUtil.<D> selectWithVariable(I.ideal.list.list, fac.nvar - 1);
798 if (p0p == null) {
799 throw new RuntimeException("no polynomial found in " + (fac.nvar - 1) + " of " + I.ideal);
800 }
801 if (logger.isInfoEnabled()) {
802 logger.info("p0 = " + p0);
803 logger.info("p0p = " + p0p);
804 }
805 int[] dep0 = p0p.degreeVector().dependencyOnVariables();
806 //System.out.println("dep0 = " + Arrays.toString(dep0));
807 if (dep0.length != 1) {
808 throw new RuntimeException("wrong number of variables " + Arrays.toString(dep0));
809 }
810
811 RingFactory<D> cfac = p0.ring.coFac;
812 ComplexRing<D> ccfac = new ComplexRing<D>(cfac);
813 GenPolynomialRing<Complex<D>> facc = new GenPolynomialRing<Complex<D>>(ccfac, p0.ring);
814 GenPolynomial<Complex<D>> p0c = PolyUtil.<D> complexFromAny(facc, p0);
815
816 GenPolynomialRing<Complex<D>> faccu = new GenPolynomialRing<Complex<D>>(ccfac, p0p.ring);
817 //GenPolynomial<Complex<D>> p0pc = PolyUtil.<D> complexFromAny(faccu, p0p);
818 //System.out.println("p0c = " + p0c);
819 //System.out.println("p0pc = " + p0pc);
820
821 List<Complex<edu.jas.application.RealAlgebraicNumber<D>>> cra = edu.jas.application.RootFactory
822 .<D> complexAlgebraicNumbersSquarefree(p0c);
823 logger.info("#roots(p0c) = " + cra.size());
824 for (Complex<edu.jas.application.RealAlgebraicNumber<D>> cr : cra) {
825 List<Complex<edu.jas.application.RealAlgebraicNumber<D>>> cl = new ArrayList<Complex<edu.jas.application.RealAlgebraicNumber<D>>>();
826 cl.add(cr);
827 can.add(cl);
828 }
829 if (fac.nvar == 1) {
830 return new IdealWithComplexAlgebraicRoots<C, D>(I, can);
831 }
832 // case i > 0:
833 for (int i = 1; i < I.upolys.size(); i++) {
834 List<List<Complex<edu.jas.application.RealAlgebraicNumber<D>>>> cn = new ArrayList<List<Complex<edu.jas.application.RealAlgebraicNumber<D>>>>();
835 GenPolynomial<D> pi = I.upolys.get(i);
836 GenPolynomial<D> pip = PolyUtil.selectWithVariable(I.ideal.list.list, fac.nvar - 1 - i);
837 if (pip == null) {
838 throw new RuntimeException("no polynomial found in " + (fac.nvar - 1 - i) + " of " + I.ideal);
839 }
840 //System.out.println("i = " + i);
841 //System.out.println("pi = " + pi);
842 if (logger.isInfoEnabled()) {
843 logger.info("pi = " + pi);
844 logger.info("pip = " + pip);
845 }
846 facc = new GenPolynomialRing<Complex<D>>(ccfac, pi.ring);
847 GenPolynomial<Complex<D>> pic = PolyUtil.<D> complexFromAny(facc, pi);
848 faccu = new GenPolynomialRing<Complex<D>>(ccfac, pip.ring);
849 //GenPolynomial<Complex<D>> pipc = PolyUtil.<D> complexFromAny(faccu, pip);
850 //System.out.println("pic = " + pic);
851 //System.out.println("pipc = " + pipc);
852 int[] depi = pip.degreeVector().dependencyOnVariables();
853 //System.out.println("depi = " + Arrays.toString(depi));
854 if (depi.length < 1 || depi.length > 2) {
855 throw new RuntimeException("wrong number of variables " + Arrays.toString(depi) + " for "
856 + pip);
857 }
858 cra = edu.jas.application.RootFactory.<D> complexAlgebraicNumbersSquarefree(pic);
859 logger.info("#roots(pic) = " + cra.size());
860 if (false && debug) {
861 boolean t = edu.jas.application.RootFactory.<D> isRoot(pic, cra);
862 if (!t) {
863 throw new RuntimeException("no roots of " + pic);
864 }
865 }
866 if (depi.length == 1) {
867 // all combinations are roots of the ideal I
868 for (Complex<edu.jas.application.RealAlgebraicNumber<D>> cr : cra) {
869 //System.out.println("cr.ring = " + cr.ring);
870 for (List<Complex<edu.jas.application.RealAlgebraicNumber<D>>> cx : can) {
871 //System.out.println("cx = " + cx);
872 List<Complex<edu.jas.application.RealAlgebraicNumber<D>>> cy = new ArrayList<Complex<edu.jas.application.RealAlgebraicNumber<D>>>();
873 cy.addAll(cx);
874 cy.add(cr);
875 cn.add(cy);
876 }
877 }
878 } else { // depi.length == 2
879 // select roots of the ideal I
880 GenPolynomial<D> pip2 = PolyUtil.removeUnusedUpperVariables(pip);
881 //System.out.println("pip2 = " + pip2.ring);
882 GenPolynomialRing<D> ufac = pip2.ring.contract(1);
883 //System.out.println("ufac = " + ufac);
884 GenPolynomialRing<Complex<D>> ucfac = new GenPolynomialRing<Complex<D>>(ccfac, ufac);
885 //System.out.println("ucfac = " + ucfac);
886 GenPolynomial<Complex<D>> pip2c = PolyUtil.<D> complexFromAny(ucfac, pip2);
887
888 GenPolynomialRing<GenPolynomial<D>> rfac = pip2.ring.recursive(1);
889 //GenPolynomial<GenPolynomial<D>> pip2r = PolyUtil.<D> recursive(rfac, pip2);
890 //System.out.println("rfac = " + rfac);
891
892 GenPolynomialRing<GenPolynomial<Complex<D>>> rcfac = new GenPolynomialRing<GenPolynomial<Complex<D>>>(
893 ucfac, rfac);
894 GenPolynomial<GenPolynomial<Complex<D>>> pip2cr = PolyUtil.<Complex<D>> recursive(rcfac,
895 pip2c);
896 //System.out.println("rcfac = " + rcfac);
897
898 int ix = fac.nvar - 1 - depi[depi.length - 1];
899 //System.out.println("ix = " + ix);
900 for (Complex<edu.jas.application.RealAlgebraicNumber<D>> cr : cra) {
901 //System.out.println("cr.ring = " + cr.ring);
902 edu.jas.application.RealAlgebraicRing<D> cring = (edu.jas.application.RealAlgebraicRing<D>) cr.ring.ring;
903 RealRootTuple<D> rroot = cring.getRoot();
904 List<RealAlgebraicNumber<D>> rlist = rroot.tuple;
905 //System.out.println("rlist = " + rlist);
906 Interval<D> vr = rlist.get(0).ring.getRoot();
907 Interval<D> vi = rlist.get(1).ring.getRoot();
908 logger.info("vr = " + vr + ", vi = " + vi);
909 Complex<D> sw = new Complex<D>(ccfac, vr.left, vi.left);
910 Complex<D> ne = new Complex<D>(ccfac, vr.right, vi.right);
911 //System.out.println("sw = " + sw + ", ne = " + ne);
912
913 GenPolynomial<Complex<D>> pip2cesw = PolyUtil.<Complex<D>> evaluateMainRecursive(ucfac,
914 pip2cr, sw);
915 GenPolynomial<Complex<D>> pip2cene = PolyUtil.<Complex<D>> evaluateMainRecursive(ucfac,
916 pip2cr, ne);
917
918 GenPolynomialRing<D> upfac = I.upolys.get(ix).ring;
919 GenPolynomialRing<Complex<D>> upcfac = new GenPolynomialRing<Complex<D>>(ccfac, upfac);
920 //System.out.println("upfac = " + upfac);
921 //System.out.println("upcfac = " + upcfac);
922 GenPolynomial<Complex<D>> pip2eswc = convertComplexComplex(upcfac, pip2cesw);
923 GenPolynomial<Complex<D>> pip2enec = convertComplexComplex(upcfac, pip2cene);
924 //System.out.println("pip2eswc = " + pip2eswc);
925 //System.out.println("pip2enec = " + pip2enec);
926 for (List<Complex<edu.jas.application.RealAlgebraicNumber<D>>> cx : can) {
927 //System.out.println("cx = " + cx);
928 ComplexRing<edu.jas.application.RealAlgebraicNumber<D>> car = cx.get(ix).ring;
929 edu.jas.application.RealAlgebraicRing<D> rar = (edu.jas.application.RealAlgebraicRing<D>) car.ring;
930 //System.out.println("car = " + car);
931 //System.out.println("rar = " + rar);
932
933 TermOrder to = new TermOrder(TermOrder.INVLEX);
934 String vvr = rar.algebraic.ring.getVars()[0];
935 String vvi = rar.algebraic.ring.getVars()[1];
936 String[] vars = new String[] { vvr, vvi };
937 GenPolynomialRing<Complex<D>> tfac = new GenPolynomialRing<Complex<D>>(ccfac, to,
938 vars);
939 //System.out.println("tfac = " + tfac);
940 GenPolynomial<Complex<D>> t = tfac.univariate(1, 1L).sum(
941 tfac.univariate(0, 1L).multiply(ccfac.getIMAG()));
942 //System.out.println("t = " + t); // t = x + i y
943 GenPolynomialRing<D> rtfac = new GenPolynomialRing<D>(cfac, tfac);
944 //System.out.println("rtfac = " + rtfac);
945
946 GenPolynomial<Complex<D>> su;
947 GenPolynomial<D> re, im;
948
949 su = PolyUtil.<Complex<D>> substituteUnivariate(pip2eswc, t);
950 //su = su.monic(); not here
951 //System.out.println("su = " + su);
952 re = PolyUtil.<D> realPartFromComplex(rtfac, su);
953 im = PolyUtil.<D> imaginaryPartFromComplex(rtfac, su);
954 //System.out.println("re = " + re);
955 //System.out.println("im = " + im);
956 edu.jas.application.RealAlgebraicNumber<D> resw = new edu.jas.application.RealAlgebraicNumber<D>(
957 rar, re);
958 //System.out.println("resw = " + resw);
959 int sswr = resw.signum();
960 edu.jas.application.RealAlgebraicNumber<D> imsw = new edu.jas.application.RealAlgebraicNumber<D>(
961 rar, im);
962 //System.out.println("imsw = " + imsw);
963 int sswi = imsw.signum();
964 //System.out.println("sswr = " + sswr + ", sswi = " + sswi);
965
966 su = PolyUtil.<Complex<D>> substituteUnivariate(pip2enec, t);
967 //su = su.monic(); not here
968 //System.out.println("su = " + su);
969 re = PolyUtil.<D> realPartFromComplex(rtfac, su);
970 im = PolyUtil.<D> imaginaryPartFromComplex(rtfac, su);
971 //System.out.println("re = " + re);
972 //System.out.println("im = " + im);
973 edu.jas.application.RealAlgebraicNumber<D> rene = new edu.jas.application.RealAlgebraicNumber<D>(
974 rar, re);
975 //System.out.println("rene = " + rene);
976 int sner = rene.signum();
977 edu.jas.application.RealAlgebraicNumber<D> imne = new edu.jas.application.RealAlgebraicNumber<D>(
978 rar, im);
979 //System.out.println("imne = " + imne);
980 int snei = imne.signum();
981 //System.out.println("sner = " + sner + ", snei = " + snei);
982
983 //System.out.println("");
984 //System.out.println("sswr = " + sswr + ", sswi = " + sswi);
985 //System.out.println("sner = " + sner + ", snei = " + snei);
986 if (sswr * sner <= 0 && sswi * snei <= 0) {
987 //System.out.println("hit, cx = " + cx + ", cr = " + cr); //.getRe().decimalMagnitude() + " i " + cr.getIm().decimalMagnitude());
988 List<Complex<edu.jas.application.RealAlgebraicNumber<D>>> cy = new ArrayList<Complex<edu.jas.application.RealAlgebraicNumber<D>>>();
989 cy.addAll(cx);
990 cy.add(cr);
991 cn.add(cy);
992 } else {
993 //System.out.println("no hit, cx = " + cx + ", cr = " + cr); //.getRe().decimalMagnitude() + " i " + cr.getIm().decimalMagnitude());
994 }
995 }
996 }
997 }
998 can = cn;
999 }
1000 IdealWithComplexAlgebraicRoots<C, D> Ic = new IdealWithComplexAlgebraicRoots<C, D>(I, can);
1001 return Ic;
1002 }
1003
1004
1005 /**
1006 * Construct complex roots for zero dimensional ideal(G).
1007 * @param I list of zero dimensional ideal with univariate irreducible
1008 * polynomials and bi-variate polynomials.
1009 * @return list of complex algebraic roots for ideal(G)
1010 */
1011 public static <C extends RingElem<C> & Rational, D extends GcdRingElem<D> & Rational> List<IdealWithComplexAlgebraicRoots<C, D>> complexAlgebraicRoots(
1012 List<IdealWithUniv<D>> I) {
1013 List<IdealWithComplexAlgebraicRoots<C, D>> lic = new ArrayList<IdealWithComplexAlgebraicRoots<C, D>>();
1014 for (IdealWithUniv<D> iu : I) {
1015 IdealWithComplexAlgebraicRoots<C, D> iuc = PolyUtilApp.<C, D> complexAlgebraicRoots(iu);
1016 //System.out.println("iuc = " + iuc);
1017 lic.add(iuc);
1018 }
1019 return lic;
1020 }
1021
1022
1023 /*
1024 * Convert to a polynomial in given ring.
1025 * @param fac result polynomial ring.
1026 * @param p polynomial.
1027 * @return polynomial in ring fac <b>Note: </b> if p can not be represented
1028 * in fac then the results are unpredictable.
1029 */
1030 static <C extends RingElem<C>> GenPolynomial<C> convert(GenPolynomialRing<C> fac, GenPolynomial<C> p) {
1031 if (fac.equals(p.factory())) {
1032 return p;
1033 }
1034 GenPolynomial<C> q = fac.parse(p.toString());
1035 if (!q.toString().equals(p.toString())) {
1036 throw new RuntimeException("convert(" + p + ") = " + q);
1037 }
1038 return q;
1039 }
1040
1041
1042 /*
1043 * Convert to a polynomial in given ring.
1044 * @param fac result polynomial ring.
1045 * @param p polynomial.
1046 * @return polynomial in ring fac <b>Note: </b> if p can not be represented
1047 * in fac then the results are unpredictable.
1048 */
1049 static <C extends RingElem<C>> GenPolynomial<Complex<C>> convertComplex(
1050 GenPolynomialRing<Complex<C>> fac, GenPolynomial<C> p) {
1051 //if (fac.equals(p.factory())) {
1052 // return p;
1053 //}
1054 GenPolynomial<Complex<C>> q = fac.parse(p.toString());
1055 if (!q.toString().equals(p.toString())) {
1056 throw new RuntimeException("convert(" + p + ") = " + q);
1057 }
1058 return q;
1059 }
1060
1061
1062 /*
1063 * Convert to a polynomial in given ring.
1064 * @param fac result polynomial ring.
1065 * @param p complex polynomial.
1066 * @return polynomial in ring fac <b>Note: </b> if p can not be represented
1067 * in fac then the results are unpredictable.
1068 */
1069 static <C extends RingElem<C>> GenPolynomial<Complex<C>> convertComplexComplex(
1070 GenPolynomialRing<Complex<C>> fac, GenPolynomial<Complex<C>> p) {
1071 //if (fac.equals(p.factory())) {
1072 // return p;
1073 //}
1074 GenPolynomial<Complex<C>> q = fac.parse(p.toString());
1075 if (!q.toString().equals(p.toString())) {
1076 throw new RuntimeException("convert(" + p + ") = " + q);
1077 }
1078 return q;
1079 }
1080
1081
1082 /**
1083 * Construct exact set of real roots for zero dimensional ideal(G).
1084 * @param I zero dimensional ideal.
1085 * @return list of coordinates of real roots for ideal(G)
1086 */
1087 public static <C extends RingElem<C> & Rational, D extends GcdRingElem<D> & Rational> List<IdealWithRealAlgebraicRoots<C, D>> realAlgebraicRoots(
1088 Ideal<D> I) {
1089 List<IdealWithUniv<D>> Ir = I.zeroDimRootDecomposition();
1090 //System.out.println("Ir = " + Ir);
1091 List<IdealWithRealAlgebraicRoots<C, D>> roots = PolyUtilApp.<C, D> realAlgebraicRoots(Ir);
1092 return roots;
1093 }
1094
1095
1096 /**
1097 * Construct primitive element for double field extension.
1098 * @param a algebraic number ring with squarefree monic minimal polynomial
1099 * @param b algebraic number ring with squarefree monic minimal polynomial
1100 * @return primitive element container with algebraic number ring c, with
1101 * Q(c) = Q(a,b)
1102 */
1103 public static <C extends GcdRingElem<C>> PrimitiveElement<C> primitiveElement(AlgebraicNumberRing<C> a,
1104 AlgebraicNumberRing<C> b) {
1105 GenPolynomial<C> ap = a.modul;
1106 GenPolynomial<C> bp = b.modul;
1107
1108 // setup bivariate polynomial ring
1109 String[] cv = new String[2];
1110 cv[0] = ap.ring.getVars()[0];
1111 cv[1] = bp.ring.getVars()[0];
1112 TermOrder to = new TermOrder(TermOrder.INVLEX);
1113 GenPolynomialRing<C> cfac = new GenPolynomialRing<C>(ap.ring.coFac, 2, to, cv);
1114 GenPolynomial<C> as = ap.extendUnivariate(cfac, 0);
1115 GenPolynomial<C> bs = bp.extendUnivariate(cfac, 1);
1116 List<GenPolynomial<C>> L = new ArrayList<GenPolynomial<C>>(2);
1117 L.add(as);
1118 L.add(bs);
1119 List<GenPolynomial<C>> Op = new ArrayList<GenPolynomial<C>>();
1120
1121 Ideal<C> id = new Ideal<C>(cfac, L);
1122 //System.out.println("id = " + id);
1123 IdealWithUniv<C> iu = id.normalPositionFor(0, 1, Op);
1124 //System.out.println("iu = " + iu);
1125
1126 // extract result polynomials
1127 List<GenPolynomial<C>> Np = iu.ideal.getList();
1128 //System.out.println("Np = " + Np);
1129 as = PolyUtil.<C> selectWithVariable(Np, 1);
1130 bs = PolyUtil.<C> selectWithVariable(Np, 0);
1131 GenPolynomial<C> cs = PolyUtil.<C> selectWithVariable(Np, 2);
1132 //System.out.println("as = " + as);
1133 //System.out.println("bs = " + bs);
1134 //System.out.println("cs = " + cs);
1135 String[] ev = new String[] { cs.ring.getVars()[0] };
1136 GenPolynomialRing<C> efac = new GenPolynomialRing<C>(ap.ring.coFac, 1, to, ev);
1137 //System.out.println("efac = " + efac);
1138 cs = cs.contractCoeff(efac);
1139 //System.out.println("cs = " + cs);
1140 as = as.reductum().contractCoeff(efac);
1141 as = as.negate();
1142 //System.out.println("as = " + as);
1143 bs = bs.reductum().contractCoeff(efac);
1144 bs = bs.negate();
1145 //System.out.println("bs = " + bs);
1146 AlgebraicNumberRing<C> c = new AlgebraicNumberRing<C>(cs);
1147 AlgebraicNumber<C> ab = new AlgebraicNumber<C>(c, as);
1148 AlgebraicNumber<C> bb = new AlgebraicNumber<C>(c, bs);
1149 PrimitiveElement<C> pe = new PrimitiveElement<C>(c, ab, bb, a, b);
1150 if (logger.isInfoEnabled()) {
1151 logger.info("primitive element = " + c);
1152 }
1153 return pe;
1154 }
1155
1156
1157 /**
1158 * Convert to primitive element ring.
1159 * @param cfac primitive element ring.
1160 * @param A algebraic number representing the generating element of a in the
1161 * new ring.
1162 * @param a algebraic number to convert.
1163 * @return a converted to the primitive element ring
1164 */
1165 public static <C extends GcdRingElem<C>> AlgebraicNumber<C> convertToPrimitiveElem(
1166 AlgebraicNumberRing<C> cfac, AlgebraicNumber<C> A, AlgebraicNumber<C> a) {
1167 GenPolynomialRing<C> aufac = a.ring.ring;
1168 GenPolynomialRing<AlgebraicNumber<C>> ar = new GenPolynomialRing<AlgebraicNumber<C>>(cfac, aufac);
1169 GenPolynomial<AlgebraicNumber<C>> aps = PolyUtil.<C> convertToAlgebraicCoefficients(ar, a.val);
1170 AlgebraicNumber<C> ac = PolyUtil.<AlgebraicNumber<C>> evaluateMain(cfac, aps, A);
1171 return ac;
1172 }
1173
1174
1175 /**
1176 * Convert coefficients to primitive element ring.
1177 * @param cfac primitive element ring.
1178 * @param A algebraic number representing the generating element of a in the
1179 * new ring.
1180 * @param a polynomial with coefficients algebraic number to convert.
1181 * @return a with coefficients converted to the primitive element ring
1182 */
1183 public static <C extends GcdRingElem<C>> GenPolynomial<AlgebraicNumber<C>> convertToPrimitiveElem(
1184 AlgebraicNumberRing<C> cfac, AlgebraicNumber<C> A, GenPolynomial<AlgebraicNumber<C>> a) {
1185 GenPolynomialRing<AlgebraicNumber<C>> cr = new GenPolynomialRing<AlgebraicNumber<C>>(cfac, a.ring);
1186 return PolyUtil.<AlgebraicNumber<C>, AlgebraicNumber<C>> map(cr, a, new CoeffConvertAlg<C>(cfac, A));
1187 }
1188
1189
1190 /**
1191 * Convert to primitive element ring.
1192 * @param cfac primitive element ring.
1193 * @param A algebraic number representing the generating element of a in the
1194 * new ring.
1195 * @param a recursive algebraic number to convert.
1196 * @return a converted to the primitive element ring
1197 */
1198 public static <C extends GcdRingElem<C>> AlgebraicNumber<C> convertToPrimitiveElem(
1199 AlgebraicNumberRing<C> cfac, AlgebraicNumber<C> A, AlgebraicNumber<C> B,
1200 AlgebraicNumber<AlgebraicNumber<C>> a) {
1201 GenPolynomial<AlgebraicNumber<C>> aps = PolyUtilApp.<C> convertToPrimitiveElem(cfac, A, a.val);
1202 AlgebraicNumber<C> ac = PolyUtil.<AlgebraicNumber<C>> evaluateMain(cfac, aps, B);
1203 return ac;
1204 }
1205
1206
1207 /**
1208 * Construct primitive element for double field extension.
1209 * @param b algebraic number ring with squarefree monic minimal polynomial
1210 * over Q(a)
1211 * @return primitive element container with algebraic number ring c, with
1212 * Q(c) = Q(a)(b)
1213 */
1214 public static <C extends GcdRingElem<C>> PrimitiveElement<C> primitiveElement(
1215 AlgebraicNumberRing<AlgebraicNumber<C>> b) {
1216 GenPolynomial<AlgebraicNumber<C>> bp = b.modul;
1217 AlgebraicNumberRing<C> a = (AlgebraicNumberRing<C>) b.ring.coFac;
1218 GenPolynomial<C> ap = a.modul;
1219
1220 // setup bivariate polynomial ring
1221 String[] cv = new String[2];
1222 cv[0] = ap.ring.getVars()[0];
1223 cv[1] = bp.ring.getVars()[0];
1224 TermOrder to = new TermOrder(TermOrder.INVLEX);
1225 GenPolynomialRing<C> cfac = new GenPolynomialRing<C>(ap.ring.coFac, 2, to, cv);
1226 GenPolynomialRing<GenPolynomial<C>> rfac = new GenPolynomialRing<GenPolynomial<C>>(a.ring, 1,
1227 bp.ring.getVars());
1228 GenPolynomial<C> as = ap.extendUnivariate(cfac, 0);
1229 GenPolynomial<GenPolynomial<C>> bss = PolyUtil.<C> fromAlgebraicCoefficients(rfac, bp);
1230 GenPolynomial<C> bs = PolyUtil.<C> distribute(cfac, bss);
1231 List<GenPolynomial<C>> L = new ArrayList<GenPolynomial<C>>(2);
1232 L.add(as);
1233 L.add(bs);
1234 List<GenPolynomial<C>> Op = new ArrayList<GenPolynomial<C>>();
1235
1236 Ideal<C> id = new Ideal<C>(cfac, L);
1237 //System.out.println("id = " + id);
1238 IdealWithUniv<C> iu = id.normalPositionFor(0, 1, Op);
1239 //System.out.println("iu = " + iu);
1240
1241 // extract result polynomials
1242 List<GenPolynomial<C>> Np = iu.ideal.getList();
1243 as = PolyUtil.<C> selectWithVariable(Np, 1);
1244 bs = PolyUtil.<C> selectWithVariable(Np, 0);
1245 GenPolynomial<C> cs = PolyUtil.<C> selectWithVariable(Np, 2);
1246 //System.out.println("as = " + as);
1247 //System.out.println("bs = " + bs);
1248 //System.out.println("cs = " + cs);
1249 String[] ev = new String[] { cs.ring.getVars()[0] };
1250 GenPolynomialRing<C> efac = new GenPolynomialRing<C>(ap.ring.coFac, 1, to, ev);
1251 // System.out.println("efac = " + efac);
1252 cs = cs.contractCoeff(efac);
1253 // System.out.println("cs = " + cs);
1254 as = as.reductum().contractCoeff(efac);
1255 as = as.negate();
1256 // System.out.println("as = " + as);
1257 bs = bs.reductum().contractCoeff(efac);
1258 bs = bs.negate();
1259 //System.out.println("bs = " + bs);
1260 AlgebraicNumberRing<C> c = new AlgebraicNumberRing<C>(cs);
1261 AlgebraicNumber<C> ab = new AlgebraicNumber<C>(c, as);
1262 AlgebraicNumber<C> bb = new AlgebraicNumber<C>(c, bs);
1263 PrimitiveElement<C> pe = new PrimitiveElement<C>(c, ab, bb); // missing ,a,b);
1264 if (logger.isInfoEnabled()) {
1265 logger.info("primitive element = " + pe);
1266 }
1267 return pe;
1268 }
1269
1270
1271 /**
1272 * Convert to primitive element ring.
1273 * @param cfac primitive element ring.
1274 * @param A algebraic number representing the generating element of a in the
1275 * new ring.
1276 * @param a polynomial with recursive algebraic number coefficients to
1277 * convert.
1278 * @return a converted to the primitive element ring
1279 */
1280 public static <C extends GcdRingElem<C>> GenPolynomial<AlgebraicNumber<C>> convertToPrimitiveElem(
1281 AlgebraicNumberRing<C> cfac, AlgebraicNumber<C> A, AlgebraicNumber<C> B,
1282 GenPolynomial<AlgebraicNumber<AlgebraicNumber<C>>> a) {
1283 GenPolynomialRing<AlgebraicNumber<C>> cr = new GenPolynomialRing<AlgebraicNumber<C>>(cfac, a.ring);
1284 return PolyUtil.<AlgebraicNumber<AlgebraicNumber<C>>, AlgebraicNumber<C>> map(cr, a,
1285 new CoeffRecConvertAlg<C>(cfac, A, B));
1286 }
1287
1288
1289 /**
1290 * Convert to RealAlgebraicNumber coefficients. Represent as polynomial with
1291 * RealAlgebraicNumber<C> coefficients from package
1292 *
1293 * <pre>
1294 * edu.jas.root
1295 * </pre>
1296 *
1297 * .
1298 * @param afac result polynomial factory.
1299 * @param A polynomial with RealAlgebraicNumber<C> coefficients to be
1300 * converted.
1301 * @return polynomial with RealAlgebraicNumber<C> coefficients.
1302 */
1303 public static <C extends GcdRingElem<C> & Rational> GenPolynomial<edu.jas.root.RealAlgebraicNumber<C>> realAlgFromRealCoefficients(
1304 GenPolynomialRing<edu.jas.root.RealAlgebraicNumber<C>> afac,
1305 GenPolynomial<edu.jas.application.RealAlgebraicNumber<C>> A) {
1306 edu.jas.root.RealAlgebraicRing<C> cfac = (edu.jas.root.RealAlgebraicRing<C>) afac.coFac;
1307 return PolyUtil.<edu.jas.application.RealAlgebraicNumber<C>, edu.jas.root.RealAlgebraicNumber<C>> map(
1308 afac, A, new ReAlgFromRealCoeff<C>(cfac));
1309 }
1310
1311
1312 /**
1313 * Convert to RealAlgebraicNumber coefficients. Represent as polynomial with
1314 * RealAlgebraicNumber<C> coefficients from package
1315 *
1316 * <pre>
1317 * edu.jas.application
1318 * </pre>
1319 *
1320 * .
1321 * @param rfac result polynomial factory.
1322 * @param A polynomial with RealAlgebraicNumber<C> coefficients to be
1323 * converted.
1324 * @return polynomial with RealAlgebraicNumber<C> coefficients.
1325 */
1326 public static <C extends GcdRingElem<C> & Rational> GenPolynomial<edu.jas.application.RealAlgebraicNumber<C>> realFromRealAlgCoefficients(
1327 GenPolynomialRing<edu.jas.application.RealAlgebraicNumber<C>> rfac,
1328 GenPolynomial<edu.jas.root.RealAlgebraicNumber<C>> A) {
1329 edu.jas.application.RealAlgebraicRing<C> cfac = (edu.jas.application.RealAlgebraicRing<C>) rfac.coFac;
1330 return PolyUtil.<edu.jas.root.RealAlgebraicNumber<C>, edu.jas.application.RealAlgebraicNumber<C>> map(
1331 rfac, A, new RealFromReAlgCoeff<C>(cfac));
1332 }
1333
1334
1335 /**
1336 * Convert to Complex<RealAlgebraicNumber> coefficients. Represent as
1337 * polynomial with Complex<RealAlgebraicNumber> coefficients, C is
1338 * e.g. BigRational.
1339 * @param pfac result polynomial factory.
1340 * @param A polynomial with Complex coefficients to be converted.
1341 * @return polynomial with Complex<RealAlgebraicNumber> coefficients.
1342 */
1343 public static <C extends GcdRingElem<C> & Rational> GenPolynomial<Complex<edu.jas.application.RealAlgebraicNumber<C>>> convertToComplexRealCoefficients(
1344 GenPolynomialRing<Complex<edu.jas.application.RealAlgebraicNumber<C>>> pfac,
1345 GenPolynomial<Complex<C>> A) {
1346 ComplexRing<edu.jas.application.RealAlgebraicNumber<C>> afac = (ComplexRing<edu.jas.application.RealAlgebraicNumber<C>>) pfac.coFac;
1347 return PolyUtil.<Complex<C>, Complex<edu.jas.application.RealAlgebraicNumber<C>>> map(pfac, A,
1348 new CoeffToComplexReal<C>(afac));
1349 }
1350
1351
1352 }
1353
1354
1355 /**
1356 * Coefficient to convert algebriac functor.
1357 */
1358 class CoeffConvertAlg<C extends GcdRingElem<C>> implements
1359 UnaryFunctor<AlgebraicNumber<C>, AlgebraicNumber<C>> {
1360
1361
1362 final protected AlgebraicNumberRing<C> afac;
1363
1364
1365 final protected AlgebraicNumber<C> A;
1366
1367
1368 public CoeffConvertAlg(AlgebraicNumberRing<C> fac, AlgebraicNumber<C> a) {
1369 if (fac == null || a == null) {
1370 throw new IllegalArgumentException("fac and a must not be null");
1371 }
1372 afac = fac;
1373 A = a;
1374 }
1375
1376
1377 public AlgebraicNumber<C> eval(AlgebraicNumber<C> c) {
1378 if (c == null) {
1379 return afac.getZERO();
1380 }
1381 return PolyUtilApp.<C> convertToPrimitiveElem(afac, A, c);
1382 }
1383 }
1384
1385
1386 /**
1387 * Coefficient recursive to convert algebriac functor.
1388 */
1389 class CoeffRecConvertAlg<C extends GcdRingElem<C>> implements
1390 UnaryFunctor<AlgebraicNumber<AlgebraicNumber<C>>, AlgebraicNumber<C>> {
1391
1392
1393 final protected AlgebraicNumberRing<C> afac;
1394
1395
1396 final protected AlgebraicNumber<C> A;
1397
1398
1399 final protected AlgebraicNumber<C> B;
1400
1401
1402 public CoeffRecConvertAlg(AlgebraicNumberRing<C> fac, AlgebraicNumber<C> a, AlgebraicNumber<C> b) {
1403 if (fac == null || a == null || b == null) {
1404 throw new IllegalArgumentException("fac, a and b must not be null");
1405 }
1406 afac = fac;
1407 A = a;
1408 B = b;
1409 }
1410
1411
1412 public AlgebraicNumber<C> eval(AlgebraicNumber<AlgebraicNumber<C>> c) {
1413 if (c == null) {
1414 return afac.getZERO();
1415 }
1416 return PolyUtilApp.<C> convertToPrimitiveElem(afac, A, B, c);
1417 }
1418 }
1419
1420
1421 /**
1422 * Coefficient to real algebriac from real algebraic functor.
1423 */
1424 class ReAlgFromRealCoeff<C extends GcdRingElem<C> & Rational> implements
1425 UnaryFunctor<edu.jas.application.RealAlgebraicNumber<C>, edu.jas.root.RealAlgebraicNumber<C>> {
1426
1427
1428 final protected edu.jas.root.RealAlgebraicRing<C> afac;
1429
1430
1431 public ReAlgFromRealCoeff(edu.jas.root.RealAlgebraicRing<C> fac) {
1432 if (fac == null) {
1433 throw new IllegalArgumentException("fac must not be null");
1434 }
1435 afac = fac;
1436 }
1437
1438
1439 public edu.jas.root.RealAlgebraicNumber<C> eval(edu.jas.application.RealAlgebraicNumber<C> c) {
1440 if (c == null) {
1441 return afac.getZERO();
1442 }
1443 return (edu.jas.root.RealAlgebraicNumber<C>) (Object) c.number; // force ignore recursion
1444 }
1445 }
1446
1447
1448 /**
1449 * Coefficient to real algebriac from algebraic functor.
1450 */
1451 class RealFromReAlgCoeff<C extends GcdRingElem<C> & Rational> implements
1452 UnaryFunctor<edu.jas.root.RealAlgebraicNumber<C>, edu.jas.application.RealAlgebraicNumber<C>> {
1453
1454
1455 final protected edu.jas.application.RealAlgebraicRing<C> rfac;
1456
1457
1458 public RealFromReAlgCoeff(edu.jas.application.RealAlgebraicRing<C> fac) {
1459 if (fac == null) {
1460 throw new IllegalArgumentException("fac must not be null");
1461 }
1462 rfac = fac;
1463 }
1464
1465
1466 public edu.jas.application.RealAlgebraicNumber<C> eval(edu.jas.root.RealAlgebraicNumber<C> c) {
1467 if (c == null) {
1468 return rfac.getZERO();
1469 }
1470 edu.jas.root.RealAlgebraicNumber<edu.jas.root.RealAlgebraicNumber<C>> rrc = (edu.jas.root.RealAlgebraicNumber<edu.jas.root.RealAlgebraicNumber<C>>) (Object) c; // force resurrect recursion
1471 return new edu.jas.application.RealAlgebraicNumber<C>(rfac, rrc);
1472 }
1473 }
1474
1475
1476 /**
1477 * Coefficient to complex real algebriac functor.
1478 */
1479 class CoeffToComplexReal<C extends GcdRingElem<C> & Rational> implements
1480 UnaryFunctor<Complex<C>, Complex<edu.jas.application.RealAlgebraicNumber<C>>> {
1481
1482
1483 final protected ComplexRing<edu.jas.application.RealAlgebraicNumber<C>> cfac;
1484
1485
1486 final edu.jas.application.RealAlgebraicRing<C> afac;
1487
1488
1489 final GenPolynomialRing<C> pfac;
1490
1491
1492 public CoeffToComplexReal(ComplexRing<edu.jas.application.RealAlgebraicNumber<C>> fac) {
1493 if (fac == null) {
1494 throw new IllegalArgumentException("fac must not be null");
1495 }
1496 cfac = fac;
1497 afac = (edu.jas.application.RealAlgebraicRing<C>) cfac.ring;
1498 pfac = afac.univs.ideal.getRing();
1499 //List<GenPolynomial<C>> gens = pfac.getGenerators();
1500 //realOne = gens.get(1);
1501 //imagOne = gens.get(2);
1502 }
1503
1504
1505 public Complex<edu.jas.application.RealAlgebraicNumber<C>> eval(Complex<C> c) {
1506 if (c == null) {
1507 return cfac.getZERO();
1508 }
1509 GenPolynomial<C> pr = new GenPolynomial<C>(pfac, c.getRe());
1510 GenPolynomial<C> pi = new GenPolynomial<C>(pfac, c.getIm());
1511 //System.out.println("pr = " + pr);
1512 //System.out.println("pi = " + pi);
1513 edu.jas.application.RealAlgebraicNumber<C> re = new edu.jas.application.RealAlgebraicNumber<C>(afac,
1514 pr);
1515 edu.jas.application.RealAlgebraicNumber<C> im = new edu.jas.application.RealAlgebraicNumber<C>(afac,
1516 pi);
1517 //System.out.println("re = " + re);
1518 //System.out.println("im = " + im);
1519 return new Complex<edu.jas.application.RealAlgebraicNumber<C>>(cfac, re, im);
1520 }
1521 }