001/* 002 * $Id: PolyUtilApp.java 6007 2020-03-29 13:34:49Z 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.SortedMap; 013import java.util.TreeMap; 014 015import org.apache.logging.log4j.Logger; 016import org.apache.logging.log4j.LogManager; 017 018import edu.jas.arith.BigDecimal; 019import edu.jas.arith.BigRational; 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 = LogManager.getLogger(PolyUtilApp.class); 062 063 064 private static final 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, BigRational 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, BigRational 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, BigRational 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, BigRational 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, BigRational 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 BigRational 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, BigRational 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, BigRational 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, BigRational 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 BigRational 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( 708 "no polynomial found in " + (fac.nvar - 1 - i) + " of " + I.ideal); 709 } 710 //System.out.println("i = " + i); 711 //System.out.println("pi = " + pi); 712 if (logger.isInfoEnabled()) { 713 logger.info("pi = " + pi); 714 logger.info("pip = " + pip); 715 } 716 int[] depi = pip.degreeVector().dependencyOnVariables(); 717 //System.out.println("depi = " + Arrays.toString(depi)); 718 if (depi.length < 1 || depi.length > 2) { 719 throw new RuntimeException("wrong number of variables " + Arrays.toString(depi)); 720 } 721 rra = RootFactory.<D> realAlgebraicNumbersIrred(pi); 722 if (logger.isInfoEnabled()) { 723 List<Interval<D>> il = new ArrayList<Interval<D>>(); 724 for (RealAlgebraicNumber<D> rr : rra) { 725 il.add(rr.ring.getRoot()); 726 } 727 logger.info("roots(pi) = " + il); 728 } 729 if (depi.length == 1) { 730 // all combinations are roots of the ideal I 731 for (RealAlgebraicNumber<D> rr : rra) { 732 //System.out.println("rr.ring = " + rr.ring); 733 for (List<RealAlgebraicNumber<D>> rx : ran) { 734 //System.out.println("rx = " + rx); 735 List<RealAlgebraicNumber<D>> ry = new ArrayList<RealAlgebraicNumber<D>>(); 736 ry.addAll(rx); 737 ry.add(rr); 738 rn.add(ry); 739 } 740 } 741 } else { // depi.length == 2 742 // select roots of the ideal I 743 GenPolynomial<D> pip2 = PolyUtil.<D> removeUnusedUpperVariables(pip); 744 //System.out.println("pip2 = " + pip2.ring); 745 GenPolynomialRing<D> ufac = pip2.ring.contract(1); 746 TermOrder to = new TermOrder(TermOrder.INVLEX); 747 GenPolynomialRing<GenPolynomial<D>> rfac = new GenPolynomialRing<GenPolynomial<D>>(ufac, 1, 748 to); // new vars 749 GenPolynomial<GenPolynomial<D>> pip2r = PolyUtil.<D> recursive(rfac, pip2); 750 int ix = fac.nvar - 1 - depi[depi.length - 1]; 751 //System.out.println("ix = " + ix); 752 for (RealAlgebraicNumber<D> rr : rra) { 753 //System.out.println("rr.ring = " + rr.ring); 754 Interval<D> rroot = rr.ring.getRoot(); 755 GenPolynomial<D> pip2el = PolyUtil.<D> evaluateMainRecursive(ufac, pip2r, rroot.left); 756 GenPolynomial<D> pip2er = PolyUtil.<D> evaluateMainRecursive(ufac, pip2r, rroot.right); 757 GenPolynomialRing<D> upfac = I.upolys.get(ix).ring; 758 GenPolynomial<D> pip2elc = convert(upfac, pip2el); 759 GenPolynomial<D> pip2erc = convert(upfac, pip2er); 760 //System.out.println("pip2elc = " + pip2elc); 761 //System.out.println("pip2erc = " + pip2erc); 762 for (List<RealAlgebraicNumber<D>> rx : ran) { 763 //System.out.println("rx = " + rx); 764 RealAlgebraicRing<D> rar = rx.get(ix).ring; 765 //System.out.println("rar = " + rar.toScript()); 766 RealAlgebraicNumber<D> rel = new RealAlgebraicNumber<D>(rar, pip2elc); 767 RealAlgebraicNumber<D> rer = new RealAlgebraicNumber<D>(rar, pip2erc); 768 int sl = rel.signum(); 769 int sr = rer.signum(); 770 //System.out.println("sl = " + sl + ", sr = " + sr + ", sl*sr = " + (sl*sr)); 771 if (sl * sr <= 0) { 772 //System.out.println("sl * sr <= 0: rar = " + rar.toScript()); 773 List<RealAlgebraicNumber<D>> ry = new ArrayList<RealAlgebraicNumber<D>>(); 774 ry.addAll(rx); 775 ry.add(rr); 776 rn.add(ry); 777 } 778 } 779 } 780 } 781 ran = rn; 782 } 783 if (logger.isInfoEnabled()) { 784 for (List<RealAlgebraicNumber<D>> rz : ran) { 785 List<Interval<D>> il = new ArrayList<Interval<D>>(); 786 for (RealAlgebraicNumber<D> rr : rz) { 787 il.add(rr.ring.getRoot()); 788 } 789 logger.info("root-tuple = " + il); 790 } 791 } 792 IdealWithRealAlgebraicRoots<D> Ir = new IdealWithRealAlgebraicRoots<D>(I, ran); 793 return Ir; 794 } 795 796 797 /** 798 * Construct real roots for zero dimensional ideal(G). 799 * @param I list of zero dimensional ideal with univariate irreducible 800 * polynomials and bi-variate polynomials. 801 * @return list of real algebraic roots for all ideal(I_i) 802 */ 803 public static <D extends GcdRingElem<D> & Rational> List<IdealWithRealAlgebraicRoots<D>> realAlgebraicRoots( 804 List<IdealWithUniv<D>> I) { 805 List<IdealWithRealAlgebraicRoots<D>> lir = new ArrayList<IdealWithRealAlgebraicRoots<D>>(I.size()); 806 for (IdealWithUniv<D> iu : I) { 807 IdealWithRealAlgebraicRoots<D> iur = PolyUtilApp.<D> realAlgebraicRoots(iu); 808 //System.out.println("iur = " + iur); 809 lir.add(iur); 810 } 811 return lir; 812 } 813 814 815 /** 816 * Construct complex roots for zero dimensional ideal(G). 817 * @param I zero dimensional ideal with univariate irreducible polynomials 818 * and bi-variate polynomials. 819 * @return complex algebraic roots for ideal(G) <b>Note:</b> implementation 820 * contains errors, do not use. 821 */ 822 public static <D extends GcdRingElem<D> & Rational> IdealWithComplexAlgebraicRoots<D> complexAlgebraicRootsWrong( // Wrong 823 IdealWithUniv<D> I) { 824 List<List<Complex<edu.jas.application.RealAlgebraicNumber<D>>>> can; 825 can = new ArrayList<List<Complex<edu.jas.application.RealAlgebraicNumber<D>>>>(); 826 if (I == null) { 827 throw new IllegalArgumentException("null ideal not permitted"); 828 } 829 if (I.ideal == null || I.upolys == null) { 830 throw new IllegalArgumentException("null ideal components not permitted " + I); 831 } 832 if (I.ideal.isZERO() || I.upolys.size() == 0) { 833 return new IdealWithComplexAlgebraicRoots<D>(I, can); 834 } 835 GenPolynomialRing<D> fac = I.ideal.list.ring; 836 if (fac.nvar == 0) { 837 return new IdealWithComplexAlgebraicRoots<D>(I, can); 838 } 839 if (fac.nvar != I.upolys.size()) { 840 throw new IllegalArgumentException("ideal not zero dimnsional: " + I); 841 } 842 // case i == 0: 843 GenPolynomial<D> p0 = I.upolys.get(0); 844 GenPolynomial<D> p0p = PolyUtil.<D> selectWithVariable(I.ideal.list.list, fac.nvar - 1); 845 if (p0p == null) { 846 throw new RuntimeException("no polynomial found in " + (fac.nvar - 1) + " of " + I.ideal); 847 } 848 if (logger.isInfoEnabled()) { 849 logger.info("p0 = " + p0); 850 logger.info("p0p = " + p0p); 851 } 852 int[] dep0 = p0p.degreeVector().dependencyOnVariables(); 853 //System.out.println("dep0 = " + Arrays.toString(dep0)); 854 if (dep0.length != 1) { 855 throw new RuntimeException("wrong number of variables " + Arrays.toString(dep0)); 856 } 857 RingFactory<D> cfac = p0.ring.coFac; 858 ComplexRing<D> ccfac = new ComplexRing<D>(cfac); 859 GenPolynomialRing<Complex<D>> facc = new GenPolynomialRing<Complex<D>>(ccfac, p0.ring); 860 GenPolynomial<Complex<D>> p0c = PolyUtil.<D> complexFromAny(facc, p0); 861 List<Complex<edu.jas.application.RealAlgebraicNumber<D>>> cra; 862 cra = edu.jas.application.RootFactoryApp.<D> complexAlgebraicNumbersSquarefree(p0c); 863 logger.info("#roots(p0c) = " + cra.size()); 864 if (debug) { 865 boolean t = edu.jas.application.RootFactoryApp.<D> isRoot(p0c, cra); 866 if (!t) { 867 throw new RuntimeException("no roots of " + p0c); 868 } 869 } 870 for (Complex<edu.jas.application.RealAlgebraicNumber<D>> cr : cra) { 871 List<Complex<edu.jas.application.RealAlgebraicNumber<D>>> cl; 872 cl = new ArrayList<Complex<edu.jas.application.RealAlgebraicNumber<D>>>(); 873 cl.add(cr); 874 can.add(cl); 875 } 876 if (fac.nvar == 1) { 877 return new IdealWithComplexAlgebraicRoots<D>(I, can); 878 } 879 // case i > 0: 880 for (int i = 1; i < I.upolys.size(); i++) { 881 List<List<Complex<edu.jas.application.RealAlgebraicNumber<D>>>> cn; 882 cn = new ArrayList<List<Complex<edu.jas.application.RealAlgebraicNumber<D>>>>(); 883 GenPolynomial<D> pi = I.upolys.get(i); 884 GenPolynomial<D> pip = PolyUtil.selectWithVariable(I.ideal.list.list, fac.nvar - 1 - i); 885 if (pip == null) { 886 throw new RuntimeException( 887 "no polynomial found in " + (fac.nvar - 1 - i) + " of " + I.ideal); 888 } 889 if (logger.isInfoEnabled()) { 890 logger.info("pi(" + i + ") = " + pi); 891 logger.info("pip = " + pip); 892 } 893 facc = new GenPolynomialRing<Complex<D>>(ccfac, pi.ring); 894 GenPolynomial<Complex<D>> pic = PolyUtil.<D> complexFromAny(facc, pi); 895 int[] depi = pip.degreeVector().dependencyOnVariables(); 896 //System.out.println("depi = " + Arrays.toString(depi)); 897 if (depi.length < 1 || depi.length > 2) { 898 throw new RuntimeException( 899 "wrong number of variables " + Arrays.toString(depi) + " for " + pip); 900 } 901 cra = edu.jas.application.RootFactoryApp.<D> complexAlgebraicNumbersSquarefree(pic); 902 logger.info("#roots(pic) = " + cra.size()); 903 if (debug) { 904 boolean t = edu.jas.application.RootFactoryApp.<D> isRoot(pic, cra); 905 if (!t) { 906 throw new RuntimeException("no roots of " + pic); 907 } 908 } 909 if (depi.length == 1) { 910 // all combinations are roots of the ideal I 911 for (Complex<edu.jas.application.RealAlgebraicNumber<D>> cr : cra) { 912 //System.out.println("cr.ring = " + cr.ring); 913 for (List<Complex<edu.jas.application.RealAlgebraicNumber<D>>> cx : can) { 914 //System.out.println("cx = " + cx); 915 List<Complex<edu.jas.application.RealAlgebraicNumber<D>>> cy; 916 cy = new ArrayList<Complex<edu.jas.application.RealAlgebraicNumber<D>>>(); 917 cy.addAll(cx); 918 cy.add(cr); 919 cn.add(cy); 920 } 921 } 922 } else { // depi.length == 2 923 // select roots of the ideal I 924 GenPolynomial<D> pip2 = PolyUtil.<D> removeUnusedUpperVariables(pip); 925 GenPolynomialRing<GenPolynomial<D>> rfac = pip2.ring.recursive(1); 926 GenPolynomialRing<D> ufac = pip2.ring.contract(1); 927 GenPolynomialRing<Complex<D>> ucfac = new GenPolynomialRing<Complex<D>>(ccfac, ufac); 928 GenPolynomialRing<Complex<D>> c2fac = new GenPolynomialRing<Complex<D>>(ccfac, pip2.ring); 929 GenPolynomial<Complex<D>> pip2c = PolyUtil.<D> complexFromAny(c2fac, pip2); 930 //System.out.println("pip2c = " + pip2c); 931 GenPolynomialRing<GenPolynomial<Complex<D>>> rcfac; 932 rcfac = new GenPolynomialRing<GenPolynomial<Complex<D>>>(ucfac, rfac); 933 GenPolynomial<GenPolynomial<Complex<D>>> pip2cr = PolyUtil.<Complex<D>> recursive(rcfac, 934 pip2c); 935 //System.out.println("pip2cr = " + pip2cr); 936 937 int ix = fac.nvar - 1 - depi[depi.length - 1]; 938 //System.out.println("ix = " + ix); 939 for (Complex<edu.jas.application.RealAlgebraicNumber<D>> cr : cra) { 940 System.out.println("cr = " + toString(cr)); // <---------------------------------- 941 edu.jas.application.RealAlgebraicRing<D> cring = (edu.jas.application.RealAlgebraicRing<D>) cr.ring.ring; 942 RealRootTuple<D> rroot = cring.getRoot(); 943 List<RealAlgebraicNumber<D>> rlist = rroot.tuple; 944 Interval<D> vr = rlist.get(0).ring.getRoot(); 945 Interval<D> vi = rlist.get(1).ring.getRoot(); 946 logger.info("vr = " + vr + ", vi = " + vi); 947 if (vr.length().isZERO()) { 948 D e = vr.left.factory().parse("1/2"); 949 D m = vr.left; //middle(); 950 vr = new Interval<D>(m.subtract(e), m.sum(e)); 951 logger.info("|vr| == 0: " + vr); 952 } 953 if (vi.length().isZERO()) { 954 D e = vi.left.factory().parse("1/2"); 955 D m = vi.left; //middle(); 956 vi = new Interval<D>(m.subtract(e), m.sum(e)); 957 logger.info("|vi| == 0: " + vi); 958 } 959 Complex<D> sw = new Complex<D>(ccfac, vr.left, vi.left); 960 Complex<D> ne = new Complex<D>(ccfac, vr.right, vi.right); 961 logger.info("sw = " + toString1(sw) + ", ne = " + toString1(ne)); 962 GenPolynomial<Complex<D>> pip2cesw, pip2cene; 963 pip2cesw = PolyUtil.<Complex<D>> evaluateMainRecursive(ucfac, pip2cr, sw); 964 pip2cene = PolyUtil.<Complex<D>> evaluateMainRecursive(ucfac, pip2cr, ne); 965 GenPolynomialRing<D> upfac = I.upolys.get(ix).ring; 966 GenPolynomialRing<Complex<D>> upcfac = new GenPolynomialRing<Complex<D>>(ccfac, upfac); 967 //System.out.println("upfac = " + upfac); 968 //System.out.println("upcfac = " + upcfac); 969 GenPolynomial<Complex<D>> pip2eswc = convertComplexComplex(upcfac, pip2cesw); 970 GenPolynomial<Complex<D>> pip2enec = convertComplexComplex(upcfac, pip2cene); 971 //System.out.println("pip2eswc = " + pip2eswc); 972 //System.out.println("pip2enec = " + pip2enec); 973 for (List<Complex<edu.jas.application.RealAlgebraicNumber<D>>> cx : can) { 974 //System.out.println("cxi = " + toString(cx.get(ix))); 975 Complex<edu.jas.application.RealAlgebraicNumber<D>> cax = cx.get(ix); 976 ComplexRing<edu.jas.application.RealAlgebraicNumber<D>> car = cax.ring; 977 edu.jas.application.RealAlgebraicRing<D> rar = (edu.jas.application.RealAlgebraicRing<D>) car.ring; 978 //System.out.println("car = " + car); 979 //System.out.println("rar = " + rar); 980 TermOrder to = new TermOrder(TermOrder.INVLEX); 981 String vvr = rar.algebraic.ring.getVars()[0]; 982 String vvi = rar.algebraic.ring.getVars()[1]; 983 String[] vars = new String[] { vvr, vvi }; 984 GenPolynomialRing<Complex<D>> tfac = new GenPolynomialRing<Complex<D>>(ccfac, to, 985 vars); 986 GenPolynomial<Complex<D>> t = tfac.univariate(1, 1L) 987 .sum(tfac.univariate(0, 1L).multiply(ccfac.getIMAG())); 988 //System.out.println("t = " + t); // t = x + i y 989 GenPolynomialRing<D> rtfac = new GenPolynomialRing<D>(cfac, tfac); 990 GenPolynomial<Complex<D>> su; 991 GenPolynomial<D> re, im; 992 su = PolyUtil.<Complex<D>> substituteUnivariate(pip2eswc, t); 993 //su = su.monic(); not here 994 re = PolyUtil.<D> realPartFromComplex(rtfac, su); 995 im = PolyUtil.<D> imaginaryPartFromComplex(rtfac, su); 996 //System.out.println("re = " + re); 997 //System.out.println("im = " + im); 998 edu.jas.application.RealAlgebraicNumber<D> resw, imsw, rene, imne; 999 resw = new edu.jas.application.RealAlgebraicNumber<D>(rar, re); 1000 //System.out.println("resw = " + resw); 1001 int sswr = resw.signum(); 1002 imsw = new edu.jas.application.RealAlgebraicNumber<D>(rar, im); 1003 //System.out.println("imsw = " + imsw); 1004 int sswi = imsw.signum(); 1005 su = PolyUtil.<Complex<D>> substituteUnivariate(pip2enec, t); 1006 //su = su.monic(); not here 1007 re = PolyUtil.<D> realPartFromComplex(rtfac, su); 1008 im = PolyUtil.<D> imaginaryPartFromComplex(rtfac, su); 1009 //System.out.println("re = " + re); 1010 //System.out.println("im = " + im); 1011 rene = new edu.jas.application.RealAlgebraicNumber<D>(rar, re); 1012 //System.out.println("rene = " + rene); 1013 int sner = rene.signum(); 1014 imne = new edu.jas.application.RealAlgebraicNumber<D>(rar, im); 1015 //System.out.println("imne = " + imne); 1016 int snei = imne.signum(); 1017 //System.out.println("sswr = " + sswr + ", sswi = " + sswi); 1018 //System.out.println("sner = " + sner + ", snei = " + snei); 1019 if ((sswr * sner <= 0 && sswi * snei <= 0)) { // wrong ! 1020 logger.info(" hit, cxi = " + toString(cx.get(ix)) + ", cr = " + toString(cr)); 1021 List<Complex<edu.jas.application.RealAlgebraicNumber<D>>> cy; 1022 cy = new ArrayList<Complex<edu.jas.application.RealAlgebraicNumber<D>>>(); 1023 cy.addAll(cx); 1024 cy.add(cr); 1025 cn.add(cy); 1026 } else { 1027 logger.info("no hit, cxi = " + toString(cx.get(ix)) + ", cr = " + toString(cr)); 1028 } 1029 } 1030 } 1031 } 1032 can = cn; 1033 } 1034 IdealWithComplexAlgebraicRoots<D> Ic = new IdealWithComplexAlgebraicRoots<D>(I, can); 1035 return Ic; 1036 } 1037 1038 1039 /** 1040 * Construct complex roots for zero dimensional ideal(G). 1041 * @param I zero dimensional ideal with univariate irreducible polynomials 1042 * and bi-variate polynomials. 1043 * @return complex algebraic roots for ideal(G) <b>Note:</b> not jet 1044 * completed oin all cases. 1045 */ 1046 public static <D extends GcdRingElem<D> & Rational> IdealWithComplexAlgebraicRoots<D> complexAlgebraicRoots( 1047 IdealWithUniv<D> I) { 1048 List<List<Complex<edu.jas.application.RealAlgebraicNumber<D>>>> can; 1049 can = new ArrayList<List<Complex<edu.jas.application.RealAlgebraicNumber<D>>>>(); 1050 if (I == null) { 1051 throw new IllegalArgumentException("null ideal not permitted"); 1052 } 1053 if (I.ideal == null || I.upolys == null) { 1054 throw new IllegalArgumentException("null ideal components not permitted " + I); 1055 } 1056 if (I.ideal.isZERO() || I.upolys.size() == 0) { 1057 return new IdealWithComplexAlgebraicRoots<D>(I, can); 1058 } 1059 GenPolynomialRing<D> fac = I.ideal.list.ring; 1060 if (fac.nvar == 0) { 1061 return new IdealWithComplexAlgebraicRoots<D>(I, can); 1062 } 1063 if (fac.nvar != I.upolys.size()) { 1064 throw new IllegalArgumentException("ideal not zero dimnsional: " + I); 1065 } 1066 // case i == 0: 1067 GenPolynomial<D> p0 = I.upolys.get(0); 1068 GenPolynomial<D> p0p = PolyUtil.<D> selectWithVariable(I.ideal.list.list, fac.nvar - 1); 1069 if (p0p == null) { 1070 throw new RuntimeException("no polynomial found in " + (fac.nvar - 1) + " of " + I.ideal); 1071 } 1072 if (logger.isInfoEnabled()) { 1073 logger.info("p0 = " + p0); 1074 logger.info("p0p = " + p0p); 1075 } 1076 int[] dep0 = p0p.degreeVector().dependencyOnVariables(); 1077 //System.out.println("dep0 = " + Arrays.toString(dep0)); 1078 if (dep0.length != 1) { 1079 throw new RuntimeException("wrong number of variables " + Arrays.toString(dep0)); 1080 } 1081 RingFactory<D> cfac = p0.ring.coFac; 1082 ComplexRing<D> ccfac = new ComplexRing<D>(cfac); 1083 GenPolynomialRing<Complex<D>> facc = new GenPolynomialRing<Complex<D>>(ccfac, p0.ring); 1084 GenPolynomial<Complex<D>> p0c = PolyUtil.<D> complexFromAny(facc, p0); 1085 List<Complex<edu.jas.application.RealAlgebraicNumber<D>>> cra; 1086 cra = edu.jas.application.RootFactoryApp.<D> complexAlgebraicNumbersSquarefree(p0c); 1087 logger.info("#roots(p0c) = " + cra.size()); 1088 for (Complex<edu.jas.application.RealAlgebraicNumber<D>> cr : cra) { 1089 List<Complex<edu.jas.application.RealAlgebraicNumber<D>>> cl; 1090 cl = new ArrayList<Complex<edu.jas.application.RealAlgebraicNumber<D>>>(); 1091 cl.add(cr); 1092 can.add(cl); 1093 } 1094 if (fac.nvar == 1) { 1095 return new IdealWithComplexAlgebraicRoots<D>(I, can); 1096 } 1097 // case i > 0: 1098 for (int i = 1; i < I.upolys.size(); i++) { 1099 List<List<Complex<edu.jas.application.RealAlgebraicNumber<D>>>> cn; 1100 cn = new ArrayList<List<Complex<edu.jas.application.RealAlgebraicNumber<D>>>>(); 1101 GenPolynomial<D> pi = I.upolys.get(i); 1102 GenPolynomial<D> pip = PolyUtil.selectWithVariable(I.ideal.list.list, fac.nvar - 1 - i); 1103 if (pip == null) { 1104 throw new RuntimeException( 1105 "no polynomial found in " + (fac.nvar - 1 - i) + " of " + I.ideal); 1106 } 1107 if (logger.isInfoEnabled()) { 1108 logger.info("pi(" + i + ") = " + pi); 1109 logger.info("pip = " + pip); 1110 } 1111 facc = new GenPolynomialRing<Complex<D>>(ccfac, pi.ring); 1112 GenPolynomial<Complex<D>> pic = PolyUtil.<D> complexFromAny(facc, pi); 1113 int[] depi = pip.degreeVector().dependencyOnVariables(); 1114 //System.out.println("depi = " + Arrays.toString(depi)); 1115 if (depi.length < 1 || depi.length > 2) { 1116 throw new RuntimeException( 1117 "wrong number of variables " + Arrays.toString(depi) + " for " + pip); 1118 } 1119 cra = edu.jas.application.RootFactoryApp.<D> complexAlgebraicNumbersSquarefree(pic); 1120 logger.info("#roots(pic) = " + cra.size()); 1121 if (depi.length == 1) { 1122 // all combinations are roots of the ideal I 1123 for (Complex<edu.jas.application.RealAlgebraicNumber<D>> cr : cra) { 1124 //System.out.println("cr.ring = " + cr.ring); 1125 for (List<Complex<edu.jas.application.RealAlgebraicNumber<D>>> cx : can) { 1126 //System.out.println("cx = " + cx); 1127 List<Complex<edu.jas.application.RealAlgebraicNumber<D>>> cy; 1128 cy = new ArrayList<Complex<edu.jas.application.RealAlgebraicNumber<D>>>(); 1129 cy.addAll(cx); 1130 cy.add(cr); 1131 cn.add(cy); 1132 } 1133 } 1134 } else { // depi.length == 2 1135 // select roots of the ideal I 1136 GenPolynomial<D> pip2 = PolyUtil.<D> removeUnusedUpperVariables(pip); 1137 pip2 = PolyUtil.<D> removeUnusedLowerVariables(pip2); 1138 pip2 = PolyUtil.<D> removeUnusedMiddleVariables(pip2); 1139 GenPolynomialRing<GenPolynomial<D>> rfac = pip2.ring.recursive(1); 1140 GenPolynomialRing<D> ufac = pip2.ring.contract(1); 1141 GenPolynomialRing<Complex<D>> ucfac = new GenPolynomialRing<Complex<D>>(ccfac, ufac); 1142 GenPolynomialRing<Complex<D>> c2fac = new GenPolynomialRing<Complex<D>>(ccfac, pip2.ring); 1143 GenPolynomial<Complex<D>> pip2c = PolyUtil.<D> complexFromAny(c2fac, pip2); 1144 GenPolynomialRing<GenPolynomial<Complex<D>>> rcfac; 1145 rcfac = new GenPolynomialRing<GenPolynomial<Complex<D>>>(ucfac, rfac); 1146 GenPolynomial<GenPolynomial<Complex<D>>> pip2cr = PolyUtil.<Complex<D>> recursive(rcfac, 1147 pip2c); 1148 //System.out.println("pip2cr = " + pip2cr); 1149 int ix = fac.nvar - 1 - depi[depi.length - 1]; 1150 //System.out.println("ix = " + ix); 1151 for (Complex<edu.jas.application.RealAlgebraicNumber<D>> cr : cra) { 1152 //System.out.println("cr = " + toString(cr)); 1153 edu.jas.application.RealAlgebraicRing<D> cring = (edu.jas.application.RealAlgebraicRing<D>) cr.ring.ring; 1154 RealRootTuple<D> rroot = cring.getRoot(); 1155 List<RealAlgebraicNumber<D>> rlist = rroot.tuple; 1156 //System.out.println("rlist = " + rlist); 1157 Interval<D> vr = rlist.get(0).ring.getRoot(); 1158 Interval<D> vi = rlist.get(1).ring.getRoot(); 1159 //logger.info("vr = " + vr + ", vi = " + vi); 1160 edu.jas.application.RealAlgebraicNumber<D> vrl, vil, vrr, vir; 1161 vrl = new edu.jas.application.RealAlgebraicNumber<D>(cring, vr.left); 1162 vil = new edu.jas.application.RealAlgebraicNumber<D>(cring, vi.left); 1163 vrr = new edu.jas.application.RealAlgebraicNumber<D>(cring, vr.right); 1164 vir = new edu.jas.application.RealAlgebraicNumber<D>(cring, vi.right); 1165 ComplexRing<edu.jas.application.RealAlgebraicNumber<D>> crr; 1166 crr = new ComplexRing<edu.jas.application.RealAlgebraicNumber<D>>(cring); 1167 Complex<edu.jas.application.RealAlgebraicNumber<D>> csw, cne; 1168 csw = new Complex<edu.jas.application.RealAlgebraicNumber<D>>(crr, vrl, vil); 1169 cne = new Complex<edu.jas.application.RealAlgebraicNumber<D>>(crr, vrr, vir); 1170 //logger.info("csw = " + toString(csw) + ", cne = " + toString(cne)); 1171 Rectangle<edu.jas.application.RealAlgebraicNumber<D>> rec; 1172 rec = new Rectangle<edu.jas.application.RealAlgebraicNumber<D>>(csw, cne); 1173 //System.out.println("rec = " + rec); 1174 for (List<Complex<edu.jas.application.RealAlgebraicNumber<D>>> cx : can) { 1175 Complex<edu.jas.application.RealAlgebraicNumber<D>> cax = cx.get(ix); 1176 //System.out.println("cax = " + toString(cax)); 1177 ComplexRing<edu.jas.application.RealAlgebraicNumber<D>> car = cax.ring; 1178 //System.out.println("car = " + car); 1179 GenPolynomialRing<Complex<edu.jas.application.RealAlgebraicNumber<D>>> pcrfac; 1180 pcrfac = new GenPolynomialRing<Complex<edu.jas.application.RealAlgebraicNumber<D>>>( 1181 car, rcfac); 1182 GenPolynomial<Complex<edu.jas.application.RealAlgebraicNumber<D>>> pcr; 1183 pcr = evaluateToComplexRealCoefficients(pcrfac, pip2cr, cax); 1184 //System.out.println("pcr = " + pcr); 1185 ComplexRoots<edu.jas.application.RealAlgebraicNumber<D>> rengine; 1186 rengine = new ComplexRootsSturm<edu.jas.application.RealAlgebraicNumber<D>>(car); 1187 long nr = 0; 1188 try { 1189 nr = rengine.complexRootCount(rec, pcr); 1190 //logger.info("rootCount = " + nr); 1191 } catch (InvalidBoundaryException e) { 1192 e.printStackTrace(); 1193 } 1194 if (nr == 1) { // one root 1195 logger.info(" hit, cxi = " + toString(cx.get(ix)) + ", cr = " + toString(cr)); 1196 List<Complex<edu.jas.application.RealAlgebraicNumber<D>>> cy; 1197 cy = new ArrayList<Complex<edu.jas.application.RealAlgebraicNumber<D>>>(); 1198 cy.addAll(cx); 1199 cy.add(cr); 1200 cn.add(cy); 1201 } else if (nr > 1) { 1202 logger.error("to many roots, cxi = " + toString(cx.get(ix)) + ", cr = " 1203 + toString(cr)); 1204 } else { // no root 1205 logger.info("no hit, cxi = " + toString(cx.get(ix)) + ", cr = " + toString(cr)); 1206 } 1207 } 1208 } 1209 } 1210 can = cn; 1211 } 1212 IdealWithComplexAlgebraicRoots<D> Ic = new IdealWithComplexAlgebraicRoots<D>(I, can); 1213 return Ic; 1214 } 1215 1216 1217 /** 1218 * String representation of a deximal approximation of a complex number. 1219 * @param c compelx number. 1220 * @return String representation of c 1221 */ 1222 public static <D extends GcdRingElem<D> & Rational> String toString( 1223 Complex<edu.jas.application.RealAlgebraicNumber<D>> c) { 1224 edu.jas.application.RealAlgebraicNumber<D> re = c.getRe(); 1225 edu.jas.application.RealAlgebraicNumber<D> im = c.getIm(); 1226 String s = re.decimalMagnitude().toString(); 1227 if (!im.isZERO()) { 1228 s = s + "i" + im.decimalMagnitude(); 1229 } 1230 return s; 1231 } 1232 1233 1234 /** 1235 * String representation of a deximal approximation of a complex number. 1236 * @param c compelx number. 1237 * @return String representation of c 1238 */ 1239 public static <D extends GcdRingElem<D> & Rational> String toString1(Complex<D> c) { 1240 D re = c.getRe(); 1241 D im = c.getIm(); 1242 String s = new BigDecimal(re.getRational()).toString(); 1243 if (!im.isZERO()) { 1244 s = s + "i" + new BigDecimal(im.getRational()); 1245 } 1246 return s; 1247 } 1248 1249 1250 /** 1251 * Construct complex roots for zero dimensional ideal(G). 1252 * @param I list of zero dimensional ideal with univariate irreducible 1253 * polynomials and bi-variate polynomials. 1254 * @return list of complex algebraic roots for ideal(G) 1255 */ 1256 public static <D extends GcdRingElem<D> & Rational> List<IdealWithComplexAlgebraicRoots<D>> complexAlgebraicRoots( 1257 List<IdealWithUniv<D>> I) { 1258 List<IdealWithComplexAlgebraicRoots<D>> lic = new ArrayList<IdealWithComplexAlgebraicRoots<D>>(); 1259 for (IdealWithUniv<D> iu : I) { 1260 IdealWithComplexAlgebraicRoots<D> iuc = PolyUtilApp.<D> complexAlgebraicRoots(iu); 1261 //System.out.println("iuc = " + iuc); 1262 lic.add(iuc); 1263 } 1264 return lic; 1265 } 1266 1267 1268 /** 1269 * Construct exact set of complex roots for zero dimensional ideal(G). 1270 * @param I zero dimensional ideal. 1271 * @return list of coordinates of complex roots for ideal(G) 1272 */ 1273 public static <D extends GcdRingElem<D> & Rational> List<IdealWithComplexAlgebraicRoots<D>> complexAlgebraicRoots( 1274 Ideal<D> I) { 1275 List<IdealWithUniv<D>> Ir = I.zeroDimRootDecomposition(); 1276 //System.out.println("Ir = " + Ir); 1277 List<IdealWithComplexAlgebraicRoots<D>> roots = PolyUtilApp.<D> complexAlgebraicRoots(Ir); 1278 return roots; 1279 } 1280 1281 1282 /* 1283 * Convert to a polynomial in given ring. 1284 * @param fac result polynomial ring. 1285 * @param p polynomial. 1286 * @return polynomial in ring fac <b>Note: </b> if p can not be represented 1287 * in fac then the results are unpredictable. 1288 */ 1289 static <C extends RingElem<C>> GenPolynomial<C> convert(GenPolynomialRing<C> fac, GenPolynomial<C> p) { 1290 if (fac.equals(p.factory())) { 1291 return p; 1292 } 1293 GenPolynomial<C> q = fac.parse(p.toString()); 1294 if (!q.toString().equals(p.toString())) { 1295 throw new RuntimeException("convert(" + p + ") = " + q); 1296 } 1297 return q; 1298 } 1299 1300 1301 /* 1302 * Convert to a polynomial in given ring. 1303 * @param fac result polynomial ring. 1304 * @param p polynomial. 1305 * @return polynomial in ring fac <b>Note: </b> if p can not be represented 1306 * in fac then the results are unpredictable. 1307 */ 1308 static <C extends RingElem<C>> GenPolynomial<Complex<C>> convertComplex(GenPolynomialRing<Complex<C>> fac, 1309 GenPolynomial<C> p) { 1310 GenPolynomial<Complex<C>> q = fac.parse(p.toString()); 1311 if (!q.toString().equals(p.toString())) { 1312 throw new RuntimeException("convert(" + p + ") = " + q); 1313 } 1314 return q; 1315 } 1316 1317 1318 /* 1319 * Convert to a polynomial in given ring. 1320 * @param fac result polynomial ring. 1321 * @param p complex polynomial. 1322 * @return polynomial in ring fac <b>Note: </b> if p can not be represented 1323 * in fac then the results are unpredictable. 1324 */ 1325 static <C extends RingElem<C>> GenPolynomial<Complex<C>> convertComplexComplex( 1326 GenPolynomialRing<Complex<C>> fac, GenPolynomial<Complex<C>> p) { 1327 if (fac.equals(p.factory())) { 1328 return p; 1329 } 1330 GenPolynomial<Complex<C>> q = fac.parse(p.toString()); 1331 if (!q.toString().equals(p.toString())) { 1332 throw new RuntimeException("convert(" + p + ") = " + q); 1333 } 1334 return q; 1335 } 1336 1337 1338 /** 1339 * Construct exact set of real roots for zero dimensional ideal(G). 1340 * @param I zero dimensional ideal. 1341 * @return list of coordinates of real roots for ideal(G) 1342 */ 1343 public static <D extends GcdRingElem<D> & Rational> List<IdealWithRealAlgebraicRoots<D>> realAlgebraicRoots( 1344 Ideal<D> I) { 1345 List<IdealWithUniv<D>> Ir = I.zeroDimRootDecomposition(); 1346 //System.out.println("Ir = " + Ir); 1347 List<IdealWithRealAlgebraicRoots<D>> roots = PolyUtilApp.<D> realAlgebraicRoots(Ir); 1348 return roots; 1349 } 1350 1351 1352 /** 1353 * Construct primitive element for double field extension. 1354 * @param a algebraic number ring with squarefree monic minimal polynomial 1355 * @param b algebraic number ring with squarefree monic minimal polynomial 1356 * @return primitive element container with algebraic number ring c, with 1357 * Q(c) = Q(a,b) 1358 */ 1359 public static <C extends GcdRingElem<C>> PrimitiveElement<C> primitiveElement(AlgebraicNumberRing<C> a, 1360 AlgebraicNumberRing<C> b) { 1361 GenPolynomial<C> ap = a.modul; 1362 GenPolynomial<C> bp = b.modul; 1363 1364 // setup bivariate polynomial ring 1365 String[] cv = new String[2]; 1366 cv[0] = ap.ring.getVars()[0]; 1367 cv[1] = bp.ring.getVars()[0]; 1368 TermOrder to = new TermOrder(TermOrder.INVLEX); 1369 GenPolynomialRing<C> cfac = new GenPolynomialRing<C>(ap.ring.coFac, 2, to, cv); 1370 GenPolynomial<C> as = ap.extendUnivariate(cfac, 0); 1371 GenPolynomial<C> bs = bp.extendUnivariate(cfac, 1); 1372 List<GenPolynomial<C>> L = new ArrayList<GenPolynomial<C>>(2); 1373 L.add(as); 1374 L.add(bs); 1375 List<GenPolynomial<C>> Op = new ArrayList<GenPolynomial<C>>(); 1376 1377 Ideal<C> id = new Ideal<C>(cfac, L); 1378 //System.out.println("id = " + id); 1379 IdealWithUniv<C> iu = id.normalPositionFor(0, 1, Op); 1380 //System.out.println("iu = " + iu); 1381 1382 // extract result polynomials 1383 List<GenPolynomial<C>> Np = iu.ideal.getList(); 1384 //System.out.println("Np = " + Np); 1385 as = PolyUtil.<C> selectWithVariable(Np, 1); 1386 bs = PolyUtil.<C> selectWithVariable(Np, 0); 1387 GenPolynomial<C> cs = PolyUtil.<C> selectWithVariable(Np, 2); 1388 //System.out.println("as = " + as); 1389 //System.out.println("bs = " + bs); 1390 //System.out.println("cs = " + cs); 1391 String[] ev = new String[] { cs.ring.getVars()[0] }; 1392 GenPolynomialRing<C> efac = new GenPolynomialRing<C>(ap.ring.coFac, 1, to, ev); 1393 //System.out.println("efac = " + efac); 1394 cs = cs.contractCoeff(efac); 1395 //System.out.println("cs = " + cs); 1396 as = as.reductum().contractCoeff(efac); 1397 as = as.negate(); 1398 //System.out.println("as = " + as); 1399 bs = bs.reductum().contractCoeff(efac); 1400 bs = bs.negate(); 1401 //System.out.println("bs = " + bs); 1402 AlgebraicNumberRing<C> c = new AlgebraicNumberRing<C>(cs); 1403 AlgebraicNumber<C> ab = new AlgebraicNumber<C>(c, as); 1404 AlgebraicNumber<C> bb = new AlgebraicNumber<C>(c, bs); 1405 PrimitiveElement<C> pe = new PrimitiveElement<C>(c, ab, bb, a, b); 1406 if (logger.isInfoEnabled()) { 1407 logger.info("primitive element = " + c); 1408 } 1409 return pe; 1410 } 1411 1412 1413 /** 1414 * Convert to primitive element ring. 1415 * @param cfac primitive element ring. 1416 * @param A algebraic number representing the generating element of a in the 1417 * new ring. 1418 * @param a algebraic number to convert. 1419 * @return a converted to the primitive element ring 1420 */ 1421 public static <C extends GcdRingElem<C>> AlgebraicNumber<C> convertToPrimitiveElem( 1422 AlgebraicNumberRing<C> cfac, AlgebraicNumber<C> A, AlgebraicNumber<C> a) { 1423 GenPolynomialRing<C> aufac = a.ring.ring; 1424 GenPolynomialRing<AlgebraicNumber<C>> ar = new GenPolynomialRing<AlgebraicNumber<C>>(cfac, aufac); 1425 GenPolynomial<AlgebraicNumber<C>> aps = PolyUtil.<C> convertToAlgebraicCoefficients(ar, a.val); 1426 AlgebraicNumber<C> ac = PolyUtil.<AlgebraicNumber<C>> evaluateMain(cfac, aps, A); 1427 return ac; 1428 } 1429 1430 1431 /** 1432 * Convert coefficients to primitive element ring. 1433 * @param cfac primitive element ring. 1434 * @param A algebraic number representing the generating element of a in the 1435 * new ring. 1436 * @param a polynomial with coefficients algebraic number to convert. 1437 * @return a with coefficients converted to the primitive element ring 1438 */ 1439 public static <C extends GcdRingElem<C>> GenPolynomial<AlgebraicNumber<C>> convertToPrimitiveElem( 1440 AlgebraicNumberRing<C> cfac, AlgebraicNumber<C> A, GenPolynomial<AlgebraicNumber<C>> a) { 1441 GenPolynomialRing<AlgebraicNumber<C>> cr = new GenPolynomialRing<AlgebraicNumber<C>>(cfac, a.ring); 1442 return PolyUtil.<AlgebraicNumber<C>, AlgebraicNumber<C>> map(cr, a, new CoeffConvertAlg<C>(cfac, A)); 1443 } 1444 1445 1446 /** 1447 * Convert to primitive element ring. 1448 * @param cfac primitive element ring. 1449 * @param A algebraic number representing the generating element of a in the 1450 * new ring. 1451 * @param a recursive algebraic number to convert. 1452 * @return a converted to the primitive element ring 1453 */ 1454 public static <C extends GcdRingElem<C>> AlgebraicNumber<C> convertToPrimitiveElem( 1455 AlgebraicNumberRing<C> cfac, AlgebraicNumber<C> A, AlgebraicNumber<C> B, 1456 AlgebraicNumber<AlgebraicNumber<C>> a) { 1457 GenPolynomial<AlgebraicNumber<C>> aps = PolyUtilApp.<C> convertToPrimitiveElem(cfac, A, a.val); 1458 AlgebraicNumber<C> ac = PolyUtil.<AlgebraicNumber<C>> evaluateMain(cfac, aps, B); 1459 return ac; 1460 } 1461 1462 1463 /** 1464 * Construct primitive element for double field extension. 1465 * @param b algebraic number ring with squarefree monic minimal polynomial 1466 * over Q(a) 1467 * @return primitive element container with algebraic number ring c, with 1468 * Q(c) = Q(a)(b) 1469 */ 1470 public static <C extends GcdRingElem<C>> PrimitiveElement<C> primitiveElement( 1471 AlgebraicNumberRing<AlgebraicNumber<C>> b) { 1472 GenPolynomial<AlgebraicNumber<C>> bp = b.modul; 1473 AlgebraicNumberRing<C> a = (AlgebraicNumberRing<C>) b.ring.coFac; 1474 GenPolynomial<C> ap = a.modul; 1475 1476 // setup bivariate polynomial ring 1477 String[] cv = new String[2]; 1478 cv[0] = ap.ring.getVars()[0]; 1479 cv[1] = bp.ring.getVars()[0]; 1480 TermOrder to = new TermOrder(TermOrder.INVLEX); 1481 GenPolynomialRing<C> cfac = new GenPolynomialRing<C>(ap.ring.coFac, 2, to, cv); 1482 GenPolynomialRing<GenPolynomial<C>> rfac = new GenPolynomialRing<GenPolynomial<C>>(a.ring, 1, 1483 bp.ring.getVars()); 1484 GenPolynomial<C> as = ap.extendUnivariate(cfac, 0); 1485 GenPolynomial<GenPolynomial<C>> bss = PolyUtil.<C> fromAlgebraicCoefficients(rfac, bp); 1486 GenPolynomial<C> bs = PolyUtil.<C> distribute(cfac, bss); 1487 List<GenPolynomial<C>> L = new ArrayList<GenPolynomial<C>>(2); 1488 L.add(as); 1489 L.add(bs); 1490 List<GenPolynomial<C>> Op = new ArrayList<GenPolynomial<C>>(); 1491 1492 Ideal<C> id = new Ideal<C>(cfac, L); 1493 //System.out.println("id = " + id); 1494 IdealWithUniv<C> iu = id.normalPositionFor(0, 1, Op); 1495 //System.out.println("iu = " + iu); 1496 1497 // extract result polynomials 1498 List<GenPolynomial<C>> Np = iu.ideal.getList(); 1499 as = PolyUtil.<C> selectWithVariable(Np, 1); 1500 bs = PolyUtil.<C> selectWithVariable(Np, 0); 1501 GenPolynomial<C> cs = PolyUtil.<C> selectWithVariable(Np, 2); 1502 //System.out.println("as = " + as); 1503 //System.out.println("bs = " + bs); 1504 //System.out.println("cs = " + cs); 1505 String[] ev = new String[] { cs.ring.getVars()[0] }; 1506 GenPolynomialRing<C> efac = new GenPolynomialRing<C>(ap.ring.coFac, 1, to, ev); 1507 // System.out.println("efac = " + efac); 1508 cs = cs.contractCoeff(efac); 1509 // System.out.println("cs = " + cs); 1510 as = as.reductum().contractCoeff(efac); 1511 as = as.negate(); 1512 // System.out.println("as = " + as); 1513 bs = bs.reductum().contractCoeff(efac); 1514 bs = bs.negate(); 1515 //System.out.println("bs = " + bs); 1516 AlgebraicNumberRing<C> c = new AlgebraicNumberRing<C>(cs); 1517 AlgebraicNumber<C> ab = new AlgebraicNumber<C>(c, as); 1518 AlgebraicNumber<C> bb = new AlgebraicNumber<C>(c, bs); 1519 PrimitiveElement<C> pe = new PrimitiveElement<C>(c, ab, bb); // missing ,a,b); 1520 if (logger.isInfoEnabled()) { 1521 logger.info("primitive element = " + pe); 1522 } 1523 return pe; 1524 } 1525 1526 1527 /** 1528 * Convert to primitive element ring. 1529 * @param cfac primitive element ring. 1530 * @param A algebraic number representing the generating element of a in the 1531 * new ring. 1532 * @param a polynomial with recursive algebraic number coefficients to 1533 * convert. 1534 * @return a converted to the primitive element ring 1535 */ 1536 public static <C extends GcdRingElem<C>> GenPolynomial<AlgebraicNumber<C>> convertToPrimitiveElem( 1537 AlgebraicNumberRing<C> cfac, AlgebraicNumber<C> A, AlgebraicNumber<C> B, 1538 GenPolynomial<AlgebraicNumber<AlgebraicNumber<C>>> a) { 1539 GenPolynomialRing<AlgebraicNumber<C>> cr = new GenPolynomialRing<AlgebraicNumber<C>>(cfac, a.ring); 1540 return PolyUtil.<AlgebraicNumber<AlgebraicNumber<C>>, AlgebraicNumber<C>> map(cr, a, 1541 new CoeffRecConvertAlg<C>(cfac, A, B)); 1542 } 1543 1544 1545 /** 1546 * Convert to RealAlgebraicNumber coefficients. Represent as polynomial with 1547 * RealAlgebraicNumber<C> coefficients from package 1548 * <code>edu.jas.root.</code> 1549 * @param afac result polynomial factory. 1550 * @param A polynomial with RealAlgebraicNumber<C> coefficients to be 1551 * converted. 1552 * @return polynomial with RealAlgebraicNumber<C> coefficients. 1553 */ 1554 public static <C extends GcdRingElem<C> & Rational> GenPolynomial<edu.jas.root.RealAlgebraicNumber<C>> realAlgFromRealCoefficients( 1555 GenPolynomialRing<edu.jas.root.RealAlgebraicNumber<C>> afac, 1556 GenPolynomial<edu.jas.application.RealAlgebraicNumber<C>> A) { 1557 edu.jas.root.RealAlgebraicRing<C> cfac = (edu.jas.root.RealAlgebraicRing<C>) afac.coFac; 1558 return PolyUtil.<edu.jas.application.RealAlgebraicNumber<C>, edu.jas.root.RealAlgebraicNumber<C>> map( 1559 afac, A, new ReAlgFromRealCoeff<C>(cfac)); 1560 } 1561 1562 1563 /** 1564 * Convert to RealAlgebraicNumber coefficients. Represent as polynomial with 1565 * RealAlgebraicNumber<C> coefficients from package 1566 * 1567 * <pre> 1568 * edu.jas.application 1569 * </pre> 1570 * 1571 * . 1572 * @param rfac result polynomial factory. 1573 * @param A polynomial with RealAlgebraicNumber<C> coefficients to be 1574 * converted. 1575 * @return polynomial with RealAlgebraicNumber<C> coefficients. 1576 */ 1577 public static <C extends GcdRingElem<C> & Rational> GenPolynomial<edu.jas.application.RealAlgebraicNumber<C>> realFromRealAlgCoefficients( 1578 GenPolynomialRing<edu.jas.application.RealAlgebraicNumber<C>> rfac, 1579 GenPolynomial<edu.jas.root.RealAlgebraicNumber<C>> A) { 1580 edu.jas.application.RealAlgebraicRing<C> cfac = (edu.jas.application.RealAlgebraicRing<C>) rfac.coFac; 1581 return PolyUtil.<edu.jas.root.RealAlgebraicNumber<C>, edu.jas.application.RealAlgebraicNumber<C>> map( 1582 rfac, A, new RealFromReAlgCoeff<C>(cfac)); 1583 } 1584 1585 1586 /** 1587 * Convert to Complex<RealAlgebraicNumber> coefficients. Represent as 1588 * polynomial with Complex<RealAlgebraicNumber> coefficients, C is 1589 * e.g. BigRational. 1590 * @param pfac result polynomial factory. 1591 * @param A polynomial with Complex coefficients to be converted. 1592 * @return polynomial with Complex<RealAlgebraicNumber> coefficients. 1593 */ 1594 public static <C extends GcdRingElem<C> & Rational> GenPolynomial<Complex<edu.jas.application.RealAlgebraicNumber<C>>> convertToComplexRealCoefficients( 1595 GenPolynomialRing<Complex<edu.jas.application.RealAlgebraicNumber<C>>> pfac, 1596 GenPolynomial<Complex<C>> A) { 1597 ComplexRing<edu.jas.application.RealAlgebraicNumber<C>> afac; 1598 afac = (ComplexRing<edu.jas.application.RealAlgebraicNumber<C>>) pfac.coFac; 1599 return PolyUtil.<Complex<C>, Complex<edu.jas.application.RealAlgebraicNumber<C>>> map(pfac, A, 1600 new CoeffToComplexReal<C>(afac)); 1601 } 1602 1603 1604 /** 1605 * Evaluate to Complex<RealAlgebraicNumber> coefficients. Represent as 1606 * polynomial with Complex<RealAlgebraicNumber> coefficients, C is 1607 * e.g. BigRational. 1608 * @param pfac result polynomial factory. 1609 * @param A = A(x,Y) a recursive polynomial with 1610 * GenPolynomial<Complex> coefficients to be converted. 1611 * @param r Complex<RealAlgebraicNumber> to be evaluated at. 1612 * @return A(r,Y), a polynomial with Complex<RealAlgebraicNumber> 1613 * coefficients. 1614 */ 1615 public static <C extends GcdRingElem<C> & Rational> GenPolynomial<Complex<edu.jas.application.RealAlgebraicNumber<C>>> evaluateToComplexRealCoefficients( 1616 GenPolynomialRing<Complex<edu.jas.application.RealAlgebraicNumber<C>>> pfac, 1617 GenPolynomial<GenPolynomial<Complex<C>>> A, 1618 Complex<edu.jas.application.RealAlgebraicNumber<C>> r) { 1619 return PolyUtil.<GenPolynomial<Complex<C>>, Complex<edu.jas.application.RealAlgebraicNumber<C>>> map( 1620 pfac, A, new EvaluateToComplexReal<C>(pfac, r)); 1621 } 1622 1623 1624} 1625 1626 1627/** 1628 * Coefficient to convert algebriac functor. 1629 */ 1630class CoeffConvertAlg<C extends GcdRingElem<C>> 1631 implements UnaryFunctor<AlgebraicNumber<C>, AlgebraicNumber<C>> { 1632 1633 1634 final protected AlgebraicNumberRing<C> afac; 1635 1636 1637 final protected AlgebraicNumber<C> A; 1638 1639 1640 public CoeffConvertAlg(AlgebraicNumberRing<C> fac, AlgebraicNumber<C> a) { 1641 if (fac == null || a == null) { 1642 throw new IllegalArgumentException("fac and a must not be null"); 1643 } 1644 afac = fac; 1645 A = a; 1646 } 1647 1648 1649 public AlgebraicNumber<C> eval(AlgebraicNumber<C> c) { 1650 if (c == null) { 1651 return afac.getZERO(); 1652 } 1653 return PolyUtilApp.<C> convertToPrimitiveElem(afac, A, c); 1654 } 1655} 1656 1657 1658/** 1659 * Coefficient recursive to convert algebriac functor. 1660 */ 1661class CoeffRecConvertAlg<C extends GcdRingElem<C>> 1662 implements UnaryFunctor<AlgebraicNumber<AlgebraicNumber<C>>, AlgebraicNumber<C>> { 1663 1664 1665 final protected AlgebraicNumberRing<C> afac; 1666 1667 1668 final protected AlgebraicNumber<C> A; 1669 1670 1671 final protected AlgebraicNumber<C> B; 1672 1673 1674 public CoeffRecConvertAlg(AlgebraicNumberRing<C> fac, AlgebraicNumber<C> a, AlgebraicNumber<C> b) { 1675 if (fac == null || a == null || b == null) { 1676 throw new IllegalArgumentException("fac, a and b must not be null"); 1677 } 1678 afac = fac; 1679 A = a; 1680 B = b; 1681 } 1682 1683 1684 public AlgebraicNumber<C> eval(AlgebraicNumber<AlgebraicNumber<C>> c) { 1685 if (c == null) { 1686 return afac.getZERO(); 1687 } 1688 return PolyUtilApp.<C> convertToPrimitiveElem(afac, A, B, c); 1689 } 1690} 1691 1692 1693/** 1694 * Coefficient to real algebriac from real algebraic functor. 1695 */ 1696class ReAlgFromRealCoeff<C extends GcdRingElem<C> & Rational> implements 1697 UnaryFunctor<edu.jas.application.RealAlgebraicNumber<C>, edu.jas.root.RealAlgebraicNumber<C>> { 1698 1699 1700 final protected edu.jas.root.RealAlgebraicRing<C> afac; 1701 1702 1703 public ReAlgFromRealCoeff(edu.jas.root.RealAlgebraicRing<C> fac) { 1704 if (fac == null) { 1705 throw new IllegalArgumentException("fac must not be null"); 1706 } 1707 afac = fac; 1708 } 1709 1710 1711 @SuppressWarnings("unchecked") 1712 public edu.jas.root.RealAlgebraicNumber<C> eval(edu.jas.application.RealAlgebraicNumber<C> c) { 1713 if (c == null) { 1714 return afac.getZERO(); 1715 } 1716 return (edu.jas.root.RealAlgebraicNumber<C>) (Object) c.number; // force ignore recursion 1717 } 1718} 1719 1720 1721/** 1722 * Coefficient to real algebriac from algebraic functor. 1723 */ 1724class RealFromReAlgCoeff<C extends GcdRingElem<C> & Rational> implements 1725 UnaryFunctor<edu.jas.root.RealAlgebraicNumber<C>, edu.jas.application.RealAlgebraicNumber<C>> { 1726 1727 1728 final protected edu.jas.application.RealAlgebraicRing<C> rfac; 1729 1730 1731 public RealFromReAlgCoeff(edu.jas.application.RealAlgebraicRing<C> fac) { 1732 if (fac == null) { 1733 throw new IllegalArgumentException("fac must not be null"); 1734 } 1735 rfac = fac; 1736 } 1737 1738 1739 @SuppressWarnings("unchecked") 1740 public edu.jas.application.RealAlgebraicNumber<C> eval(edu.jas.root.RealAlgebraicNumber<C> c) { 1741 if (c == null) { 1742 return rfac.getZERO(); 1743 } 1744 edu.jas.root.RealAlgebraicNumber<edu.jas.root.RealAlgebraicNumber<C>> rrc = (edu.jas.root.RealAlgebraicNumber<edu.jas.root.RealAlgebraicNumber<C>>) (Object) c; // force resurrect recursion 1745 return new edu.jas.application.RealAlgebraicNumber<C>(rfac, rrc); 1746 } 1747} 1748 1749 1750/** 1751 * Coefficient to complex real algebriac functor. 1752 */ 1753class CoeffToComplexReal<C extends GcdRingElem<C> & Rational> 1754 implements UnaryFunctor<Complex<C>, Complex<edu.jas.application.RealAlgebraicNumber<C>>> { 1755 1756 1757 final protected ComplexRing<edu.jas.application.RealAlgebraicNumber<C>> cfac; 1758 1759 1760 final edu.jas.application.RealAlgebraicRing<C> afac; 1761 1762 1763 final GenPolynomialRing<C> pfac; 1764 1765 1766 public CoeffToComplexReal(ComplexRing<edu.jas.application.RealAlgebraicNumber<C>> fac) { 1767 if (fac == null) { 1768 throw new IllegalArgumentException("fac must not be null"); 1769 } 1770 cfac = fac; 1771 afac = (edu.jas.application.RealAlgebraicRing<C>) cfac.ring; 1772 pfac = afac.univs.ideal.getRing(); 1773 } 1774 1775 1776 public Complex<edu.jas.application.RealAlgebraicNumber<C>> eval(Complex<C> c) { 1777 if (c == null) { 1778 return cfac.getZERO(); 1779 } 1780 GenPolynomial<C> pr, pi; 1781 pr = new GenPolynomial<C>(pfac, c.getRe()); 1782 pi = new GenPolynomial<C>(pfac, c.getIm()); 1783 //System.out.println("pr = " + pr); 1784 //System.out.println("pi = " + pi); 1785 edu.jas.application.RealAlgebraicNumber<C> re, im; 1786 re = new edu.jas.application.RealAlgebraicNumber<C>(afac, pr); 1787 im = new edu.jas.application.RealAlgebraicNumber<C>(afac, pi); 1788 //System.out.println("re = " + re); 1789 //System.out.println("im = " + im); 1790 return new Complex<edu.jas.application.RealAlgebraicNumber<C>>(cfac, re, im); 1791 } 1792} 1793 1794 1795/** 1796 * Polynomial coefficient to complex real algebriac evaluation functor. 1797 */ 1798class EvaluateToComplexReal<C extends GcdRingElem<C> & Rational> implements 1799 UnaryFunctor<GenPolynomial<Complex<C>>, Complex<edu.jas.application.RealAlgebraicNumber<C>>> { 1800 1801 1802 final protected GenPolynomialRing<Complex<edu.jas.application.RealAlgebraicNumber<C>>> pfac; 1803 1804 1805 final protected ComplexRing<edu.jas.application.RealAlgebraicNumber<C>> cfac; 1806 1807 1808 final protected Complex<edu.jas.application.RealAlgebraicNumber<C>> root; 1809 1810 1811 public EvaluateToComplexReal(GenPolynomialRing<Complex<edu.jas.application.RealAlgebraicNumber<C>>> fac, 1812 Complex<edu.jas.application.RealAlgebraicNumber<C>> r) { 1813 if (fac == null) { 1814 throw new IllegalArgumentException("fac must not be null"); 1815 } 1816 if (r == null) { 1817 throw new IllegalArgumentException("r must not be null"); 1818 } 1819 pfac = fac; 1820 cfac = (ComplexRing<edu.jas.application.RealAlgebraicNumber<C>>) fac.coFac; 1821 root = r; 1822 //System.out.println("cfac = " + cfac); 1823 //System.out.println("root = " + root); 1824 } 1825 1826 1827 public Complex<edu.jas.application.RealAlgebraicNumber<C>> eval(GenPolynomial<Complex<C>> c) { 1828 if (c == null) { 1829 return cfac.getZERO(); 1830 } 1831 //System.out.println("c = " + c); 1832 GenPolynomial<Complex<edu.jas.application.RealAlgebraicNumber<C>>> cp; 1833 cp = PolyUtilApp.<C> convertToComplexRealCoefficients(pfac, c); 1834 Complex<edu.jas.application.RealAlgebraicNumber<C>> cr; 1835 cr = PolyUtil.<Complex<edu.jas.application.RealAlgebraicNumber<C>>> evaluateMain(cfac, cp, root); 1836 return cr; 1837 } 1838}