# Program to calculate the real roots
# of a cubic equation
# 
# It based on 
# A new approach to solving the cubic Cardan's
# solution revealed
# R.W.D. Nickalls
# The mathematical gazette, 1993, vol 77, pp. 354-359
#
# Fraga  Jun 6, 2010

use strict;

die "Args: a b c d\n" if $#ARGV != 3;

my $a = $ARGV[0]; 
my $b = $ARGV[1]; 
my $c = $ARGV[2]; 
my $d = $ARGV[3]; 

my $xN = -$b/(3.0*$a);
print "xN = $xN\n";

my $yN = $xN*( $xN*($a*$xN + $b) + $c ) + $d;
print "yN = $yN\n";
my $yN2 = $yN*$yN;

my $delta2 = ($b*$b - 3.0*$a*$c)/(9.0*$a*$a); 
my $h2 = 4.0*$a*$a*$delta2*$delta2*$delta2;

print "yN^2=$yN2, h^2=$h2\n";

if ( abs($yN2 - $h2) < 1e-8 ) {
	if ( abs( $yN2 ) < 1e-8 ) {
		print "Tres raices reales iguales en\n";
		print "$xN\n";

	}
	else {
		print "Dos reales iguales y la tercera distinta\n";
		my $del = cubic_root( $yN/(2.0*$a) );
		my $r1 = $xN + $del;
		my $r3 = $xN - 2.0*$del;
		print "Dos raices en $r1\n";
		print "Una raiz en $r3\n";
	}
}
elsif ( $yN2 > $h2 ) {
	print "Una raiz real\n";
	my $v = sqrt( $yN2 - $h2 );
	my $v1 = cubic_root( (-$yN + $v)/(2*$a) );
	my $v2 = cubic_root( (-$yN - $v)/(2*$a) );

	my $r1 = $xN + $v1 + $v2;
	print "Raiz en $r1\n";
}
else {
	print "Tres raices reales distintas\n";

	my $delta = sqrt( $delta2 );
	my $c = -$yN/( 2.0*$a*$delta2*$delta );
	my $s = sqrt( 1.0 - $c*$c );
	my $theta = atan2( $s, $c )/3.0;
	print "Theta = $theta\n";
	$delta *= 2.0;
	my $r1 = $xN + $delta*cos( $theta );
	my $r2 = $xN + $delta*cos( $theta + 2.09439510266 );
	my $r3 = $xN + $delta*cos( $theta + 4.18879020533 );
	print "Raices: $r1 $r2 $r3\n";
}

sub cubic_root( $ ) 
{
	my $v = shift;

	my $signo = 1;
	if( $v < 0.0 ) {
		$signo = -1;
		$v *= -1;
	}
	my $del = $signo*exp( log($v)/3.0 );

	return( $del );
} 
