AntiDerivative.java

package jscl.math;

import jscl.math.function.*;
import jscl.math.polynomial.Monomial;
import jscl.math.polynomial.Polynomial;
import jscl.math.polynomial.UnivariatePolynomial;

import javax.annotation.Nonnull;

public class AntiDerivative {
    UnivariatePolynomial factory;
    PolynomialWithSyzygy syzygy;
    Generic result;

    AntiDerivative(Variable variable) {
        factory = (UnivariatePolynomial) Polynomial.factory(variable);
        syzygy = (PolynomialWithSyzygy) PolynomialWithSyzygy.factory(variable);
    }

    public static Generic compute(Fraction fraction, Variable variable) {
        AntiDerivative s = new AntiDerivative(variable);
        s.compute(fraction);
        return s.getValue();
    }

    public static Generic compute(Root root, Variable variable) throws NotIntegrableException {
        int d = root.degree();
        Generic a[] = root.getParameters();
        boolean b = d > 0;
        b = b && a[0].negate().isIdentity(variable);
        for (int i = 1; i < d; i++) b = b && a[i].signum() == 0;
        b = b && a[d].compareTo(JsclInteger.valueOf(1)) == 0;
        if (b) {
            return new Pow(
                    a[0].negate(),
                    new Inverse(JsclInteger.valueOf(d)).selfExpand()
            ).antiDerivative(0);
        } else {
            throw new NotIntegrableException();
        }
    }

    void compute(Fraction fraction) {
        Debug.println("antiDerivative");
        Debug.increment();
        Generic g[] = fraction.getParameters();
        Generic r[] = reduce(g[0], g[1]);
        r = divideAndRemainder(r[0], r[1]);
        Generic s = new Inverse(r[2]).selfExpand();
        Generic p = r[0].multiply(s);
        Generic a = r[1].multiply(s);
        result = p.antiDerivative(factory.variable()).add(hermite(a, g[1]));
        Debug.decrement();
    }

    Generic[] reduce(Generic n, Generic d) {
        Debug.println("reduce(" + n + ", " + d + ")");
        Polynomial pn = factory.valueOf(n);
        Polynomial pd = factory.valueOf(d);
        Polynomial gcd = pn.gcd(pd);
        return new Generic[]{
                pn.divide(gcd).genericValue(),
                pd.divide(gcd).genericValue()
        };
    }

    Generic[] divideAndRemainder(Generic n, Generic d) {
        Debug.println("divideAndRemainder(" + n + ", " + d + ")");
        Polynomial pn = syzygy.valueof(n, 0);
        Polynomial pd = syzygy.valueof(d, 1);
        PolynomialWithSyzygy pr = (PolynomialWithSyzygy) pn.remainderUpToCoefficient(pd);
        return new Generic[]{
                pr.syzygy[1].genericValue().negate(),
                pr.genericValue(),
                pr.syzygy[0].genericValue()
        };
    }

    Generic[] bezout(Generic a, Generic b) {
        Debug.println("bezout(" + a + ", " + b + ")");
        Polynomial pa = syzygy.valueof(a, 0);
        Polynomial pb = syzygy.valueof(b, 1);
        PolynomialWithSyzygy gcd = (PolynomialWithSyzygy) pa.gcd(pb);
        return new Generic[]{
                gcd.syzygy[0].genericValue(),
                gcd.syzygy[1].genericValue(),
                gcd.genericValue()
        };
    }

    Generic hermite(Generic a, Generic d) {
        Debug.println("hermite(" + a + ", " + d + ")");
        UnivariatePolynomial sd[] = ((UnivariatePolynomial) factory.valueOf(d)).squarefreeDecomposition();
        int m = sd.length - 1;
        if (m < 2) return trager(a, d);
        else {
            Generic u = sd[0].genericValue();
            for (int i = 1; i < m; i++) {
                u = u.multiply(sd[i].genericValue().pow(i));
            }
            Generic v = sd[m].genericValue();
            Generic vprime = sd[m].derivative().genericValue();
            Generic uvprime = u.multiply(vprime);
            Generic r[] = bezout(uvprime, v);
            Generic b = r[0].multiply(a);
            Generic c = r[1].multiply(a);
            Generic s = r[2];
            r = divideAndRemainder(b, v);
            b = r[1];
            c = c.multiply(r[2]).add(r[0].multiply(uvprime));
            s = new Inverse(s.multiply(r[2]).multiply(JsclInteger.valueOf(1 - m))).selfExpand();
            b = b.multiply(s);
            c = c.multiply(s);
            Generic bprime = ((UnivariatePolynomial) factory.valueOf(b)).derivative().genericValue();
            return new Fraction(b, v.pow(m - 1)).selfExpand().add(hermite(JsclInteger.valueOf(1 - m).multiply(c).subtract(u.multiply(bprime)), u.multiply(v.pow(m - 1))));
        }
    }

