001 /*
002 * $Id: RootFactory.java 3665 2011-06-13 11:08:16Z kredel $
003 */
004
005 package edu.jas.application;
006
007
008 import java.util.ArrayList;
009 import java.util.List;
010 import java.util.Set;
011 import java.util.Map;
012
013 import org.apache.log4j.Logger;
014
015 import edu.jas.arith.Rational;
016 import edu.jas.poly.Complex;
017 import edu.jas.poly.ComplexRing;
018 import edu.jas.poly.GenPolynomial;
019 import edu.jas.poly.GenPolynomialRing;
020 import edu.jas.poly.PolyUtil;
021 import edu.jas.poly.TermOrder;
022 import edu.jas.root.ComplexRoots;
023 import edu.jas.root.ComplexRootsSturm;
024 import edu.jas.root.RealRootTuple;
025 import edu.jas.root.Interval;
026 import edu.jas.structure.GcdRingElem;
027 import edu.jas.structure.RingFactory;
028 import edu.jas.ufd.SquarefreeAbstract;
029 import edu.jas.ufd.SquarefreeFactory;
030
031
032 /**
033 * Roots factory.
034 * @author Heinz Kredel
035 */
036 public class RootFactory {
037
038
039 private static final Logger logger = Logger.getLogger(RootFactory.class);
040
041
042 //private static boolean debug = logger.isDebugEnabled();
043
044
045 /**
046 * Is complex algebraic number a root of a polynomial.
047 * @param f univariate polynomial.
048 * @param r complex algebraic number.
049 * @return true, if f(r) == 0, else false;
050 */
051 public static <C extends GcdRingElem<C> & Rational>
052 boolean isRoot(GenPolynomial<Complex<C>> f, Complex<RealAlgebraicNumber<C>> r) {
053 ComplexRing<RealAlgebraicNumber<C>> cr = r.factory();
054 GenPolynomialRing<Complex<RealAlgebraicNumber<C>>> cfac
055 = new GenPolynomialRing<Complex<RealAlgebraicNumber<C>>>(cr,f.factory());
056 GenPolynomial<Complex<RealAlgebraicNumber<C>>> p;
057 p = PolyUtilApp.<C> convertToComplexRealCoefficients(cfac,f);
058 //System.out.println("p = " + p);
059 // test algebraic part
060 Complex<RealAlgebraicNumber<C>> a = PolyUtil.<Complex<RealAlgebraicNumber<C>>> evaluateMain(cr,p,r);
061 boolean t = a.isZERO();
062 if ( !t ) {
063 logger.info("p(r) = " + a + ", p = " + f + ", r = " + r);
064 return t;
065 }
066 // test real part
067 RealAlgebraicRing<C> rring = (RealAlgebraicRing<C>)cr.ring;
068 RealRootTuple<C> rroot = rring.getRoot();
069 List<edu.jas.root.RealAlgebraicNumber<C>> rlist = rroot.tuple;
070 //System.out.println("rlist = " + rlist);
071 Interval<C> vr = rlist.get(0).ring.getRoot();
072 Interval<C> vi = rlist.get(1).ring.getRoot();
073 ComplexRing<C> ccfac = new ComplexRing<C>((RingFactory<C>)vr.left.factory());
074 Complex<C> sw = new Complex<C>(ccfac,vr.left,vi.left);
075 Complex<C> ne = new Complex<C>(ccfac,vr.right,vi.right);
076 Complex<C> epsw = PolyUtil.<Complex<C>> evaluateMain(ccfac, f, sw);
077 Complex<C> epne = PolyUtil.<Complex<C>> evaluateMain(ccfac, f, ne);
078 int rootre = (epsw.getRe().signum()*epne.getRe().signum());
079 int rootim = (epsw.getIm().signum()*epne.getIm().signum());
080 t = (rootre <= 0 && rootim <= 0);
081 if ( !t ) {
082 logger.debug("vr = " + vr + ", vi = " + vi);
083 logger.info("sw = " + sw + ", ne = " + ne);
084 //System.out.println("root(re) = " + rootre + ", root(im) = " + rootim);
085 logger.info("p(root): p = " + f + ", epsw = " + epsw + ", epne = " + epne);
086 return t;
087 }
088 //System.out.println("r = " + r.getRe().magnitude() + " i " + r.getIm().magnitude());
089 return t;
090 }
091
092
093 /**
094 * Is complex algebraic number a root of a polynomial.
095 * @param f univariate polynomial.
096 * @param R list of complex algebraic numbers.
097 * @return true, if f(r) == 0 for all r in R, else false;
098 */
099 public static <C extends GcdRingElem<C> & Rational>
100 boolean isRoot(GenPolynomial<Complex<C>> f, List<Complex<RealAlgebraicNumber<C>>> R) {
101 for ( Complex<RealAlgebraicNumber<C>> r : R ) {
102 boolean t = isRoot(f,r);
103 if ( !t ) {
104 return false;
105 }
106 }
107 return true;
108 }
109
110
111 /**
112 * Complex algebraic number roots.
113 * @param f univariate polynomial.
114 * @return a list of different complex algebraic numbers, with f(c) == 0 for c in roots.
115 */
116 public static <C extends GcdRingElem<C> & Rational>
117 List<Complex<RealAlgebraicNumber<C>>> complexAlgebraicNumbersComplex(GenPolynomial<Complex<C>> f) {
118 GenPolynomialRing<Complex<C>> pfac = f.factory();
119 if (pfac.nvar != 1) {
120 throw new IllegalArgumentException("only for univariate polynomials");
121 }
122 ComplexRing<C> cfac = (ComplexRing<C>) pfac.coFac;
123 ComplexRoots<C> cr = new ComplexRootsSturm<C>(cfac);
124 SquarefreeAbstract<Complex<C>> engine = SquarefreeFactory.<Complex<C>> getImplementation(cfac);
125 Map<GenPolynomial<Complex<C>>,Long> F = engine.squarefreeFactors(f.monic());
126 Set<GenPolynomial<Complex<C>>> S = F.keySet();
127 //System.out.println("S = " + S);
128 List<Complex<RealAlgebraicNumber<C>>> list = new ArrayList<Complex<RealAlgebraicNumber<C>>>();
129 for (GenPolynomial<Complex<C>> sp : S) {
130 if (sp.isConstant() || sp.isZERO()) {
131 continue;
132 }
133 List<Complex<RealAlgebraicNumber<C>>> ls = RootFactory.<C> complexAlgebraicNumbersSquarefree(sp);
134 long m = F.get(sp);
135 for (long i = 0L; i < m; i++) {
136 list.addAll(ls);
137 }
138 }
139 return list;
140 }
141
142
143 /**
144 * Complex algebraic number roots.
145 * @param f univariate squarefree polynomial.
146 * @return a list of different complex algebraic numbers, with f(c) == 0 for c in roots.
147 */
148 public static <C extends GcdRingElem<C> & Rational>
149 List<Complex<RealAlgebraicNumber<C>>> complexAlgebraicNumbersSquarefree(GenPolynomial<Complex<C>> f) {
150 GenPolynomialRing<Complex<C>> pfac = f.factory();
151 if (pfac.nvar != 1) {
152 throw new IllegalArgumentException("only for univariate polynomials");
153 }
154 ComplexRing<C> cfac = (ComplexRing<C>) pfac.coFac;
155 ComplexRoots<C> cr = new ComplexRootsSturm<C>(cfac);
156 TermOrder to = new TermOrder(TermOrder.INVLEX);
157 GenPolynomialRing<Complex<C>> tfac = new GenPolynomialRing<Complex<C>>(cfac, 2, to); //,vars); //tord?
158 //System.out.println("tfac = " + tfac);
159 GenPolynomial<Complex<C>> t = tfac.univariate(1, 1L).sum(
160 tfac.univariate(0, 1L).multiply(cfac.getIMAG()));
161 //System.out.println("t = " + t); // t = x + i y
162
163 GenPolynomialRing<C> rfac = new GenPolynomialRing<C>(cfac.ring, tfac); //tord?
164 //System.out.println("rfac = " + rfac);
165
166 List<Complex<RealAlgebraicNumber<C>>> list = new ArrayList<Complex<RealAlgebraicNumber<C>>>();
167 GenPolynomial<Complex<C>> sp = f;
168 if (sp.isConstant() || sp.isZERO()) {
169 return list;
170 }
171 GenPolynomial<Complex<C>> su = PolyUtil.<Complex<C>> substituteUnivariate(sp, t);
172 //System.out.println("su = " + su);
173 su = su.monic();
174 //System.out.println("su = " + su);
175 GenPolynomial<C> re = PolyUtil.<C> realPartFromComplex(rfac, su);
176 GenPolynomial<C> im = PolyUtil.<C> imaginaryPartFromComplex(rfac, su);
177 if ( logger.isInfoEnabled() ) {
178 logger.debug("rfac = " + rfac.toScript());
179 logger.info("t = " + t);
180 logger.info("re = " + re.toScript());
181 logger.info("im = " + im.toScript());
182 }
183 List<GenPolynomial<C>> li = new ArrayList<GenPolynomial<C>>(2);
184 li.add(re);
185 li.add(im);
186 Ideal<C> id = new Ideal<C>(rfac, li);
187 //System.out.println("id = " + id);
188
189 List<IdealWithUniv<C>> idul = id.zeroDimRootDecomposition();
190
191 IdealWithRealAlgebraicRoots<C, C> idr;
192 for (IdealWithUniv<C> idu : idul) {
193 //System.out.println("---idu = " + idu);
194 idr = PolyUtilApp.<C, C> realAlgebraicRoots(idu);
195 //System.out.println("---idr = " + idr);
196 for (List<edu.jas.root.RealAlgebraicNumber<C>> crr : idr.ran) {
197 //System.out.println("crr = " + crr);
198 RealRootTuple<C> root = new RealRootTuple<C>(crr);
199 //System.out.println("root = " + root);
200 RealAlgebraicRing<C> car = new RealAlgebraicRing<C>(idu, root);
201 //System.out.println("car = " + car);
202 List<RealAlgebraicNumber<C>> gens = car.generators();
203 //System.out.println("gens = " + gens);
204 int sg = gens.size();
205 RealAlgebraicNumber<C> rre = gens.get(sg-2);
206 RealAlgebraicNumber<C> rim = gens.get(sg-1);
207 ComplexRing<RealAlgebraicNumber<C>> cring = new ComplexRing<RealAlgebraicNumber<C>>(car);
208 Complex<RealAlgebraicNumber<C>> crn = new Complex<RealAlgebraicNumber<C>>(cring,rre,rim);
209 //System.out.println("crn = " + crn.toScript());
210
211 boolean it;
212 int count = 0;
213 do { // refine intervals if necessary
214 Interval<C> vr = crr.get(0).ring.getRoot();
215 Interval<C> vi = crr.get(1).ring.getRoot();
216 ComplexRing<C> ccfac = new ComplexRing<C>((RingFactory<C>)vr.left.factory());
217 Complex<C> sw = new Complex<C>(ccfac,vr.left,vi.left);
218 Complex<C> ne = new Complex<C>(ccfac,vr.right,vi.right);
219 Complex<C> epsw = PolyUtil.<Complex<C>> evaluateMain(ccfac, f, sw);
220 Complex<C> epne = PolyUtil.<Complex<C>> evaluateMain(ccfac, f, ne);
221 int rootre = (epsw.getRe().signum()*epne.getRe().signum());
222 int rootim = (epsw.getIm().signum()*epne.getIm().signum());
223 it = (rootre <= 0 && rootim <= 0);
224 if ( !it ) {
225 logger.info("refine intervals: vr = " + vr + ", vi = " + vi);
226 //System.out.println("crn = " + crn.getRe().magnitude() + " i " + crn.getIm().magnitude());
227 //System.out.println("sw = " + sw + ", ne = " + ne);
228 //System.out.println("epsw = " + epsw + ", epne = " + epne);
229 //System.out.println("root(re) = " + rootre + ", root(im) = " + rootim);
230 crn.getRe().ring.realRing.halfInterval();
231 //System.out.println("crn.re = " + crn.getRe().ring.realRing);
232 crn.getIm().ring.realRing.halfInterval();
233 //System.out.println("crn.im = " + crn.getIm().ring.realRing);
234 //edu.jas.root.RealAlgebraicRing<C> rrr;
235 //rrr = (edu.jas.root.RealAlgebraicRing<C>) crn.getRe().ring.realRing.algebraic.ring.coFac;
236 //rrr.halfInterval();
237 //System.out.println("rrr.re = " + rrr);
238 //rrr = (edu.jas.root.RealAlgebraicRing<C>) crn.getIm().ring.realRing.algebraic.ring.coFac;
239 //rrr.halfInterval();
240 //System.out.println("rrr.im = " + rrr);
241 if ( count++ > 2 ) {
242 //throw new RuntimeException("no roots of " + f);
243 logger.info("break in root refinement of " + f);
244 it = true;
245 }
246 }
247 } while ( !it );
248 list.add(crn);
249 }
250 }
251 return list;
252 }
253
254
255 /* todo
256 * Complex algebraic numbers.
257 * @param f univariate polynomial.
258 * @param eps rational precision.
259 * @return a list of different complex algebraic numbers.
260 public static <C extends GcdRingElem<C> & Rational>
261 List<ComplexAlgebraicNumber<C>> complexAlgebraicNumbersComplex(GenPolynomial<Complex<C>> f, BigRational eps) {
262 }
263 */
264
265 }