CubicRealPolynomial.js

import DeveloperError from "./DeveloperError.js";
import QuadraticRealPolynomial from "./QuadraticRealPolynomial.js";

/**
 * Defines functions for 3rd order polynomial functions of one variable with only real coefficients.
 *
 * @namespace CubicRealPolynomial
 */
const CubicRealPolynomial = {};

/**
 * Provides the discriminant of the cubic equation from the supplied coefficients.
 *
 * @param {Number} a The coefficient of the 3rd order monomial.
 * @param {Number} b The coefficient of the 2nd order monomial.
 * @param {Number} c The coefficient of the 1st order monomial.
 * @param {Number} d The coefficient of the 0th order monomial.
 * @returns {Number} The value of the discriminant.
 */
CubicRealPolynomial.computeDiscriminant = function (a, b, c, d) {
  //>>includeStart('debug', pragmas.debug);
  if (typeof a !== "number") {
    throw new DeveloperError("a is a required number.");
  }
  if (typeof b !== "number") {
    throw new DeveloperError("b is a required number.");
  }
  if (typeof c !== "number") {
    throw new DeveloperError("c is a required number.");
  }
  if (typeof d !== "number") {
    throw new DeveloperError("d is a required number.");
  }
  //>>includeEnd('debug');

  const a2 = a * a;
  const b2 = b * b;
  const c2 = c * c;
  const d2 = d * d;

  const discriminant =
    18.0 * a * b * c * d +
    b2 * c2 -
    27.0 * a2 * d2 -
    4.0 * (a * c2 * c + b2 * b * d);
  return discriminant;
};

function computeRealRoots(a, b, c, d) {
  const A = a;
  const B = b / 3.0;
  const C = c / 3.0;
  const D = d;

  const AC = A * C;
  const BD = B * D;
  const B2 = B * B;
  const C2 = C * C;
  const delta1 = A * C - B2;
  const delta2 = A * D - B * C;
  const delta3 = B * D - C2;

  const discriminant = 4.0 * delta1 * delta3 - delta2 * delta2;
  let temp;
  let temp1;

  if (discriminant < 0.0) {
    let ABar;
    let CBar;
    let DBar;

    if (B2 * BD >= AC * C2) {
      ABar = A;
      CBar = delta1;
      DBar = -2.0 * B * delta1 + A * delta2;
    } else {
      ABar = D;
      CBar = delta3;
      DBar = -D * delta2 + 2.0 * C * delta3;
    }

    const s = DBar < 0.0 ? -1.0 : 1.0; // This is not Math.Sign()!
    const temp0 = -s * Math.abs(ABar) * Math.sqrt(-discriminant);
    temp1 = -DBar + temp0;

    const x = temp1 / 2.0;
    const p = x < 0.0 ? -Math.pow(-x, 1.0 / 3.0) : Math.pow(x, 1.0 / 3.0);
    const q = temp1 === temp0 ? -p : -CBar / p;

    temp = CBar <= 0.0 ? p + q : -DBar / (p * p + q * q + CBar);

    if (B2 * BD >= AC * C2) {
      return [(temp - B) / A];
    }

    return [-D / (temp + C)];
  }

  const CBarA = delta1;
  const DBarA = -2.0 * B * delta1 + A * delta2;

  const CBarD = delta3;
  const DBarD = -D * delta2 + 2.0 * C * delta3;

  const squareRootOfDiscriminant = Math.sqrt(discriminant);
  const halfSquareRootOf3 = Math.sqrt(3.0) / 2.0;

  let theta = Math.abs(Math.atan2(A * squareRootOfDiscriminant, -DBarA) / 3.0);
  temp = 2.0 * Math.sqrt(-CBarA);
  let cosine = Math.cos(theta);
  temp1 = temp * cosine;
  let temp3 = temp * (-cosine / 2.0 - halfSquareRootOf3 * Math.sin(theta));

  const numeratorLarge = temp1 + temp3 > 2.0 * B ? temp1 - B : temp3 - B;
  const denominatorLarge = A;

  const root1 = numeratorLarge / denominatorLarge;

  theta = Math.abs(Math.atan2(D * squareRootOfDiscriminant, -DBarD) / 3.0);
  temp = 2.0 * Math.sqrt(-CBarD);
  cosine = Math.cos(theta);
  temp1 = temp * cosine;
  temp3 = temp * (-cosine / 2.0 - halfSquareRootOf3 * Math.sin(theta));

  const numeratorSmall = -D;
  const denominatorSmall = temp1 + temp3 < 2.0 * C ? temp1 + C : temp3 + C;

  const root3 = numeratorSmall / denominatorSmall;

  const E = denominatorLarge * denominatorSmall;
  const F =
    -numeratorLarge * denominatorSmall - denominatorLarge * numeratorSmall;
  const G = numeratorLarge * numeratorSmall;

  const root2 = (C * F - B * G) / (-B * F + C * E);

  if (root1 <= root2) {
    if (root1 <= root3) {
      if (root2 <= root3) {
        return [root1, root2, root3];
      }
      return [root1, root3, root2];
    }
    return [root3, root1, root2];
  }
  if (root1 <= root3) {
    return [root2, root1, root3];
  }
  if (root2 <= root3) {
    return [root2, root3, root1];
  }
  return [root3, root2, root1];
}

