#!/usr/bin/perl -w

use Math::Complex;
use constant two_pi => 6.28318530717959;

unshift @ARGV, (0) x (3 - $#ARGV);
@roots = cubic(@ARGV);
print "@roots\n";

# linear($a, $b) solves the equation ax + b = 0, returning x.
#
sub linear {
    my ($a, $b) = @_;
    return unless $a;
    return -$b / $a;
}

# quadratic($a, $b, $c) solves this equation for x:
#
#   2
# ax  + bx + c = 0
#
# It returns a list of the two values of x.  Unlike
# the quadratic() shown earlier, the coefficients a, b, and c
# are allowed to be complex.
#
sub quadratic {
    my ($a, $b, $c) = @_;

    return linear($b, $c) unless $a;
    my ($sgn) = 1;
    $sgn = $b/abs($b) if $b;
    if (ref($a) || ref($b) || ref($c)) {
        my ($tmp) = Math::Complex->new(0, 0);
        $tmp = ref($b) ? ~$b : $b;
        $tmp *= sqrt($b * $b - 4 * $a * $c);
        $sgn = -1 if (ref($tmp) && $tmp->Re < 0) or $tmp < 0;
    }
    my ($tmp) = -0.5 * ($b + $sgn * sqrt($b * $b - 4 * $a * $c));
    return ($tmp / $a, $c / $tmp);
}

# cubic($a, $b, $c, $d) solves this equation for x:
#
#   3     2
# ax  + bx  + cx  + d = 0
#
# It returns a list of the three values of x.
#
# Derived from Numerical Recipes in C (Press et al.)
#
sub cubic {
    my ($a, $b, $c, $d) = @_;
    return quadratic($b, $c, $d) unless $a;
    ($a, $b, $c) = ($b / $a, $c / $a, $d / $a);
    my ($q) = ($a ** 2 - (3 * $b)) / 9;
    my ($r) = ((2 * ($a ** 3)) - (9 * $a * $b) + (27 * $c)) / 54;
    if (!ref($q) && !ref($r) && ($r ** 2) < ($q ** 3)) {
        my ($theta) = acos($r / ($q ** 1.5));
        my ($gain) = -2 * sqrt($q);
        my ($bias) = $a / 3;
        return ($gain * cos($theta / 3) - $bias,
                $gain * cos(($theta + two_pi) / 3) - $bias,
                $gain * cos(($theta - two_pi) / 3) - $bias);
    } else {
        my ($sgn) = 1;
        my ($tmp) = sqrt($r ** 2 - $q ** 3);
        my ($rconj) = $r;
        ref($rconj) && ($rconj = ~$rconj);
        $rconj *= $tmp;
        $sgn = -1 if (ref($rconj) && $rconj->Re < 0) or $rconj < 0;
        $s = Math::Complex->new($sgn, 0);
        $s = $s * $tmp + $r;
        $s **= 1/3;
        $s = -$s;
        $t = ($s ? ($q / $s) : 0);
        return ($s + $t - $a / 3,
                -0.5 * ($s+$t) + sqrt(-1) * sqrt(3)/2 * ($s-$t) - ($a/3),
                -0.5 * ($s+$t) - sqrt(-1) * sqrt(3)/2 * ($s-$t) - ($a/3));
    }
}