Noble Gas Xenon Noble Gas Xenon - 1 month ago 33
C++ Question

Cube root of a number, C++

I would like to ask a very short question, and it is as follows: in finding the cube root of a number (both neg. and pos.) in C++, how does one restrict the output to real solutions only?
I am currently writing a program to solve a cubic with Cardano's formula, and one of the intermediate variables I am using randomly outputs the complex and real cube roots - and I only need the real roots.

(E.g. in evaluating the cube root of -0.0127378, the three roots would be 0.11677095+0.202253218i, −0.2335419, 0.11677095−0.202253218i - I wish to ignore the complex ones for substitution into a later formula)

Thank you!

EDIT: Solved it! :) I created a signum function and tweaked the sign after taking the power of the absolute value of SPrime and TPrime, so now it carries forward only the real cube root.

/* ... */

#include <iostream>
#include <cmath>
#include <complex>
#include <cstdio>
#include <cassert>

using namespace std;

int signum(std::complex<double> z)
{
if (z.real() < 0 || z.imag() < 0) return -1;
else if (z.real() >= 0 || z.imag() >= 0) return 1;
}

// POST: The function is intended to solve a cubic equation with coefficients a, b, c and d., such that
// ax^3 + bx^2 + cx + d = 0. If there exist infinitely many solutions, we output -1, i.e. if a=b=c=d=0
// (trivial solution).

void solve(std::complex<double> a, std::complex<double> b, std::complex<double> c, std::complex<double> d, std::complex<double>& x1, std::complex<double>& x2, std::complex<double>& x3)
{
complex<double> i = complex<double> (0, 1.0);
// Consider implementing Cardano's method for obtaining the solution of a degree 3 polynomial, as suggested
// We must hence define the discriminant D of such an equation through complex doubles Q and R
std::complex<double> Q;
Q = (3.0*a*c - pow(b, 2)) / (9.0*pow(a, 2));
cout << "Q=" << Q << endl;

std::complex<double> R;
R = (9.0*a*b*c - 27.0*d*pow(a, 2) - 2.0*pow(b, 3)) / (54.0*pow(a, 3));
cout << "R=" << R << endl;

std::complex<double> D;
D = pow(Q, 3) + pow(R, 2);

// Possible types of output for discriminant
if (abs(D) < 0.0)
{
cout << "The cubic has three distinct, real roots." << endl;
}

else if (abs(D) == 0.0)
{
cout << "The cubic has three real roots, at least two of which are equal." << endl;
}

else if (abs(D) > 0.0)
{
cout << "The cubic has one real root and two complex conjugate roots." << endl;
}

// Defining two further complex double variables S and T, which are required to obtain the final solution for x1, x2 and x3
std::complex<double> S;
std::complex<double> SPrime;
SPrime = R+sqrt(Q*Q*Q + R*R);
cout << "SPrime=" << SPrime << endl;

if (signum(SPrime) == -1)
{
S = (-1)*pow(abs(SPrime), 0.3333333333333);
}

else if (signum(SPrime) == 1)
{
S = pow(abs(SPrime), 0.3333333333333);
}

cout << "S=" << S << endl;

std::complex<double> T;
std::complex<double> TPrime;
TPrime = (R-sqrt(Q*Q*Q + R*R));

if (signum(TPrime) == -1)
{
T = (-1)*pow(abs(TPrime), 0.3333333333333);
}

else if (signum(TPrime) == 1)
{
T = pow(abs(TPrime), 0.3333333333333);
}

cout << "T=" << T << endl;
cout << "TPrime= " << TPrime << endl;

// Expressions for the solutions
x1 = S + T - (b/(3.0*a));
x2 = (-0.5)*(S + T) - (b/(3.0*a)) + (sqrt(3.0)*0.5)*(S - T)*i;
x3 = conj(x2);

if (abs(x1) < 0.000000000001)
{
x1 = 0;
}
}

// Driver code
int main ()
{
// Taking user input for a, b, c and d
std::complex<double> a, b, c, d, x1, x2, x3;
cout << "Please enter the coefficients of the polynomial in successive order." << endl;
cin >> a >> b >> c >> d;

solve (a, b, c, d, x1, x2, x3);
cout << x1 << ", " << x2 << ", " << x3 << "." << endl;

return 0;
}

Answer

The problem as you're stating it can be solved trivially (with real numbers the cubic root of -x is the opposite of the cubic root of x):

double cuberoot(double x) {
    if (x < 0) {
        return -pow(-x, 1.0/3.0);
    } else if (x > 0) {
        return pow(x, 1.0/3.0);
    } else {
        return 0;
    }
 }

If the input is instead in general complex and you're looking for the "most real" (principal) cubic root then you can just use std::pow with complex numbers:

std::complex<double> cuberoot(std::complex<double> z) {
    return std::pow(z, 1.0/3);
}

because it's defined already to return that result.

To compute it explicitly you can divide the argument by three and take the cubic root of the modulo:

std::complex<double> cuberoot(std::complex<double> z) {
    return std::polar(pow(std::abs(z), 1.0/3),
                      std::arg(z)/3);
}