Fix T87779: Asymmetric vertex positions in circles primitives

Add sin_cos_from_fraction which ensures each quadrant has matching
values when their sign is flipped.

3 files changed, 83 insertions, 14 deletions

diff --git a/source/blender/blenlib/BLI_math_rotation.h b/source/blender/blenlib/BLI_math_rotation.h
index 192ad482a69..fef51fa780e 100644
--- a/source/blender/blenlib/BLI_math_rotation.h
+++ b/source/blender/blenlib/BLI_math_rotation.h
@@ -176,6 +176,26 @@ void mat3_to_quat_is_ok(float q[4], const float mat[3][3]);
 
 /* Other. */
 
+/**
+ * Utility function that performs `sinf` & `cosf` where the quadrants of the circle
+ * will have exactly matching values when their sign is flipped.
+ * This works as long as the denominator can be divided by 2 or 4,
+ * otherwise `sinf` & `cosf` are used without any additional logic.
+ *
+ * Besides adjustments to precision, this function is the equivalent of:
+ * \code {.c}
+ * float phi = (2 * M_PI) * (float)i / (float)denominator;
+ * *r_sin = sinf(phi);
+ * *r_cos = cosf(phi);
+ * \endcode
+ *
+ * \param numerator: An integer factor in [0..denominator] (inclusive).
+ * \param denominator: The faction denominator (typically the number of segments of the circle).
+ * \param r_sin: The resulting sine.
+ * \param r_cos: The resulting cosine.
+ */
+void sin_cos_from_fraction(const int numerator, const int denominator, float *r_sin, float *r_cos);
+
 void print_qt(const char *str, const float q[4]);
 
 #define print_qt_id(q) print_qt(STRINGIFY(q), q)
diff --git a/source/blender/blenlib/intern/math_rotation.c b/source/blender/blenlib/intern/math_rotation.c
index 92223bdf1d5..f0bfc7c21e1 100644
--- a/source/blender/blenlib/intern/math_rotation.c
+++ b/source/blender/blenlib/intern/math_rotation.c
@@ -915,6 +915,55 @@ float tri_to_quat(float q[4], const float a[3], const float b[3], const float c[
   return len;
 }
 
+void sin_cos_from_fraction(const int numerator, const int denominator, float *r_sin, float *r_cos)
+{
+  BLI_assert((numerator <= denominator) && (denominator > 0));
+  if ((denominator & 3) == 0) {
+    const int denominator_4 = denominator / 4;
+    if (numerator <= denominator_4) {
+      /* Fall through. */
+    }
+    else {
+      if (numerator <= denominator_4 * 2) {
+        const float phi = (float)(2.0 * M_PI) *
+                          ((float)(numerator - denominator_4) / (float)denominator);
+        *r_sin = cosf(phi);
+        *r_cos = -sinf(phi);
+      }
+      else if (numerator <= denominator_4 * 3) {
+        const float phi = (float)(2.0 * M_PI) *
+                          ((float)(numerator - (denominator_4 * 2)) / (float)denominator);
+        *r_sin = -sinf(phi);
+        *r_cos = -cosf(phi);
+      }
+      else {
+        const float phi = (float)(2.0 * M_PI) *
+                          ((float)(numerator - (denominator_4 * 3)) / (float)denominator);
+        *r_cos = sinf(phi);
+        *r_sin = -cosf(phi);
+      }
+      return;
+    }
+  }
+  else if ((denominator & 1) == 0) {
+    const int denominator_2 = denominator / 2;
+    if (numerator <= denominator_2) {
+      /* Fall through. */
+    }
+    else {
+      const float phi = (float)(2.0 * M_PI) *
+                        ((float)(numerator - denominator_2) / (float)denominator);
+      *r_sin = -sinf(phi);
+      *r_cos = -cosf(phi);
+      return;
+    }
+  }
+
+  const float phi = (float)(2.0 * M_PI) * ((float)numerator / (float)denominator);
+  *r_sin = sinf(phi);
+  *r_cos = cosf(phi);
+}
+
 void print_qt(const char *str, const float q[4])
 {
   printf("%s: %.3f %.3f %.3f %.3f\n", str, q[0], q[1], q[2], q[3]);
diff --git a/source/blender/bmesh/operators/bmo_primitive.c b/source/blender/bmesh/operators/bmo_primitive.c
index 432c7590f3c..2ed0964d735 100644
--- a/source/blender/bmesh/operators/bmo_primitive.c
+++ b/source/blender/bmesh/operators/bmo_primitive.c
@@ -855,12 +855,12 @@ void bmo_create_uvsphere_exec(BMesh *bm, BMOperator *op)
   /* one segment first */
   for (a = 0; a <= tot; a++) {
     /* Going in this direction, then edge extruding, makes normals face outward */
-    /* Calculate with doubles for higher precision, see: T87779. */
-    const float phi = M_PI * ((double)a / (double)tot);
+    float sin_phi, cos_phi;
+    sin_cos_from_fraction(a, tot, &sin_phi, &cos_phi);
 
     vec[0] = 0.0;
-    vec[1] = rad * sinf(phi);
-    vec[2] = rad * cosf(phi);
+    vec[1] = rad * sin_phi;
+    vec[2] = rad * cos_phi;
     eve = BM_vert_create(bm, vec, NULL, BM_CREATE_NOP);
     BMO_vert_flag_enable(bm, eve, VERT_MARK);
 
@@ -1262,11 +1262,9 @@ void bmo_create_circle_exec(BMesh *bm, BMOperator *op)
 
   for (a = 0; a < segs; a++) {
     /* Going this way ends up with normal(s) upward */
-
-    /* Calculate with doubles for higher precision, see: T87779. */
-    const float phi = (2.0 * M_PI) * ((double)a / (double)segs);
-    vec[0] = -radius * sinf(phi);
-    vec[1] = radius * cosf(phi);
+    sin_cos_from_fraction(a, segs, &vec[0], &vec[1]);
+    vec[0] *= -radius;
+    vec[1] *= radius;
     vec[2] = 0.0f;
     mul_m4_v3(mat, vec);
     v1 = BM_vert_create(bm, vec, NULL, BM_CREATE_NOP);
@@ -1394,15 +1392,17 @@ void bmo_create_cone_exec(BMesh *bm, BMOperator *op)
 
   for (int i = 0; i < segs; i++) {
     /* Calculate with doubles for higher precision, see: T87779. */
-    const float phi = (2.0 * M_PI) * ((double)i / (double)segs);
-    vec[0] = rad1 * sinf(phi);
-    vec[1] = rad1 * cosf(phi);
+    float sin_phi, cos_phi;
+    sin_cos_from_fraction(i, segs, &sin_phi, &cos_phi);
+
+    vec[0] = rad1 * sin_phi;
+    vec[1] = rad1 * cos_phi;
     vec[2] = -depth_half;
     mul_m4_v3(mat, vec);
     v1 = BM_vert_create(bm, vec, NULL, BM_CREATE_NOP);
 
-    vec[0] = rad2 * sinf(phi);
-    vec[1] = rad2 * cosf(phi);
+    vec[0] = rad2 * sin_phi;
+    vec[1] = rad2 * cos_phi;
     vec[2] = depth_half;
     mul_m4_v3(mat, vec);
     v2 = BM_vert_create(bm, vec, NULL, BM_CREATE_NOP);