com.opengamma.analytics.math.rootfinding.CubicRootFinder.java Source code

Introduction

Here is the source code for com.opengamma.analytics.math.rootfinding.CubicRootFinder.java

Source

/**
 * Copyright (C) 2009 - present by OpenGamma Inc. and the OpenGamma group of companies
 *
 * Please see distribution for license.
 */
package com.opengamma.analytics.math.rootfinding;

import org.apache.commons.lang.Validate;

import com.opengamma.analytics.math.function.RealPolynomialFunction1D;
import com.opengamma.analytics.math.number.ComplexNumber;
import com.opengamma.util.CompareUtils;

/**
 * Class that calculates the roots of a cubic equation. 
 * <p>
 * As the polynomial has real coefficients, the roots of the cubic can be found using the method described
 * <a href="http://mathworld.wolfram.com/CubicFormula.html">here</a>.
 */
public class CubicRootFinder implements Polynomial1DRootFinder<ComplexNumber> {
    private static final double TWO_PI = 2 * Math.PI;

    /**
     * {@inheritDoc}
     * @throws IllegalArgumentException If the function is not cubic
     */
    @Override
    public ComplexNumber[] getRoots(final RealPolynomialFunction1D function) {
        Validate.notNull(function, "function");
        final double[] coefficients = function.getCoefficients();
        Validate.isTrue(coefficients.length == 4, "Function is not a cubic");
        final double divisor = coefficients[3];
        final double a = coefficients[2] / divisor;
        final double b = coefficients[1] / divisor;
        final double c = coefficients[0] / divisor;
        final double aSq = a * a;
        final double q = (aSq - 3 * b) / 9;
        final double r = (2 * a * aSq - 9 * a * b + 27 * c) / 54;
        final double rSq = r * r;
        final double qCb = q * q * q;
        final double constant = a / 3;
        if (rSq < qCb) {
            final double mult = -2 * Math.sqrt(q);
            final double theta = Math.acos(r / Math.sqrt(qCb));
            return new ComplexNumber[] { new ComplexNumber(mult * Math.cos(theta / 3) - constant, 0),
                    new ComplexNumber(mult * Math.cos((theta + TWO_PI) / 3) - constant, 0),
                    new ComplexNumber(mult * Math.cos((theta - TWO_PI) / 3) - constant, 0) };
        }
        final double s = -Math.signum(r) * Math.cbrt(Math.abs(r) + Math.sqrt(rSq - qCb));
        final double t = CompareUtils.closeEquals(s, 0, 1e-16) ? 0 : q / s;
        final double sum = s + t;
        final double real = -0.5 * sum - constant;
        final double imaginary = Math.sqrt(3) * (s - t) / 2;
        return new ComplexNumber[] { new ComplexNumber(sum - constant, 0), new ComplexNumber(real, imaginary),
                new ComplexNumber(real, -imaginary) };
    }
}