PolynomialDomain.java

package algebra.domains; 
 
import java.io.Serializable; 
 
import algebra.Utilities; 
import algebra.data.Polynomial; 
import lang.Functions.Function2; 
 
public class PolynomialDomain<T extends Serializable, S extends Ring<T>> extends CompositeDomain<Polynomial<T>, T, S> { 
 
    // Helper functions for creating domains by extension.  
 
    // Fully typed versions, for extensions in up to three indeterminates.   
 
    public static <T extends Serializable, S extends Ring<T>> PolynomialDomain<T, S> extend(S subDomain, String varName) { 
        return new PolynomialDomain<T, S>(varName, subDomain); 
    } 
 
    public static <T extends Serializable, S extends Ring<T>> PolynomialDomain<Polynomial<T>, PolynomialDomain<T, S>> extend(S subDomain, String v1, String v2) { 
        return extend(extend(subDomain, v2), v1); 
    } 
 
    public static <T extends Serializable, S extends Ring<T>> PolynomialDomain<Polynomial<Polynomial<T>>, PolynomialDomain<Polynomial<T>, PolynomialDomain<T, S>>> extend(S subDomain, String v1, String v2, String v3) { 
        return extend(extend(subDomain, v2, v3), v1); 
    } 
 
    // General extension, returning a polynomial domain (with unspecified sub-domains)  
    // when there is at least one indeterminate.   
 
    public static PolynomialDomain extend(Ring subDomain, String var1, String... vars) { 
        Ring result = subDomain; 
        for (int i = vars.length - 1; i >= 0; i--) { 
            result = extend(result, vars[i]); 
        } 
        return extend(result, var1); 
    } 
 
    public final String varName; 
 
    public PolynomialDomain(String varName, S subDomain) { 
        super((Class) Polynomial.class, subDomain); 
        this.varName = varName; 
    } 
 
    protected boolean objectEquals(Polynomial<T> a, Polynomial<T> b) { 
        return a.deg == b.deg && subDomain.equals(a.lc, b.lc) && equals(a.red, b.red); 
    } 
 
    // Note, the canonical form for zero is null; the restriction that leading  
    // coefficients are non-zero (because the rep is sparse) dictates that zero cannot  
    // consistently be represented by an instance of the Polynomial class.  
    public boolean isZero(Polynomial<T> a) { 
        return a == zero(); 
    } 
 
    public Polynomial<T> polynomial(T lc, int deg, Polynomial<T> red) { 
        if (isZero(lc)) { 
            return red; 
        } 
        return new Polynomial<T>(lc, deg, red); 
    } 
 
    public Polynomial<T> polynomial(T lc, int deg) { 
        return polynomial(lc, deg, Polynomial.ZERO); 
    } 
 
    public Polynomial<T> polynomial(T lc) { 
        return polynomial(lc, 0); 
    } 
 
    public Polynomial<T> polynomial(T[] a) { 
        Polynomial<T> result = zero(); 
        for (int i = 0; i < a.length; i++) { 
            result = polynomial(a[i], i, result); 
        } 
        return result; 
    } 
 
    public Polynomial<T> promote(T a) { 
        return polynomial(a); 
    } 
 
    public Polynomial<T> generator() { 
        return polynomial(subDomain.one(), 1); 
    } 
 
    public Polynomial<T> neg(Polynomial<T> a) { 
        if (a == zero()) { 
            return a; 
        } 
        return polynomial(neg(a.lc), a.deg, neg(a.red)); 
    } 
 
    public Polynomial<T> sum(Polynomial<T> a, Polynomial<T> b) { 
        if (a == zero()) { 
            return b; 
        } 
        if (b == zero()) { 
            return a; 
        } 
        int n1 = a.deg; 
        int n2 = b.deg; 
        if (n1 < n2) { 
            return polynomial(b.lc, b.deg, sum(a, b.red)); 
        } 
        if (n2 < n1) { 
            return polynomial(a.lc, a.deg, sum(a.red, b)); 
        } 
        return polynomial(sum(a.lc, b.lc), a.deg, sum(a.red, b.red)); 
    } 
 
