Source code

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