Re: solve cubic equations
by vrk (Chaplain) on May 03, 2017 at 08:02 UTC

It's always fun to implement numerical algorithms, especially when a closed form solution exists. You're missing some roots, though... The fundamental theorem of algebra shows that a polynomial of degree n has exactly n roots, although some (or all) of them can be complex. Your cubic solver finds only the real roots.
By the way, for finding the roots of polynomials of higher degree, there's an iterative solver on CPAN (Math::Polynomial::Solve). You'll need it if you go higher than degree 5.
 [reply] 

Quite so. If you want the complex roots, use $t multiplied by the three cube roots of unity (1 and 1/2 ± i*sqrt(3)/2) in the first if branch.
 [reply] 
Re: solve cubic equations
by no_slogan (Deacon) on May 03, 2017 at 16:33 UTC

BTW, the some of the solutions returned by the second if branch may be the same. Mathematicians call these "multiple roots." If $u == 0, there is a triple root at $x = $s. If $t == 0, there is a single root at $x = $s  $u and a double root at x = $s + $u/2 (because cos(2/3*pi) == cos(2/3*pi) == 1/2).  [reply] 
Re: solve cubic equations
by bduggan (Pilgrim) on May 05, 2017 at 19:30 UTC

I prefer the solution involving radicals.
What's cool is that you can use square roots to solve quadratics, third roots to solve cubics, fourth roots to solve quartics, but you can't use just radicals from 5 and up.
Perl 6's builtin roots function comes in handy.
This equation comes right from the wikipedia entry, and uses the names of the variables there.
#!/usr/bin/env perl6
sub cubic(\a,\b,\c,\d) {
my \Δ0 = b²  3 × a × c;
# note: special case when Δ0 == 0
my \Δ1 = 2 * b³  9 × a × b × c + 27 × a² × d;
my \C = ( ( Δ1 + sqrt( Δ1²  4 × Δ0³ + 0i) ) / 2 ).roots(3)[0];
my \ς = 1.roots(3); # cubic roots of unity
return [0,1,2].map: > \k {
( 1 / ( 3 × a ) ) × ( b + ς[k] × C + Δ0 / ( C × ς[k] ) )
}
}
my @vals = cubic(1,10,10,10);
# test
use Test;
plan 3;
my $f = > \x { x³ + 10 * x² + 10 * x  10 };
isapprox $f( @vals[0] ), 0, 'first value';
isapprox $f( @vals[1] ), 0, 'second value';
isapprox $f( @vals[2] ), 0, 'third value';
I had a hard time doing this with a code block on perlmonks, so I made a gist instead.  [reply] 

my \C = ( ( Δ1 + sqrt( Δ1²  4 × Δ0³ + 0i) ) / 2 ).roots(3)[0];
Mathematically, your code is equivalent. (The arctan and cos are hidden inside roots.) Numerically, yours is potentially unstable. The first + in the line above might lose accuracy due to subtracting nearlyequal numbers. See Numerical Recipes, chapter 5.6. Unfortunately, sometimes we have to choose between code that looks nice and code that works.
The trouble with Perl 6 is that it's full of people who think it's a good idea to have a class called "Cool" which contains a hodgepodge of numerical and stringmanipulation methods.
 [reply] 

 Yes, it's totally possible to optimize for accuracy, and yes, some situations call for that.
 Allomorphic typing is cool.
 [reply] 

Re: solve cubic equations (Python)
by Anonymous Monk on May 03, 2017 at 15:58 UTC

This is an area where I actually prefer to use Python instead of Perl, anytime I implement MATHS!
import math
def cubic(a,b,c):
q = a*a/9  b/3
r = (a*a/27  b/6)*a + c/2
s = a/3
d = r*r  q*q*q
if d > 0:
t = (math.sqrt(d) + abs(r)) ** (1/3)
u = (t + q/t)
return s  u if r > 0 else s + u
else:
t = math.atan2(math.sqrt(d), r) / 3
u = 2 * math.sqrt(q)
return (
s  u * math.cos(t),
s  u * math.cos(t + 2/3*math.pi),
s  u * math.cos(t  2/3*math.pi)
)
print( cubic(10,10,10) )
 [reply] [d/l] 

Not sure what the advantage of Python is in this case. Your code is basically the same, except that it doesn't have any dollar signs or semicolons. I never liked Python's "x if cond else y" as a replacement for "cond ? x : y", it just doesn't seem right to put the condition in the middle of the expression. But I have to admit that Python's builtin multiple precision ints are a huge advantage for certain kinds of math programming.
 [reply] 

 [reply] 




 [reply] 

A reply falls below the community's threshold of quality. You may see it by logging in.


This code does not produce the right answer.
It produces
(9.433363736780528, 9.433363736780528, 9.433363736780528)
Adding ".0" to all the numbers in the code seems to help.
 [reply] [d/l] 

Yep that's Python, if the first assignment (declaration) is an integer the variable is typed to be an integer and all operators will only perform integer operations.
Cheers Rolf
_{(addicted to the Perl Programming Language and ☆☆☆☆ :)
Je suis Charlie!
}
 [reply] 





Thanks. The point was simply to show how using the right tool for the job leads to a cleaner solution. For all things string related, I would pick Python as a last resort.
 [reply] 