Source code

001/*
002 * $Id: PartialFraction.java 5940 2018-10-19 08:53:13Z kredel $
003 */
004
005package edu.jas.ufd;
006
007
008import java.io.Serializable;
009import java.util.ArrayList;
010import java.util.HashMap;
011import java.util.HashSet;
012import java.util.List;
013import java.util.Map;
014import java.util.Set;
015
016import org.apache.logging.log4j.Logger;
017import org.apache.logging.log4j.LogManager; 
018
019import edu.jas.poly.AlgebraicNumber;
020import edu.jas.poly.AlgebraicNumberRing;
021import edu.jas.poly.GenPolynomial;
022import edu.jas.poly.GenPolynomialRing;
023import edu.jas.poly.PolyUtil;
024import edu.jas.structure.GcdRingElem;
025
026
027/**
028 * Container for the partial fraction decomposition of a squarefree denominator.
029 * num/den = sum( a_i / d_i )
030 * @author Heinz Kredel
031 * @param <C> coefficient type
032 */
033
034public class PartialFraction<C extends GcdRingElem<C>> implements Serializable {
035
036
037    private static final Logger logger = LogManager.getLogger(PartialFraction.class);
038
039
040    /**
041     * Original numerator polynomial coefficients from C and deg(num) &lt;
042     * deg(den).
043     */
044    public final GenPolynomial<C> num;
045
046
047    /**
048     * Original (irreducible) denominator polynomial coefficients from C.
049     */
050    public final GenPolynomial<C> den;
051
052
053    /**
054     * List of numbers from C.
055     */
056    public final List<C> cfactors;
057
058
059    /**
060     * List of linear factors of the denominator with coefficients from C.
061     */
062    public final List<GenPolynomial<C>> cdenom;
063
064
065    /**
066     * List of algebraic numbers of an algebraic field extension over C.
067     */
068    public final List<AlgebraicNumber<C>> afactors;
069
070
071    /**
072     * List of factors of the denominator with coefficients from an
073     * AlgebraicNumberRing&lt;C&gt;.
074     */
075    public final List<GenPolynomial<AlgebraicNumber<C>>> adenom;
076
077
078    /**
079     * Constructor.
080     * @param n numerator GenPolynomial over C.
081     * @param d irreducible denominator GenPolynomial over C.
082     * @param cf list of elements a_i.
083     * @param cd list of linear factors d_i of d.
084     * @param af list of algebraic elements a_i.
085     * @param ad list of linear (irreducible) factors d_i of d with algebraic
086     *            coefficients. n/d = sum( a_i / d_i )
087     */
088    public PartialFraction(GenPolynomial<C> n, GenPolynomial<C> d, List<C> cf, List<GenPolynomial<C>> cd,
089                    List<AlgebraicNumber<C>> af, List<GenPolynomial<AlgebraicNumber<C>>> ad) {
090        num = n;
091        den = d;
092        cfactors = cf;
093        cdenom = cd;
094        afactors = af;
095        adenom = ad;
096        for (GenPolynomial<C> p : cdenom) {
097            if (p.degree(0) > 1) {
098                throw new IllegalArgumentException("polynomial not linear, p = " + p);
099            }
100        }
101        for (GenPolynomial<AlgebraicNumber<C>> a : adenom) {
102            if (a.degree(0) > 1) {
103                throw new IllegalArgumentException("polynomial not linear, a = " + a);
104            }
105        }
106    }
107
108
109    /**
110     * Get the String representation.
111     * @see java.lang.Object#toString()
112     */
113    @Override
114    public String toString() {
115        StringBuffer sb = new StringBuffer();
116        sb.append("(" + num.toString() + ")");
117        sb.append(" / ");
118        sb.append("(" + den.toString() + ")");
119        sb.append(" =\n");
120        boolean first = true;
121        for (int i = 0; i < cfactors.size(); i++) {
122            C cp = cfactors.get(i);
123            if (first) {
124                first = false;
125            } else {
126                sb.append(" + ");
127            }
128            sb.append("(" + cp.toString() + ")");
129            GenPolynomial<C> p = cdenom.get(i);
130            sb.append(" / (" + p.toString() + ")");
131        }
132        if (!first && afactors.size() > 0) {
133            sb.append(" + ");
134        }
135        first = true;
136        for (int i = 0; i < afactors.size(); i++) {
137            if (first) {
138                first = false;
139            } else {
140                sb.append(" + ");
141            }
142            AlgebraicNumber<C> ap = afactors.get(i);
143            AlgebraicNumberRing<C> ar = ap.factory();
144            GenPolynomial<AlgebraicNumber<C>> p = adenom.get(i);
145            if (p.degree(0) < ar.modul.degree(0) && ar.modul.degree(0) > 2) {
146                sb.append("sum_(" + ar.getGenerator() + " in ");
147                sb.append("rootOf(" + ar.modul + ") ) ");
148            } else {
149                //sb.append("sum_("+ar+") ");
150            }
151            sb.append("(" + ap.toString() + ")");
152            sb.append(" / (" + p.toString() + ")");
153            //sb.append(" ## over " + ap.factory() + "\n");
154        }
155        return sb.toString();
156    }
157
158
159    /**
160     * Get a scripting compatible string representation.
161     * @return script compatible representation for this container.
162     * @see edu.jas.structure.ElemFactory#toScript()
163     */
164    public String toScript() {
165        // Python case
166        StringBuffer sb = new StringBuffer();
167        sb.append(num.toScript());
168        sb.append(" / ");
169        sb.append(den.toScript());
170        sb.append(" = ");
171        boolean first = true;
172        int i = 0;
173        for (C cp : cfactors) {
174            if (first) {
175                first = false;
176            } else {
177                sb.append(" + ");
178            }
179            sb.append(cp.toScript());
180            GenPolynomial<C> p = cdenom.get(i);
181            sb.append(" / " + p.toScript());
182        }
183        if (!first) {
184            sb.append(" + ");
185        }
186        first = true;
187        i = 0;
188        for (AlgebraicNumber<C> ap : afactors) {
189            if (first) {
190                first = false;
191            } else {
192                sb.append(" + ");
193            }
194            AlgebraicNumberRing<C> ar = ap.factory();
195            GenPolynomial<AlgebraicNumber<C>> p = adenom.get(i);
196            if (p.degree(0) < ar.modul.degree(0) && ar.modul.degree(0) > 2) {
197                sb.append("sum_(" + ar.getGenerator().toScript() + " in ");
198                sb.append("rootOf(" + ar.modul.toScript() + ") ) ");
199            } else {
200                //sb.append("sum_("+ar+") ");
201            }
202            sb.append(ap.toScript());
203            sb.append(" / " + p.toScript());
204            //sb.append(" ## over " + ap.toScriptFactory() + "\n");
205        }
206        return sb.toString();
207    }
208
209
210    /**
211     * Hash code for this Factors.
212     * @see java.lang.Object#hashCode()
213     */
214    @Override
215    public int hashCode() {
216        int h = num.hashCode();
217        h = h * 37 + den.hashCode();
218        h = h * 37 + cfactors.hashCode();
219        h = h * 37 + cdenom.hashCode();
220        h = h * 37 + afactors.hashCode();
221        h = h * 37 + adenom.hashCode();
222        return h;
223    }
224
225
226    /**
227     * Comparison with any other object.
228     * @see java.lang.Object#equals(java.lang.Object)
229     */
230    @Override
231    @SuppressWarnings("unchecked")
232    public boolean equals(Object B) {
233        if (B == null) {
234            return false;
235        }
236        if (!(B instanceof PartialFraction)) {
237            return false;
238        }
239        PartialFraction<C> a = (PartialFraction<C>) B;
240        boolean t = num.equals(a.num) && den.equals(a.den);
241        if (!t) {
242            return t;
243        }
244        t = cfactors.equals(a.cfactors);
245        if (!t) {
246            return t;
247        }
248        t = cdenom.equals(a.cdenom);
249        if (!t) {
250            return t;
251        }
252        t = afactors.equals(a.afactors);
253        if (!t) {
254            return t;
255        }
256        t = adenom.equals(a.adenom);
257        return t;
258    }
259
260
261    /**
262     * Test if correct partial fraction. num/den = sum( a_i / d_i )
263     */
264    @SuppressWarnings("unchecked")
265    public boolean isPartialFraction() {
266        QuotientRing<C> qfac = new QuotientRing<C>(num.ring);
267        // num / den
268        Quotient<C> q = new Quotient<C>(qfac, num, den);
269        //System.out.println("q = " + q);
270        Quotient<C> qs = qfac.getZERO();
271        int i = 0;
272        for (C c : cfactors) {
273            GenPolynomial<C> cp = cdenom.get(i++);
274            // plus c / cp
275            GenPolynomial<C> cd = num.ring.getONE().multiply(c);
276            Quotient<C> qq = new Quotient<C>(qfac, cd, cp);
277            qs = qs.sum(qq);
278        }
279        //System.out.println("qs = " + qs);
280        if (afactors.isEmpty()) {
281            return q.compareTo(qs) == 0;
282        }
283
284        // sort by extension field
285        Set<AlgebraicNumberRing<C>> fields = new HashSet<AlgebraicNumberRing<C>>();
286        for (AlgebraicNumber<C> ap : afactors) {
287            if (ap.ring.depth() > 1) {
288                logger.warn("extension field depth to high"); // todo
289            }
290            fields.add(ap.ring);
291        }
292        //System.out.println("fields = " + fields);
293        Map<AlgebraicNumberRing<C>, List<AlgebraicNumber<C>>> facs = new HashMap<AlgebraicNumberRing<C>, List<AlgebraicNumber<C>>>();
294        for (AlgebraicNumber<C> ap : afactors) {
295            List<AlgebraicNumber<C>> cf = facs.get(ap.ring);
296            if (cf == null) {
297                cf = new ArrayList<AlgebraicNumber<C>>();
298            }
299            cf.add(ap);
300            facs.put(ap.ring, cf);
301        }
302        //System.out.println("facs = " + facs);
303        Map<AlgebraicNumberRing<C>, List<GenPolynomial<AlgebraicNumber<C>>>> pfacs = new HashMap<AlgebraicNumberRing<C>, List<GenPolynomial<AlgebraicNumber<C>>>>();
304        for (GenPolynomial<AlgebraicNumber<C>> ap : adenom) {
305            AlgebraicNumberRing<C> ar = (AlgebraicNumberRing<C>) ap.ring.coFac;
306            List<GenPolynomial<AlgebraicNumber<C>>> cf = pfacs.get(ar);
307            if (cf == null) {
308                cf = new ArrayList<GenPolynomial<AlgebraicNumber<C>>>();
309            }
310            cf.add(ap);
311            pfacs.put(ar, cf);
312        }
313        //System.out.println("pfacs = " + pfacs);
314
315        // check algebraic parts 
316        boolean sumMissing = false;
317        for (AlgebraicNumberRing<C> ar : fields) {
318            if (ar.modul.degree(0) > 2) { //&& p.degree(0) < ar.modul.degree(0) ?
319                sumMissing = true;
320            }
321            List<AlgebraicNumber<C>> cf = facs.get(ar);
322            List<GenPolynomial<AlgebraicNumber<C>>> cfp = pfacs.get(ar);
323            GenPolynomialRing<AlgebraicNumber<C>> apfac = cfp.get(0).ring;
324            QuotientRing<AlgebraicNumber<C>> aqfac = new QuotientRing<AlgebraicNumber<C>>(apfac);
325            Quotient<AlgebraicNumber<C>> aq = aqfac.getZERO();
326            i = 0;
327            for (AlgebraicNumber<C> c : cf) {
328                GenPolynomial<AlgebraicNumber<C>> cp = cfp.get(i++);
329                // plus c / cp
330                GenPolynomial<AlgebraicNumber<C>> cd = apfac.getONE().multiply(c);
331                Quotient<AlgebraicNumber<C>> qq = new Quotient<AlgebraicNumber<C>>(aqfac, cd, cp);
332                //System.out.println("qq = " + qq);
333                aq = aq.sum(qq);
334            }
335            //System.out.println("aq = " + aq);
336            GenPolynomialRing<C> cfac = ar.ring;
337            GenPolynomialRing<GenPolynomial<C>> prfac = new GenPolynomialRing<GenPolynomial<C>>(cfac, apfac);
338            GenPolynomial<GenPolynomial<C>> pqnum = PolyUtil.<C> fromAlgebraicCoefficients(prfac, aq.num);
339            GenPolynomial<GenPolynomial<C>> pqden = PolyUtil.<C> fromAlgebraicCoefficients(prfac, aq.den);
340            //System.out.println("pq = (" + pqnum + ") / (" + pqden + ")");
341
342            C one = cfac.coFac.getONE(); // varaible should no more occur in coefficient
343            GenPolynomialRing<C> pfac = new GenPolynomialRing<C>(cfac.coFac, prfac);
344            GenPolynomial<C> pnum = PolyUtil.<C> evaluateFirstRec(cfac, pfac, pqnum, one);
345            GenPolynomial<C> pden = PolyUtil.<C> evaluateFirstRec(cfac, pfac, pqden, one);
346            //System.out.println("p = (" + pnum + ") / (" + pden + ")");
347
348            // iterate if multiple field extensions
349            while (cfac.coFac instanceof AlgebraicNumberRing) {
350                //System.out.println("cfac.coFac = " + cfac.coFac.toScript());
351                AlgebraicNumberRing<C> ar2 = (AlgebraicNumberRing<C>) cfac.coFac;
352                cfac = ar2.ring;
353                prfac = new GenPolynomialRing<GenPolynomial<C>>(cfac, apfac);
354                GenPolynomial<AlgebraicNumber<C>> prnum = (GenPolynomial<AlgebraicNumber<C>>) pnum;
355                GenPolynomial<AlgebraicNumber<C>> prden = (GenPolynomial<AlgebraicNumber<C>>) pden;
356                pqnum = PolyUtil.<C> fromAlgebraicCoefficients(prfac, prnum);
357                pqden = PolyUtil.<C> fromAlgebraicCoefficients(prfac, prden);
358                one = cfac.coFac.getONE(); // varaible should no more occur in coefficient
359                pfac = new GenPolynomialRing<C>(cfac.coFac, prfac);
360                pnum = PolyUtil.<C> evaluateFirstRec(cfac, pfac, pqnum, one);
361                pden = PolyUtil.<C> evaluateFirstRec(cfac, pfac, pqden, one);
362            }
363
364            Quotient<C> qq = new Quotient<C>(qfac, pnum, pden);
365            //System.out.println("qq = " + qq);
366            qs = qs.sum(qq);
367        }
368        boolean cmp = q.compareTo(qs) == 0;
369        if (!cmp) {
370            System.out.println("q != qs: " + q + " != " + qs);
371        }
372        return cmp || sumMissing;
373    }
374
375}