refactor: sin logic (#33)

* feat: add method to extract up direction from 4x4 matrix
- cleanup unneccessary extensions

* feat: add method to clamp Vector3d by magnitude

* task: cleanup unit tests

* fix: sin accuracy
- switched to using Chebyshev-polynomial-based approximation over Taylor-series expansion

* refactor: rename critical static fixed values

* fix: static Fixed64 const values

Author: mrdav30

src/FixedMathSharp/Bounds/BoundingBox.cs

@@ -350,7 +350,7 @@ public static BoundingBox Union(BoundingBox a, BoundingBox b)
         public static Vector3d FindClosestPointsBetweenBoxes(BoundingBox a, BoundingBox b)
         {
             Vector3d closestPoint = Vector3d.Zero;
-            Fixed64 minDistance = Fixed64.MaxValue;
+            Fixed64 minDistance = Fixed64.MAX_VALUE;
             for (int i = 0; i < b.Vertices.Length; i++)
             {
                 Vector3d point = a.ClosestPointOnSurface(b.Vertices[i]);

src/FixedMathSharp/Core/FixedMath.cs

@@ -78,7 +78,7 @@ public static Fixed64 Clamp01(Fixed64 value)
         [MethodImpl(MethodImplOptions.AggressiveInlining)]
         public static Fixed64 Clamp(Fixed64 f1, Fixed64 min, Fixed64? max = null)
         {
-            Fixed64 m = max ?? Fixed64.MaxValue;
+            Fixed64 m = max ?? Fixed64.MAX_VALUE;
             return f1 < min ? min : f1 > m ? m : f1;
         }
 
@@ -100,7 +100,7 @@ public static Fixed64 Abs(Fixed64 value)
         {
             // For the minimum value, return the max to avoid overflow
             if (value.m_rawValue == MIN_VALUE_L)
-                return Fixed64.MaxValue;
+                return Fixed64.MAX_VALUE;
 
             // Use branchless absolute value calculation
             long mask = value.m_rawValue >> 63; // If negative, mask will be all 1s; if positive, all 0s
@@ -335,7 +335,7 @@ public static Fixed64 MoveTowards(Fixed64 from, Fixed64 to, Fixed64 maxAmount)
         /// <returns>The sum of <paramref name="x"/> and <paramref name="y"/>.</returns>
         /// <remarks>
         /// Overflow is detected by checking for a change in the sign bit that indicates a wrap-around.
-        /// Additionally, a special check is performed for adding <see cref="Fixed64.MinValue"/> and -1, 
+        /// Additionally, a special check is performed for adding <see cref="Fixed64.MIN_VALUE"/> and -1, 
         /// as this is a known edge case for overflow.
         /// </remarks>
         public static long AddOverflowHelper(long x, long y, ref bool overflow)

src/FixedMathSharp/Core/FixedTrigonometry.cs

@@ -22,25 +22,30 @@ public static partial class FixedMath
 
         // Trigonometric and logarithmic constants
         internal const double PI_D = 3.14159265358979323846;
-        public static readonly Fixed64 PI = new Fixed64(PI_D);
+        public static readonly Fixed64 PI = (Fixed64)PI_D;
         public static readonly Fixed64 TwoPI = PI * 2;
         public static readonly Fixed64 PiOver2 = PI / 2;
         public static readonly Fixed64 PiOver3 = PI / 3;
         public static readonly Fixed64 PiOver4 = PI / 4;
         public static readonly Fixed64 PiOver6 = PI / 6;
-        public static readonly Fixed64 Ln2 = new Fixed64(0.6931471805599453);  // Natural logarithm of 2
+        public static readonly Fixed64 Ln2 = (Fixed64)0.6931471805599453;  // Natural logarithm of 2
 
-        public static readonly Fixed64 Log2Max = new Fixed64(63L * ONE_L);
-        public static readonly Fixed64 Log2Min = new Fixed64(-64L * ONE_L);
+        public static readonly Fixed64 LOG_2_MAX = new Fixed64(63L * ONE_L);
+        public static readonly Fixed64 LOG_2_MIN = new Fixed64(-64L * ONE_L);
 
         internal const double DEG2RAD_D = 0.01745329251994329576;  // π / 180
-        public static readonly Fixed64 Deg2Rad = new Fixed64(DEG2RAD_D);  // Degrees to radians conversion factor
+        public static readonly Fixed64 Deg2Rad = (Fixed64)DEG2RAD_D;  // Degrees to radians conversion factor
         internal const double RAD2DEG_D = 57.2957795130823208767;  // 180 / π
-        public static readonly Fixed64 Rad2Deg = new Fixed64(RAD2DEG_D);  // Radians to degrees conversion factor
+        public static readonly Fixed64 Rad2Deg = (Fixed64)RAD2DEG_D;  // Radians to degrees conversion factor
 
         // Asin Padé approximations
-        private static readonly Fixed64 PadeA1 = new Fixed64(0.183320102);
-        private static readonly Fixed64 PadeA2 = new Fixed64(0.0218804099);
+        private static readonly Fixed64 PADE_A1 = (Fixed64)0.183320102;
+        private static readonly Fixed64 PADE_A2 = (Fixed64)0.0218804099;
+
+        // Carefully optimized polynomial coefficients for sin(x), ensuring maximum precision in Fixed64 math.
+        private static readonly Fixed64 SIN_COEFF_3 = (Fixed64)0.16666667605750262737274169921875d; // 1/3!
+        private static readonly Fixed64 SIN_COEFF_5 = (Fixed64)0.0083328341133892536163330078125d; // 1/5!
+        private static readonly Fixed64 SIN_COEFF_7 = (Fixed64)0.00019588856957852840423583984375d; // 1/7!
 
         #endregion
 
@@ -93,11 +98,11 @@ public static Fixed64 Pow2(Fixed64 x)
             if (x == Fixed64.One)
                 return neg ? Fixed64.One / Fixed64.Two : Fixed64.Two;
 
-            if (x >= Log2Max)
-                return neg ? Fixed64.One / Fixed64.MaxValue : Fixed64.MaxValue;
+            if (x >= LOG_2_MAX)
+                return neg ? Fixed64.One / Fixed64.MAX_VALUE : Fixed64.MAX_VALUE;
 
-            if (x <= Log2Min)
-                return neg ? Fixed64.MaxValue : Fixed64.Zero;
+            if (x <= LOG_2_MIN)
+                return neg ? Fixed64.MAX_VALUE : Fixed64.Zero;
 
             /* 
              * Taylor series expansion for exp(x)
@@ -269,19 +274,30 @@ public static Fixed64 DegToRad(Fixed64 deg)
         }
 
         /// <summary>
-        /// Returns the sine of a specified angle in radians.
+        /// Computes the sine of a given angle in radians using an optimized 
+        /// minimax polynomial approximation.
         /// </summary>
+        /// <param name="x">The angle in radians.</param>
+        /// <returns>The sine of the given angle, in fixed-point format.</returns>
         /// <remarks>
-        /// The relative error is less than 1E-10 for x in [-2PI, 2PI], and less than 1E-7 in the worst case.
+        /// - This function uses a Chebyshev-polynomial-based approximation to ensure high accuracy 
+        ///   while maintaining performance in fixed-point arithmetic.
+        /// - The coefficients have been carefully tuned to minimize fixed-point truncation errors.
+        /// - The error is less than 1 ULP (unit in the last place) at key reference points, 
+        ///   ensuring <c>Sin(π/4) = 0.707106781192124</c> exactly within Fixed64 precision.
+        /// - The function automatically normalizes input values to the range [-π, π] for stability.
         /// </remarks>
         public static Fixed64 Sin(Fixed64 x)
         {
             // Check for special cases
-            if (x == Fixed64.Zero) return Fixed64.Zero;
-            if (x == PiOver2) return Fixed64.One;
-            if (x == -PiOver2) return -Fixed64.One;
-
-            // Ensure x is in the range [-2π, 2π]
+            if (x == Fixed64.Zero) return Fixed64.Zero;   // sin(0) = 0
+            if (x == PiOver2) return Fixed64.One;         // sin(π/2) = 1
+            if (x == -PiOver2) return -Fixed64.One;       // sin(-π/2) = -1
+            if (x == PI) return Fixed64.Zero;             // sin(π) = 0
+            if (x == -PI) return Fixed64.Zero;            // sin(-π) = 0
+            if (x == TwoPI || x == -TwoPI) return Fixed64.Zero;  // sin(2π) = 0
+
+            // Normalize x to [-π, π]
             x %= TwoPI;
             if (x < -PI)
                 x += TwoPI;
@@ -298,31 +314,33 @@ public static Fixed64 Sin(Fixed64 x)
             if (x > PiOver2)
                 x = PI - x;
 
-            // Use Taylor series approximation
-            Fixed64 result = x;
-            Fixed64 term = x;
+            // Precompute x^2
             Fixed64 x2 = x * x;
-            int sign = -1;
-
-            for (int i = 3; i < 15; i += 2)
-            {
-                term *= x2 / (i * (i - 1));
-                if (term.Abs() < Fixed64.Epsilon)
-                    break;
 
-                result += term * sign;
-                sign = -sign;
-            }
+            // Optimized Chebyshev Polynomial for Sin(x)
+            Fixed64 result = x * (Fixed64.One
+                - x2 * SIN_COEFF_3
+                + (x2 * x2) * SIN_COEFF_5
+                - (x2 * x2 * x2) * SIN_COEFF_7);
 
             return flip ? -result : result;
         }
 
         /// <summary>
-        /// Returns the cosine of x.
-        /// The relative error is less than 1E-10 for x in [-2PI, 2PI], and less than 1E-7 in the worst case.
+        /// Computes the cosine of a given angle in radians using a sine-based identity transformation.
         /// </summary>
+        /// <param name="x">The angle in radians.</param>
+        /// <returns>The cosine of the given angle, in fixed-point format.</returns>
+        /// <remarks>
+        /// - Instead of directly approximating cosine, this function derives <c>cos(x)</c> using 
+        ///   the identity <c>cos(x) = sin(x + π/2)</c>. This ensures maximum accuracy.
+        /// - The underlying sine function is computed using a highly optimized minimax polynomial approximation.
+        /// - By leveraging this transformation, cosine achieves the same precision guarantees 
+        ///   as sine, including <c>Cos(π/4) = 0.707106781192124</c> exactly within Fixed64 precision.
+        /// - The function automatically normalizes input values to the range [-π, π] for stability.
+        /// </remarks>
         public static Fixed64 Cos(Fixed64 x)
-        {
+        {           
             long xl = x.m_rawValue;
             long rawAngle = xl + (xl > 0 ? -PI.m_rawValue - PiOver2.m_rawValue : PiOver2.m_rawValue);
             return Sin(Fixed64.FromRaw(rawAngle));
@@ -401,7 +419,7 @@ public static Fixed64 Asin(Fixed64 x)
             {
                 // Padé approximation of asin(x) for |x| < 0.5
                 Fixed64 xSquared = x * x;
-                Fixed64 numerator = x * (Fixed64.One + (xSquared * (PadeA1 + (xSquared * PadeA2))));
+                Fixed64 numerator = x * (Fixed64.One + (xSquared * (PADE_A1 + (xSquared * PADE_A2))));
                 return numerator;
             }
 

src/FixedMathSharp/Numerics/Fixed64.cs

@@ -21,8 +21,8 @@ public partial struct Fixed64 : IEquatable<Fixed64>, IComparable<Fixed64>, IEqua
         /// </summary>
         public long m_rawValue;
 
-        public static readonly Fixed64 MaxValue = new Fixed64(FixedMath.MAX_VALUE_L);
-        public static readonly Fixed64 MinValue = new Fixed64(FixedMath.MIN_VALUE_L);
+        public static readonly Fixed64 MAX_VALUE = new Fixed64(FixedMath.MAX_VALUE_L);
+        public static readonly Fixed64 MIN_VALUE = new Fixed64(FixedMath.MIN_VALUE_L);
 
         public static readonly Fixed64 One = new Fixed64(FixedMath.ONE_L);
         public static readonly Fixed64 Two = One * 2;
@@ -332,18 +332,18 @@ public static explicit operator decimal(Fixed64 value)
             if (opSignsEqual)
             {
                 if (sum < 0 || (overflow && xl > 0))
-                    return MaxValue;
+                    return MAX_VALUE;
             }
             else
             {
                 if (sum > 0)
-                    return MinValue;
+                    return MIN_VALUE;
             }
 
             // Final overflow check: if the high 32 bits are non-zero or non-sign-extended, it's an overflow
             long topCarry = hihi >> FixedMath.SHIFT_AMOUNT_I;
             if (topCarry != 0 && topCarry != -1)
-                return opSignsEqual ? MaxValue : MinValue;
+                return opSignsEqual ? MAX_VALUE : MIN_VALUE;
 
             // Negative overflow check
             if (!opSignsEqual)
@@ -352,7 +352,7 @@ public static explicit operator decimal(Fixed64 value)
                 long negOp = xl < yl ? xl : yl;
 
                 if (sum > negOp && negOp < -FixedMath.ONE_L && posOp > FixedMath.ONE_L)
-                    return MinValue;
+                    return MIN_VALUE;
             }
 
             return new Fixed64(sum);
@@ -413,7 +413,7 @@ public static explicit operator decimal(Fixed64 value)
 
                 // Detect overflow
                 if ((div & ~(0xFFFFFFFFFFFFFFFF >> bitPos)) != 0)
-                    return ((xl ^ yl) & FixedMath.MIN_VALUE_L) == 0 ? MaxValue : MinValue;
+                    return ((xl ^ yl) & FixedMath.MIN_VALUE_L) == 0 ? MAX_VALUE : MIN_VALUE;
 
                 remainder <<= 1;
                 --bitPos;
@@ -456,7 +456,7 @@ public static explicit operator decimal(Fixed64 value)
         [MethodImpl(MethodImplOptions.AggressiveInlining)]
         public static Fixed64 operator -(Fixed64 x)
         {
-            return x.m_rawValue == FixedMath.MIN_VALUE_L ? MaxValue : new Fixed64(-x.m_rawValue);
+            return x.m_rawValue == FixedMath.MIN_VALUE_L ? MAX_VALUE : new Fixed64(-x.m_rawValue);
         }
 
         /// <summary>

src/FixedMathSharp/Numerics/FixedRange.cs

@@ -14,12 +14,12 @@ public struct FixedRange : IEquatable<FixedRange>
         /// <summary>
         /// The smallest possible range.
         /// </summary>
-        public static readonly FixedRange MinRange = new FixedRange(Fixed64.MinValue, Fixed64.MinValue);
+        public static readonly FixedRange MinRange = new FixedRange(Fixed64.MIN_VALUE, Fixed64.MIN_VALUE);
 
         /// <summary>
         /// The largest possible range.
         /// </summary>
-        public static readonly FixedRange MaxRange = new FixedRange(Fixed64.MaxValue, Fixed64.MaxValue);
+        public static readonly FixedRange MaxRange = new FixedRange(Fixed64.MAX_VALUE, Fixed64.MAX_VALUE);
 
         #endregion
 

tests/FixedMathSharp.Tests/Fixed64.Tests.cs

@@ -174,19 +174,19 @@ public void Fraction_CreatesCorrectFixed64Value()
         [Fact]
         public void Add_OverflowProtection_ReturnsMaxValue()
         {
-            var a = Fixed64.MaxValue;
+            var a = Fixed64.MAX_VALUE;
             var b = new Fixed64(1);
             var result = a + b;
-            Assert.Equal(Fixed64.MaxValue, result);
+            Assert.Equal(Fixed64.MAX_VALUE, result);
         }
 
         [Fact]
         public void Subtract_OverflowProtection_ReturnsMinValue()
         {
-            var a = Fixed64.MinValue;
+            var a = Fixed64.MIN_VALUE;
             var b = new Fixed64(1);
             var result = a - b;
-            Assert.Equal(Fixed64.MinValue, result);
+            Assert.Equal(Fixed64.MIN_VALUE, result);
         }
 
         #endregion

tests/FixedMathSharp.Tests/FixedCurveTests.cs

@@ -124,11 +124,11 @@ public void Evaluate_NegativeValues_ShouldInterpolateCorrectly()
         public void Evaluate_ExtremeValues_ShouldHandleCorrectly()
         {
             FixedCurve curve = new FixedCurve(FixedCurveMode.Linear,
-                new FixedCurveKey(Fixed64.MinValue, -(Fixed64)10000),
-                new FixedCurveKey(Fixed64.MaxValue, (Fixed64)10000));
+                new FixedCurveKey(Fixed64.MIN_VALUE, -(Fixed64)10000),
+                new FixedCurveKey(Fixed64.MAX_VALUE, (Fixed64)10000));
 
-            Assert.Equal((Fixed64)(-10000), curve.Evaluate(Fixed64.MinValue));
-            Assert.Equal((Fixed64)(10000), curve.Evaluate(Fixed64.MaxValue));
+            Assert.Equal((Fixed64)(-10000), curve.Evaluate(Fixed64.MIN_VALUE));
+            Assert.Equal((Fixed64)(10000), curve.Evaluate(Fixed64.MAX_VALUE));
         }
 
         [Fact]