001 /*
002 * $Id: ComplexRootsSturm.java 3613 2011-04-28 08:45:32Z kredel $
003 */
004
005 package edu.jas.root;
006
007
008 import java.util.ArrayList;
009 //import java.util.Arrays;
010 import java.util.List;
011
012 import org.apache.log4j.Logger;
013
014 import edu.jas.arith.Rational;
015 import edu.jas.arith.BigRational;
016 import edu.jas.poly.Complex;
017 import edu.jas.poly.ComplexRing;
018 import edu.jas.poly.GenPolynomial;
019 import edu.jas.poly.PolyUtil;
020 import edu.jas.structure.RingElem;
021 import edu.jas.structure.RingFactory;
022 import edu.jas.util.ArrayUtil;
023
024
025 /**
026 * Complex roots implemented by Sturm sequences. Algorithms use exact method
027 * derived from Wilf's numeric Routh-Hurwitz method.
028 * @param <C> coefficient type.
029 * @author Heinz Kredel
030 */
031 public class ComplexRootsSturm<C extends RingElem<C> & Rational> extends ComplexRootsAbstract<C> {
032
033
034 private static final Logger logger = Logger.getLogger(ComplexRootsSturm.class);
035
036
037 private final boolean debug = logger.isDebugEnabled();
038
039
040 /**
041 * Constructor.
042 * @param cf coefficient factory.
043 */
044 public ComplexRootsSturm(RingFactory<Complex<C>> cf) {
045 super(cf);
046 //ufd = GCDFactory.<Complex<C>> getImplementation(cf);
047 }
048
049
050 /**
051 * Cauchy index of rational function f/g on interval.
052 * @param a interval bound for I = [a,b].
053 * @param b interval bound for I = [a,b].
054 * @param f univariate polynomial.
055 * @param g univariate polynomial.
056 * @return winding number of f/g in I.
057 */
058 public long indexOfCauchy(C a, C b, GenPolynomial<C> f, GenPolynomial<C> g) {
059 List<GenPolynomial<C>> S = sturmSequence(g, f);
060 //System.out.println("S = " + S);
061 if (debug) {
062 logger.info("sturmSeq = " + S);
063 }
064 RingFactory<C> cfac = f.ring.coFac;
065 List<C> l = PolyUtil.<C> evaluateMain(cfac, S, a);
066 List<C> r = PolyUtil.<C> evaluateMain(cfac, S, b);
067 long v = RootUtil.<C> signVar(l) - RootUtil.<C> signVar(r);
068 //System.out.println("v = " + v);
069 // if (v < 0L) {
070 // v = -v;
071 // }
072 return v;
073 }
074
075
076 /**
077 * Routh index of complex function f + i g on interval.
078 * @param a interval bound for I = [a,b].
079 * @param b interval bound for I = [a,b].
080 * @param f univariate polynomial.
081 * @param g univariate polynomial != 0.
082 * @return index number of f + i g.
083 */
084 public long[] indexOfRouth(C a, C b, GenPolynomial<C> f, GenPolynomial<C> g) {
085 List<GenPolynomial<C>> S = sturmSequence(f, g);
086 //System.out.println("S = " + S);
087 RingFactory<C> cfac = f.ring.coFac;
088 List<C> l = PolyUtil.<C> evaluateMain(cfac, S, a);
089 List<C> r = PolyUtil.<C> evaluateMain(cfac, S, b);
090 long v = RootUtil.<C> signVar(l) - RootUtil.<C> signVar(r);
091 //System.out.println("v = " + v);
092
093 long d = f.degree(0);
094 if (d < g.degree(0)) {
095 d = g.degree(0);
096 }
097 //System.out.println("d = " + d);
098 long ui = (d - v) / 2;
099 long li = (d + v) / 2;
100 //System.out.println("upper = " + ui);
101 //System.out.println("lower = " + li);
102 return new long[] { ui, li };
103 }
104
105
106 /**
107 * Sturm sequence.
108 * @param f univariate polynomial.
109 * @param g univariate polynomial.
110 * @return a Sturm sequence for f and g.
111 */
112 public List<GenPolynomial<C>> sturmSequence(GenPolynomial<C> f, GenPolynomial<C> g) {
113 List<GenPolynomial<C>> S = new ArrayList<GenPolynomial<C>>();
114 if (f == null || f.isZERO()) {
115 return S;
116 }
117 if (f.isConstant()) {
118 S.add(f.monic());
119 return S;
120 }
121 GenPolynomial<C> F = f;
122 S.add(F);
123 GenPolynomial<C> G = g; //PolyUtil.<C> baseDeriviative(f);
124 while (!G.isZERO()) {
125 GenPolynomial<C> r = F.remainder(G);
126 F = G;
127 G = r.negate();
128 S.add(F/*.monic()*/);
129 }
130 //System.out.println("F = " + F);
131 if (F.isConstant()) {
132 return S;
133 }
134 // make squarefree
135 List<GenPolynomial<C>> Sp = new ArrayList<GenPolynomial<C>>(S.size());
136 for (GenPolynomial<C> p : S) {
137 p = p.divide(F);
138 Sp.add(p);
139 }
140 return Sp;
141 }
142
143
144 /**
145 * Complex root count of complex polynomial on rectangle.
146 * @param rect rectangle.
147 * @param a univariate complex polynomial.
148 * @return root count of a in rectangle.
149 */
150 @Override
151 public long complexRootCount(Rectangle<C> rect, GenPolynomial<Complex<C>> a)
152 throws InvalidBoundaryException {
153 return windingNumber(rect, a);
154 }
155
156
157 /**
158 * Winding number of complex function A on rectangle.
159 * @param rect rectangle.
160 * @param A univariate complex polynomial.
161 * @return winding number of A arround rect.
162 */
163 public long windingNumber(Rectangle<C> rect, GenPolynomial<Complex<C>> A) throws InvalidBoundaryException {
164 Boundary<C> bound = new Boundary<C>(rect, A); // throws InvalidBoundaryException
165 ComplexRing<C> cr = (ComplexRing<C>) A.ring.coFac;
166 RingFactory<C> cf = cr.ring;
167 C zero = cf.getZERO();
168 C one = cf.getONE();
169 long ix = 0L;
170 for (int i = 0; i < 4; i++) {
171 long ci = indexOfCauchy(zero, one, bound.getRealPart(i), bound.getImagPart(i));
172 //System.out.println("ci["+i+","+(i+1)+"] = " + ci);
173 ix += ci;
174 }
175 if (ix % 2L != 0) {
176 throw new InvalidBoundaryException("odd winding number " + ix);
177 }
178 return ix / 2L;
179 }
180
181
182 /**
183 * List of complex roots of complex polynomial a on rectangle.
184 * @param rect rectangle.
185 * @param a univariate squarefree complex polynomial.
186 * @return list of complex roots.
187 */
188 @Override
189 public List<Rectangle<C>> complexRoots(Rectangle<C> rect, GenPolynomial<Complex<C>> a)
190 throws InvalidBoundaryException {
191 ComplexRing<C> cr = (ComplexRing<C>) a.ring.coFac;
192 List<Rectangle<C>> roots = new ArrayList<Rectangle<C>>();
193 if ( a.isConstant() || a.isZERO() ) {
194 return roots;
195 }
196 //System.out.println("rect = " + rect);
197 long n = windingNumber(rect, a);
198 if (n < 0) { // can this happen?
199 throw new RuntimeException("negative winding number " + n);
200 //System.out.println("negative winding number " + n);
201 //return roots;
202 }
203 if (n == 0) {
204 return roots;
205 }
206 if (n == 1) {
207 roots.add(rect);
208 return roots;
209 }
210 Complex<C> eps = cr.fromInteger(1);
211 eps = eps.divide(cr.fromInteger(1000)); // 1/1000
212 //System.out.println("eps = " + eps);
213 //System.out.println("rect = " + rect);
214 // construct new center
215 Complex<C> delta = rect.corners[3].subtract(rect.corners[1]);
216 delta = delta.divide(cr.fromInteger(2));
217 //System.out.println("delta = " + delta);
218 boolean work = true;
219 while (work) {
220 Complex<C> center = rect.corners[1].sum(delta);
221 //System.out.println("center = " + toDecimal(center));
222 if (debug) {
223 logger.info("new center = " + center);
224 }
225 try {
226 Complex<C>[] cp = (Complex<C>[]) copyOfComplex(rect.corners, 4);
227 // (Complex<C>[]) new Complex[4]; cp[0] = rect.corners[0];
228 // cp[0] fix
229 cp[1] = new Complex<C>(cr, cp[1].getRe(), center.getIm());
230 cp[2] = center;
231 cp[3] = new Complex<C>(cr, center.getRe(), cp[3].getIm());
232 Rectangle<C> nw = new Rectangle<C>(cp);
233 //System.out.println("nw = " + nw);
234 List<Rectangle<C>> nwr = complexRoots(nw, a);
235 //System.out.println("#nwr = " + nwr.size());
236 roots.addAll(nwr);
237 if ( roots.size() == a.degree(0) ) {
238 work = false;
239 break;
240 }
241
242 cp = (Complex<C>[]) copyOfComplex(rect.corners, 4);
243 cp[0] = new Complex<C>(cr, cp[0].getRe(), center.getIm());
244 // cp[1] fix
245 cp[2] = new Complex<C>(cr, center.getRe(), cp[2].getIm());
246 cp[3] = center;
247 Rectangle<C> sw = new Rectangle<C>(cp);
248 //System.out.println("sw = " + sw);
249 List<Rectangle<C>> swr = complexRoots(sw, a);
250 //System.out.println("#swr = " + swr.size());
251 roots.addAll(swr);
252 if ( roots.size() == a.degree(0) ) {
253 work = false;
254 break;
255 }
256
257 cp = (Complex<C>[]) copyOfComplex(rect.corners, 4);
258 cp[0] = center;
259 cp[1] = new Complex<C>(cr, center.getRe(), cp[1].getIm());
260 // cp[2] fix
261 cp[3] = new Complex<C>(cr, cp[3].getRe(), center.getIm());
262 Rectangle<C> se = new Rectangle<C>(cp);
263 //System.out.println("se = " + se);
264 List<Rectangle<C>> ser = complexRoots(se, a);
265 //System.out.println("#ser = " + ser.size());
266 roots.addAll(ser);
267 if ( roots.size() == a.degree(0) ) {
268 work = false;
269 break;
270 }
271
272 cp = (Complex<C>[]) copyOfComplex(rect.corners, 4);
273 cp[0] = new Complex<C>(cr, center.getRe(), cp[0].getIm());
274 cp[1] = center;
275 cp[2] = new Complex<C>(cr, cp[2].getRe(), center.getIm());
276 // cp[3] fix
277 Rectangle<C> ne = new Rectangle<C>(cp);
278 //System.out.println("ne = " + ne);
279 List<Rectangle<C>> ner = complexRoots(ne, a);
280 //System.out.println("#ner = " + ner.size());
281 roots.addAll(ner);
282 work = false;
283 } catch (InvalidBoundaryException e) {
284 // repeat with new center
285 delta = delta.sum(delta.multiply(eps)); // distort
286 //System.out.println("new delta = " + toDecimal(delta));
287 eps = eps.sum(eps.multiply(cr.getIMAG()));
288 }
289 }
290 return roots;
291 }
292
293
294 /**
295 * Invariant rectangle for algebraic number.
296 * @param rect root isolating rectangle for f which contains exactly one root.
297 * @param f univariate polynomial, non-zero.
298 * @param g univariate polynomial, gcd(f,g) == 1.
299 * @return v a new rectangle contained in rect such that g(w) != 0 for w in v.
300 */
301 @Override
302 public Rectangle<C> invariantRectangle(Rectangle<C> rect,
303 GenPolynomial<Complex<C>> f,
304 GenPolynomial<Complex<C>> g)
305 throws InvalidBoundaryException {
306 Rectangle<C> v = rect;
307 if (g == null || g.isZERO()) {
308 return v;
309 }
310 if (g.isConstant()) {
311 return v;
312 }
313 if (f == null || f.isZERO() || f.isConstant()) { // ?
314 return v;
315 }
316 BigRational len = v.rationalLength();
317 BigRational half = new BigRational(1,2);
318 while (true) {
319 long n = windingNumber(v, g);
320 //System.out.println("n = " + n);
321 if (n < 0) { // can this happen?
322 throw new RuntimeException("negative winding number " + n);
323 }
324 if (n == 0) {
325 return v;
326 }
327 len = len.multiply(half);
328 Rectangle<C> v1 = v;
329 v = complexRootRefinement(v,f,len);
330 if ( v.equals(v1) ) {
331 //System.out.println("len = " + len);
332 if ( !f.gcd(g).isONE() ) {
333 System.out.println("f.gcd(g) = " + f.gcd(g));
334 throw new RuntimeException("no convergence " + v);
335 }
336 //break; // no convergence
337 }
338 }
339 //return v;
340 }
341
342 }