| | | 1 | | using System; |
| | | 2 | | using System.Runtime.CompilerServices; |
| | | 3 | | |
| | | 4 | | namespace FixedMathSharp |
| | | 5 | | { |
| | | 6 | | public static partial class FixedMath |
| | | 7 | | { |
| | | 8 | | #region Fields and Constants |
| | | 9 | | |
| | | 10 | | /// <summary> |
| | | 11 | | /// Provides a lookup table of integer powers of 10 from 10^0 to 10^9. |
| | | 12 | | /// </summary> |
| | | 13 | | /// <remarks> |
| | | 14 | | /// This array can be used to efficiently retrieve the value of 10 raised to an integer |
| | | 15 | | /// exponent within the supported range, avoiding repeated calculations. |
| | | 16 | | /// The index corresponds to the exponent. |
| | | 17 | | /// </remarks> |
| | 1 | 18 | | public static readonly int[] Pow10Lookup = { |
| | 1 | 19 | | 1, // 10^0 = 1 |
| | 1 | 20 | | 10, // 10^1 = 10 |
| | 1 | 21 | | 100, // 10^2 = 100 |
| | 1 | 22 | | 1000, // 10^3 = 1000 |
| | 1 | 23 | | 10000, // 10^4 = 10000 |
| | 1 | 24 | | 100000, // 10^5 = 100000 |
| | 1 | 25 | | 1000000, // 10^6 = 1000000 |
| | 1 | 26 | | 10000000, // 10^7 = 1000000 |
| | 1 | 27 | | 100000000, // 10^8 = 1000000 |
| | 1 | 28 | | 1000000000, // 10^9 = 1000000 |
| | 1 | 29 | | }; |
| | | 30 | | |
| | | 31 | | // Trigonometric and logarithmic constants |
| | | 32 | | |
| | | 33 | | /// <summary> |
| | | 34 | | /// Represents the mathematical constant π (pi). |
| | | 35 | | /// </summary> |
| | | 36 | | /// <remarks>The value is approximately 3.14159265358979323846.</remarks> |
| | 1 | 37 | | public static readonly Fixed64 PI = (Fixed64)3.14159265358979323846; |
| | | 38 | | |
| | | 39 | | /// <summary> |
| | | 40 | | /// Represents the mathematical constant 2π as a fixed-point value. |
| | | 41 | | /// </summary> |
| | | 42 | | /// <remarks>This value is commonly used to represent a full rotation in radians.</remarks> |
| | 1 | 43 | | public static readonly Fixed64 TwoPI = PI * Fixed64.Two; |
| | | 44 | | /// <summary> |
| | | 45 | | /// Represents the mathematical constant π divided by 2 as a fixed-point value. |
| | | 46 | | /// </summary> |
| | | 47 | | /// <remarks>This value is used where π/2 is required.</remarks> |
| | 1 | 48 | | public static readonly Fixed64 PiOver2 = PI / new Fixed64(2); |
| | | 49 | | /// <summary> |
| | | 50 | | /// Represents the mathematical constant π divided by 3 as a fixed-point value. |
| | | 51 | | /// </summary> |
| | 1 | 52 | | public static readonly Fixed64 PiOver3 = PI / new Fixed64(3); |
| | | 53 | | /// <summary> |
| | | 54 | | /// Represents the mathematical constant π divided by 4 as a fixed-point value. |
| | | 55 | | /// </summary> |
| | 1 | 56 | | public static readonly Fixed64 PiOver4 = PI / new Fixed64(4); |
| | | 57 | | /// <summary> |
| | | 58 | | /// Represents the mathematical constant π divided by 6 as a fixed-point value. |
| | | 59 | | /// </summary> |
| | 1 | 60 | | public static readonly Fixed64 PiOver6 = PI / new Fixed64(6); |
| | | 61 | | /// <summary> |
| | | 62 | | /// Represents the multiplicative inverse of the mathematical constant π (pi) as a fixed-point value. |
| | | 63 | | /// </summary> |
| | 1 | 64 | | public static readonly Fixed64 InvertedPI = Fixed64.One / PI; |
| | | 65 | | |
| | | 66 | | /// <summary> |
| | | 67 | | /// Represents the fixed-point value for 180. |
| | | 68 | | /// </summary> |
| | 1 | 69 | | public static readonly Fixed64 OneEighty = new(180); |
| | | 70 | | |
| | | 71 | | /// <summary> |
| | | 72 | | /// Degrees to radians conversion factor (π / 180) |
| | | 73 | | /// </summary> |
| | 1 | 74 | | public static readonly Fixed64 Deg2Rad = PI / OneEighty; |
| | | 75 | | /// <summary> |
| | | 76 | | /// Radians to degrees conversion factor (180 / π) |
| | | 77 | | /// </summary> |
| | 1 | 78 | | public static readonly Fixed64 Rad2Deg = OneEighty / PI; |
| | | 79 | | |
| | | 80 | | /// <summary> |
| | | 81 | | /// Represents the natural logarithm of 2 as a fixed-point value. |
| | | 82 | | /// </summary> |
| | 1 | 83 | | public static readonly Fixed64 Ln2 = new(0.6931471805599453); |
| | | 84 | | /// <summary> |
| | | 85 | | /// Represents the maximum value for the base-2 logarithm that can be represented by a Fixed64 instance. |
| | | 86 | | /// </summary> |
| | | 87 | | /// <remarks> |
| | | 88 | | /// This constant is useful when performing logarithmic calculations to ensure results |
| | | 89 | | /// do not exceed the representable range of the Fixed64 type. |
| | | 90 | | /// </remarks> |
| | 1 | 91 | | public static readonly Fixed64 LOG_2_MAX = new(63L * ONE_L); |
| | | 92 | | /// <summary> |
| | | 93 | | /// Represents the minimum base-2 logarithm value supported by the Fixed64 type. |
| | | 94 | | /// </summary> |
| | | 95 | | /// <remarks> |
| | | 96 | | /// This constant can be used as a lower bound when performing logarithmic calculations |
| | | 97 | | /// with Fixed64 values to prevent underflow or invalid results. |
| | | 98 | | /// </remarks> |
| | 1 | 99 | | public static readonly Fixed64 LOG_2_MIN = new(-64L * ONE_L); |
| | | 100 | | |
| | | 101 | | // Asin Padé approximations |
| | 1 | 102 | | private static readonly Fixed64 PADE_A1 = (Fixed64)0.183320102; |
| | 1 | 103 | | private static readonly Fixed64 PADE_A2 = (Fixed64)0.0218804099; |
| | | 104 | | |
| | | 105 | | // Carefully optimized polynomial coefficients for sin(x), ensuring maximum precision in Fixed64 math. |
| | 1 | 106 | | private static readonly Fixed64 SIN_COEFF_3 = (Fixed64)0.16666667605750262737274169921875d; // 1/3! |
| | 1 | 107 | | private static readonly Fixed64 SIN_COEFF_5 = (Fixed64)0.0083328341133892536163330078125d; // 1/5! |
| | 1 | 108 | | private static readonly Fixed64 SIN_COEFF_7 = (Fixed64)0.00019588856957852840423583984375d; // 1/7! |
| | | 109 | | |
| | | 110 | | #endregion |
| | | 111 | | |
| | | 112 | | #region FixedTrigonometry Operations |
| | | 113 | | |
| | | 114 | | /// <summary> |
| | | 115 | | /// Raises the base number b to the power of exp. |
| | | 116 | | /// Uses logarithms to compute power efficiently for fixed-point values. |
| | | 117 | | /// </summary> |
| | | 118 | | /// <exception cref="DivideByZeroException"> |
| | | 119 | | /// The base was Fixed64.Zero, with a negative expFixed64.Onent |
| | | 120 | | /// </exception> |
| | | 121 | | /// <exception cref="ArgumentOutOfRangeException"> |
| | | 122 | | /// The base was negative, with a non-Fixed64.Zero expFixed64.Onent |
| | | 123 | | /// </exception> |
| | | 124 | | public static Fixed64 Pow(Fixed64 b, Fixed64 exp) |
| | 6 | 125 | | { |
| | 6 | 126 | | if (b == Fixed64.One) |
| | 1 | 127 | | return Fixed64.One; |
| | | 128 | | |
| | 5 | 129 | | if (exp.m_rawValue == 0) |
| | 1 | 130 | | return Fixed64.One; |
| | | 131 | | |
| | 4 | 132 | | if (b.m_rawValue == 0) |
| | 2 | 133 | | { |
| | 2 | 134 | | if (exp.m_rawValue < 0) |
| | 1 | 135 | | throw new DivideByZeroException("Cannot raise 0 to a negative power."); |
| | | 136 | | |
| | 1 | 137 | | return Fixed64.Zero; |
| | | 138 | | } |
| | | 139 | | |
| | 2 | 140 | | Fixed64 log2 = Log2(b); // Calculate logarithm base 2 |
| | 2 | 141 | | return Pow2(exp * log2); // Raise 2 to the power of log2 result |
| | 5 | 142 | | } |
| | | 143 | | |
| | | 144 | | /// <summary> |
| | | 145 | | /// Raises 2 to the power of x. |
| | | 146 | | /// Provides high accuracy for small values of x. |
| | | 147 | | /// </summary> |
| | | 148 | | public static Fixed64 Pow2(Fixed64 x) |
| | 8 | 149 | | { |
| | 8 | 150 | | if (x.m_rawValue == 0) |
| | 1 | 151 | | return Fixed64.One; |
| | | 152 | | |
| | | 153 | | // Handle negative expFixed64.Onents by using the reciprocal |
| | 7 | 154 | | bool neg = x.m_rawValue < 0; |
| | 7 | 155 | | if (neg) |
| | 4 | 156 | | x = -x; |
| | | 157 | | |
| | 7 | 158 | | if (x == Fixed64.One) |
| | 3 | 159 | | return neg ? Fixed64.One / Fixed64.Two : Fixed64.Two; |
| | | 160 | | |
| | 4 | 161 | | if (x >= LOG_2_MAX) |
| | 2 | 162 | | return neg ? Fixed64.One / new Fixed64(MAX_VALUE_L) : new Fixed64(MAX_VALUE_L); |
| | | 163 | | |
| | | 164 | | /* |
| | | 165 | | * Taylor series expansion for exp(x) |
| | | 166 | | * From term n, we get term n+1 by multiplying with x/n. |
| | | 167 | | * When the sum term drops to Fixed64.Zero, we can stop summing. |
| | | 168 | | */ |
| | 2 | 169 | | int integerPart = (int)x.Floor(); |
| | 2 | 170 | | x = Fixed64.FromRaw(x.m_rawValue & MAX_SHIFTED_AMOUNT_UI); // Fractional part |
| | | 171 | | |
| | 2 | 172 | | var result = Fixed64.One; |
| | 2 | 173 | | var term = Fixed64.One; |
| | 2 | 174 | | int i = 1; |
| | 13 | 175 | | while (term.m_rawValue != 0) |
| | 11 | 176 | | { |
| | 11 | 177 | | term = FastMul(FastMul(x, term), Ln2) / (Fixed64)i; |
| | 11 | 178 | | result += term; |
| | 11 | 179 | | i++; |
| | 11 | 180 | | } |
| | | 181 | | |
| | 2 | 182 | | result = Fixed64.FromRaw(result.m_rawValue << integerPart); |
| | 2 | 183 | | if (neg) |
| | 1 | 184 | | result = Fixed64.One / result; |
| | | 185 | | |
| | 2 | 186 | | return result; |
| | 8 | 187 | | } |
| | | 188 | | |
| | | 189 | | /// <summary> |
| | | 190 | | /// Returns the base-2 logarithm of a specified number. |
| | | 191 | | /// Provides at least 9 decimals of accuracy. |
| | | 192 | | /// </summary> |
| | | 193 | | /// <remarks> |
| | | 194 | | /// This implementation is based on Clay. S. Turner's fast binary logarithm algorithm |
| | | 195 | | /// (C. S. Turner, "A Fast Binary Logarithm Algorithm", IEEE Signal Processing Mag., pp. 124,140, Sep. 2010.) |
| | | 196 | | /// </remarks> |
| | | 197 | | public static Fixed64 Log2(Fixed64 x) |
| | 6 | 198 | | { |
| | 6 | 199 | | if (x.m_rawValue <= 0) |
| | 1 | 200 | | throw new ArgumentOutOfRangeException(nameof(x), "Cannot compute logarithm of non-positive number."); |
| | | 201 | | |
| | 5 | 202 | | long b = 1U << (SHIFT_AMOUNT_I - 1); // Initial value for binary logarithm |
| | 5 | 203 | | long y = 0; // Result accumulator |
| | 5 | 204 | | long rawX = x.m_rawValue; |
| | | 205 | | |
| | | 206 | | // Adjust rawX to the correct range [1, 2) |
| | 6 | 207 | | while (rawX < ONE_L) |
| | 1 | 208 | | { |
| | 1 | 209 | | rawX <<= 1; |
| | 1 | 210 | | y -= ONE_L; |
| | 1 | 211 | | } |
| | | 212 | | |
| | 11 | 213 | | while (rawX >= (ONE_L << 1)) |
| | 6 | 214 | | { |
| | 6 | 215 | | rawX >>= 1; |
| | 6 | 216 | | y += ONE_L; |
| | 6 | 217 | | } |
| | | 218 | | |
| | 5 | 219 | | Fixed64 z = Fixed64.FromRaw(rawX); // Remaining fraction |
| | | 220 | | |
| | 330 | 221 | | for (int i = 0; i < SHIFT_AMOUNT_I; i++) |
| | 160 | 222 | | { |
| | 160 | 223 | | z = FastMul(z, z); |
| | 160 | 224 | | if (z.m_rawValue >= (ONE_L << 1)) |
| | 15 | 225 | | { |
| | 15 | 226 | | z = Fixed64.FromRaw(z.m_rawValue >> 1); |
| | 15 | 227 | | y += b; |
| | 15 | 228 | | } |
| | 160 | 229 | | b >>= 1; |
| | 160 | 230 | | } |
| | | 231 | | |
| | 5 | 232 | | return Fixed64.FromRaw(y); |
| | 5 | 233 | | } |
| | | 234 | | |
| | | 235 | | /// <summary> |
| | | 236 | | /// Returns the natural logarithm of a specified fixed-point number. |
| | | 237 | | /// Provides at least 7 decimals of accuracy. |
| | | 238 | | /// </summary> |
| | | 239 | | public static Fixed64 Ln(Fixed64 x) |
| | 2 | 240 | | { |
| | 2 | 241 | | if (x.m_rawValue <= 0) |
| | 1 | 242 | | throw new ArgumentOutOfRangeException(nameof(x), "Cannot compute logarithm of non-positive number."); |
| | | 243 | | |
| | 1 | 244 | | return FastMul(Log2(x), Ln2).Round(); |
| | 1 | 245 | | } |
| | | 246 | | |
| | | 247 | | /// <summary> |
| | | 248 | | /// Returns the square root of a specified fixed-point number. |
| | | 249 | | /// </summary> |
| | | 250 | | public static Fixed64 Sqrt(Fixed64 x) |
| | 807 | 251 | | { |
| | 807 | 252 | | if (x.m_rawValue < 0) |
| | 1 | 253 | | throw new ArgumentOutOfRangeException(nameof(x), "Cannot compute square root of a negative number."); |
| | | 254 | | |
| | 806 | 255 | | ulong num = (ulong)x.m_rawValue; |
| | 806 | 256 | | ulong result = 0UL; |
| | 806 | 257 | | ulong bit = 1UL << (sizeof(long) * 8) - 2; // second-to-top bit of a 64-bit integer |
| | | 258 | | |
| | | 259 | | // Adjust the bit position to a suitable starting point |
| | 11694 | 260 | | while (bit > num) |
| | 10888 | 261 | | bit >>= 2; |
| | | 262 | | |
| | | 263 | | // Perform the square root calculation using bitwise shifts |
| | 4836 | 264 | | for (int i = 0; i < 2; ++i) |
| | 1612 | 265 | | { |
| | | 266 | | // Calculate the top bits of the square root result |
| | 29412 | 267 | | while (bit != 0) |
| | 27800 | 268 | | { |
| | 27800 | 269 | | if (num >= result + bit) |
| | 10339 | 270 | | { |
| | 10339 | 271 | | num -= result + bit; |
| | 10339 | 272 | | result = (result >> 1) + bit; |
| | 10339 | 273 | | } |
| | | 274 | | else |
| | 17461 | 275 | | { |
| | 17461 | 276 | | result >>= 1; |
| | 17461 | 277 | | } |
| | | 278 | | |
| | 27800 | 279 | | bit >>= 2; |
| | 27800 | 280 | | } |
| | | 281 | | |
| | 1612 | 282 | | if (i == 0) |
| | 806 | 283 | | { |
| | | 284 | | // Process it again to get the remaining bits |
| | 806 | 285 | | if (num > ((1UL << SHIFT_AMOUNT_I) - 1)) |
| | 1 | 286 | | { |
| | | 287 | | // Handle large remainders by adjusting the result |
| | 1 | 288 | | num -= result; |
| | 1 | 289 | | num = (num << SHIFT_AMOUNT_I) - (ulong)Fixed64.Half.m_rawValue; |
| | 1 | 290 | | result = (result << SHIFT_AMOUNT_I) + (ulong)Fixed64.Half.m_rawValue; |
| | 1 | 291 | | } |
| | | 292 | | else |
| | 805 | 293 | | { |
| | 805 | 294 | | num <<= SHIFT_AMOUNT_I; |
| | 805 | 295 | | result <<= SHIFT_AMOUNT_I; |
| | 805 | 296 | | } |
| | | 297 | | |
| | 806 | 298 | | bit = 1UL << (SHIFT_AMOUNT_I - 2); |
| | 806 | 299 | | } |
| | 1612 | 300 | | } |
| | | 301 | | |
| | | 302 | | // Rounding: round up if necessary |
| | 806 | 303 | | if (num > result && (num - result) > (result >> 1)) |
| | 131 | 304 | | ++result; |
| | | 305 | | |
| | 806 | 306 | | return Fixed64.FromRaw((long)result); |
| | 806 | 307 | | } |
| | | 308 | | |
| | | 309 | | /// <summary> |
| | | 310 | | /// Converts a value in radians to degrees. |
| | | 311 | | /// </summary> |
| | | 312 | | [MethodImpl(MethodImplOptions.AggressiveInlining)] |
| | | 313 | | public static Fixed64 RadToDeg(Fixed64 rad) |
| | 20 | 314 | | { |
| | 20 | 315 | | return (rad * OneEighty) / PI; |
| | 20 | 316 | | } |
| | | 317 | | |
| | | 318 | | /// <summary> |
| | | 319 | | /// Converts a value in degrees to radians. |
| | | 320 | | /// </summary> |
| | | 321 | | [MethodImpl(MethodImplOptions.AggressiveInlining)] |
| | | 322 | | public static Fixed64 DegToRad(Fixed64 deg) |
| | 95 | 323 | | { |
| | 95 | 324 | | return (deg * PI) / OneEighty; |
| | 95 | 325 | | } |
| | | 326 | | |
| | | 327 | | /// <summary> |
| | | 328 | | /// Computes the sine of a given angle in radians using an optimized |
| | | 329 | | /// minimax polynomial approximation. |
| | | 330 | | /// </summary> |
| | | 331 | | /// <param name="x">The angle in radians.</param> |
| | | 332 | | /// <returns>The sine of the given angle, in fixed-point format.</returns> |
| | | 333 | | /// <remarks> |
| | | 334 | | /// - This function uses a Chebyshev-polynomial-based approximation to ensure high accuracy |
| | | 335 | | /// while maintaining performance in fixed-point arithmetic. |
| | | 336 | | /// - The coefficients have been carefully tuned to minimize fixed-point truncation errors. |
| | | 337 | | /// - The error is less than 1 ULP (unit in the last place) at key reference points, |
| | | 338 | | /// ensuring <c>Sin(π/4) = 0.707106781192124</c> exactly within Fixed64 precision. |
| | | 339 | | /// - The function automatically normalizes input values to the range [-π, π] for stability. |
| | | 340 | | /// </remarks> |
| | | 341 | | public static Fixed64 Sin(Fixed64 x) |
| | 355 | 342 | | { |
| | | 343 | | // Check for special cases |
| | 413 | 344 | | if (x == Fixed64.Zero) return Fixed64.Zero; // sin(0) = 0 |
| | 379 | 345 | | if (x == PiOver2) return Fixed64.One; // sin(π/2) = 1 |
| | 217 | 346 | | if (x == -PiOver2) return -Fixed64.One; // sin(-π/2) = -1 |
| | 214 | 347 | | if (x == PI) return Fixed64.Zero; // sin(π) = 0 |
| | 236 | 348 | | if (x == -PI) return Fixed64.Zero; // sin(-π) = 0 |
| | 190 | 349 | | if (x == TwoPI || x == -TwoPI) return Fixed64.Zero; // sin(2π) = 0 |
| | | 350 | | |
| | | 351 | | // Normalize x to [-π, π] |
| | 186 | 352 | | x %= TwoPI; |
| | 186 | 353 | | if (x < -PI) |
| | 89 | 354 | | x += TwoPI; |
| | 97 | 355 | | else if (x > PI) |
| | 1 | 356 | | x -= TwoPI; |
| | | 357 | | |
| | 186 | 358 | | bool flip = false; |
| | 186 | 359 | | if (x < Fixed64.Zero) |
| | 4 | 360 | | { |
| | 4 | 361 | | x = -x; |
| | 4 | 362 | | flip = true; |
| | 4 | 363 | | } |
| | | 364 | | |
| | 186 | 365 | | if (x > PiOver2) |
| | 90 | 366 | | x = PI - x; |
| | | 367 | | |
| | | 368 | | // Precompute x^2 |
| | 186 | 369 | | Fixed64 x2 = x * x; |
| | | 370 | | |
| | | 371 | | // Optimized Chebyshev Polynomial for Sin(x) |
| | 186 | 372 | | Fixed64 result = x * (Fixed64.One |
| | 186 | 373 | | - x2 * SIN_COEFF_3 |
| | 186 | 374 | | + (x2 * x2) * SIN_COEFF_5 |
| | 186 | 375 | | - (x2 * x2 * x2) * SIN_COEFF_7); |
| | | 376 | | |
| | 186 | 377 | | return flip ? -result : result; |
| | 355 | 378 | | } |
| | | 379 | | |
| | | 380 | | /// <summary> |
| | | 381 | | /// Computes the cosine of a given angle in radians using a sine-based identity transformation. |
| | | 382 | | /// </summary> |
| | | 383 | | /// <param name="x">The angle in radians.</param> |
| | | 384 | | /// <returns>The cosine of the given angle, in fixed-point format.</returns> |
| | | 385 | | /// <remarks> |
| | | 386 | | /// - Instead of directly approximating cosine, this function derives <c>cos(x)</c> using |
| | | 387 | | /// the identity <c>cos(x) = sin(x + π/2)</c>. This ensures maximum accuracy. |
| | | 388 | | /// - The underlying sine function is computed using a highly optimized minimax polynomial approximation. |
| | | 389 | | /// - By leveraging this transformation, cosine achieves the same precision guarantees |
| | | 390 | | /// as sine, including <c>Cos(π/4) = 0.707106781192124</c> exactly within Fixed64 precision. |
| | | 391 | | /// - The function automatically normalizes input values to the range [-π, π] for stability. |
| | | 392 | | /// </remarks> |
| | | 393 | | public static Fixed64 Cos(Fixed64 x) |
| | 174 | 394 | | { |
| | 174 | 395 | | long xl = x.m_rawValue; |
| | 174 | 396 | | long rawAngle = xl + (xl > 0 ? -PI.m_rawValue - PiOver2.m_rawValue : PiOver2.m_rawValue); |
| | 174 | 397 | | return Sin(Fixed64.FromRaw(rawAngle)); |
| | 174 | 398 | | } |
| | | 399 | | |
| | | 400 | | /// <summary> |
| | | 401 | | /// Calculates the cosine value corresponding to a given sine value, assuming the angle is in the first or |
| | | 402 | | /// second quadrant. |
| | | 403 | | /// </summary> |
| | | 404 | | /// <remarks> |
| | | 405 | | /// This method returns the principal (non-negative) value of the cosine. |
| | | 406 | | /// If the input is outside the valid range for sine values, the result may not be meaningful. |
| | | 407 | | /// </remarks> |
| | | 408 | | /// <param name="sin">The sine of the angle. Must be in the range [-1, 1].</param> |
| | | 409 | | /// <returns>The cosine of the angle, computed as the positive square root of (1 - sin²).</returns> |
| | | 410 | | public static Fixed64 SinToCos(Fixed64 sin) |
| | 1 | 411 | | { |
| | 1 | 412 | | return Sqrt(Fixed64.One - sin * sin); |
| | 1 | 413 | | } |
| | | 414 | | |
| | | 415 | | /// <summary> |
| | | 416 | | /// Returns the tangent of x. |
| | | 417 | | /// </summary> |
| | | 418 | | /// <remarks> |
| | | 419 | | /// This function is not well-tested. It may be wildly inaccurate. |
| | | 420 | | /// </remarks> |
| | | 421 | | public static Fixed64 Tan(Fixed64 x) |
| | 12 | 422 | | { |
| | | 423 | | // Check for special cases |
| | 13 | 424 | | if (x == Fixed64.Zero) return Fixed64.Zero; |
| | 15 | 425 | | if (x == PiOver4) return Fixed64.One; |
| | 8 | 426 | | if (x == -PiOver4) return -Fixed64.One; |
| | | 427 | | |
| | | 428 | | // Normalize x to [-π/2, π/2] |
| | 6 | 429 | | x %= PI; |
| | 6 | 430 | | if (x < -PiOver2) |
| | 1 | 431 | | x += PI; |
| | 5 | 432 | | else if (x > PiOver2) |
| | 1 | 433 | | x -= PI; |
| | | 434 | | |
| | | 435 | | // Use continued fraction to approximate tan(x) |
| | 6 | 436 | | Fixed64 x2 = x * x; |
| | 6 | 437 | | Fixed64 numerator = x; |
| | 6 | 438 | | Fixed64 denominator = Fixed64.One; |
| | | 439 | | |
| | | 440 | | // Iterate over the continued fraction terms |
| | 6 | 441 | | Fixed64 prevDenominator = denominator; |
| | 6 | 442 | | int start = x.Abs() > PiOver6 ? 19 : 13; |
| | 114 | 443 | | for (int i = start; i >= 1; i -= 2) |
| | 51 | 444 | | { |
| | 51 | 445 | | denominator = (Fixed64)i - (x2 / denominator); |
| | 51 | 446 | | if ((denominator - prevDenominator).Abs() < Fixed64.Epsilon) |
| | 0 | 447 | | break; |
| | 51 | 448 | | prevDenominator = denominator; |
| | 51 | 449 | | } |
| | | 450 | | |
| | 6 | 451 | | return numerator / denominator; |
| | 12 | 452 | | } |
| | | 453 | | |
| | | 454 | | /// <summary> |
| | | 455 | | /// Returns the arc-sine of a fixed-point number x, which is the angle in radians |
| | | 456 | | /// whose sine is x, using a combination of a Taylor series expansion and trigonometric identities. |
| | | 457 | | /// |
| | | 458 | | /// For values of x near ±1, the identity asin(x) = π/2 - acos(x) is used for stability. |
| | | 459 | | /// For values of x near 0, a Taylor series expansion is used. |
| | | 460 | | /// </summary> |
| | | 461 | | /// <param name="x">The input value (sine) whose arcsine is to be computed. Should be in the range [-1, 1].</par |
| | | 462 | | /// <returns>The arc-sine of x in radians.</returns> |
| | | 463 | | /// <exception cref="ArithmeticException">Thrown if x is outside the domain [-1, 1].</exception> |
| | | 464 | | public static Fixed64 Asin(Fixed64 x) |
| | 11 | 465 | | { |
| | | 466 | | // Ensure x is within the domain [-1, 1] |
| | 11 | 467 | | if (x < -Fixed64.One || x > Fixed64.One) |
| | 2 | 468 | | throw new ArithmeticException("Input out of domain for Asin: " + x); |
| | | 469 | | |
| | | 470 | | // Handle boundary cases for -1 and 1 |
| | 10 | 471 | | if (x == Fixed64.One) return PiOver2; // asin(1) = π/2 |
| | 9 | 472 | | if (x == -Fixed64.One) return -PiOver2; // asin(-1) = -π/2 |
| | | 473 | | |
| | | 474 | | // Special case handling for asin(0.5) -> π/6 and asin(-0.5) -> -π/6 |
| | 8 | 475 | | if (x == Fixed64.Half) return PiOver6; |
| | 7 | 476 | | if (x == -Fixed64.Half) return -PiOver6; |
| | | 477 | | |
| | | 478 | | // For values close to 0, use a Padé approximation for better precision |
| | 5 | 479 | | if (x.Abs() < Fixed64.Half) |
| | 3 | 480 | | { |
| | | 481 | | // Padé approximation of asin(x) for |x| < 0.5 |
| | 3 | 482 | | Fixed64 xSquared = x * x; |
| | 3 | 483 | | Fixed64 numerator = x * (Fixed64.One + (xSquared * (PADE_A1 + (xSquared * PADE_A2)))); |
| | 3 | 484 | | return numerator; |
| | | 485 | | } |
| | | 486 | | |
| | | 487 | | // For values closer to ±1, use the identity: asin(x) = π/2 - acos(x) for stability |
| | 2 | 488 | | return x > Fixed64.Zero |
| | 2 | 489 | | ? PiOver2 - Acos(x) |
| | 2 | 490 | | : -PiOver2 + Acos(-x); |
| | 9 | 491 | | } |
| | | 492 | | |
| | | 493 | | /// <summary> |
| | | 494 | | /// Returns the arccosine of the specified number x, calculated using a combination of the atan and sqrt functio |
| | | 495 | | /// </summary> |
| | | 496 | | /// <param name="x">The input value whose arccosine is to be computed. Should be in the range [-1, 1].</param> |
| | | 497 | | /// <returns>The arccosine of x in radians.</returns> |
| | | 498 | | /// <exception cref="ArgumentOutOfRangeException">Thrown if x is outside the domain [-1, 1].</exception> |
| | | 499 | | public static Fixed64 Acos(Fixed64 x) |
| | 20 | 500 | | { |
| | 20 | 501 | | if (x < -Fixed64.One || x > Fixed64.One) |
| | 2 | 502 | | throw new ArgumentOutOfRangeException(nameof(x), "Input out of domain for Acos: " + x); |
| | | 503 | | |
| | | 504 | | // For values near 1 or -1, the result is directly known. |
| | 20 | 505 | | if (x == Fixed64.One) return Fixed64.Zero; // acos(1) = 0 |
| | 17 | 506 | | if (x == -Fixed64.One) return PI; // acos(-1) = π |
| | 21 | 507 | | if (x == Fixed64.Zero) return PiOver2; // acos(0) = π/2 |
| | | 508 | | |
| | | 509 | | // Compute using the relationship acos(x) = atan(sqrt(1 - x^2) / x) + π/2 when x is negative |
| | 9 | 510 | | var sqrtTerm = Sqrt(Fixed64.One - x * x); // sqrt(1 - x^2) |
| | 9 | 511 | | var atanTerm = Atan(sqrtTerm / x); |
| | | 512 | | |
| | 9 | 513 | | return x < Fixed64.Zero |
| | 9 | 514 | | ? atanTerm + PI // acos(-x) = atan(...) + π |
| | 9 | 515 | | : atanTerm; // Otherwise, return just atan(sqrt(...)) |
| | 18 | 516 | | } |
| | | 517 | | |
| | | 518 | | /// <summary> |
| | | 519 | | /// Returns the arctangent of the specified number, using a more accurate approximation for larger values. |
| | | 520 | | /// This function has at least 7 decimals of accuracy. |
| | | 521 | | /// </summary> |
| | | 522 | | public static Fixed64 Atan(Fixed64 z) |
| | 78 | 523 | | { |
| | 85 | 524 | | if (z == Fixed64.Zero) return Fixed64.Zero; |
| | 85 | 525 | | if (z == Fixed64.One) return PiOver4; |
| | 63 | 526 | | if (z == -Fixed64.One) return -PiOver4; |
| | | 527 | | |
| | 51 | 528 | | bool neg = z < Fixed64.Zero; |
| | 67 | 529 | | if (neg) z = -z; |
| | | 530 | | |
| | | 531 | | // Adjust series for z > 0.5 using the identity. |
| | | 532 | | Fixed64 adjustedResult; |
| | 51 | 533 | | if (z > Fixed64.Half) |
| | 24 | 534 | | { |
| | | 535 | | // Apply the identity: atan(z) = π/4 - atan((1 - z) / (1 + z)) |
| | 24 | 536 | | Fixed64 transformedZ = (Fixed64.One - z) / (Fixed64.One + z); |
| | 24 | 537 | | adjustedResult = PiOver4 - Atan(transformedZ); |
| | 24 | 538 | | } |
| | | 539 | | else |
| | 27 | 540 | | { |
| | | 541 | | // Use extended Taylor series directly for better precision on small z. |
| | 27 | 542 | | Fixed64 zSq = z * z; |
| | | 543 | | |
| | 27 | 544 | | Fixed64 result = z; |
| | 27 | 545 | | Fixed64 term = z; |
| | 27 | 546 | | int sign = -1; |
| | | 547 | | |
| | 184 | 548 | | for (int i = 3; i < 15; i += 2) |
| | 85 | 549 | | { |
| | 85 | 550 | | term *= zSq; |
| | 85 | 551 | | Fixed64 nextTerm = term / i; |
| | 85 | 552 | | if (nextTerm.Abs() < Fixed64.Epsilon) |
| | 20 | 553 | | break; |
| | | 554 | | |
| | 65 | 555 | | result += nextTerm * sign; |
| | 65 | 556 | | sign = -sign; |
| | 65 | 557 | | } |
| | | 558 | | |
| | 27 | 559 | | adjustedResult = result; |
| | 27 | 560 | | } |
| | | 561 | | |
| | 51 | 562 | | return neg ? -adjustedResult : adjustedResult; |
| | 78 | 563 | | } |
| | | 564 | | |
| | | 565 | | /// <summary> |
| | | 566 | | /// Computes the angle whose tangent is the quotient of two specified numbers. |
| | | 567 | | /// </summary> |
| | | 568 | | /// <remarks> |
| | | 569 | | /// Uses a fixed-point arithmetic approximation for the arc tangent function, which is more efficient than using |
| | | 570 | | /// especially on systems where floating-point operations are expensive. |
| | | 571 | | /// </remarks> |
| | | 572 | | /// <param name="y">The y-coordinate of the point to which the angle is measured.</param> |
| | | 573 | | /// <param name="x">The x-coordinate of the point to which the angle is measured.</param> |
| | | 574 | | /// <returns>An angle, θ, measured in radians, such that -π ≤ θ ≤ π, and tan(θ) = y / x, |
| | | 575 | | /// taking into account the quadrants of the inputs to determine the sign of the result.</returns> |
| | | 576 | | public static Fixed64 Atan2(Fixed64 y, Fixed64 x) |
| | 27 | 577 | | { |
| | 27 | 578 | | if (x == Fixed64.Zero) |
| | 3 | 579 | | { |
| | 3 | 580 | | if (y > Fixed64.Zero) |
| | 1 | 581 | | return PiOver2; |
| | 2 | 582 | | if (y == Fixed64.Zero) |
| | 1 | 583 | | return Fixed64.Zero; |
| | 1 | 584 | | return -PiOver2; |
| | | 585 | | } |
| | | 586 | | |
| | 24 | 587 | | Fixed64 atan = Atan(y / x); |
| | | 588 | | |
| | | 589 | | // Adjust based on the quadrant |
| | 24 | 590 | | if (x < Fixed64.Zero) |
| | 8 | 591 | | { |
| | 8 | 592 | | if (y >= Fixed64.Zero) |
| | 4 | 593 | | { |
| | | 594 | | // Second quadrant |
| | 4 | 595 | | return atan + PI; |
| | | 596 | | } |
| | | 597 | | else |
| | 4 | 598 | | { |
| | | 599 | | // Third quadrant |
| | 4 | 600 | | return atan - PI; |
| | | 601 | | } |
| | | 602 | | } |
| | | 603 | | |
| | | 604 | | // First or fourth quadrant |
| | 16 | 605 | | return atan; |
| | 27 | 606 | | } |
| | | 607 | | |
| | | 608 | | #endregion |
| | | 609 | | } |
| | | 610 | | } |