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