Solving a Cubic Function

I am trying to write a function where I can pass four integers (A, B, C, D) and get the zeros of the cubic function Ax3 + Bx2 + Cx + D.

If this were quadratic, it would be simple, just type the quadratic formula in the code. For cubic, however, I can’t seem to recall there being such an equation.

Is there a way to have the arduino do this (maybe a library?) or an arithmetic equation that would spit out the answer similar to quadratic eq?

There's a similar formula for cubics: Cubic equation - Wikipedia

Pieter

electricviolin:
I am trying to write a function where I can pass four integers (A, B, C, D) and get the zeros of the cubic function Ax3 + Bx2 + Cx + D.

If this were quadratic, it would be simple, just type the quadratic formula in the code. For cubic, however, I can’t seem to recall there being such an equation.

Is there a way to have the arduino do this (maybe a library?) or an arithmetic equation that would spit out the answer similar to quadratic eq?

something like this maybe:

#define pi 3.141593

struct real_solutions_t {
  uint8_t number_of_solutions;
  float solutions[3];
};

uint8_t solve_cubic(float a, float b, float c, float d, struct real_solutions_t *soln)
{
  float _a, _b, _c;
  float half = 1.0 / 2;
  float third = 1.0 / 3;

  if (a != 0) {
    _a = b / a;
    _b = c / a;
    _c = d / a;
  }
  else {
    return 0;  //not a cubic funtions
  }

  float p = _b - ((_a * _a) / 3);
  float q = ((2 * _a * _a * _a) / 27) - ((_a * _b) / 3) + _c;
  float disciminant = ((q * q) / 4) + ((p * p * p) / 27);

  if (disciminant > 0.000001) { // not using 0 to account for float accuracy 
    soln->number_of_solutions = 1; //one real solution
    soln->solutions[0] = pow((pow(disciminant, half) - (q / 2)), third) + pow((-1 * (pow(disciminant, half)) - (q / 2)), third) - (_a / 3);
  }
  else if (disciminant < 0 ) {
    p = pow((-1 * p), half);
    float f = (2 / pow(3, half)) * p;
    float g = (3 / 2) * pow(3, half);
    soln->number_of_solutions = 3; //3 real solutions
    soln->solutions[0] = f * sin(third * asin((g * q) / (p * p * p))) - (_a / 3);
    soln->solutions[1] = f * sin(third * asin(((g * q) / (p * p * p))) + (pi / 3)) - (_a / 3);
    soln->solutions[2] = f * cos(third * asin(((g * q) / (p * p * p))) + (pi / 6)) - (_a / 3);
  }
  else { //assumes disciminant = 0 
    soln->number_of_solutions = 2; //2 real solutions (x2 = x3, repeated root)
    if (q < 0) {
      q *= -1;
      soln->solutions[0] = (2 * pow((q / 2), third)) - (_a / 3);
      soln->solutions[1] = (-1 * pow((q / 2), third)) - (_a / 3);
    }
    else {
      soln->solutions[0] = (-2 * pow((q / 2), third)) - (_a / 3);
      soln->solutions[1] = pow((q / 2), third) - (_a / 3);
    }
  }

  return 1;

}

void setup() {
  struct real_solutions_t real_solutions;

  Serial.begin(115200);

  uint8_t i = solve_cubic(1, -4, 5, -2, &real_solutions);

  if (i) {
    Serial.print(F("Number of Real Solutions: "));
    Serial.println(real_solutions.number_of_solutions);

    if (real_solutions.number_of_solutions == 1) {
      Serial.print(F("x = "));
      Serial.println(real_solutions.solutions[0]);
    }
    else if (real_solutions.number_of_solutions == 2) {
      Serial.print(F("x1 = "));
      Serial.println(real_solutions.solutions[0]);
      Serial.print(F("x2 = x3 = "));
      Serial.println(real_solutions.solutions[1]);
    }
    else {
      Serial.print(F("x1 = "));
      Serial.println(real_solutions.solutions[0]);
      Serial.print(F("x2 = "));
      Serial.println(real_solutions.solutions[1]);
      Serial.print(F("x3 = "));
      Serial.println(real_solutions.solutions[2]);
    }
  }
  else {
    Serial.println(F("NOT A CUBIC EQUATION!"));
  }
}

void loop() {
  //nothing here
}

output:
Number of Real Solutions: 2
x1 = 2.00
x2 = x3 = 1.00

I used this material to implement the code:
http://www2.trinity.unimelb.edu.au/~rbroekst/MathX/Cubic%20Formula.pdf