    Generic trager(Generic a, Generic d) {
        Debug.println("trager(" + a + ", " + d + ")");
        Variable t = new TechnicalVariable("t");
        UnivariatePolynomial pd = (UnivariatePolynomial) factory.valueOf(d);
        UnivariatePolynomial pa = (UnivariatePolynomial) factory.valueOf(a).subtract(pd.derivative().multiply(t.expressionValue()));
        UnivariatePolynomial rs[] = pd.remainderSequence(pa);
        Polynomial fact = UnivariatePolynomial.factory(t);
        for (int i = 0; i < rs.length; i++)
            if (rs[i] != null)
                rs[i] = (UnivariatePolynomial) fact.valueOf((i > 0 ? rs[i].normalize() : rs[i]).genericValue());
        UnivariatePolynomial q[] = rs[0].squarefreeDecomposition();
        int m = q.length - 1;
        Generic s = JsclInteger.valueOf(0);
        for (int i = 1; i <= m; i++) {
            for (int j = 0; j < q[i].degree(); j++) {
                Generic a2 = new Root(q[i], j).selfExpand();
                s = s.add(a2.multiply(new Ln(i == pd.degree() ? d : rs[i].substitute(a2)).selfExpand()));
            }
        }
        return s;
    }

    Generic getValue() {
        return result;
    }
}

class PolynomialWithSyzygy extends UnivariatePolynomial {
    Polynomial syzygy[] = new Polynomial[2];

    PolynomialWithSyzygy(Variable variable) {
        super(variable);
    }

    public static Polynomial factory(Variable variable) {
        return new PolynomialWithSyzygy(variable);
    }

    @Nonnull
    public Polynomial subtract(@Nonnull Polynomial that) {
        PolynomialWithSyzygy p2 = (PolynomialWithSyzygy) that;
        PolynomialWithSyzygy p = (PolynomialWithSyzygy) super.subtract(p2);
        for (int i = 0; i < syzygy.length; i++) p.syzygy[i] = syzygy[i].subtract(p2.syzygy[i]);
        return p;
    }

    public Polynomial multiply(Generic generic) {
        PolynomialWithSyzygy p = (PolynomialWithSyzygy) super.multiply(generic);
        for (int i = 0; i < syzygy.length; i++) p.syzygy[i] = syzygy[i].multiply(generic);
        return p;
    }

    public Polynomial multiply(Monomial monomial, Generic generic) {
        PolynomialWithSyzygy p = (PolynomialWithSyzygy) super.multiply(monomial, generic);
        for (int i = 0; i < syzygy.length; i++) p.syzygy[i] = syzygy[i].multiply(monomial).multiply(generic);
        return p;
    }

    public Polynomial divide(Generic generic) throws ArithmeticException {
        PolynomialWithSyzygy p = (PolynomialWithSyzygy) super.divide(generic);
        for (int i = 0; i < syzygy.length; i++) p.syzygy[i] = syzygy[i].divide(generic);
        return p;
    }

    public Polynomial remainderUpToCoefficient(Polynomial polynomial) throws ArithmeticException {
        PolynomialWithSyzygy p = this;
        PolynomialWithSyzygy q = (PolynomialWithSyzygy) polynomial;
        if (p.signum() == 0) return p;
        int d = p.degree();
        int d2 = q.degree();
        for (int i = d - d2; i >= 0; i--) {
            Generic c1 = p.get(i + d2);
            Generic c2 = q.get(d2);
            Generic c = c1.gcd(c2);
            c1 = c1.divide(c);
            c2 = c2.divide(c);
            p = (PolynomialWithSyzygy) p.multiply(c2).subtract(q.multiply(monomial(Literal.valueOf(variable, i)), c1)).normalize();
        }
        return p;
    }

    public Polynomial gcd(Polynomial polynomial) {
        Polynomial p = this;
        Polynomial q = polynomial;
        while (q.signum() != 0) {
            Polynomial r = p.remainderUpToCoefficient(q);
            p = q;
            q = r;
        }
        return p;
    }

    public Generic gcd() {
        Generic a = super.gcd();
        for (int i = 0; i < syzygy.length; i++) a = a.gcd(syzygy[i].gcd());
        return a.signum() == signum() ? a : a.negate();
    }

    public PolynomialWithSyzygy valueof(Generic generic, int n) {
        PolynomialWithSyzygy p = (PolynomialWithSyzygy) newinstance();
        p.init(generic, n);
        return p;
    }

    void init(Generic generic, int n) {
        init(generic);
        for (int i = 0; i < syzygy.length; i++)
            syzygy[i] = Polynomial.factory(variable).valueOf(JsclInteger.valueOf(i == n ? 1 : 0));
    }

    protected UnivariatePolynomial newinstance() {
        return new PolynomialWithSyzygy(variable);
    }
}