123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238 |
- 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
- */
- var 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');
- var a2 = a * a;
- var b2 = b * b;
- var c2 = c * c;
- var d2 = d * d;
- var 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) {
- var A = a;
- var B = b / 3.0;
- var C = c / 3.0;
- var D = d;
- var AC = A * C;
- var BD = B * D;
- var B2 = B * B;
- var C2 = C * C;
- var delta1 = A * C - B2;
- var delta2 = A * D - B * C;
- var delta3 = B * D - C2;
- var discriminant = 4.0 * delta1 * delta3 - delta2 * delta2;
- var temp;
- var temp1;
- if (discriminant < 0.0) {
- var ABar;
- var CBar;
- var 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;
- }
- var s = DBar < 0.0 ? -1.0 : 1.0; // This is not Math.Sign()!
- var temp0 = -s * Math.abs(ABar) * Math.sqrt(-discriminant);
- temp1 = -DBar + temp0;
- var x = temp1 / 2.0;
- var p = x < 0.0 ? -Math.pow(-x, 1.0 / 3.0) : Math.pow(x, 1.0 / 3.0);
- var 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)];
- }
- var CBarA = delta1;
- var DBarA = -2.0 * B * delta1 + A * delta2;
- var CBarD = delta3;
- var DBarD = -D * delta2 + 2.0 * C * delta3;
- var squareRootOfDiscriminant = Math.sqrt(discriminant);
- var halfSquareRootOf3 = Math.sqrt(3.0) / 2.0;
- var theta = Math.abs(Math.atan2(A * squareRootOfDiscriminant, -DBarA) / 3.0);
- temp = 2.0 * Math.sqrt(-CBarA);
- var cosine = Math.cos(theta);
- temp1 = temp * cosine;
- var temp3 = temp * (-cosine / 2.0 - halfSquareRootOf3 * Math.sin(theta));
- var numeratorLarge = temp1 + temp3 > 2.0 * B ? temp1 - B : temp3 - B;
- var denominatorLarge = A;
- var 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));
- var numeratorSmall = -D;
- var denominatorSmall = temp1 + temp3 < 2.0 * C ? temp1 + C : temp3 + C;
- var root3 = numeratorSmall / denominatorSmall;
- var E = denominatorLarge * denominatorSmall;
- var F =
- -numeratorLarge * denominatorSmall - denominatorLarge * numeratorSmall;
- var G = numeratorLarge * numeratorSmall;
- var 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');
- var roots;
- var 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;
- var 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;
|