zerawiel::Poly3.java

package zerawiel;

import org.apache.commons.math3.complex.Complex;

public class Poly3 {
    public Complex a,b,c,d;
    public Complex x0,x1,x2;
    public Complex p,q;

    /**
     * Create a polynomial given coefficinents.
     * @throws Exception
     *
     */
    public Poly3(Complex a, Complex b, Complex c, Complex d) throws Exception {
        this.a = a;
        this.b = b;
        this.c = c;
        this.d = d;
        if (a.equals(Complex.ZERO)) throw new Exception();
    }

    // p = c/a-b*b/(3*a*a);
    // q = 2*b*b*b/(27*a*a*a)+d/a-b*c/(3*a*a);
    private void canonicalForm()
    {
        p = c.divide(a).subtract(b.multiply(b).divide(a.multiply(a).multiply(3)));
        q = (b.multiply(b).multiply(b).multiply(2).divide(a.multiply(a).multiply(a).multiply(27)))
            .add(d.divide(a)).subtract(b.multiply(c).divide(a.multiply(a).multiply(3)));
    }

    //x = y - b/3a
    private Complex ytox(Complex y)
    {
        return y.subtract(b.divide(a.multiply(3)));
    }

    int find_m_from_eps(Complex z, Complex[] eps) throws Exception
    {
        for (int i=0; i<3; i++)
        {
            Complex r = z.subtract(eps[i]);
            if (ComplexEx.getR(r)<1e-12)
                return i;
        }
        throw new Exception(); //it is not possible
    }

    int find_l_for_m(int m) throws Exception
    {
        switch (m)
        {
        case 0: return 0;
        case 1: return 2;
        case 2: return 1;
        default: throw new Exception(); //it is not possible
        }
    }

    public void solve() throws Exception
    {
        canonicalForm();
        if (ComplexEx.getR(p)<1e-12)
        {
            System.out.println("p=0");
            Complex roots[] = ComplexEx.rootn(3, Complex.ZERO.subtract(q));
            x0 = ytox(roots[0]);
            x1 = ytox(roots[1]);
            x2 = ytox(roots[2]);
        }
        else
        {
            System.out.println("p!=0");
            Complex temp = (q.multiply(q).add(p.multiply(p).multiply(p).multiply((double)4/27))).sqrt();
            Complex z0 = Complex.ZERO.subtract(q).add(temp);
            z0 = z0.divide(2);
            Complex []roots = ComplexEx.rootn(3, z0);
            Complex v0 = roots[0];
            Complex []roots1 = ComplexEx.rootn(3, Complex.ZERO.subtract(q).subtract(z0));
            Complex uaster = roots1[2];

            temp = v0.multiply(uaster).divide(p).multiply(-3);
            Complex[] eps = new Complex[3];
            eps[0] = new Complex(1,0);
            eps[1] = new Complex(-1.0/2,Math.sqrt(3)/2);
            eps[2] = new Complex(-1.0/2,-Math.sqrt(3)/2);
            int m = find_m_from_eps(temp,eps);
            int n = find_l_for_m(m);
            Complex u0 = eps[n].multiply(uaster);
            Complex y0 = v0.add(u0);
            Complex y1 = (eps[1].multiply(v0)).add(eps[1].multiply(eps[1]).multiply(u0));
            Complex y2 = (eps[1].multiply(eps[1]).multiply(v0)).add(eps[1].multiply(u0));

            x0 = ytox(y0);
            x1 = ytox(y1);
            x2 = ytox(y2);
        }
    }

    public void print()
    {
        System.out.format("(%f,%f)x^3 + (%f,%f)x^2 + (%f,%f)x +(%f,%f)\n",
                a.getReal(),a.getImaginary(),
                b.getReal(),b.getImaginary(),
                c.getReal(),c.getImaginary(),
                d.getReal(),d.getImaginary());
        System.out.format("p = (%f,%f), q = (%f,%f)\n",
                p.getReal(),p.getImaginary(),
                q.getReal(),q.getImaginary());
        System.out.format("x0 = (%f,%f), x1 = (%f,%f), x2 = (%f,%f)\n",
                x0.getReal(),x0.getImaginary(),
                x1.getReal(),x1.getImaginary(),
                x2.getReal(),x2.getImaginary());

    }
}