/**
 * Provides the real valued roots of the cubic polynomial with the provided coefficients.
 *
 * @param {Number} a The coefficient of the 3rd order monomial.
 * @param {Number} b The coefficient of the 2nd order monomial.
 * @param {Number} c The coefficient of the 1st order monomial.
 * @param {Number} d The coefficient of the 0th order monomial.
 * @returns {Number[]} The real valued roots.
 */
CubicRealPolynomial.computeRealRoots = function (a, b, c, d) {
  //>>includeStart('debug', pragmas.debug);
  if (typeof a !== "number") {
    throw new DeveloperError("a is a required number.");
  }
  if (typeof b !== "number") {
    throw new DeveloperError("b is a required number.");
  }
  if (typeof c !== "number") {
    throw new DeveloperError("c is a required number.");
  }
  if (typeof d !== "number") {
    throw new DeveloperError("d is a required number.");
  }
  //>>includeEnd('debug');

  let roots;
  let ratio;
  if (a === 0.0) {
    // Quadratic function: b * x^2 + c * x + d = 0.
    return QuadraticRealPolynomial.computeRealRoots(b, c, d);
  } else if (b === 0.0) {
    if (c === 0.0) {
      if (d === 0.0) {
        // 3rd order monomial: a * x^3 = 0.
        return [0.0, 0.0, 0.0];
      }

      // a * x^3 + d = 0
      ratio = -d / a;
      const root =
        ratio < 0.0 ? -Math.pow(-ratio, 1.0 / 3.0) : Math.pow(ratio, 1.0 / 3.0);
      return [root, root, root];
    } else if (d === 0.0) {
      // x * (a * x^2 + c) = 0.
      roots = QuadraticRealPolynomial.computeRealRoots(a, 0, c);

      // Return the roots in ascending order.
      if (roots.Length === 0) {
        return [0.0];
      }
      return [roots[0], 0.0, roots[1]];
    }

    // Deflated cubic polynomial: a * x^3 + c * x + d= 0.
    return computeRealRoots(a, 0, c, d);
  } else if (c === 0.0) {
    if (d === 0.0) {
      // x^2 * (a * x + b) = 0.
      ratio = -b / a;
      if (ratio < 0.0) {
        return [ratio, 0.0, 0.0];
      }
      return [0.0, 0.0, ratio];
    }
    // a * x^3 + b * x^2 + d = 0.
    return computeRealRoots(a, b, 0, d);
  } else if (d === 0.0) {
    // x * (a * x^2 + b * x + c) = 0
    roots = QuadraticRealPolynomial.computeRealRoots(a, b, c);

    // Return the roots in ascending order.
    if (roots.length === 0) {
      return [0.0];
    } else if (roots[1] <= 0.0) {
      return [roots[0], roots[1], 0.0];
    } else if (roots[0] >= 0.0) {
      return [0.0, roots[0], roots[1]];
    }
    return [roots[0], 0.0, roots[1]];
  }

  return computeRealRoots(a, b, c, d);
};
export default CubicRealPolynomial;