Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

feat: Added utility function solveCubic() #1696

Merged
merged 6 commits into from
Jun 5, 2022
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
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) {
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I wonder if this should be in utils or if we should have a dedicated math section?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

"Math sees itself as a dreamy poet. Science sees it as a supplier of specialized technical equipment. And herein we find one of the greatest paradoxes of human inquiry: These two views, both valid, are hard to reconcile. If math is an equipment supplier, why is its equipment so strangely poetic? And if math is a poet, then why is its poetry so unexpectedly useful?"

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Haha, that can be the top level doc for the math section 😄
I agree that we can have a math directory and a barrel file for it, I'm guessing we'll have quite a lot of more things going in there later.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

That can be done on a separate PR though since it will move the other file too

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));
}
}