Source code

001/*
002 * $Id$
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.logging.log4j.Logger;
014import org.apache.logging.log4j.LogManager; 
015
016import edu.jas.arith.Rational;
017import edu.jas.poly.Complex;
018import edu.jas.poly.ComplexRing;
019import edu.jas.poly.GenPolynomial;
020import edu.jas.poly.GenPolynomialRing;
021import edu.jas.poly.AlgebraicNumberRing;
022import edu.jas.poly.AlgebraicNumber;
023import edu.jas.poly.PolyUtil;
024import edu.jas.poly.TermOrder;
025import edu.jas.root.Interval;
026import edu.jas.root.RealRootTuple;
027import edu.jas.root.AlgebraicRoots;
028import edu.jas.structure.GcdRingElem;
029import edu.jas.structure.RingFactory;
030import edu.jas.ufd.SquarefreeAbstract;
031import edu.jas.ufd.SquarefreeFactory;
032
033
034/**
035 * Roots factory.
036 * @author Heinz Kredel
037 */
038public class RootFactoryApp {
039
040
041    private static final Logger logger = LogManager.getLogger(RootFactoryApp.class);
042
043
044    private static final boolean debug = logger.isDebugEnabled();
045
046
047    /**
048     * Is complex algebraic number a root of a polynomial.
049     * @param f univariate polynomial.
050     * @param r complex algebraic number.
051     * @return true, if f(r) == 0, else false;
052     */
053    public static <C extends GcdRingElem<C> & Rational> boolean isRootRealCoeff(GenPolynomial<C> f,
054                    Complex<RealAlgebraicNumber<C>> r) {
055        RingFactory<C> cfac = f.ring.coFac;
056        ComplexRing<C> ccfac = new ComplexRing<C>(cfac);
057        GenPolynomialRing<Complex<C>> facc = new GenPolynomialRing<Complex<C>>(ccfac, f.ring);
058        GenPolynomial<Complex<C>> fc = PolyUtil.<C> complexFromAny(facc, f);
059        return isRoot(fc, r);
060    }
061
062
063    /**
064     * Is complex algebraic number a root of a polynomial.
065     * @param f univariate polynomial.
066     * @param r complex algebraic number.
067     * @return true, if f(r) == 0, else false;
068     */
069    public static <C extends GcdRingElem<C> & Rational> boolean isRoot(GenPolynomial<Complex<C>> f,
070                    Complex<RealAlgebraicNumber<C>> r) {
071        ComplexRing<RealAlgebraicNumber<C>> cr = r.factory();
072        GenPolynomialRing<Complex<RealAlgebraicNumber<C>>> cfac = new GenPolynomialRing<Complex<RealAlgebraicNumber<C>>>(
073                        cr, f.factory());
074        GenPolynomial<Complex<RealAlgebraicNumber<C>>> p;
075        p = PolyUtilApp.<C> convertToComplexRealCoefficients(cfac, f);
076        // test algebraic part
077        Complex<RealAlgebraicNumber<C>> a = PolyUtil.<Complex<RealAlgebraicNumber<C>>> evaluateMain(cr, p, r);
078        boolean t = a.isZERO();
079        if (!t) {
080            logger.info("f(r) = {}, f = {}, r  = {}", a, f, r);
081            return t;
082        }
083        // test approximation? not working
084        return true;
085    }
086
087
088    /**
089     * Is complex algebraic number a root of a polynomial.
090     * @param f univariate polynomial.
091     * @param R list of complex algebraic numbers.
092     * @return true, if f(r) == 0 for all r in R, else false;
093     */
094    public static <C extends GcdRingElem<C> & Rational> boolean isRoot(GenPolynomial<Complex<C>> f,
095                    List<Complex<RealAlgebraicNumber<C>>> R) {
096        for (Complex<RealAlgebraicNumber<C>> r : R) {
097            boolean t = isRoot(f, r);
098            if (!t) {
099                return false;
100            }
101        }
102        return true;
103    }
104
105
106    /**
107     * Complex algebraic number roots.
108     * @param f univariate polynomial.
109     * @return a list of different complex algebraic numbers, with f(c) == 0 for
110     *         c in roots.
111     */
112    public static <C extends GcdRingElem<C> & Rational> List<Complex<RealAlgebraicNumber<C>>> complexAlgebraicNumbersComplex(
113                    GenPolynomial<Complex<C>> f) {
114        GenPolynomialRing<Complex<C>> pfac = f.factory();
115        if (pfac.nvar != 1) {
116            throw new IllegalArgumentException("only for univariate polynomials");
117        }
118        ComplexRing<C> cfac = (ComplexRing<C>) pfac.coFac;
119        SquarefreeAbstract<Complex<C>> engine = SquarefreeFactory.<Complex<C>> getImplementation(cfac);
120        Map<GenPolynomial<Complex<C>>, Long> F = engine.squarefreeFactors(f.monic());
121        //System.out.println("S = " + F.keySet());
122        List<Complex<RealAlgebraicNumber<C>>> list = new ArrayList<Complex<RealAlgebraicNumber<C>>>();
123        for (Map.Entry<GenPolynomial<Complex<C>>,Long> me : F.entrySet()) {
124            GenPolynomial<Complex<C>> sp = me.getKey();
125            if (sp.isConstant() || sp.isZERO()) {
126                continue;
127            }
128            List<Complex<RealAlgebraicNumber<C>>> ls = RootFactoryApp.<C> complexAlgebraicNumbersSquarefree(sp);
129            long m = me.getValue();
130            for (long i = 0L; i < m; i++) {
131                list.addAll(ls);
132            }
133        }
134        return list;
135    }
136
137
138    /**
139     * Complex algebraic number roots.
140     * @param f univariate squarefree polynomial.
141     * @return a list of different complex algebraic numbers, with f(c) == 0 for
142     *         c in roots.
143     */
144    public static <C extends GcdRingElem<C> & Rational> 
145      List<Complex<RealAlgebraicNumber<C>>> complexAlgebraicNumbersSquarefree(
146                   GenPolynomial<Complex<C>> f) {
147        GenPolynomialRing<Complex<C>> pfac = f.factory();
148        if (pfac.nvar != 1) {
149            throw new IllegalArgumentException("only for univariate polynomials");
150        }
151        ComplexRing<C> cfac = (ComplexRing<C>) pfac.coFac;
152        TermOrder to = new TermOrder(TermOrder.INVLEX);
153        GenPolynomialRing<Complex<C>> tfac = new GenPolynomialRing<Complex<C>>(cfac, 2, to); //,new vars); //tord?
154        //System.out.println("tfac = " + tfac);
155        GenPolynomial<Complex<C>> t = tfac.univariate(1, 1L).sum(
156                        tfac.univariate(0, 1L).multiply(cfac.getIMAG()));
157        //System.out.println("t = " + t); // t = x + i y
158        GenPolynomialRing<C> rfac = new GenPolynomialRing<C>(cfac.ring, tfac); //tord?
159        //System.out.println("rfac = " + rfac);
160        List<Complex<RealAlgebraicNumber<C>>> list = new ArrayList<Complex<RealAlgebraicNumber<C>>>();
161        GenPolynomial<Complex<C>> sp = f;
162        if (sp.isConstant() || sp.isZERO()) {
163            return list;
164        }
165        // substitute t = x + i y
166        GenPolynomial<Complex<C>> su = PolyUtil.<Complex<C>> substituteUnivariate(sp, t);
167        //System.out.println("su = " + su);
168        su = su.monic();
169        //System.out.println("su = " + su);
170        GenPolynomial<C> re = PolyUtil.<C> realPartFromComplex(rfac, su);
171        GenPolynomial<C> im = PolyUtil.<C> imaginaryPartFromComplex(rfac, su);
172        if (debug) {
173            logger.debug("rfac = {}", () -> rfac.toScript());
174            logger.debug("t  = {}, re = {}, im = {}", t, re.toScript(), im.toScript()); // ()-> not possible
175        }
176        List<GenPolynomial<C>> li = new ArrayList<GenPolynomial<C>>(2);
177        li.add(re);
178        li.add(im);
179        Ideal<C> id = new Ideal<C>(rfac, li);
180        //System.out.println("id = " + id);
181        List<IdealWithUniv<C>> idul = id.zeroDimRootDecomposition();
182
183        IdealWithRealAlgebraicRoots<C> idr;
184        for (IdealWithUniv<C> idu : idul) {
185            //System.out.println("---idu = " + idu);
186            idr = PolyUtilApp.<C> realAlgebraicRoots(idu);
187            //System.out.println("---idr = " + idr);
188            for (List<edu.jas.root.RealAlgebraicNumber<C>> crr : idr.ran) {
189                //System.out.println("crr = " + crr);
190                RealRootTuple<C> root = new RealRootTuple<C>(crr);
191                //System.out.println("root = " + root);
192                RealAlgebraicRing<C> car = new RealAlgebraicRing<C>(idu, root);
193                //System.out.println("car = " + car);
194                List<RealAlgebraicNumber<C>> gens = car.generators();
195                //System.out.println("gens = " + gens);
196                int sg = gens.size();
197                RealAlgebraicNumber<C> rre = gens.get(sg - 2);
198                RealAlgebraicNumber<C> rim = gens.get(sg - 1);
199                ComplexRing<RealAlgebraicNumber<C>> cring = new ComplexRing<RealAlgebraicNumber<C>>(car);
200                Complex<RealAlgebraicNumber<C>> crn = new Complex<RealAlgebraicNumber<C>>(cring, rre, rim);
201                //System.out.println("crn = " + crn + " in " + crn.ring);
202                // refine intervals if necessary, not meaningful
203                list.add(crn);
204            }
205        }
206        return list;
207    }
208
209
210    /* approximation ?
211    List<ComplexAlgebraicNumber<C>> complexAlgebraicNumbersComplex(GenPolynomial<Complex<C>> f, BigRational eps)
212    */
213
214
215    /**
216     * Root reduce of real and complex algebraic numbers.
217     * @param a container of real and complex algebraic numbers.
218     * @param b container of real and complex algebraic numbers.
219     * @return container of real and complex algebraic numbers 
220     *         of the primitive element of a and b.
221     */
222    public static <C extends GcdRingElem<C> & Rational> 
223           AlgebraicRootsPrimElem<C> rootReduce(AlgebraicRoots<C> a, AlgebraicRoots<C> b) {
224        return rootReduce(a.getAlgebraicRing(), b.getAlgebraicRing());
225    }
226
227
228    /**
229     * Root reduce of real and complex algebraic numbers.
230     * @param a polynomial.
231     * @param b polynomial.
232     * @return container of real and complex algebraic numbers 
233     *         of the primitive element of a and b.
234     */
235    public static <C extends GcdRingElem<C> & Rational> 
236           AlgebraicRootsPrimElem<C> rootReduce(GenPolynomial<C> a, GenPolynomial<C> b) {
237        AlgebraicNumberRing<C> anr = new AlgebraicNumberRing<C>(a);
238        AlgebraicNumberRing<C> bnr = new AlgebraicNumberRing<C>(b);
239        return rootReduce(anr, bnr);
240    }
241
242
243    /**
244     * Root reduce of real and complex algebraic numbers.
245     * @param a algebraic number ring.
246     * @param b algebraic number ring.
247     * @return container of real and complex algebraic numbers 
248     *         of the primitive element of a and b.
249     */
250    public static <C extends GcdRingElem<C> & Rational> 
251           AlgebraicRootsPrimElem<C> rootReduce(AlgebraicNumberRing<C> a, AlgebraicNumberRing<C> b) {
252        PrimitiveElement<C> pe = PolyUtilApp.<C>primitiveElement(a, b);
253        AlgebraicRoots<C> ar = edu.jas.root.RootFactory.<C>algebraicRoots(pe.primitiveElem.modul);
254        return new AlgebraicRootsPrimElem<C>(ar, pe);
255    }
256
257
258    /**
259     * Roots of unity of real and complex algebraic numbers.
260     * @param ar container of real and complex algebraic numbers with primitive element.
261     * @return container of real and complex algebraic numbers which are roots
262     *         of unity.
263     */
264    public static <C extends GcdRingElem<C> & Rational> 
265           AlgebraicRootsPrimElem<C> rootsOfUnity(AlgebraicRootsPrimElem<C> ar) {
266        AlgebraicRoots<C> ur = edu.jas.root.RootFactory.rootsOfUnity(ar);
267        if (ar.pelem == null) {
268            return new AlgebraicRootsPrimElem<C>(ur, ar.pelem);
269        }
270        List<AlgebraicNumber<C>> al = new ArrayList<AlgebraicNumber<C>>();
271        long d = ar.pelem.primitiveElem.modul.degree();
272        AlgebraicNumber<C> c = ar.pelem.A;
273        AlgebraicNumber<C> m = c.ring.getONE();
274        for (long i = 1; i <= d; i++) {
275            m = m.multiply(c);
276            if (m.isRootOfUnity()) {
277                if (!al.contains(m)) {
278                    al.add(m);
279                }
280            }
281        }
282        c = ar.pelem.B;
283        m = c.ring.getONE();
284        for (long i = 1; i <= d; i++) {
285            m = m.multiply(c);
286            if (m.isRootOfUnity()) {
287                if (!al.contains(m)) {
288                    al.add(m);
289                }
290            }
291        }
292        return new AlgebraicRootsPrimElem<C>(ur, ar.pelem, al);
293    }
294
295}