src/libraries/System.Private.CoreLib/src/System/MathF.cs

   1 // Licensed to the .NET Foundation under one or more agreements.
   2 // The .NET Foundation licenses this file to you under the MIT license.
   3
   4 // ===================================================================================================
   5 // Portions of the code implemented below are based on the 'Berkeley SoftFloat Release 3e' algorithms.
   6 // ===================================================================================================
   7
   8 /*============================================================
   9 **
  10 ** Purpose: Some single-precision floating-point math operations
  11 **
  12 ===========================================================*/
  13
  14 using System.Diagnostics;
  15 using System.Numerics;
  16 using System.Runtime.CompilerServices;
  17 using System.Runtime.Intrinsics;
  18 using System.Runtime.Intrinsics.X86;
  19 using System.Runtime.Intrinsics.Arm;
  20
  21 namespace System
  22 {
  23     public static partial class MathF
  24     {
  25         public const float E = 2.71828183f;
  26
  27         public const float PI = 3.14159265f;
  28
  29         public const float Tau = 6.283185307f;
  30
  31         private const int maxRoundingDigits = 6;
  32
  33         // This table is required for the Round function which can specify the number of digits to round to
  34         private static readonly float[] roundPower10Single = new float[] {
  35             1e0f, 1e1f, 1e2f, 1e3f, 1e4f, 1e5f, 1e6f
  36         };
  37
  38         private const float singleRoundLimit = 1e8f;
  39
  40         private const float SCALEB_C1 = 1.7014118E+38f; // 0x1p127f
  41
  42         private const float SCALEB_C2 = 1.1754944E-38f; // 0x1p-126f
  43
  44         private const float SCALEB_C3 = 16777216f; // 0x1p24f
  45
  46         [MethodImpl(MethodImplOptions.AggressiveInlining)]
  47         public static float Abs(float x)
  48         {
  49             return Math.Abs(x);
  50         }
  51
  52         public static float BitDecrement(float x)
  53         {
  54             int bits = BitConverter.SingleToInt32Bits(x);
  55
  56             if ((bits & 0x7F800000) >= 0x7F800000)
  57             {
  58                 // NaN returns NaN
  59                 // -Infinity returns -Infinity
  60                 // +Infinity returns float.MaxValue
  61                 return (bits == 0x7F800000) ? float.MaxValue : x;
  62             }
  63
  64             if (bits == 0x00000000)
  65             {
  66                 // +0.0 returns -float.Epsilon
  67                 return -float.Epsilon;
  68             }
  69
  70             // Negative values need to be incremented
  71             // Positive values need to be decremented
  72
  73             bits += ((bits < 0) ? +1 : -1);
  74             return BitConverter.Int32BitsToSingle(bits);
  75         }
  76
  77         public static float BitIncrement(float x)
  78         {
  79             int bits = BitConverter.SingleToInt32Bits(x);
  80
  81             if ((bits & 0x7F800000) >= 0x7F800000)
  82             {
  83                 // NaN returns NaN
  84                 // -Infinity returns float.MinValue
  85                 // +Infinity returns +Infinity
  86                 return (bits == unchecked((int)(0xFF800000))) ? float.MinValue : x;
  87             }
  88
  89             if (bits == unchecked((int)(0x80000000)))
  90             {
  91                 // -0.0 returns float.Epsilon
  92                 return float.Epsilon;
  93             }
  94
  95             // Negative values need to be decremented
  96             // Positive values need to be incremented
  97
  98             bits += ((bits < 0) ? -1 : +1);
  99             return BitConverter.Int32BitsToSingle(bits);
 100         }
 101
 102         [MethodImpl(MethodImplOptions.AggressiveInlining)]
 103         public static float CopySign(float x, float y)
 104         {
 105             if (Sse.IsSupported || AdvSimd.IsSupported)
 106             {
 107                 return VectorMath.ConditionalSelectBitwise(Vector128.CreateScalarUnsafe(-0.0f), Vector128.CreateScalarUnsafe(y), Vector128.CreateScalarUnsafe(x)).ToScalar();
 108             }
 109             else
 110             {
 111                 return SoftwareFallback(x, y);
 112             }
 113
 114             static float SoftwareFallback(float x, float y)
 115             {
 116                 const int signMask = 1 << 31;
 117
 118                 // This method is required to work for all inputs,
 119                 // including NaN, so we operate on the raw bits.
 120                 int xbits = BitConverter.SingleToInt32Bits(x);
 121                 int ybits = BitConverter.SingleToInt32Bits(y);
 122
 123                 // Remove the sign from x, and remove everything but the sign from y
 124                 xbits &= ~signMask;
 125                 ybits &= signMask;
 126
 127                 // Simply OR them to get the correct sign
 128                 return BitConverter.Int32BitsToSingle(xbits | ybits);
 129             }
 130         }
 131
 132         public static float IEEERemainder(float x, float y)
 133         {
 134             if (float.IsNaN(x))
 135             {
 136                 return x; // IEEE 754-2008: NaN payload must be preserved
 137             }
 138
 139             if (float.IsNaN(y))
 140             {
 141                 return y; // IEEE 754-2008: NaN payload must be preserved
 142             }
 143
 144             float regularMod = x % y;
 145
 146             if (float.IsNaN(regularMod))
 147             {
 148                 return float.NaN;
 149             }
 150
 151             if ((regularMod == 0) && float.IsNegative(x))
 152             {
 153                 return float.NegativeZero;
 154             }
 155
 156             float alternativeResult = (regularMod - (Abs(y) * Sign(x)));
 157
 158             if (Abs(alternativeResult) == Abs(regularMod))
 159             {
 160                 float divisionResult = x / y;
 161                 float roundedResult = Round(divisionResult);
 162
 163                 if (Abs(roundedResult) > Abs(divisionResult))
 164                 {
 165                     return alternativeResult;
 166                 }
 167                 else
 168                 {
 169                     return regularMod;
 170                 }
 171             }
 172
 173             if (Abs(alternativeResult) < Abs(regularMod))
 174             {
 175                 return alternativeResult;
 176             }
 177             else
 178             {
 179                 return regularMod;
 180             }
 181         }
 182
 183         public static float Log(float x, float y)
 184         {
 185             if (float.IsNaN(x))
 186             {
 187                 return x; // IEEE 754-2008: NaN payload must be preserved
 188             }
 189
 190             if (float.IsNaN(y))
 191             {
 192                 return y; // IEEE 754-2008: NaN payload must be preserved
 193             }
 194
 195             if (y == 1)
 196             {
 197                 return float.NaN;
 198             }
 199
 200             if ((x != 1) && ((y == 0) || float.IsPositiveInfinity(y)))
 201             {
 202                 return float.NaN;
 203             }
 204
 205             return Log(x) / Log(y);
 206         }
 207
 208         [MethodImpl(MethodImplOptions.AggressiveInlining)]
 209         public static float Max(float x, float y)
 210         {
 211             return Math.Max(x, y);
 212         }
 213
 214         public static float MaxMagnitude(float x, float y)
 215         {
 216             // This matches the IEEE 754:2019 `maximumMagnitude` function
 217             //
 218             // It propagates NaN inputs back to the caller and
 219             // otherwise returns the input with a larger magnitude.
 220             // It treats +0 as larger than -0 as per the specification.
 221
 222             float ax = Abs(x);
 223             float ay = Abs(y);
 224
 225             if ((ax > ay) || float.IsNaN(ax))
 226             {
 227                 return x;
 228             }
 229
 230             if (ax == ay)
 231             {
 232                 return float.IsNegative(x) ? y : x;
 233             }
 234
 235             return y;
 236         }
 237
 238         [MethodImpl(MethodImplOptions.AggressiveInlining)]
 239         public static float Min(float x, float y)
 240         {
 241             return Math.Min(x, y);
 242         }
 243
 244         public static float MinMagnitude(float x, float y)
 245         {
 246             // This matches the IEEE 754:2019 `minimumMagnitude` function
 247             //
 248             // It propagates NaN inputs back to the caller and
 249             // otherwise returns the input with a larger magnitude.
 250             // It treats +0 as larger than -0 as per the specification.
 251
 252             float ax = Abs(x);
 253             float ay = Abs(y);
 254
 255             if ((ax < ay) || float.IsNaN(ax))
 256             {
 257                 return x;
 258             }
 259
 260             if (ax == ay)
 261             {
 262                 return float.IsNegative(x) ? x : y;
 263             }
 264
 265             return y;
 266         }
 267
 268         /// <summary>Returns an estimate of the reciprocal of a specified number.</summary>
 269         /// <param name="x">The number whose reciprocal is to be estimated.</param>
 270         /// <returns>An estimate of the reciprocal of <paramref name="x" />.</returns>
 271         /// <remarks>
 272         ///    <para>On x86/x64 hardware this may use the <c>RCPSS</c> instruction which has a maximum relative error of <c>1.5 * 2^-12</c>.</para>
 273         ///    <para>On ARM64 hardware this may use the <c>FRECPE</c> instruction which performs a single Newton-Raphson iteration.</para>
 274         ///    <para>On hardware without specialized support, this may just return <c>1.0 / x</c>.</para>
 275         /// </remarks>
 276         [MethodImpl(MethodImplOptions.AggressiveInlining)]
 277         public static float ReciprocalEstimate(float x)
 278         {
 279             if (Sse.IsSupported)
 280             {
 281                 return Sse.ReciprocalScalar(Vector128.CreateScalarUnsafe(x)).ToScalar();
 282             }
 283             else if (AdvSimd.Arm64.IsSupported)
 284             {
 285                 return AdvSimd.Arm64.ReciprocalEstimateScalar(Vector64.CreateScalarUnsafe(x)).ToScalar();
 286             }
 287             else
 288             {
 289                 return 1.0f / x;
 290             }
 291         }
 292
 293         /// <summary>Returns an estimate of the reciprocal square root of a specified number.</summary>
 294         /// <param name="x">The number whose reciprocal square root is to be estimated.</param>
 295         /// <returns>An estimate of the reciprocal square root <paramref name="x" />.</returns>
 296         /// <remarks>
 297         ///    <para>On x86/x64 hardware this may use the <c>RSQRTSS</c> instruction which has a maximum relative error of <c>1.5 * 2^-12</c>.</para>
 298         ///    <para>On ARM64 hardware this may use the <c>FRSQRTE</c> instruction which performs a single Newton-Raphson iteration.</para>
 299         ///    <para>On hardware without specialized support, this may just return <c>1.0 / Sqrt(x)</c>.</para>
 300         /// </remarks>
 301         [MethodImpl(MethodImplOptions.AggressiveInlining)]
 302         public static float ReciprocalSqrtEstimate(float x)
 303         {
 304             if (Sse.IsSupported)
 305             {
 306                 return Sse.ReciprocalSqrtScalar(Vector128.CreateScalarUnsafe(x)).ToScalar();
 307             }
 308             else if (AdvSimd.Arm64.IsSupported)
 309             {
 310                 return AdvSimd.Arm64.ReciprocalSquareRootEstimateScalar(Vector64.CreateScalarUnsafe(x)).ToScalar();
 311             }
 312             else
 313             {
 314                 return 1.0f / Sqrt(x);
 315             }
 316         }
 317
 318         [Intrinsic]
 319         public static float Round(float x)
 320         {
 321             // ************************************************************************************
 322             // IMPORTANT: Do not change this implementation without also updating MathF.Round(float),
 323             //            FloatingPointUtils::round(double), and FloatingPointUtils::round(float)
 324             // ************************************************************************************
 325
 326             // This is based on the 'Berkeley SoftFloat Release 3e' algorithm
 327
 328             uint bits = (uint)BitConverter.SingleToInt32Bits(x);
 329             int exponent = float.ExtractExponentFromBits(bits);
 330
 331             if (exponent <= 0x7E)
 332             {
 333                 if ((bits << 1) == 0)
 334                 {
 335                     // Exactly +/- zero should return the original value
 336                     return x;
 337                 }
 338
 339                 // Any value less than or equal to 0.5 will always round to exactly zero
 340                 // and any value greater than 0.5 will always round to exactly one. However,
 341                 // we need to preserve the original sign for IEEE compliance.
 342
 343                 float result = ((exponent == 0x7E) && (float.ExtractSignificandFromBits(bits) != 0)) ? 1.0f : 0.0f;
 344                 return CopySign(result, x);
 345             }
 346
 347             if (exponent >= 0x96)
 348             {
 349                 // Any value greater than or equal to 2^23 cannot have a fractional part,
 350                 // So it will always round to exactly itself.
 351
 352                 return x;
 353             }
 354
 355             // The absolute value should be greater than or equal to 1.0 and less than 2^23
 356             Debug.Assert((0x7F <= exponent) && (exponent <= 0x95));
 357
 358             // Determine the last bit that represents the integral portion of the value
 359             // and the bits representing the fractional portion
 360
 361             uint lastBitMask = 1U << (0x96 - exponent);
 362             uint roundBitsMask = lastBitMask - 1;
 363
 364             // Increment the first fractional bit, which represents the midpoint between
 365             // two integral values in the current window.
 366
 367             bits += lastBitMask >> 1;
 368
 369             if ((bits & roundBitsMask) == 0)
 370             {
 371                 // If that overflowed and the rest of the fractional bits are zero
 372                 // then we were exactly x.5 and we want to round to the even result
 373
 374                 bits &= ~lastBitMask;
 375             }
 376             else
 377             {
 378                 // Otherwise, we just want to strip the fractional bits off, truncating
 379                 // to the current integer value.
 380
 381                 bits &= ~roundBitsMask;
 382             }
 383
 384             return BitConverter.Int32BitsToSingle((int)bits);
 385         }
 386
 387         [MethodImpl(MethodImplOptions.AggressiveInlining)]
 388         public static float Round(float x, int digits)
 389         {
 390             return Round(x, digits, MidpointRounding.ToEven);
 391         }
 392
 393         [MethodImpl(MethodImplOptions.AggressiveInlining)]
 394         public static float Round(float x, MidpointRounding mode)
 395         {
 396             return Round(x, 0, mode);
 397         }
 398
 399         public static unsafe float Round(float x, int digits, MidpointRounding mode)
 400         {
 401             if ((digits < 0) || (digits > maxRoundingDigits))
 402             {
 403                 throw new ArgumentOutOfRangeException(nameof(digits), SR.ArgumentOutOfRange_RoundingDigits_MathF);
 404             }
 405
 406             if (mode < MidpointRounding.ToEven || mode > MidpointRounding.ToPositiveInfinity)
 407             {
 408                 throw new ArgumentException(SR.Format(SR.Argument_InvalidEnumValue, mode, nameof(MidpointRounding)), nameof(mode));
 409             }
 410
 411             if (Abs(x) < singleRoundLimit)
 412             {
 413                 float power10 = roundPower10Single[digits];
 414
 415                 x *= power10;
 416
 417                 switch (mode)
 418                 {
 419                     // Rounds to the nearest value; if the number falls midway,
 420                     // it is rounded to the nearest value with an even least significant digit
 421                     case MidpointRounding.ToEven:
 422                     {
 423                         x = Round(x);
 424                         break;
 425                     }
 426                     // Rounds to the nearest value; if the number falls midway,
 427                     // it is rounded to the nearest value above (for positive numbers) or below (for negative numbers)
 428                     case MidpointRounding.AwayFromZero:
 429                     {
 430                         float fraction = ModF(x, &x);
 431
 432                         if (Abs(fraction) >= 0.5f)
 433                         {
 434                             x += Sign(fraction);
 435                         }
 436
 437                         break;
 438                     }
 439                     // Directed rounding: Round to the nearest value, toward to zero
 440                     case MidpointRounding.ToZero:
 441                     {
 442                         x = Truncate(x);
 443                         break;
 444                     }
 445                     // Directed Rounding: Round down to the next value, toward negative infinity
 446                     case MidpointRounding.ToNegativeInfinity:
 447                     {
 448                         x = Floor(x);
 449                         break;
 450                     }
 451                     // Directed rounding: Round up to the next value, toward positive infinity
 452                     case MidpointRounding.ToPositiveInfinity:
 453                     {
 454                         x = Ceiling(x);
 455                         break;
 456                     }
 457                     default:
 458                     {
 459                         throw new ArgumentException(SR.Format(SR.Argument_InvalidEnumValue, mode, nameof(MidpointRounding)), nameof(mode));
 460                     }
 461                 }
 462
 463                 x /= power10;
 464             }
 465
 466             return x;
 467         }
 468
 469         [MethodImpl(MethodImplOptions.AggressiveInlining)]
 470         public static int Sign(float x)
 471         {
 472             return Math.Sign(x);
 473         }
 474
 475         public static unsafe float Truncate(float x)
 476         {
 477             ModF(x, &x);
 478             return x;
 479         }
 480
 481         public static float ScaleB(float x, int n)
 482         {
 483             // Implementation based on https://git.musl-libc.org/cgit/musl/tree/src/math/scalblnf.c
 484             //
 485             // Performs the calculation x * 2^n efficiently. It constructs a float from 2^n by building
 486             // the correct biased exponent. If n is greater than the maximum exponent (127) or less than
 487             // the minimum exponent (-126), adjust x and n to compute correct result.
 488
 489             float y = x;
 490             if (n > 127)
 491             {
 492                 y *= SCALEB_C1;
 493                 n -= 127;
 494                 if (n > 127)
 495                 {
 496                     y *= SCALEB_C1;
 497                     n -= 127;
 498                     if (n > 127)
 499                     {
 500                         n = 127;
 501                     }
 502                 }
 503             }
 504             else if (n < -126)
 505             {
 506                 y *= SCALEB_C2 * SCALEB_C3;
 507                 n += 126 - 24;
 508                 if (n < -126)
 509                 {
 510                     y *= SCALEB_C2 * SCALEB_C3;
 511                     n += 126 - 24;
 512                     if (n < -126)
 513                     {
 514                         n = -126;
 515                     }
 516                 }
 517             }
 518
 519             float u = BitConverter.Int32BitsToSingle(((int)(0x7f + n) << 23));
 520             return y * u;
 521         }
 522     }
 523 }