Source code

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