Source code

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