From 31784ca0b05082042003f847be2b4004da83edb6 Mon Sep 17 00:00:00 2001 From: Pasha Stetsenko Date: Sun, 5 Jun 2022 04:41:46 -0700 Subject: [PATCH] feat: Added utility function solveCubic() (#1696) The function solves a cubic equation, and can be used for a variety of purposes, such as handling the Bezier curves. --- packages/flame/lib/src/utils/solve_cubic.dart | 59 ++++++++++ .../flame/test/utils/solve_cubic_test.dart | 109 ++++++++++++++++++ 2 files changed, 168 insertions(+) create mode 100644 packages/flame/lib/src/utils/solve_cubic.dart create mode 100644 packages/flame/test/utils/solve_cubic_test.dart diff --git a/packages/flame/lib/src/utils/solve_cubic.dart b/packages/flame/lib/src/utils/solve_cubic.dart new file mode 100644 index 00000000000..6d40bdd970e --- /dev/null +++ b/packages/flame/lib/src/utils/solve_cubic.dart @@ -0,0 +1,59 @@ +import 'dart:math'; + +import 'package:flame/src/utils/solve_quadratic.dart'; + +/// Solves cubic equation `ax³ + bx² + cx + d == 0`. +/// +/// Depending on the coefficients, either 1 or 3 real solutions may be returned. +/// In degenerate cases, when some of the roots of the equation coincide, we +/// return all such roots without deduplication. +/// +/// If coefficient [a] is equal to zero, then we solve the equation as a +/// quadratic one (see [solveQuadratic]). +List solveCubic(double a, double b, double c, double d) { + if (a == 0) { + return solveQuadratic(b, c, d); + } + if (b == 0) { + return _solveDepressedCubic(c / a, d / a); + } else { + final s = b / (3 * a); + final p = c / a - 3 * s * s; + final q = d / a - (p + s * s) * s; + return _solveDepressedCubic(p, q).map((t) => t - s).toList(); + } +} + +/// Solves cubic equation `x³ + px + q == 0`. +List _solveDepressedCubic(double p, double q) { + final discriminant = q * q / 4 + p * p * p / 27; + // If the discriminant is very close to zero, then we will treat this as if + // it was equal to zero. + if (discriminant.abs() < discriminantEpsilon) { + final x1 = _cubicRoot(q / 2); + final x2 = -2 * x1; + return [x1, x1, x2]; + } else if (discriminant > 0) { + final w = _cubicRoot(q.abs() / 2 + sqrt(discriminant)); + return [(p / (3 * w) - w) * q.sign]; + } else { + final f = 2 * sqrt(-p / 3); + final v = acos(3 * q / (f * p)) / 3; + final x0 = f * cos(v); + final x1 = f * cos(v - 1 / 3 * tau); + final x2 = f * cos(v - 2 / 3 * tau); + return [x0, x1, x2]; + } +} + +double _cubicRoot(double x) { + // Note: `pow(x, y)` function cannot handle negative values of `x` + if (x >= 0) { + return pow(x, 1 / 3).toDouble(); + } else { + return -pow(-x, 1 / 3).toDouble(); + } +} + +const discriminantEpsilon = 1e-15; +const tau = 2 * pi; diff --git a/packages/flame/test/utils/solve_cubic_test.dart b/packages/flame/test/utils/solve_cubic_test.dart new file mode 100644 index 00000000000..b4a72e5bda8 --- /dev/null +++ b/packages/flame/test/utils/solve_cubic_test.dart @@ -0,0 +1,109 @@ +import 'package:flame/src/utils/solve_cubic.dart'; +import 'package:flame_test/flame_test.dart'; +import 'package:test/test.dart'; + +void main() { + group('solveCubic', () { + const repeatCount = 3; + + testRandom( + 'solve equation with 3 roots', + (rnd) { + final x1 = rnd.nextDouble() * 5 - 1; + final x2 = rnd.nextDouble() * 0.3 - 0.2; + final x3 = rnd.nextDouble() * 2 - 1; + // a(x - x1)(x - x2)(x - x3) == 0 + final a = rnd.nextDouble() + 1e-6; + final b = -(x1 + x2 + x3) * a; + final c = (x1 * x2 + x2 * x3 + x3 * x1) * a; + final d = -x1 * x2 * x3 * a; + final solutions = solveCubic(a, b, c, d); + if (((x1 - x2) * (x2 - x3) * (x3 - x1)).abs() > 1e-5) { + check(solutions, [x1, x2, x3]); + } + }, + repeatCount: repeatCount, + ); + + testRandom( + 'solve equation with 2 roots', + (rnd) { + final x1 = rnd.nextDouble() * 5 - 1; + final x2 = rnd.nextDouble() * 0.3 - 0.2; + // a(x - x1)(x - x2)² == 0 + final a = rnd.nextDouble(); + final b = -(x1 + x2 + x2) * a; + final c = (2 * x1 * x2 + x2 * x2) * a; + final d = -x1 * x2 * x2 * a; + final solutions = solveCubic(a, b, c, d); + if (solutions.length == 1) { + check(solutions, [x1]); + } else { + check(solutions, [x1, x2, x2]); + } + }, + repeatCount: repeatCount, + ); + + test('solve equation with 1 triple root', () { + check(solveCubic(1, -3, 3, -1), [1, 1, 1]); + check(solveCubic(10, 30, 30, 10), [-1, -1, -1]); + const x = 2.78; + check(solveCubic(1, -3 * x, 3 * x * x, -x * x * x), [x, x, x]); + }); + + testRandom( + 'solve equation with 1 real root', + (rnd) { + final x1 = rnd.nextDouble() * 5 - 1; + final x2 = rnd.nextDouble() * 0.3 - 0.2; + // a(x - x1)((x - x2)² + 0.5) == 0 + final a = rnd.nextDouble(); + final b = -(x1 + 2 * x2) * a; + final c = (2 * x1 * x2 + x2 * x2 + 0.5) * a; + final d = -x1 * (x2 * x2 + 0.5) * a; + final solutions = solveCubic(a, b, c, d); + check(solutions, [x1]); + }, + repeatCount: repeatCount, + ); + + test('solve degenerate equation', () { + check(solveCubic(0, 1, 2, 1), [-1, -1]); + check(solveCubic(0, 0, 1, 3), [-3]); + }); + + test('solve depressed equation', () { + check(solveCubic(1, 0, -7, 6), [1, 2, -3]); + check(solveCubic(0.1, 0, -0.7, -0.6), [-1, -2, 3]); + }); + + testRandom( + 'solve random equation', + (rnd) { + final a = rnd.nextDouble(); + final b = rnd.nextDouble() * 2; + final c = rnd.nextDouble() * 4; + final d = rnd.nextDouble() * 6; + final solutions = solveCubic(a, b, c, d); + for (final x in solutions) { + expect(a * x * x * x + b * x * x + c * x + d, closeTo(0, 1e-6)); + } + }, + repeatCount: repeatCount, + ); + }); +} + +void check(List list1, List list2) { + expect( + list1.length, + equals(list2.length), + reason: 'solutions are: $list1 vs $list2', + ); + list1.sort(); + list2.sort(); + for (var i = 0; i < list1.length; i++) { + expect(list1[i], closeTo(list2[i], 1e-6)); + } +}