Skip to content

Commit

Permalink
feat: Added utility function solveCubic() (#1696)
Browse files Browse the repository at this point in the history
The function solves a cubic equation, and can be used for a variety of purposes, such as handling the Bezier curves.
  • Loading branch information
st-pasha authored Jun 5, 2022
1 parent 4b19a1b commit 31784ca
Show file tree
Hide file tree
Showing 2 changed files with 168 additions and 0 deletions.
59 changes: 59 additions & 0 deletions packages/flame/lib/src/utils/solve_cubic.dart
Original file line number Diff line number Diff line change
@@ -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<double> 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<double> _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;
109 changes: 109 additions & 0 deletions packages/flame/test/utils/solve_cubic_test.dart
Original file line number Diff line number Diff line change
@@ -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<double> list1, List<double> 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));
}
}

0 comments on commit 31784ca

Please sign in to comment.