polynomialRoot.js 5.3 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122
  1. import { factory } from '../../utils/factory.js';
  2. var name = 'polynomialRoot';
  3. var dependencies = ['typed', 'isZero', 'equalScalar', 'add', 'subtract', 'multiply', 'divide', 'sqrt', 'unaryMinus', 'cbrt', 'typeOf', 'im', 're'];
  4. export var createPolynomialRoot = /* #__PURE__ */factory(name, dependencies, _ref => {
  5. var {
  6. typed,
  7. isZero,
  8. equalScalar,
  9. add,
  10. subtract,
  11. multiply,
  12. divide,
  13. sqrt,
  14. unaryMinus,
  15. cbrt,
  16. typeOf,
  17. im,
  18. re
  19. } = _ref;
  20. /**
  21. * Finds the numerical values of the distinct roots of a polynomial with real or complex coefficients.
  22. * Currently operates only on linear, quadratic, and cubic polynomials using the standard
  23. * formulas for the roots.
  24. *
  25. * Syntax:
  26. *
  27. * polynomialRoot(constant, linearCoeff, quadraticCoeff, cubicCoeff)
  28. *
  29. * Examples:
  30. * // linear
  31. * math.polynomialRoot(6, 3) // [-2]
  32. * math.polynomialRoot(math.complex(6,3), 3) // [-2 - i]
  33. * math.polynomialRoot(math.complex(6,3), math.complex(2,1)) // [-3 + 0i]
  34. * // quadratic
  35. * math.polynomialRoot(2, -3, 1) // [2, 1]
  36. * math.polynomialRoot(8, 8, 2) // [-2]
  37. * math.polynomialRoot(-2, 0, 1) // [1.4142135623730951, -1.4142135623730951]
  38. * math.polynomialRoot(2, -2, 1) // [1 + i, 1 - i]
  39. * math.polynomialRoot(math.complex(1,3), math.complex(-3, -2), 1) // [2 + i, 1 + i]
  40. * // cubic
  41. * math.polynomialRoot(-6, 11, -6, 1) // [1, 3, 2]
  42. * math.polynomialRoot(-8, 0, 0, 1) // [-1 - 1.7320508075688774i, 2, -1 + 1.7320508075688774i]
  43. * math.polynomialRoot(0, 8, 8, 2) // [0, -2]
  44. * math.polynomialRoot(1, 1, 1, 1) // [-1 + 0i, 0 - i, 0 + i]
  45. *
  46. * See also:
  47. * cbrt, sqrt
  48. *
  49. * @param {... number | Complex} coeffs
  50. * The coefficients of the polynomial, starting with with the constant coefficent, followed
  51. * by the linear coefficient and subsequent coefficients of increasing powers.
  52. * @return {Array} The distinct roots of the polynomial
  53. */
  54. return typed(name, {
  55. 'number|Complex, ...number|Complex': (constant, restCoeffs) => {
  56. var coeffs = [constant, ...restCoeffs];
  57. while (coeffs.length > 0 && isZero(coeffs[coeffs.length - 1])) {
  58. coeffs.pop();
  59. }
  60. if (coeffs.length < 2) {
  61. throw new RangeError("Polynomial [".concat(constant, ", ").concat(restCoeffs, "] must have a non-zero non-constant coefficient"));
  62. }
  63. switch (coeffs.length) {
  64. case 2:
  65. // linear
  66. return [unaryMinus(divide(coeffs[0], coeffs[1]))];
  67. case 3:
  68. {
  69. // quadratic
  70. var [c, b, a] = coeffs;
  71. var denom = multiply(2, a);
  72. var d1 = multiply(b, b);
  73. var d2 = multiply(4, a, c);
  74. if (equalScalar(d1, d2)) return [divide(unaryMinus(b), denom)];
  75. var discriminant = sqrt(subtract(d1, d2));
  76. return [divide(subtract(discriminant, b), denom), divide(subtract(unaryMinus(discriminant), b), denom)];
  77. }
  78. case 4:
  79. {
  80. // cubic, cf. https://en.wikipedia.org/wiki/Cubic_equation
  81. var [d, _c, _b, _a] = coeffs;
  82. var _denom = unaryMinus(multiply(3, _a));
  83. var D0_1 = multiply(_b, _b);
  84. var D0_2 = multiply(3, _a, _c);
  85. var D1_1 = add(multiply(2, _b, _b, _b), multiply(27, _a, _a, d));
  86. var D1_2 = multiply(9, _a, _b, _c);
  87. if (equalScalar(D0_1, D0_2) && equalScalar(D1_1, D1_2)) {
  88. return [divide(_b, _denom)];
  89. }
  90. var Delta0 = subtract(D0_1, D0_2);
  91. var Delta1 = subtract(D1_1, D1_2);
  92. var discriminant1 = add(multiply(18, _a, _b, _c, d), multiply(_b, _b, _c, _c));
  93. var discriminant2 = add(multiply(4, _b, _b, _b, d), multiply(4, _a, _c, _c, _c), multiply(27, _a, _a, d, d));
  94. if (equalScalar(discriminant1, discriminant2)) {
  95. return [divide(subtract(multiply(4, _a, _b, _c), add(multiply(9, _a, _a, d), multiply(_b, _b, _b))), multiply(_a, Delta0)),
  96. // simple root
  97. divide(subtract(multiply(9, _a, d), multiply(_b, _c)), multiply(2, Delta0)) // double root
  98. ];
  99. }
  100. // OK, we have three distinct roots
  101. var Ccubed;
  102. if (equalScalar(D0_1, D0_2)) {
  103. Ccubed = Delta1;
  104. } else {
  105. Ccubed = divide(add(Delta1, sqrt(subtract(multiply(Delta1, Delta1), multiply(4, Delta0, Delta0, Delta0)))), 2);
  106. }
  107. var allRoots = true;
  108. var rawRoots = cbrt(Ccubed, allRoots).toArray().map(C => divide(add(_b, C, divide(Delta0, C)), _denom));
  109. return rawRoots.map(r => {
  110. if (typeOf(r) === 'Complex' && equalScalar(re(r), re(r) + im(r))) {
  111. return re(r);
  112. }
  113. return r;
  114. });
  115. }
  116. default:
  117. throw new RangeError("only implemented for cubic or lower-order polynomials, not ".concat(coeffs));
  118. }
  119. }
  120. });
  121. });