001/*
002 * JScience - Java(TM) Tools and Libraries for the Advancement of Sciences.
003 * Copyright (C) 2006 - JScience (http://jscience.org/)
004 * All rights reserved.
005 * 
006 * Permission to use, copy, modify, and distribute this software is
007 * freely granted, provided that this notice is preserved.
008 */
009package org.jscience.mathematics.function;
010
011import java.util.List;
012import java.util.Map;
013import java.util.Set;
014
015import org.jscience.mathematics.structure.GroupAdditive;
016import org.jscience.mathematics.structure.GroupMultiplicative;
017import org.jscience.mathematics.structure.Ring;
018
019import javolution.util.FastMap;
020import javolution.util.FastTable;
021import javolution.context.ObjectFactory;
022import javolution.text.Text;
023import javolution.text.TextBuilder;
024
025/**
026 * <p> This class represents a mathematical expression involving a sum of powers
027 *     in one or more {@link Variable variables} multiplied by 
028 *     coefficients (such as <code>x² + x·y + 3y²</code>).</p>
029 *     
030 * <p> Polynomials are characterized by the type of variable they operate 
031 *     upon. For example:[code]
032 *           Variable<Amount<?>> varX = new Variable.Local<Amount<?>>("x");
033 *           Polynomial<Amount<?>> x = Polynomial.valueOf(Amount.valueOf(1, SI.METER), varX);
034 *     and
035 *           Variable<Complex> varX = new Variable.Local<Complex>("x");
036 *           Polynomial<Complex> x = Polynomial.valueOf(Complex.ONE, varX);[/code]
037 *     are two different polynomials, the first one operates on physical 
038 *     {@link org.jscience.physics.amount.Amount measures},
039 *     whereas the second operates on 
040 *     {@link org.jscience.mathematics.number.Complex complex} numbers.</p>
041 *     
042 * <p> Terms (others than {@link Term#ONE ONE}) having zero (additive identity) 
043 *     for coefficient are automatically removed.</p>    
044 * 
045 * @author  <a href="mailto:jean-marie@dautelle.com">Jean-Marie Dautelle</a>
046 * @version 3.1, April 1, 2006
047 */
048public class Polynomial<R extends Ring<R>> extends Function<R, R> implements 
049         Ring<Polynomial<R>> {
050
051    /**
052     * Holds the terms to coefficients mapping 
053     * (never empty, holds Term.ONE when constant)
054     */
055    final FastMap<Term, R> _termToCoef = new FastMap<Term, R>();
056
057    /**
058     * Default constructor.
059     */
060    Polynomial() {
061    }
062
063    /**
064     * Returns an univariate polynomial of degree one with the specified 
065     * coefficient multiplier.
066     * 
067     * @param coefficient the coefficient for the variable of degree 1.  
068     * @param variable the variable for this polynomial.
069     * @return <code>valueOf(coefficient, Term.valueOf(variable, 1))</code>
070     */
071    public static <R extends Ring<R>> Polynomial<R> valueOf(R coefficient,
072            Variable<R> variable) {
073        return valueOf(coefficient, Term.valueOf(variable, 1));
074    }
075
076    /**
077     * Returns a polynomial corresponding to the specified {@link Term term}
078     * with the specified coefficient multiplier.
079     * 
080     * @param coefficient the coefficient multiplier.  
081     * @param term the term multiplicand.
082     * @return <code>coefficient * term</code>
083     */
084    public static <R extends Ring<R>> Polynomial<R> valueOf(R coefficient,
085            Term term) {
086        if (term.equals(Term.ONE)) return Constant.valueOf(coefficient);
087        if (isZero(coefficient)) return Constant.valueOf(coefficient);
088        Polynomial<R> p = Polynomial.newInstance();
089        p._termToCoef.put(term, coefficient);
090        return p;
091    }
092
093    private static boolean isZero(GroupAdditive<?> coefficient) {
094        return coefficient.equals(coefficient.opposite());
095    }
096    
097    /**
098     * Returns the terms of this polynomial.
099     * 
100     * @return this polynomial's terms.
101     */
102    public Set<Term> getTerms() {
103        return _termToCoef.unmodifiable().keySet();
104    }
105
106    /**
107     * Returns the coefficient for the specified term.
108     * 
109     * @param term the term for which the coefficient is returned.
110     * @return the coefficient for the specified term or <code>null</code>
111     *         if this polynomial does not contain the specified term.
112     */
113    public final R getCoefficient(Term term) {
114        return _termToCoef.get(term);
115    }
116
117    /**
118     * Returns the order of this polynomial for the specified variable.
119     * 
120     * @return the polynomial order relative to the specified variable.
121     */
122    public int getOrder(Variable<R> v) {
123        int order = 0;
124        for (Term term : _termToCoef.keySet()) {
125            int power = term.getPower(v);
126            if (power > order) {
127                order = power;
128            }
129        }
130        return order;
131    }
132
133    /**
134     * Returns the sum of this polynomial with a constant polynomial 
135     * having the specified value (convenience method).
136     * 
137     * @param constantValue the value of the constant polynomial to add.
138     * @return <code>this + Constant.valueOf(constantValue)</code>
139     */
140    public Polynomial<R> plus(R constantValue) {
141        return this.plus(Constant.valueOf(constantValue));
142    }
143
144    /**
145     * Returns the product of this polynomial with a constant polynomial 
146     * having the specified value (convenience method).
147     * 
148     * @param constantValue the value of the constant polynomial to multiply.
149     * @return <code>this · Constant.valueOf(constantValue)</code>
150     */
151    public Polynomial<R> times(R constantValue) {
152        return this.times(Constant.valueOf(constantValue));
153    }
154
155    /**
156     * Returns the sum of two polynomials.
157     * 
158     * @param that the polynomial being added.
159     * @return <code>this + that</code>
160     */
161    public Polynomial<R> plus(Polynomial<R> that) {
162        Polynomial<R> result = Polynomial.newInstance();
163        R zero = null;
164        result._termToCoef.putAll(this._termToCoef);
165        result._termToCoef.putAll(that._termToCoef);
166        for (FastMap.Entry<Term, R> e = result._termToCoef.head(), 
167                end = result._termToCoef.tail(); (e = e.getNext()) != end;) {
168            Term term = e.getKey();
169            R thisCoef = this._termToCoef.get(term);
170            R thatCoef = that._termToCoef.get(term);
171            if ((thisCoef != null) && (thatCoef != null)) {
172                R sum = thisCoef.plus(thatCoef);
173                if (isZero(sum)) { // Remove entry (be careful iterating)
174                    FastMap.Entry<Term, R> prev = e.getPrevious();
175                    result._termToCoef.remove(term);
176                    e = prev; // Move back to previous entry.
177                    zero = sum; // To be used if constant polynomial.
178                } else {
179                    result._termToCoef.put(term, sum);
180                }
181            } // Else the current coefficient is correct.
182        }
183        if (result._termToCoef.size() == 0) return Constant.valueOf(zero);
184        return result;
185    }
186
187    /**
188     * Returns the opposite of this polynomial.
189     * 
190     * @return <code> - this</code>
191     */
192    public Polynomial<R> opposite() {
193        Polynomial<R> result = Polynomial.newInstance();
194        for (FastMap.Entry<Term, R> e = _termToCoef.head(), 
195                end = _termToCoef.tail(); (e = e.getNext()) != end;) {
196            result._termToCoef.put(e.getKey(), e.getValue().opposite());
197        }
198        return result;
199    }
200    
201    /**
202     * Returns the difference of two polynomials.
203     * 
204     * @param that the polynomial being subtracted.
205     * @return <code>this - that</code>
206     */
207    public Polynomial<R> minus(Polynomial<R> that) {
208        return this.plus(that.opposite());
209    }
210
211    /**
212     * Returns the product of two polynomials.
213     * 
214     * @param that the polynomial multiplier.
215     * @return <code>this · that</code>
216     */
217    public Polynomial<R> times(Polynomial<R> that) {
218        Polynomial<R> result = Polynomial.newInstance();
219        R zero = null;
220        for (Map.Entry<Term, R> entry1 : this._termToCoef.entrySet()) {
221            Term t1 = entry1.getKey();
222            R c1 = entry1.getValue();
223            for (Map.Entry<Term, R> entry2 : that._termToCoef.entrySet()) {
224                Term t2 = entry2.getKey();
225                R c2 = entry2.getValue();
226                Term t = t1.times(t2);
227                R c = c1.times(c2);
228                R prev = result.getCoefficient(t);
229                R coef = (prev != null) ? prev.plus(c) : c;
230                if (isZero(coef)) {
231                    zero = coef;
232                } else {
233                    result._termToCoef.put(t, coef);
234                }
235            }
236        }
237        if (result._termToCoef.size() == 0) return Constant.valueOf(zero);
238        return result;
239    }
240
241    /**
242     * Returns the composition of this polynomial with the one specified.
243     *
244     * @param  that the polynomial for which the return value is passed as
245     *         argument to this function.
246     * @return the polynomial <code>(this o that)</code>
247     * @throws FunctionException if this function is not univariate.
248     */
249    public Polynomial<R> compose(Polynomial<R> that) {
250        List<Variable<R>> variables = getVariables();
251        if (getVariables().size() != 1)
252            throw new FunctionException("This polynomial is not monovariate");
253        Variable<R> v = variables.get(0);
254        Polynomial<R> result = null;
255        for (Map.Entry<Term, R> entry : this._termToCoef.entrySet()) {
256            Term term = entry.getKey();
257            Constant<R> cst = Constant.valueOf(entry.getValue());
258            int power = term.getPower(v);
259            if (power > 0) {
260                Polynomial<R> fn = that.pow(power);
261                result = (result != null) ? result.plus(cst.times(fn)) : cst
262                        .times(fn);
263            } else { // power = 0
264                result = (result != null) ? result.plus(cst) : cst;
265            }
266        }
267        return result;
268    }
269
270    ////////////////////////////////////////////////////////////////
271    // Overrides parent method potentially returning polynomials  //
272    ////////////////////////////////////////////////////////////////
273    
274    @SuppressWarnings("unchecked")
275    @Override
276    public <Z> Function<Z, R> compose(Function<Z, R> that) {
277      return (Function<Z, R>) ((that instanceof Polynomial) ?
278      compose((Polynomial)that) : super.compose(that)); 
279    }
280
281    @Override
282    public Polynomial<R> differentiate(Variable<R> v) {
283        if (this.getOrder(v) > 0) {
284            Polynomial<R> result = null;
285            for (Map.Entry<Term, R> entry : this._termToCoef.entrySet()) {
286                Term term = entry.getKey();
287                R coef = entry.getValue();
288                int power = term.getPower(v);
289                if (power > 0) {
290                    R newCoef = multiply(coef, power);
291                    Term newTerm = term.divide(Term.valueOf(v, 1));
292                    Polynomial<R> p = valueOf(newCoef, newTerm);
293                    result = (result != null) ? result.plus(p) : p;
294                }
295            }
296            return result;
297        } else { // Returns zero.
298            R coef = _termToCoef.values().iterator().next();
299            return Constant.valueOf(coef.plus(coef.opposite()));
300        }
301    }
302
303    private static <R extends Ring<R>> R multiply(R o, int n) {
304        if (n <= 0)
305            throw new IllegalArgumentException("n: " + n
306                    + " zero or negative values not allowed");
307        R shift2 = o;
308        R result = null;
309        while (n >= 1) { // Iteration.
310            if ((n & 1) == 1) {
311                result = (result == null) ? shift2 : result.plus(shift2);
312            }
313            shift2 = shift2.plus(shift2);
314            n >>>= 1;
315        }
316        return result;
317    }
318
319    @SuppressWarnings("unchecked")
320    @Override
321    public Polynomial<R> integrate(Variable<R> v) {
322        Polynomial<R> result = null;
323        for (Map.Entry<Term, R> entry : this._termToCoef.entrySet()) {
324            Term term = entry.getKey();
325            R coef = entry.getValue();
326            int power = term.getPower(v);
327            R newCoef = (R)((GroupMultiplicative)multiply((R)((GroupMultiplicative)coef).inverse(), power + 1)).inverse();
328            Term newTerm = term.times(Term.valueOf(v, 1));
329            Polynomial<R> p = valueOf(newCoef, newTerm);
330            result = (result != null) ? result.plus(p) : p;
331        }
332        return result;
333    }
334
335    @SuppressWarnings("unchecked")
336    @Override
337    public Function<R, R> plus(Function<R, R> that) {
338        return (that instanceof Polynomial) ?
339                this.plus((Polynomial<R>)that) : super.plus(that);       
340    }
341
342    @SuppressWarnings("unchecked")
343    @Override
344    public Function<R, R> minus(Function<R, R> that) {
345        return (that instanceof Polynomial) ?
346                this.minus((Polynomial<R>)that) : super.minus(that);       
347    }
348
349    @SuppressWarnings("unchecked")
350    @Override
351    public Function<R, R> times(Function<R, R> that) {
352        return (that instanceof Polynomial) ?
353                this.times((Polynomial<R>)that) : super.times(that);       
354    }
355    
356    @SuppressWarnings("unchecked")
357    @Override
358    public Polynomial<R> pow(int n) {
359        return (Polynomial<R>) super.pow(n);
360    }
361
362    @SuppressWarnings("unchecked")
363    @Override
364    public List<Variable<R>> getVariables() {
365        // We multiply all terms togethers, the resulting product 
366        // will hold all variabgles (powers are always positive).
367        Term product = _termToCoef.head().getNext().getKey();
368        for (FastMap.Entry<Term, R> e = _termToCoef.head().getNext(), 
369                end = _termToCoef.tail(); (e = e.getNext()) != end;) {
370            product = product.times(e.getKey());
371        }
372        FastTable vars = FastTable.newInstance();
373        for (int i=0, n = product.size(); i < n; i++) {
374            vars.add(product.getVariable(i));
375        }
376        return vars;
377    }
378
379    @Override
380    @SuppressWarnings("unchecked")
381    public R evaluate() {
382        R sum = null;
383        for (Map.Entry<Term, R> entry : _termToCoef.entrySet()) {
384            Term term = entry.getKey();
385            R coef = entry.getValue();
386            R termValue = (R) term.evaluate();
387            R value = (termValue != null) ? coef.times(termValue) : coef;
388            sum = (sum == null) ? value : sum.plus(value);
389        }
390        return sum;
391    }
392
393    @Override
394    public boolean equals(Object obj) {
395        if (!(obj instanceof Polynomial))
396            return false;
397        Polynomial<?> that = (Polynomial<?>) obj;
398        return this._termToCoef.equals(that._termToCoef);
399    }
400
401    @Override
402    public int hashCode() {
403        return _termToCoef.hashCode();
404    }
405
406    @Override
407    public Text toText() {
408        FastTable<Term> terms = FastTable.newInstance();
409        terms.addAll(_termToCoef.keySet());
410        terms.sort();
411        TextBuilder tb = TextBuilder.newInstance();
412        for (int i=0, n = terms.size(); i < n; i++) {
413            if (i != 0) {
414                tb.append(" + ");
415            }
416            tb.append('[').append(_termToCoef.get(terms.get(i)));
417            tb.append(']').append(terms.get(i));
418        }
419        return tb.toText();
420    }
421    
422
423    /**
424     * Returns a copy of this polynomial 
425     * {@link javolution.context.AllocatorContext allocated} 
426     * by the calling thread (possibly on the stack).
427     *     
428     * @return an identical and independant copy of this polynomial.
429     */
430    public Polynomial<R> copy() {
431        Polynomial<R> p = Polynomial.newInstance();
432        for (Map.Entry<Term, R> entry : _termToCoef.entrySet()) {
433            p._termToCoef.put(entry.getKey().copy(), entry.getValue());
434        }    
435        return p;
436    }
437
438    @SuppressWarnings("unchecked")
439    private static <R extends Ring<R>> Polynomial<R> newInstance() {
440        Polynomial p = FACTORY.object();
441        p._termToCoef.clear();
442        return p;
443    }
444    
445    @SuppressWarnings("unchecked")
446    private static final ObjectFactory<Polynomial> FACTORY = new ObjectFactory<Polynomial>() {
447        protected Polynomial create() {
448            return new Polynomial();
449        }
450    };
451
452    private static final long serialVersionUID = 1L;
453
454}