    // Multiply by an element of the coeficient domain.  
    public Polynomial<T> prd(T u, Polynomial<T> v) { 
        if (v == zero()) { 
            return v; 
        } 
        return polynomial(prd(u, v.lc), v.deg, prd(u, v.red)); 
    } 
 
    public Polynomial<T> prd(Polynomial<T> u, Polynomial<T> v) { 
        if (u == zero()) { 
            return u; 
        } 
        if (v == zero()) { 
            return v; 
        } 
        return polynomial(prd(u.lc, v.lc), u.deg + v.deg, sum(prd(u.red, v), prd(polynomial(u.lc, u.deg), v.red))); 
    } 
 
    // Divide by an element of the coeficient domain.  
    public Polynomial<T> div(Polynomial<T> u, T v) { 
        if (u == zero()) { 
            return u; 
        } 
        return polynomial(div(u.lc, v), u.deg, div(u.red, v)); 
    } 
 
    public <R> R divide(Function2<Polynomial<T>, Polynomial<T>, R> f, Polynomial<T> u, Polynomial<T> v) { 
        if (v == zero()) { 
            throw new ArithmeticException("Divide by zero"); 
        } 
        if (u == zero() || u.deg < v.deg) { 
            return f.of(zero(), u); 
        } 
        int delta = u.deg - v.deg; 
        T r = div(u.lc, v.lc); 
        QuoAndRem<Polynomial<T>> tmp = divide(QUO_AND_REM, sub(u.red, prd(polynomial(r, delta), v.red)), v); 
        return f.of(polynomial(r, delta, tmp.quo), tmp.rem); 
    } 
 
    public Polynomial<T> quo(Polynomial<T> a, Polynomial<T> b) { 
        return divide(QUO, a, b); 
    } 
 
    public Polynomial<T> rem(Polynomial<T> a, Polynomial<T> b) { 
        return divide(REM, a, b); 
    } 
 
    // This will normally throw when the coef domain is not a field (as in the case of  
    // multivatrate polynomials).  
    public Polynomial<T> div(Polynomial<T> a, Polynomial<T> b) throws ArithmeticException { 
        return divide(DIVIDE_OR_FAIL, a, b); 
    } 
 
    // u.deg >= v.deg 
    private T subresultant_rec(Polynomial<T> u, Polynomial<T> v, T g, T h) { 
        if (v == zero()) { 
            return h; 
        } 
        int delta = u.deg - v.deg; 
        T nh = div(pow(v.lc, delta), pow(h, delta - 1)); 
//        for now, leave out out the optimisation below         
//        if (v.deg == 0) { return nh; } 
        Polynomial<T> pseudoRem = rem(prd(pow(v.lc, delta + 1), u), v); 
        return subresultant_rec(v, div(pseudoRem, prd(g, pow(h, delta))), v.lc, nh); 
    } 
 
    public T subresultant(Polynomial<T> u, Polynomial<T> v) { 
        T one = subDomain.one(); 
        return degree(u) > degree(v) ? subresultant_rec(u, v, one, one) : subresultant_rec(v, u, one, one); 
    } 
 
    public Polynomial<T> parse(String s) { 
        return Utilities.polynomial(this, s); 
    } 
 
    public int degree(Polynomial<T> p) { 
        return p == zero() ? 0 : p.deg; 
    } 
 
    public T[] coeficients(Polynomial<T> p) { 
        T[] result = new VectorDomain<T>(subDomain, degree(p) + 1).zero(); // need to fill with zeros as rep is sparse.  
        for (Polynomial<T> r = p; r != zero(); r = r.red) { 
            result[r.deg] = r.lc; 
        } 
        return result; 
    } 
 
    // inefficient 
    public String toString(Polynomial<T> a) { 
        if (a == zero()) { 
            return "0"; 
        } 
        String coefs = toString(a.lc); 
        String coef = coefs.equals("1") && a.deg != 0 ? "" : subDomain instanceof CompositeDomain ? "(" + coefs + ")" : coefs; 
        String pow = a.deg == 0 ? "" : a.deg == 1 ? varName : varName + "^" + a.deg; 
        String mul = coef != "" && pow != "" ? "*" : ""; 
        String rest = a.red == zero() ? "" : " + " + toString(a.red); 
        return coef + mul + pow + rest; 
    } 
}