summaryrefslogtreecommitdiffstats
path: root/libdimension
diff options
context:
space:
mode:
Diffstat (limited to 'libdimension')
-rw-r--r--libdimension/polynomial.c6
-rw-r--r--libdimension/tests/polynomial.c192
2 files changed, 163 insertions, 35 deletions
diff --git a/libdimension/polynomial.c b/libdimension/polynomial.c
index 9491861..3f64091 100644
--- a/libdimension/polynomial.c
+++ b/libdimension/polynomial.c
@@ -358,9 +358,9 @@ dmnsn_solve_cubic(double poly[4], double x[3])
x[0] = -bdiv3;
} else {
/* Two real roots; one duplicate */
- double t1 = 3.0*q/p, t2 = -t1/2.0;
- x[0] = dmnsn_max(t1, t2);
- x[1] = dmnsn_min(t1, t2);
+ double t1 = -(3.0*q)/(2.0*p), t2 = -2.0*t1;
+ x[0] = dmnsn_max(t1, t2) - bdiv3;
+ x[1] = dmnsn_min(t1, t2) - bdiv3;
if (x[1] >= dmnsn_epsilon)
return 2;
}
diff --git a/libdimension/tests/polynomial.c b/libdimension/tests/polynomial.c
index 335fd27..b734cc2 100644
--- a/libdimension/tests/polynomial.c
+++ b/libdimension/tests/polynomial.c
@@ -24,59 +24,187 @@
#include "tests.h"
#include "../polynomial.c"
+#include <stdarg.h>
-/* poly[] = 2*(x + 1)*(x - 1.2345)*(x - 2.3456)*(x - 5)*(x - 100) */
-static const double poly[6] = {
- [5] = 2.0,
- [4] = -215.1602,
- [3] = 1540.4520864,
- [2] = -2430.5727856,
- [1] = -1292.541872,
- [0] = 2895.6432,
-};
+#define DMNSN_CLOSE_ENOUGH 1.0e-6
-static double roots[5];
-static size_t nroots;
+static void
+dmnsn_assert_roots(const double poly[], size_t degree, size_t nroots_ex, ...)
+{
+ double roots[degree - 1];
+ size_t nroots = dmnsn_polynomial_solve(poly, degree, roots);
+ ck_assert_int_eq(nroots, nroots_ex);
+
+ va_list ap;
+ va_start(ap, nroots_ex);
+ for (size_t i = 0; i < nroots; ++i) {
+ double root_ex = va_arg(ap, double);
+ bool found = false;
+ for (size_t j = 0; j < nroots; ++j) {
+ if (fabs(root_ex - roots[j]) >= dmnsn_epsilon) {
+ continue;
+ }
+ found = true;
+
+ double evroot = dmnsn_polynomial_evaluate(poly, degree, roots[j]);
+ ck_assert(fabs(evroot) < DMNSN_CLOSE_ENOUGH);
+
+ double evmin = dmnsn_polynomial_evaluate(poly, degree, roots[j] - dmnsn_epsilon);
+ double evmax = dmnsn_polynomial_evaluate(poly, degree, roots[j] + dmnsn_epsilon);
+ ck_assert(fabs(evroot) <= fabs(evmin) && fabs(evroot) <= fabs(evmax));
+
+ break;
+ }
+
+ if (!found) {
+ for (size_t j = 0; j < nroots; ++j) {
+ fprintf(stderr, "roots[%zu] == %.17g\n", j, roots[j]);
+ }
+ ck_abort_msg("Expected root %.17g not found", root_ex);
+ }
+ }
+ va_end(ap);
+}
-#define DMNSN_CLOSE_ENOUGH 1.0e-6
-DMNSN_TEST_SETUP(polynomial)
+DMNSN_TEST(linear, no_positive_roots)
{
- nroots = dmnsn_polynomial_solve(poly, 5, roots);
+ /* poly[] = x + 1 */
+ static const double poly[] = {
+ [1] = 1.0,
+ [0] = 1.0,
+ };
+ dmnsn_assert_roots(poly, 1, 0);
}
-DMNSN_TEST(polynomial, finds_positive_roots)
+DMNSN_TEST(linear, one_root)
{
- ck_assert_int_eq(nroots, 4);
+ /* poly[] = x - 1 */
+ static const double poly[] = {
+ [1] = 1.0,
+ [0] = -1.0,
+ };
+ dmnsn_assert_roots(poly, 1, 1, 1.0);
}
-DMNSN_TEST(polynomial, local_min_roots)
+
+DMNSN_TEST(quadratic, no_roots)
{
- for (size_t i = 0; i < nroots; ++i) {
- double evmin = dmnsn_polynomial_evaluate(poly, 5, roots[i] - dmnsn_epsilon);
- double ev = dmnsn_polynomial_evaluate(poly, 5, roots[i]);
- double evmax = dmnsn_polynomial_evaluate(poly, 5, roots[i] + dmnsn_epsilon);
- ck_assert(fabs(ev) < fabs(evmin) && fabs(ev) < fabs(evmax));
- }
+ /* poly[] = x^2 + 1 */
+ static const double poly[] = {
+ [2] = 1.0,
+ [1] = 0.0,
+ [0] = 1.0,
+ };
+ dmnsn_assert_roots(poly, 2, 0);
}
-DMNSN_TEST(polynomial, accurate_roots)
+DMNSN_TEST(quadratic, no_positive_roots)
{
- for (size_t i = 0; i < nroots; ++i) {
- double ev = dmnsn_polynomial_evaluate(poly, 5, roots[i]);
- ck_assert(fabs(ev) < DMNSN_CLOSE_ENOUGH);
- }
+ /* poly[] = (x + 1)^2 */
+ static const double poly[] = {
+ [2] = 1.0,
+ [1] = 2.0,
+ [0] = 1.0,
+ };
+ dmnsn_assert_roots(poly, 2, 0);
+}
+
+DMNSN_TEST(quadratic, one_positive_root)
+{
+ /* poly[] = (x + 1)*(x - 1) */
+ static const double poly[] = {
+ [2] = 1.0,
+ [1] = 0.0,
+ [0] = -1.0,
+ };
+ dmnsn_assert_roots(poly, 2, 1, 1.0);
+}
+
+DMNSN_TEST(quadratic, two_roots)
+{
+ /* poly[] = (x - 1.2345)*(x - 2.3456) */
+ static const double poly[] = {
+ [2] = 1.0,
+ [1] = -3.5801,
+ [0] = 2.8956432,
+ };
+ dmnsn_assert_roots(poly, 2, 2, 1.2345, 2.3456);
+}
+
+
+DMNSN_TEST(cubic, no_positive_roots)
+{
+ /* poly[] = x^3 + 1 */
+ static const double poly[] = {
+ [3] = 1.0,
+ [2] = 0.0,
+ [1] = 0.0,
+ [0] = 1.0,
+ };
+ dmnsn_assert_roots(poly, 3, 0);
+}
+
+DMNSN_TEST(cubic, one_root)
+{
+ /* poly[] = x^3 - 1 */
+ static const double poly[] = {
+ [3] = 1.0,
+ [2] = 0.0,
+ [1] = 0.0,
+ [0] = -1.0,
+ };
+ dmnsn_assert_roots(poly, 3, 1, 1.0);
+}
+
+DMNSN_TEST(cubic, two_roots)
+{
+ /* poly[] = (x - 1)*(x - 4)^2 **/
+ static const double poly[] = {
+ [3] = 1.0,
+ [2] = -9.0,
+ [1] = 24.0,
+ [0] = -16.0,
+ };
+ dmnsn_assert_roots(poly, 3, 2, 1.0, 4.0);
+}
+
+DMNSN_TEST(cubic, three_roots)
+{
+ /* poly[] = (x - 1.2345)*(x - 2.3456)*(x - 100) */
+ static const double poly[] = {
+ [3] = 1.0,
+ [2] = -103.5801,
+ [1] = 360.9056432,
+ [0] = -289.56432,
+ };
+ dmnsn_assert_roots(poly, 3, 3, 1.2345, 2.3456, 100.0);
+}
+
+
+DMNSN_TEST(quintic, four_roots)
+{
+ /* poly[] = 2*(x + 1)*(x - 1.2345)*(x - 2.3456)*(x - 5)*(x - 100) */
+ static const double poly[] = {
+ [5] = 2.0,
+ [4] = -215.1602,
+ [3] = 1540.4520864,
+ [2] = -2430.5727856,
+ [1] = -1292.541872,
+ [0] = 2895.6432,
+ };
+ dmnsn_assert_roots(poly, 5, 4, 1.2345, 2.3456, 5.0, 100.0);
}
/* repeated_root[] = (x - 1)^6 */
static const double repeated_root[7] = {
- [6] = 1.0,
- [5] = -6.0,
+ [6] = 1.0,
+ [5] = -6.0,
[4] = 15.0,
[3] = -20.0,
[2] = 15.0,
- [1] = -6.0,
- [0] = 1.0,
+ [1] = -6.0,
+ [0] = 1.0,
};
DMNSN_TEST(stability, equal_bounds)
generated by cgit v1.2.3 (git 2.39.1) at 2024-11-22 04:34:47 +0000