Fix T83023: incorrect shape of cyclic F-Curve with only two points.

The equation solver didn't handle the one unknown case correctly.

3 files changed, 87 insertions, 1 deletions

diff --git a/source/blender/blenlib/CMakeLists.txt b/source/blender/blenlib/CMakeLists.txt
index 473cc4d65f3..bcc8c181af4 100644
--- a/source/blender/blenlib/CMakeLists.txt
+++ b/source/blender/blenlib/CMakeLists.txt
@@ -400,6 +400,7 @@ if(WITH_GTESTS)
     tests/BLI_math_color_test.cc
     tests/BLI_math_geom_test.cc
     tests/BLI_math_matrix_test.cc
+    tests/BLI_math_solvers_test.cc
     tests/BLI_math_vector_test.cc
     tests/BLI_memiter_test.cc
     tests/BLI_memory_utils_test.cc
diff --git a/source/blender/blenlib/intern/math_solvers.c b/source/blender/blenlib/intern/math_solvers.c
index cda3d9b66a2..e366d834fc4 100644
--- a/source/blender/blenlib/intern/math_solvers.c
+++ b/source/blender/blenlib/intern/math_solvers.c
@@ -137,9 +137,24 @@ bool BLI_tridiagonal_solve_cyclic(
     return false;
   }
 
+  /* Degenerate case not handled correctly by the generic formula. */
+  if (count == 1) {
+    r_x[0] = d[0] / (a[0] + b[0] + c[0]);
+
+    return isfinite(r_x[0]);
+  }
+
+  /* Degenerate case that works but can be simplified. */
+  if (count == 2) {
+    float a2[2] = {0, a[1] + c[1]};
+    float c2[2] = {a[0] + c[0], 0};
+
+    return BLI_tridiagonal_solve(a2, b, c2, d, r_x, count);
+  }
+
+  /* If not really cyclic, fall back to the simple solver. */
   float a0 = a[0], cN = c[count - 1];
 
-  /* if not really cyclic, fall back to the simple solver */
   if (a0 == 0.0f && cN == 0.0f) {
     return BLI_tridiagonal_solve(a, b, c, d, r_x, count);
   }
diff --git a/source/blender/blenlib/tests/BLI_math_solvers_test.cc b/source/blender/blenlib/tests/BLI_math_solvers_test.cc
new file mode 100644
index 00000000000..08cb8571490
--- /dev/null
+++ b/source/blender/blenlib/tests/BLI_math_solvers_test.cc
@@ -0,0 +1,70 @@
+/* Apache License, Version 2.0 */
+
+#include "testing/testing.h"
+
+#include "BLI_math_solvers.h"
+
+TEST(math_solvers, Tridiagonal1)
+{
+  float a[1] = {1};  // ignored
+  float b[1] = {2};
+  float c[1] = {1};  // ignored
+  float d[1] = {4};
+  float x[1];
+
+  EXPECT_TRUE(BLI_tridiagonal_solve(a, b, c, d, x, 1));
+  EXPECT_FLOAT_EQ(x[0], 2);
+}
+
+TEST(math_solvers, Tridiagonal3)
+{
+  float a[3] = {1, 2, 3};  // 1 ignored
+  float b[3] = {4, 5, 6};
+  float c[3] = {7, 8, 9};  // 9 ignored
+  float d[3] = {18, 36, 24};
+  float x[3];
+
+  EXPECT_TRUE(BLI_tridiagonal_solve(a, b, c, d, x, 3));
+  EXPECT_FLOAT_EQ(x[0], 1);
+  EXPECT_FLOAT_EQ(x[1], 2);
+  EXPECT_FLOAT_EQ(x[2], 3);
+}
+
+TEST(math_solvers, CyclicTridiagonal1)
+{
+  float a[1] = {1};
+  float b[1] = {2};
+  float c[1] = {1};
+  float d[1] = {4};
+  float x[1];
+
+  EXPECT_TRUE(BLI_tridiagonal_solve_cyclic(a, b, c, d, x, 1));
+  EXPECT_FLOAT_EQ(x[0], 1);
+}
+
+TEST(math_solvers, CyclicTridiagonal2)
+{
+  float a[2] = {1, 2};
+  float b[2] = {3, 4};
+  float c[2] = {5, 6};
+  float d[2] = {15, 16};
+  float x[2];
+
+  EXPECT_TRUE(BLI_tridiagonal_solve_cyclic(a, b, c, d, x, 2));
+  EXPECT_FLOAT_EQ(x[0], 1);
+  EXPECT_FLOAT_EQ(x[1], 2);
+}
+
+TEST(math_solvers, CyclicTridiagonal3)
+{
+  float a[3] = {1, 2, 3};
+  float b[3] = {4, 5, 6};
+  float c[3] = {7, 8, 9};
+  float d[3] = {21, 36, 33};
+  float x[3];
+
+  EXPECT_TRUE(BLI_tridiagonal_solve_cyclic(a, b, c, d, x, 3));
+  EXPECT_FLOAT_EQ(x[0], 1);
+  EXPECT_FLOAT_EQ(x[1], 2);
+  EXPECT_FLOAT_EQ(x[2], 3);
+}