diff --git a/src/libraries/System.Private.CoreLib/src/System/Double.cs b/src/libraries/System.Private.CoreLib/src/System/Double.cs
index d4ac8530645ee1..6ea753ba60eadb 100644
--- a/src/libraries/System.Private.CoreLib/src/System/Double.cs
+++ b/src/libraries/System.Private.CoreLib/src/System/Double.cs
@@ -8,6 +8,8 @@
using System.Numerics;
using System.Runtime.CompilerServices;
using System.Runtime.InteropServices;
+using System.Runtime.Intrinsics.Arm;
+using System.Runtime.Intrinsics.X86;
using System.Runtime.Versioning;
namespace System
@@ -1972,22 +1974,242 @@ public static double CosPi(double x)
return result;
}
+ // `pi / 180` and `180 / pi` are not exactly representable, so each is carried as a
+ // `head + mid + tail` triple giving about 159 significant bits, where each term is the
+ // correctly rounded remainder left by the ones above it.
+ //
+ // Two limbs are not always enough. Their accumulated error is about `2^-53` ulp of the
+ // result, while enumerating every significand whose exact product lands nearest a rounding
+ // midpoint gives worst cases only `2^-55.58` ulp away for degrees-to-radians and `2^-55.99`
+ // ulp for the reverse. Those sit well inside that error, so two limbs cannot decide them.
+ // The third limb takes the error to about `2^-106` ulp, which clears the worst case by 50
+ // bits.
+ internal const double DegreesToRadiansHead = 0.017453292519943295; // 0x1.1DF46A2529D39p-6
+ internal const double DegreesToRadiansMid = 2.9486522708701687E-19; // 0x1.5C1D8BECDD291p-62
+ internal const double DegreesToRadiansTail = -1.3427726813345382E-35; // -0x1.1D937FA428858p-116
+ internal const double RadiansToDegreesHead = 57.29577951308232; // 0x1.CA5DC1A63C1F8p+5
+ internal const double RadiansToDegreesMid = -1.9878495670576283E-15; // -0x1.1E7AB456405F9p-49
+ internal const double RadiansToDegreesTail = -1.6833394980391744E-31; // -0x1.B505196FABB41p-103
+
+ // Recovering the roundoff of `value * mid` is what both the two limb and the three limb form
+ // rest on, and below these that roundoff falls onto the subnormal grid, losing up to one
+ // smallest subnormal. That loss has to stay under the closest an input brings the exact
+ // product to a rounding boundary, which is the `2^-55.58` and `2^-55.99` ulp above. Writing
+ // the smallest subnormal as `2^(emin - p + 1)` and an ulp of the result as
+ // `2^(e + e_head - p + 1)`, the `p` cancels and the floor on the input exponent `e` is just
+ // `emin - e_head + w`, giving `-960.42` and `-971.01`.
+ internal const double DegreesToRadiansMin = 1.0261342003245941E-289; // 0x1p-960
+ internal const double RadiansToDegreesMin = 5.010420900022432E-293; // 0x1p-971
+
+ // At or above these `value * head` overflows, which the triple cannot recover from since the
+ // roundoff of an infinite product is `inf - inf`. Paired with the minimums above they bound
+ // the domain MultiplyWide covers, which is what the vector paths range test on. `pi / 180`
+ // is under one, so nothing finite overflows there and the bound only excludes the non-finite;
+ // the single unsigned range test covers that the same as any other value.
+ internal const double DegreesToRadiansMax = PositiveInfinity;
+ internal const double RadiansToDegreesMax = 3.137566414384587E+306; // 0x1.1DF46A2529D39p+1018
+
+ // The cold path lifts `value` by `ConversionScaleUp` so the whole computation runs as it
+ // would in the normal range. Scaling by a power of two is exact in both directions, so any
+ // exponent in the valid window gives bit-identical results; only the window itself matters.
+ //
+ // The floor is the larger of two conditions, evaluated for the smallest input `2^(emin-p+1)`
+ // and taken over both directions, where `k` is the scale exponent being solved for,
+ // `p` is the precision in bits including the implicit leading bit, `emin` is the minimum
+ // normal exponent, and `e_head`/`e_mid` are the constants' own exponents. The `p - 1` term
+ // is the gap from the smallest subnormal to the smallest normal, so it is the lift needed to
+ // bring any input into the normal range:
+ //
+ // k >= (p - 1) - e_mid so `scaled * mid` stays normal
+ // k >= 2 * (p - 1) - e_head so the smallest Veltkamp sub-product in the no-FMA
+ // MultiplyRoundoff stays normal, since it is `2^-(p-1)`
+ // relative to `scaled * head`
+ //
+ // For `double` those give 114 and 110 for degrees-to-radians, 101 and 99 for the reverse,
+ // so the floor is `2^114`. The ceiling is `2^1022`, keeping `ConversionScaleDown` normal.
+ //
+ // `tail` gets no condition. Losing it outright costs `2^-110` relative, and since a result
+ // is under `2^p` ulp that is under `2^-57` ulp, which the worst case above already clears.
+ //
+ // Both conditions are sufficient rather than tight, since an input small enough to violate
+ // either produces a result with too few significand bits for what is lost to reach the
+ // rounding position.
+ private const double ConversionScaleUp = 2.076918743413931E+34; // 0x1p+114
+ private const double ConversionScaleDown = 4.81482486096809E-35; // 0x1p-114
+
+ // Keeps the top 26 significand bits, which is the Veltkamp split MultiplyRoundoff uses.
+ internal const ulong ConversionSplitMask = 0xFFFF_FFFF_F800_0000;
+
+ // Bounds the error of the two limb form relative to `value * head`, which is what lets that
+ // form tell whether it resolved the rounding. Three things are dropped or rounded away
+ // there, none of them reaching `2^-105`:
+ //
+ // `value * (C - head - mid)`, which `tail` puts under `2^-108`
+ // the roundoff of `value * mid`, so `2^-53` of a term already under `2^-54`
+ // the roundoff of their sum, so `2^-53` of a sum already under `2^-52`
+ //
+ // `2^-100` clears the total by more than 16x, leaving room for rounding the ends of the
+ // interval it forms as well. A wider bound stays correct and only sends more inputs to the
+ // three limb form; at this width that is about one in `2^46`.
+ internal const double ConversionErrorScale = 7.888609052210118E-31; // 0x1p-100
+
+ // The roundoff of `value * head` is itself exactly representable, and recovering it is
+ // what lets each limb be spent without losing what fell off the end. A fused multiply-add
+ // gives it directly; otherwise splitting both operands into 26 and 27 bit halves makes
+ // every sub-product exact, which recovers the same value using only multiplies and adds.
+ [MethodImpl(MethodImplOptions.AggressiveInlining)]
+ private static double MultiplyRoundoff(double value, double head, double product)
+ {
+ if (Fma.IsSupported || AdvSimd.Arm64.IsSupported)
+ {
+ return FusedMultiplyAdd(value, head, -product);
+ }
+
+ double headHi = BitConverter.UInt64BitsToDouble(BitConverter.DoubleToUInt64Bits(head) & ConversionSplitMask);
+ double headLo = head - headHi;
+
+ double valueHi = BitConverter.UInt64BitsToDouble(BitConverter.DoubleToUInt64Bits(value) & ConversionSplitMask);
+ double valueLo = value - valueHi;
+
+ return (((((valueHi * headHi) - product) + (valueHi * headLo)) + (valueLo * headHi)) + (valueLo * headLo));
+ }
+
+ // Computes `value * (head + mid + tail)` as the unevaluated pair `result + residual`, where
+ // `result` is the correctly rounded product. `value * head` must be known finite and
+ // non-zero, which is what lets `head` be the nearest representable value rather than being
+ // biased to keep a sign.
+ //
+ // The three limb products are accumulated exactly apart from the two terms already below
+ // `2^-106` relative, then folded back onto `product` in decreasing order of magnitude. The
+ // final fold rounds twice, so the correction is first rounded to odd, which is the standard
+ // way of making the second rounding agree with rounding the exact value once: the odd
+ // neighbour is never a midpoint, so nothing that decides the outcome can be discarded.
+ [MethodImpl(MethodImplOptions.NoInlining)]
+ private static double MultiplyWide(double value, double head, double mid, double tail, out double residual)
+ {
+ double product = value * head;
+ double productError = MultiplyRoundoff(value, head, product); // exact
+ double middle = value * mid;
+ double middleError = MultiplyRoundoff(value, mid, middle); // exact but for the
+ // subnormal loss the
+ // minimums bound
+ double bottom = middleError + (value * tail);
+
+ // The magnitudes of `productError` and `middle` are not ordered, so the sum needs the
+ // form that does not assume it
+ double sum = productError + middle;
+ double bias = sum - productError;
+ double sumError = (productError - (sum - bias)) + (middle - bias);
+ double lower = sumError + bottom;
+
+ double upper = product + sum;
+ double upperError = sum - (upper - product); // exact, |product| >= |sum|
+
+ double correction = upperError + lower;
+ double lost = (upperError - correction) + lower; // exact
+
+ if ((lost != 0) && ((BitConverter.DoubleToUInt64Bits(correction) & 1) == 0))
+ {
+ // Consecutive values alternate in significand parity, so exactly one of the two
+ // bracketing values is odd and it lies in the direction of what was lost
+ correction = (lost > 0) ? BitIncrement(correction) : BitDecrement(correction);
+ }
+
+ double result = upper + correction;
+ residual = (upper - result) + correction; // exact
+ return result;
+ }
+
+ // Computes `value * (head + mid + tail)` as a single rounding, over the domain the minimums
+ // and maximums above bound. Two limbs are enough for all but the rare input, and whether
+ // they were is itself decidable: their sum is known to within `bound`, so if both ends of
+ // that interval round to the same value then so does the exact product.
+ [MethodImpl(MethodImplOptions.AggressiveInlining)]
+ private static double MultiplyWide(double value, double head, double mid, double tail)
+ {
+ double product = value * head;
+
+ if (!IsFinite(product) || (product == 0))
+ {
+ // The roundoff of an infinite product is `inf - inf`, and a zero product would
+ // come back from the sum having lost its sign
+ return product;
+ }
+
+ double sum = MultiplyRoundoff(value, head, product) + (value * mid);
+ double bound = Abs(product) * ConversionErrorScale;
+ double result = product + (sum - bound);
+
+ if (result == (product + (sum + bound)))
+ {
+ return result;
+ }
+
+ return MultiplyWide(value, head, mid, tail, out _);
+ }
+
+ // Computes `value * (head + mid + tail)` for the inputs too small for MultiplyWide, namely
+ // those below the minimums above. Scaling up moves the product back into the range where
+ // the triple is resolved exactly; scaling back down then rounds a second time, with
+ // `residual` breaking any tie that rounding would otherwise resolve to even.
+ [MethodImpl(MethodImplOptions.NoInlining)]
+ private static double ScaleAndMultiplyWide(double value, double head, double mid, double tail)
+ {
+ double scaled = value * ConversionScaleUp;
+
+ if (scaled == 0)
+ {
+ // Multiplying keeps the sign of zero, which the sum would lose
+ return scaled * head;
+ }
+
+ double sum = MultiplyWide(scaled, head, mid, tail, out double residual);
+ double result = sum * ConversionScaleDown;
+
+ if (residual == 0)
+ {
+ return result;
+ }
+
+ double dropped = sum - (result * ConversionScaleUp); // exact
+
+ if (dropped == 0)
+ {
+ // Scaling back down lost nothing, so it was already the single rounding
+ return result;
+ }
+
+ double next = (dropped > 0) ? BitIncrement(result) : BitDecrement(result);
+
+ if ((Abs(next - result) * ConversionScaleUp) != (2.0 * Abs(dropped)))
+ {
+ return result;
+ }
+
+ // `sum` sat exactly between `result` and `next`, so the true value did not
+ return ((residual > 0) == (dropped > 0)) ? next : result;
+ }
+
///
public static double DegreesToRadians(double degrees)
{
- // NOTE: Don't change the algorithm without consulting the DIM
- // which elaborates on why this implementation was chosen
+ if (Abs(degrees) < DegreesToRadiansMin)
+ {
+ return ScaleAndMultiplyWide(degrees, DegreesToRadiansHead, DegreesToRadiansMid, DegreesToRadiansTail);
+ }
- return (degrees * Pi) / 180.0;
+ return MultiplyWide(degrees, DegreesToRadiansHead, DegreesToRadiansMid, DegreesToRadiansTail);
}
///
public static double RadiansToDegrees(double radians)
{
- // NOTE: Don't change the algorithm without consulting the DIM
- // which elaborates on why this implementation was chosen
+ if (Abs(radians) < RadiansToDegreesMin)
+ {
+ return ScaleAndMultiplyWide(radians, RadiansToDegreesHead, RadiansToDegreesMid, RadiansToDegreesTail);
+ }
- return (radians * 180.0) / Pi;
+ return MultiplyWide(radians, RadiansToDegreesHead, RadiansToDegreesMid, RadiansToDegreesTail);
}
///
diff --git a/src/libraries/System.Private.CoreLib/src/System/Half.cs b/src/libraries/System.Private.CoreLib/src/System/Half.cs
index 37fdc90a016d46..6ec122ca0b9a31 100644
--- a/src/libraries/System.Private.CoreLib/src/System/Half.cs
+++ b/src/libraries/System.Private.CoreLib/src/System/Half.cs
@@ -2324,22 +2324,10 @@ public static bool TryParsePartial(ReadOnlySpan utf8Text, NumberStyles sty
public static Half CosPi(Half x) => (Half)float.CosPi((float)x);
///
- public static Half DegreesToRadians(Half degrees)
- {
- // NOTE: Don't change the algorithm without consulting the DIM
- // which elaborates on why this implementation was chosen
-
- return (Half)float.DegreesToRadians((float)degrees);
- }
+ public static Half DegreesToRadians(Half degrees) => (Half)float.DegreesToRadians((float)degrees);
///
- public static Half RadiansToDegrees(Half radians)
- {
- // NOTE: Don't change the algorithm without consulting the DIM
- // which elaborates on why this implementation was chosen
-
- return (Half)float.RadiansToDegrees((float)radians);
- }
+ public static Half RadiansToDegrees(Half radians) => (Half)float.RadiansToDegrees((float)radians);
///
public static Half Sin(Half x) => (Half)MathF.Sin((float)x);
diff --git a/src/libraries/System.Private.CoreLib/src/System/Number.DecimalIeee754.cs b/src/libraries/System.Private.CoreLib/src/System/Number.DecimalIeee754.cs
index 033777187b62fc..eaa08464a24a67 100644
--- a/src/libraries/System.Private.CoreLib/src/System/Number.DecimalIeee754.cs
+++ b/src/libraries/System.Private.CoreLib/src/System/Number.DecimalIeee754.cs
@@ -1188,6 +1188,64 @@ internal static TValue MultiplyDecimalIeee754(TValue left, TVa
return NumberToDecimalIeee754BitsFromWide(resultSign, productHigh, productLow, productExponent, sticky: false);
}
+ ///
+ /// Multiplies the IEEE 754 decimal value represented by by the positive constant
+ /// (( × 10^Precision) + ) ×
+ /// 10^, rounds the exact product once (round-to-nearest,
+ /// ties-to-even), and returns its bit pattern.
+ ///
+ ///
+ /// Splitting the constant at the format precision is what keeps this within the double-width machinery
+ /// the arithmetic operators already use: both halves fit a single limb, so each partial product is one
+ /// , and the low product is exactly Precision digits below the
+ /// high one, so shifting it into place is a digit drop whose remainder can only ever be sticky.
+ ///
+ /// How many digits a constant needs is a property of that constant, not of the format: for an irrational
+ /// constant the closest a Precision digit significand drives the exact product to a midpoint has
+ /// no closed form, so it has to be enumerated against the truncation error once per constant rather than
+ /// carried over from a narrower format. See Decimal128 for the existing pair.
+ ///
+ internal static TValue MultiplyByWideConstantDecimalIeee754(TValue bits, TValue constantHead, TValue constantTail, int constantExponent)
+ where TDecimal : unmanaged, IDecimalIeee754ParseAndFormatInfo
+ where TValue : unmanaged, IBinaryInteger
+ {
+ if (TDecimal.IsNaN(bits))
+ {
+ return PropagateNaN(bits, bits);
+ }
+
+ if (TDecimal.IsInfinity(bits))
+ {
+ // The constant is positive and finite, so the infinity carries through with its sign
+ return TDecimal.IsNegative(bits) ? TDecimal.NegativeInfinity : TDecimal.PositiveInfinity;
+ }
+
+ DecodedDecimalIeee754 value = UnpackDecimalIeee754(bits);
+
+ int precision = TDecimal.Precision;
+
+ if (TValue.IsZero(value.Significand))
+ {
+ // The IEEE 754 preferred exponent for a product is the sum of the operand exponents, and the
+ // constant's own quantum is 10^constantExponent
+ return DecimalIeee754FiniteNumberBinaryEncoding(value.Signed, TValue.Zero, value.UnbiasedExponent + constantExponent);
+ }
+
+ // `significand × constantHead` spans at most `2p + 2` digits and `significand × constantTail` at most
+ // `2p`, so both fit the limb pair, and everything the alignment drop removes is far enough below the
+ // result's last digit to only ever reach it through `sticky`.
+ WideMultiply(value.Significand, constantHead, out TValue high, out TValue low);
+ WideMultiply(value.Significand, constantTail, out TValue tailHigh, out TValue tailLow);
+
+ bool sticky = false;
+ DropDigits(ref tailHigh, ref tailLow, precision, ref sticky, out int roundDigit);
+ sticky |= roundDigit != 0;
+
+ WideAdd(high, low, tailHigh, tailLow, out high, out low);
+
+ return NumberToDecimalIeee754BitsFromWide(value.Signed, high, low, value.UnbiasedExponent + constantExponent + precision, sticky);
+ }
+
///
/// Computes (left × right) + addend for three IEEE 754 decimal values represented by their raw bit
/// patterns, rounds the exact result once (round-to-nearest, ties-to-even), and returns its bit pattern.
@@ -2467,14 +2525,20 @@ private static void WideScaleByPow10(ref TValue high, ref TVal
private static void WideDropLowDigits(ref TValue high, ref TValue low, int dropCount, ref bool sticky)
where TValue : unmanaged, IBinaryInteger
{
- for (int i = 0; i < dropCount; i++)
+ int chunkPower = SinglePassPow10(out TValue chunkDivisor);
+
+ while (dropCount > 0)
{
if (TValue.IsZero(high) && TValue.IsZero(low))
{
break;
}
- sticky |= WideDivideByTen(ref high, ref low) != 0;
+ int step = Math.Min(dropCount, chunkPower);
+ TValue divisor = (step == chunkPower) ? chunkDivisor : WidePow10(step);
+
+ sticky |= !TValue.IsZero(WideDivideByPow10(ref high, ref low, divisor));
+ dropCount -= step;
}
}
@@ -2498,6 +2562,23 @@ private static void AlignWideToCommonExponent(ref TValue high,
}
}
+ ///
+ /// Adds the double-width magnitudes (, ) and
+ /// (, ), which the caller guarantees cannot carry
+ /// out of the high limb.
+ ///
+ private static void WideAdd(TValue leftHigh, TValue leftLow, TValue rightHigh, TValue rightLow, out TValue high, out TValue low)
+ where TValue : unmanaged, IBinaryInteger
+ {
+ low = leftLow + rightLow;
+ high = leftHigh + rightHigh;
+
+ if (low < leftLow)
+ {
+ high += TValue.One;
+ }
+ }
+
///
/// Compares two double-width magnitudes, returning a negative value, zero, or a positive value according to
/// whether the first is less than, equal to, or greater than the second.
@@ -2527,27 +2608,71 @@ private static void WideSubtract(TValue leftHigh, TValue leftLow, TValue
}
///
- /// Divides the double-width value in (, ) by ten in place,
- /// returning the discarded least-significant decimal digit. Used to strip low-order digits during rounding.
+ /// Gets the largest power of ten that a single pass can strip,
+ /// returning its exponent and, in , its value.
///
///
- /// Only the 128-bit format reaches this helper: the 32-bit and 64-bit formats widen the limb pair to a
- /// single native integer and divide directly (see and
- /// ). The Intel reference implementation avoids hardware
- /// division here by multiplying with precomputed reciprocals of powers of ten (e.g.
- /// bid_reciprocals10_64) and shifting; this helper instead uses direct integer division for
- /// simplicity, and adopting the reciprocal-multiply tables is a possible future performance optimization.
+ /// The bound is 10^k <= 2^(bits / 2), which is what keeps the running remainder times the
+ /// half-width base inside : 10^4 for a 32-bit limb, 10^9 for a
+ /// 64-bit limb, and 10^19 for a 128-bit limb.
///
- private static int WideDivideByTen(ref TValue high, ref TValue low)
+ private static int SinglePassPow10(out TValue divisor)
+ where TValue : unmanaged, IBinaryInteger
+ {
+ TValue ten = TValue.CreateTruncating(10);
+ TValue limit = (TValue.One << ((TValue.Zero.GetByteCount() * 8) / 2)) / ten;
+
+ divisor = TValue.One;
+ int power = 0;
+
+ while (divisor <= limit)
+ {
+ divisor *= ten;
+ power++;
+ }
+
+ return power;
+ }
+
+ ///
+ /// Computes 10^, which the caller guarantees fits a single limb.
+ ///
+ private static TValue WidePow10(int power)
where TValue : unmanaged, IBinaryInteger
{
TValue ten = TValue.CreateTruncating(10);
+ TValue result = TValue.One;
+ for (int i = 0; i < power; i++)
+ {
+ result *= ten;
+ }
+
+ return result;
+ }
+
+ ///
+ /// Divides the double-width value in (, ) by
+ /// in place, returning the discarded low-order digits. Used to strip low-order
+ /// digits during rounding, a whole chunk of them at a time.
+ ///
+ ///
+ /// The rounding paths only reach this helper with a two-limb value for the 128-bit format: the 32-bit and
+ /// 64-bit formats widen the limb pair to a single native integer and divide directly (see
+ /// and ). Addition
+ /// reaches it for every format through . The Intel reference
+ /// implementation avoids hardware division here by multiplying with precomputed reciprocals of powers of ten
+ /// (e.g. bid_reciprocals10_64) and shifting; this helper instead uses direct integer division for
+ /// simplicity, and adopting the reciprocal-multiply tables is a possible future performance optimization.
+ ///
+ private static TValue WideDivideByPow10(ref TValue high, ref TValue low, TValue divisor)
+ where TValue : unmanaged, IBinaryInteger
+ {
if (TValue.IsZero(high))
{
- (TValue quotient, TValue remainder) = TValue.DivRem(low, ten);
+ (TValue quotient, TValue remainder) = TValue.DivRem(low, divisor);
low = quotient;
- return int.CreateTruncating(remainder);
+ return remainder;
}
int bits = TValue.Zero.GetByteCount() * 8;
@@ -2555,17 +2680,18 @@ private static int WideDivideByTen(ref TValue high, ref TValue low)
TValue lowMask = (TValue.One << half) - TValue.One;
TValue baseValue = TValue.One << half;
- // Long division of the four half-width limbs (most significant first) by ten. Because the running
- // remainder stays below ten, `remainder * baseValue + limb` never exceeds the integer width.
+ // Long division of the four half-width limbs (most significant first). Because the running remainder
+ // stays below the divisor, which `SinglePassPow10` caps at `baseValue`, the intermediate
+ // `remainder * baseValue + limb` never exceeds the integer width.
TValue rem = TValue.Zero;
- (TValue q3, rem) = TValue.DivRem((rem * baseValue) + (high >> half), ten);
- (TValue q2, rem) = TValue.DivRem((rem * baseValue) + (high & lowMask), ten);
- (TValue q1, rem) = TValue.DivRem((rem * baseValue) + (low >> half), ten);
- (TValue q0, rem) = TValue.DivRem((rem * baseValue) + (low & lowMask), ten);
+ (TValue q3, rem) = TValue.DivRem((rem * baseValue) + (high >> half), divisor);
+ (TValue q2, rem) = TValue.DivRem((rem * baseValue) + (high & lowMask), divisor);
+ (TValue q1, rem) = TValue.DivRem((rem * baseValue) + (low >> half), divisor);
+ (TValue q0, rem) = TValue.DivRem((rem * baseValue) + (low & lowMask), divisor);
high = (q3 << half) | q2;
low = (q1 << half) | q0;
- return int.CreateTruncating(rem);
+ return rem;
}
///
@@ -2599,11 +2725,14 @@ private static int WideDigitCount(TValue high, TValue low)
}
int count = 0;
+ int chunkPower = SinglePassPow10(out TValue chunkDivisor);
+ // A non-zero high limb means the value spans at least the full limb width, which is more digits than a
+ // single pass strips, so the count stays additive across each pass.
while (!TValue.IsZero(high))
{
- WideDivideByTen(ref high, ref low);
- count++;
+ WideDivideByPow10(ref high, ref low, chunkDivisor);
+ count += chunkPower;
}
return count + TDecimal.CountDigits(low);
@@ -2735,16 +2864,22 @@ private static TValue DropDigits(ref TValue high, ref TValue l
return low;
}
- for (int i = 0; i < dropCount; i++)
+ // Every removed digit below the most-significant one only contributes stickiness, so they come off a
+ // whole single-pass chunk at a time; the rounding digit is then stripped on its own.
+ int chunkPower = SinglePassPow10(out TValue chunkDivisor);
+ int remaining = dropCount - 1;
+
+ while (remaining > 0)
{
- if (i > 0)
- {
- sticky |= roundDigit != 0;
- }
+ int step = Math.Min(remaining, chunkPower);
+ TValue divisor = (step == chunkPower) ? chunkDivisor : WidePow10(step);
- roundDigit = WideDivideByTen(ref high, ref low);
+ sticky |= !TValue.IsZero(WideDivideByPow10(ref high, ref low, divisor));
+ remaining -= step;
}
+ roundDigit = int.CreateTruncating(WideDivideByPow10(ref high, ref low, TValue.CreateTruncating(10)));
+
Debug.Assert(TValue.IsZero(high));
return low;
}
diff --git a/src/libraries/System.Private.CoreLib/src/System/Numerics/BFloat16.cs b/src/libraries/System.Private.CoreLib/src/System/Numerics/BFloat16.cs
index d97bda482c571a..ae6b65813be1f8 100644
--- a/src/libraries/System.Private.CoreLib/src/System/Numerics/BFloat16.cs
+++ b/src/libraries/System.Private.CoreLib/src/System/Numerics/BFloat16.cs
@@ -2045,22 +2045,10 @@ public static bool TryParsePartial(ReadOnlySpan utf8Text, NumberStyles sty
public static BFloat16 CosPi(BFloat16 x) => (BFloat16)float.CosPi((float)x);
///
- public static BFloat16 DegreesToRadians(BFloat16 degrees)
- {
- // NOTE: Don't change the algorithm without consulting the DIM
- // which elaborates on why this implementation was chosen
-
- return (BFloat16)float.DegreesToRadians((float)degrees);
- }
+ public static BFloat16 DegreesToRadians(BFloat16 degrees) => (BFloat16)float.DegreesToRadians((float)degrees);
///
- public static BFloat16 RadiansToDegrees(BFloat16 radians)
- {
- // NOTE: Don't change the algorithm without consulting the DIM
- // which elaborates on why this implementation was chosen
-
- return (BFloat16)float.RadiansToDegrees((float)radians);
- }
+ public static BFloat16 RadiansToDegrees(BFloat16 radians) => (BFloat16)float.RadiansToDegrees((float)radians);
///
public static BFloat16 Sin(BFloat16 x) => (BFloat16)float.Sin((float)x);
diff --git a/src/libraries/System.Private.CoreLib/src/System/Numerics/Decimal128.cs b/src/libraries/System.Private.CoreLib/src/System/Numerics/Decimal128.cs
index b91149c6148b97..e3c27363a55b0d 100644
--- a/src/libraries/System.Private.CoreLib/src/System/Numerics/Decimal128.cs
+++ b/src/libraries/System.Private.CoreLib/src/System/Numerics/Decimal128.cs
@@ -61,6 +61,25 @@ public readonly struct Decimal128
private const ulong SNaNMaskUpper = 0x7E00_0000_0000_0000;
private const ulong InfinityMaskUpper = 0x7800_0000_0000_0000;
+ // `pi / 180` and `180 / pi` cannot be written exactly in decimal, so each is kept as two
+ // integers and an exponent, `(Head * 10^34 + Tail) * 10^Exponent`. That is just the digits of
+ // the value with the point moved: `pi / 180 = 1.745_329...e-2`, `180 / pi = 5.729_577...e+1`,
+ // and since the pair is a `2p + 2` digit integer, `Exponent` is `-(2p + 1)` plus that `-2` or `+1`.
+ //
+ // `Tail` is `p` digits and `Head` is `p + 2`, since `Head` is bounded by the `UInt128` it is
+ // stored in rather than by the format. Those two extra digits are the accuracy: the product
+ // is exact and rounded only once, so the digits dropped off the constant are the only error,
+ // and they have to stay below the closest a `p` digit input can push the product to a
+ // rounding tie. Every format stops at two, `Decimal32` being the one with no room for more,
+ // since `180 / pi` at `p + 3` digits no longer fits its `uint`.
+ private static UInt128 DegreesToRadiansHead => new UInt128(0x0021_9D23_38C9_2382, 0xBF86_F9E5_ACDE_A057); // 174_532_925_199_432_957_692_369_076_848_861_271
+ private static UInt128 DegreesToRadiansTail => new UInt128(0x0000_A9BF_270D_B49E, 0x82D0_06FB_7D2C_4001); // 3_442_871_888_541_725_456_097_191_440_171_009
+ private const int DegreesToRadiansExponent = -71;
+
+ private static UInt128 RadiansToDegreesHead => new UInt128(0x006E_5900_1427_58B2, 0xDB41_A116_5AA0_7F37); // 572_957_795_130_823_208_767_981_548_141_051_703
+ private static UInt128 RadiansToDegreesTail => new UInt128(0x0000_9FC5_74E0_7F0D, 0x2D46_04B9_0C4C_EA5D); // 3_240_547_246_656_432_154_916_024_386_120_285
+ private const int RadiansToDegreesExponent = -68;
+
/// Gets a value that represents positive infinity.
public static Decimal128 PositiveInfinity => new Decimal128(PositiveInfinityValue);
@@ -954,6 +973,9 @@ bool IFloatingPoint.TryWriteSignificandLittleEndian(Span desti
///
public static Decimal128 Cosh(Decimal128 x) => new Decimal128(Number.CoshDecimalIeee754(new UInt128(x._upper, x._lower)));
+ ///
+ public static Decimal128 DegreesToRadians(Decimal128 degrees) => new Decimal128(Number.MultiplyByWideConstantDecimalIeee754(new UInt128(degrees._upper, degrees._lower), DegreesToRadiansHead, DegreesToRadiansTail, DegreesToRadiansExponent));
+
///
public static Decimal128 Exp(Decimal128 x) => new Decimal128(Number.ExpDecimalIeee754(new UInt128(x._upper, x._lower)));
@@ -1011,6 +1033,9 @@ bool IFloatingPoint.TryWriteSignificandLittleEndian(Span desti
///
public static Decimal128 Pow(Decimal128 x, Decimal128 y) => new Decimal128(Number.PowDecimalIeee754(new UInt128(x._upper, x._lower), new UInt128(y._upper, y._lower)));
+ ///
+ public static Decimal128 RadiansToDegrees(Decimal128 radians) => new Decimal128(Number.MultiplyByWideConstantDecimalIeee754(new UInt128(radians._upper, radians._lower), RadiansToDegreesHead, RadiansToDegreesTail, RadiansToDegreesExponent));
+
///
public static Decimal128 ReciprocalEstimate(Decimal128 x) => One / x;
diff --git a/src/libraries/System.Private.CoreLib/src/System/Numerics/Decimal32.cs b/src/libraries/System.Private.CoreLib/src/System/Numerics/Decimal32.cs
index b46866b0525ac3..3f165d8ade6857 100644
--- a/src/libraries/System.Private.CoreLib/src/System/Numerics/Decimal32.cs
+++ b/src/libraries/System.Private.CoreLib/src/System/Numerics/Decimal32.cs
@@ -64,6 +64,15 @@ internal Decimal32(uint value)
private const uint MaxInternalValue = 0x77F8_967F; // +9.999_999 * 10^96; aka +9_999_999 * 10^90
private const uint MinInternalValue = 0xF7F8_967F; // -9.999_999 * 10^96; aka -9_999_999 * 10^90
+ // See `Decimal128` for how these are laid out and why they carry the digits they do.
+ private const uint DegreesToRadiansHead = 174_532_925;
+ private const uint DegreesToRadiansTail = 1_994_330;
+ private const int DegreesToRadiansExponent = -17;
+
+ private const uint RadiansToDegreesHead = 572_957_795;
+ private const uint RadiansToDegreesTail = 1_308_232;
+ private const int RadiansToDegreesExponent = -14;
+
/// Gets a value that represents positive infinity.
public static Decimal32 PositiveInfinity => new Decimal32(PositiveInfinityValue);
@@ -975,6 +984,9 @@ bool IFloatingPoint.TryWriteSignificandLittleEndian(Span destin
///
public static Decimal32 Cosh(Decimal32 x) => new Decimal32(Number.CoshDecimalIeee754(x._value));
+ ///
+ public static Decimal32 DegreesToRadians(Decimal32 degrees) => new Decimal32(Number.MultiplyByWideConstantDecimalIeee754(degrees._value, DegreesToRadiansHead, DegreesToRadiansTail, DegreesToRadiansExponent));
+
///
public static Decimal32 Exp(Decimal32 x) => new Decimal32(Number.ExpDecimalIeee754(x._value));
@@ -1032,6 +1044,9 @@ bool IFloatingPoint.TryWriteSignificandLittleEndian(Span destin
///
public static Decimal32 Pow(Decimal32 x, Decimal32 y) => new Decimal32(Number.PowDecimalIeee754(x._value, y._value));
+ ///
+ public static Decimal32 RadiansToDegrees(Decimal32 radians) => new Decimal32(Number.MultiplyByWideConstantDecimalIeee754(radians._value, RadiansToDegreesHead, RadiansToDegreesTail, RadiansToDegreesExponent));
+
///
public static Decimal32 ReciprocalEstimate(Decimal32 x) => One / x;
diff --git a/src/libraries/System.Private.CoreLib/src/System/Numerics/Decimal64.cs b/src/libraries/System.Private.CoreLib/src/System/Numerics/Decimal64.cs
index 56fde15e2c7d9c..fbc675d8d688bd 100644
--- a/src/libraries/System.Private.CoreLib/src/System/Numerics/Decimal64.cs
+++ b/src/libraries/System.Private.CoreLib/src/System/Numerics/Decimal64.cs
@@ -59,6 +59,15 @@ public readonly struct Decimal64
private const ulong MaxInternalValue = 0x77FB_86F2_6FC0_FFFF; // 9.999_999_999_999_999 * 10^384; aka 9_999_999_999_999_999 * 10^369
private const ulong MinInternalValue = 0xF7FB_86F2_6FC0_FFFF; // -9.999_999_999_999_999 * 10^384; aka -9_999_999_999_999_999 * 10^369
+ // See `Decimal128` for how these are laid out and why they carry the digits they do.
+ private const ulong DegreesToRadiansHead = 174_532_925_199_432_957;
+ private const ulong DegreesToRadiansTail = 6_923_690_768_488_613;
+ private const int DegreesToRadiansExponent = -35;
+
+ private const ulong RadiansToDegreesHead = 572_957_795_130_823_208;
+ private const ulong RadiansToDegreesTail = 7_679_815_481_410_517;
+ private const int RadiansToDegreesExponent = -32;
+
/// Gets a value that represents positive infinity.
public static Decimal64 PositiveInfinity => new Decimal64(PositiveInfinityValue);
@@ -968,6 +977,9 @@ bool IFloatingPoint.TryWriteSignificandLittleEndian(Span destin
///
public static Decimal64 Cosh(Decimal64 x) => new Decimal64(Number.CoshDecimalIeee754(x._value));
+ ///
+ public static Decimal64 DegreesToRadians(Decimal64 degrees) => new Decimal64(Number.MultiplyByWideConstantDecimalIeee754(degrees._value, DegreesToRadiansHead, DegreesToRadiansTail, DegreesToRadiansExponent));
+
///
public static Decimal64 Exp(Decimal64 x) => new Decimal64(Number.ExpDecimalIeee754(x._value));
@@ -1025,6 +1037,9 @@ bool IFloatingPoint.TryWriteSignificandLittleEndian(Span destin
///
public static Decimal64 Pow(Decimal64 x, Decimal64 y) => new Decimal64(Number.PowDecimalIeee754(x._value, y._value));
+ ///
+ public static Decimal64 RadiansToDegrees(Decimal64 radians) => new Decimal64(Number.MultiplyByWideConstantDecimalIeee754(radians._value, RadiansToDegreesHead, RadiansToDegreesTail, RadiansToDegreesExponent));
+
///
public static Decimal64 ReciprocalEstimate(Decimal64 x) => One / x;
diff --git a/src/libraries/System.Private.CoreLib/src/System/Numerics/ITrigonometricFunctions.cs b/src/libraries/System.Private.CoreLib/src/System/Numerics/ITrigonometricFunctions.cs
index eedb7952551548..7a1a1f7926fc5b 100644
--- a/src/libraries/System.Private.CoreLib/src/System/Numerics/ITrigonometricFunctions.cs
+++ b/src/libraries/System.Private.CoreLib/src/System/Numerics/ITrigonometricFunctions.cs
@@ -109,6 +109,17 @@ static virtual TSelf DegreesToRadians(TSelf degrees)
// This is particularly true for `RadiansToDegrees` where 180 is exactly representable
// and so allows an exactly representable intermediate value to be computed when overflow
// doesn't occur.
+ //
+ // Which of the two forms wins is decided by where the nearest representable `pi / 180`
+ // falls relative to the infinitely precise value, which lands more or less uniformly for
+ // each number of significand bits. Folding therefore ranges from far worse to far better
+ // depending on the type, while this form stays between 55% and 75% correctly rounded
+ // everywhere. A DIM cannot pick the better half of that, since it can only reach
+ // `TSelf.Pi` and so must round twice to get there.
+ //
+ // A concrete type has neither constraint. It can name the constant directly and knows
+ // what its hardware and any wider type give it, so it can be correctly rounded over the
+ // whole domain, which every type here does and any other implementation is free to do.
return (degrees * TSelf.Pi) / TSelf.CreateChecked(180);
}
diff --git a/src/libraries/System.Private.CoreLib/src/System/Numerics/Vector.cs b/src/libraries/System.Private.CoreLib/src/System/Numerics/Vector.cs
index 636e38112b3169..6dad48f30f9526 100644
--- a/src/libraries/System.Private.CoreLib/src/System/Numerics/Vector.cs
+++ b/src/libraries/System.Private.CoreLib/src/System/Numerics/Vector.cs
@@ -1180,7 +1180,7 @@ public static Vector DegreesToRadians(Vector degrees)
{
if (IsHardwareAccelerated)
{
- return VectorMath.DegreesToRadians, double>(degrees);
+ return VectorMath.DegreesToRadiansDouble, Vector>(degrees);
}
else
{
@@ -1195,7 +1195,7 @@ public static Vector DegreesToRadians(Vector degrees)
{
if (IsHardwareAccelerated)
{
- return VectorMath.DegreesToRadians, float>(degrees);
+ return VectorMath.DegreesToRadiansSingle, Vector>(degrees);
}
else
{
@@ -2848,7 +2848,7 @@ public static Vector RadiansToDegrees(Vector radians)
{
if (IsHardwareAccelerated)
{
- return VectorMath.RadiansToDegrees, double>(radians);
+ return VectorMath.RadiansToDegreesDouble, Vector>(radians);
}
else
{
@@ -2862,7 +2862,7 @@ public static Vector RadiansToDegrees(Vector radians)
{
if (IsHardwareAccelerated)
{
- return VectorMath.RadiansToDegrees, float>(radians);
+ return VectorMath.RadiansToDegreesSingle, Vector>(radians);
}
else
{
diff --git a/src/libraries/System.Private.CoreLib/src/System/Runtime/Intrinsics/Vector128.cs b/src/libraries/System.Private.CoreLib/src/System/Runtime/Intrinsics/Vector128.cs
index 829430f245d40f..bdce7f65bb07ac 100644
--- a/src/libraries/System.Private.CoreLib/src/System/Runtime/Intrinsics/Vector128.cs
+++ b/src/libraries/System.Private.CoreLib/src/System/Runtime/Intrinsics/Vector128.cs
@@ -1736,7 +1736,7 @@ public static Vector128 DegreesToRadians(Vector128 degrees)
{
if (IsHardwareAccelerated)
{
- return VectorMath.DegreesToRadians, double>(degrees);
+ return VectorMath.DegreesToRadiansDouble, Vector128>(degrees);
}
else
{
@@ -1754,7 +1754,14 @@ public static Vector128 DegreesToRadians(Vector128 degrees)
{
if (IsHardwareAccelerated)
{
- return VectorMath.DegreesToRadians, float>(degrees);
+ if (Vector256.IsHardwareAccelerated)
+ {
+ return VectorMath.DegreesToRadiansSingle, Vector256>(degrees);
+ }
+ else
+ {
+ return VectorMath.DegreesToRadiansSingle, Vector128>(degrees);
+ }
}
else
{
@@ -3110,7 +3117,7 @@ public static Vector128 RadiansToDegrees(Vector128 radians)
{
if (IsHardwareAccelerated)
{
- return VectorMath.RadiansToDegrees, double>(radians);
+ return VectorMath.RadiansToDegreesDouble, Vector128>(radians);
}
else
{
@@ -3128,7 +3135,14 @@ public static Vector128 RadiansToDegrees(Vector128 radians)
{
if (IsHardwareAccelerated)
{
- return VectorMath.RadiansToDegrees, float>(radians);
+ if (Vector256.IsHardwareAccelerated)
+ {
+ return VectorMath.RadiansToDegreesSingle, Vector256>(radians);
+ }
+ else
+ {
+ return VectorMath.RadiansToDegreesSingle, Vector128>(radians);
+ }
}
else
{
diff --git a/src/libraries/System.Private.CoreLib/src/System/Runtime/Intrinsics/Vector256.cs b/src/libraries/System.Private.CoreLib/src/System/Runtime/Intrinsics/Vector256.cs
index 9134f25f408861..171f95fa258dc4 100644
--- a/src/libraries/System.Private.CoreLib/src/System/Runtime/Intrinsics/Vector256.cs
+++ b/src/libraries/System.Private.CoreLib/src/System/Runtime/Intrinsics/Vector256.cs
@@ -1830,7 +1830,7 @@ public static Vector256 DegreesToRadians(Vector256 degrees)
{
if (IsHardwareAccelerated)
{
- return VectorMath.DegreesToRadians, double>(degrees);
+ return VectorMath.DegreesToRadiansDouble, Vector256>(degrees);
}
else
{
@@ -1848,7 +1848,14 @@ public static Vector256 DegreesToRadians(Vector256 degrees)
{
if (IsHardwareAccelerated)
{
- return VectorMath.DegreesToRadians, float>(degrees);
+ if (Vector512.IsHardwareAccelerated)
+ {
+ return VectorMath.DegreesToRadiansSingle, Vector512>(degrees);
+ }
+ else
+ {
+ return VectorMath.DegreesToRadiansSingle, Vector256>(degrees);
+ }
}
else
{
@@ -3153,7 +3160,7 @@ public static Vector256 RadiansToDegrees(Vector256 radians)
{
if (IsHardwareAccelerated)
{
- return VectorMath.RadiansToDegrees, double>(radians);
+ return VectorMath.RadiansToDegreesDouble, Vector256>(radians);
}
else
{
@@ -3171,7 +3178,14 @@ public static Vector256 RadiansToDegrees(Vector256 radians)
{
if (IsHardwareAccelerated)
{
- return VectorMath.RadiansToDegrees, float>(radians);
+ if (Vector512.IsHardwareAccelerated)
+ {
+ return VectorMath.RadiansToDegreesSingle, Vector512>(radians);
+ }
+ else
+ {
+ return VectorMath.RadiansToDegreesSingle, Vector256>(radians);
+ }
}
else
{
diff --git a/src/libraries/System.Private.CoreLib/src/System/Runtime/Intrinsics/Vector512.cs b/src/libraries/System.Private.CoreLib/src/System/Runtime/Intrinsics/Vector512.cs
index 5ab9fb26a39d43..ccf5ffd1e9c8bd 100644
--- a/src/libraries/System.Private.CoreLib/src/System/Runtime/Intrinsics/Vector512.cs
+++ b/src/libraries/System.Private.CoreLib/src/System/Runtime/Intrinsics/Vector512.cs
@@ -1885,7 +1885,7 @@ public static Vector512 DegreesToRadians(Vector512 degrees)
{
if (IsHardwareAccelerated)
{
- return VectorMath.DegreesToRadians, double>(degrees);
+ return VectorMath.DegreesToRadiansDouble, Vector512>(degrees);
}
else
{
@@ -1903,7 +1903,7 @@ public static Vector512 DegreesToRadians(Vector512 degrees)
{
if (IsHardwareAccelerated)
{
- return VectorMath.DegreesToRadians, float>(degrees);
+ return VectorMath.DegreesToRadiansSingle, Vector512>(degrees);
}
else
{
@@ -3194,7 +3194,7 @@ public static Vector512 RadiansToDegrees(Vector512 radians)
{
if (IsHardwareAccelerated)
{
- return VectorMath.RadiansToDegrees, double>(radians);
+ return VectorMath.RadiansToDegreesDouble, Vector512>(radians);
}
else
{
@@ -3212,7 +3212,7 @@ public static Vector512 RadiansToDegrees(Vector512 radians)
{
if (IsHardwareAccelerated)
{
- return VectorMath.RadiansToDegrees, float>(radians);
+ return VectorMath.RadiansToDegreesSingle, Vector512>(radians);
}
else
{
diff --git a/src/libraries/System.Private.CoreLib/src/System/Runtime/Intrinsics/Vector64.cs b/src/libraries/System.Private.CoreLib/src/System/Runtime/Intrinsics/Vector64.cs
index 3639b3ae6551e6..f4b04d311cfc29 100644
--- a/src/libraries/System.Private.CoreLib/src/System/Runtime/Intrinsics/Vector64.cs
+++ b/src/libraries/System.Private.CoreLib/src/System/Runtime/Intrinsics/Vector64.cs
@@ -1679,7 +1679,7 @@ public static Vector64 DegreesToRadians(Vector64 degrees)
{
if (IsHardwareAccelerated)
{
- return VectorMath.DegreesToRadians, double>(degrees);
+ return VectorMath.DegreesToRadiansDouble, Vector64>(degrees);
}
else
{
@@ -1696,7 +1696,14 @@ public static Vector64 DegreesToRadians(Vector64 degrees)
{
if (IsHardwareAccelerated)
{
- return VectorMath.DegreesToRadians, float>(degrees);
+ if (Vector128.IsHardwareAccelerated)
+ {
+ return VectorMath.DegreesToRadiansSingle, Vector128>(degrees);
+ }
+ else
+ {
+ return VectorMath.DegreesToRadiansSingle, Vector64>(degrees);
+ }
}
else
{
@@ -3277,7 +3284,7 @@ public static Vector64 RadiansToDegrees(Vector64 radians)
{
if (IsHardwareAccelerated)
{
- return VectorMath.RadiansToDegrees, double>(radians);
+ return VectorMath.RadiansToDegreesDouble, Vector64>(radians);
}
else
{
@@ -3294,7 +3301,14 @@ public static Vector64 RadiansToDegrees(Vector64 radians)
{
if (IsHardwareAccelerated)
{
- return VectorMath.RadiansToDegrees, float>(radians);
+ if (Vector128.IsHardwareAccelerated)
+ {
+ return VectorMath.RadiansToDegreesSingle, Vector128>(radians);
+ }
+ else
+ {
+ return VectorMath.RadiansToDegreesSingle, Vector64>(radians);
+ }
}
else
{
diff --git a/src/libraries/System.Private.CoreLib/src/System/Runtime/Intrinsics/VectorMath.cs b/src/libraries/System.Private.CoreLib/src/System/Runtime/Intrinsics/VectorMath.cs
index bdbe668efa326b..510206a9ad3d1c 100644
--- a/src/libraries/System.Private.CoreLib/src/System/Runtime/Intrinsics/VectorMath.cs
+++ b/src/libraries/System.Private.CoreLib/src/System/Runtime/Intrinsics/VectorMath.cs
@@ -338,15 +338,126 @@ public static TVector CopySign(TVector value, TVector sign)
}
}
+ // Computes `value * (head + mid)`, and reports per element whether that resolved the
+ // rounding. See `double.MultiplyWide` for the interval this tests.
[MethodImpl(MethodImplOptions.AggressiveInlining)]
- public static TVector DegreesToRadians(TVector degrees)
- where TVector : unmanaged, ISimdVector
- where T : IFloatingPointIeee754
+ private static TVectorDouble MultiplyWideDouble(TVectorDouble value, double head, double mid, out TVectorUInt64 resolved)
+ where TVectorDouble : unmanaged, ISimdVector
+ where TVectorUInt64 : unmanaged, ISimdVector
+ {
+ TVectorDouble vhead = Create(head);
+
+ TVectorDouble product = value * vhead;
+ TVectorDouble sum = MultiplyRoundoffDouble(value, head, vhead, product)
+ + (value * Create(mid));
+
+ TVectorDouble bound = TVectorDouble.Abs(product) * Create(double.ConversionErrorScale);
+ TVectorDouble result = product + (sum - bound);
+
+ // Handed back as a mask rather than reduced here, so the caller can fold it into the
+ // range test and branch once, and as a `ulong` mask because reducing the same all ones
+ // pattern as `double` is several times more expensive for no difference in meaning.
+ resolved = Unsafe.BitCast(TVectorDouble.Equals(result, product + (sum + bound)));
+ return result;
+ }
+
+ // The vector form of `double.MultiplyRoundoff`, and it has to be that rather than
+ // `MultiplyAddEstimate`: that one is free to evaluate as a separate multiply and add, and
+ // then it returns `product - product`, silently dropping the term the two limb form is
+ // built on. `head` is a constant, so only `value` is split at runtime.
+ [MethodImpl(MethodImplOptions.AggressiveInlining)]
+ private static TVectorDouble MultiplyRoundoffDouble(TVectorDouble value, double head, TVectorDouble vhead, TVectorDouble product)
+ where TVectorDouble : unmanaged, ISimdVector
+ where TVectorUInt64 : unmanaged, ISimdVector
+ {
+ if (Fma.IsSupported || AdvSimd.Arm64.IsSupported)
+ {
+ return FusedMultiplyAdd(value, vhead, -product);
+ }
+
+ double headHi = BitConverter.UInt64BitsToDouble(BitConverter.DoubleToUInt64Bits(head) & double.ConversionSplitMask);
+
+ TVectorDouble vheadHi = Create(headHi);
+ TVectorDouble vheadLo = Create(head - headHi);
+
+ TVectorDouble valueHi = Unsafe.BitCast(
+ Unsafe.BitCast(value) & TVectorUInt64.Create(double.ConversionSplitMask)
+ );
+ TVectorDouble valueLo = value - valueHi;
+
+ return ((((valueHi * vheadHi) - product) + (valueHi * vheadLo)) + (valueLo * vheadHi)) + (valueLo * vheadLo);
+ }
+
+ [MethodImpl(MethodImplOptions.AggressiveInlining)]
+ public static TVectorDouble DegreesToRadiansDouble(TVectorDouble degrees)
+ where TVectorDouble : unmanaged, ISimdVector
+ where TVectorUInt64 : unmanaged, ISimdVector
+ {
+ // See `double.DegreesToRadians` for why the constant is carried as a `head + mid + tail`
+ // triple, and `double.DegreesToRadiansMin` and `Max` for the window the two limb form
+ // used here is valid over.
+ //
+ // That window is a range test rather than a pair of compares: bit patterns of
+ // non-negative values increase with magnitude, so subtracting the low end makes one
+ // unsigned compare reject both ends at once, zero and the non-finite values included.
+ // An element outside it, or one the two limb form cannot resolve, sends the whole
+ // vector to the scalar path rather than growing a third limb across it. Both are
+ // elementwise masks, so they reduce to one branch together.
+
+ TVectorUInt64 ux = Unsafe.BitCast(TVectorDouble.Abs(degrees))
+ - TVectorUInt64.Create(BitConverter.DoubleToUInt64Bits(double.DegreesToRadiansMin));
+
+ TVectorUInt64 window = TVectorUInt64.Create(BitConverter.DoubleToUInt64Bits(double.DegreesToRadiansMax)
+ - BitConverter.DoubleToUInt64Bits(double.DegreesToRadiansMin));
+
+ TVectorDouble result = MultiplyWideDouble(
+ degrees,
+ double.DegreesToRadiansHead,
+ double.DegreesToRadiansMid,
+ out TVectorUInt64 resolved
+ );
+
+ if (TVectorUInt64.EqualsAll(resolved & TVectorUInt64.LessThan(ux, window), TVectorUInt64.AllBitsSet))
+ {
+ return result;
+ }
+
+ return ScalarFallback(degrees);
+
+ static TVectorDouble ScalarFallback(TVectorDouble x)
+ {
+ TVectorDouble result = TVectorDouble.Zero;
+
+ for (int i = 0; i < TVectorDouble.ElementCount; i++)
+ {
+ double scalar = double.DegreesToRadians(x[i]);
+ result = result.WithElement(i, scalar);
+ }
+
+ return result;
+ }
+ }
+
+ [MethodImpl(MethodImplOptions.AggressiveInlining)]
+ public static TVectorSingle DegreesToRadiansSingle(TVectorSingle degrees)
+ where TVectorSingle : unmanaged, ISimdVector
+ where TVectorDouble : unmanaged, ISimdVector
{
- // NOTE: Don't change the algorithm without consulting the DIM
- // which elaborates on why this implementation was chosen
+ // See `float.DegreesToRadians` for why widening is correctly rounded here.
- return (degrees * TVector.Create(T.Pi)) / TVector.Create(T.CreateTruncating(180));
+ TVectorDouble scale = Create(double.DegreesToRadiansHead);
+
+ if (TVectorSingle.ElementCount == TVectorDouble.ElementCount)
+ {
+ return Narrow(
+ Widen(degrees) * scale
+ );
+ }
+
+ return Narrow(
+ WidenLower(degrees) * scale,
+ WidenUpper(degrees) * scale
+ );
}
public static TVectorDouble ExpDouble(TVectorDouble x)
@@ -1681,14 +1792,66 @@ public static TVector MinNumber(TVector x, TVector y)
}
[MethodImpl(MethodImplOptions.AggressiveInlining)]
- public static TVector RadiansToDegrees(TVector radians)
- where TVector : unmanaged, ISimdVector
- where T : IFloatingPointIeee754
+ public static TVectorDouble RadiansToDegreesDouble(TVectorDouble radians)
+ where TVectorDouble : unmanaged, ISimdVector
+ where TVectorUInt64 : unmanaged, ISimdVector
{
- // NOTE: Don't change the algorithm without consulting the DIM
- // which elaborates on why this implementation was chosen
+ // See `DegreesToRadiansDouble` for the shape of this.
+
+ TVectorUInt64 ux = Unsafe.BitCast(TVectorDouble.Abs(radians))
+ - TVectorUInt64.Create(BitConverter.DoubleToUInt64Bits(double.RadiansToDegreesMin));
- return (radians * TVector.Create(T.CreateTruncating(180))) / TVector.Create(T.Pi);
+ TVectorUInt64 window = TVectorUInt64.Create(BitConverter.DoubleToUInt64Bits(double.RadiansToDegreesMax)
+ - BitConverter.DoubleToUInt64Bits(double.RadiansToDegreesMin));
+
+ TVectorDouble result = MultiplyWideDouble(
+ radians,
+ double.RadiansToDegreesHead,
+ double.RadiansToDegreesMid,
+ out TVectorUInt64 resolved
+ );
+
+ if (TVectorUInt64.EqualsAll(resolved & TVectorUInt64.LessThan(ux, window), TVectorUInt64.AllBitsSet))
+ {
+ return result;
+ }
+
+ return ScalarFallback(radians);
+
+ static TVectorDouble ScalarFallback(TVectorDouble x)
+ {
+ TVectorDouble result = TVectorDouble.Zero;
+
+ for (int i = 0; i < TVectorDouble.ElementCount; i++)
+ {
+ double scalar = double.RadiansToDegrees(x[i]);
+ result = result.WithElement(i, scalar);
+ }
+
+ return result;
+ }
+ }
+
+ [MethodImpl(MethodImplOptions.AggressiveInlining)]
+ public static TVectorSingle RadiansToDegreesSingle(TVectorSingle radians)
+ where TVectorSingle : unmanaged, ISimdVector
+ where TVectorDouble : unmanaged, ISimdVector
+ {
+ // See `DegreesToRadiansSingle` for the shape of this.
+
+ TVectorDouble scale = Create(double.RadiansToDegreesHead);
+
+ if (TVectorSingle.ElementCount == TVectorDouble.ElementCount)
+ {
+ return Narrow(
+ Widen(radians) * scale
+ );
+ }
+
+ return Narrow(
+ WidenLower(radians) * scale,
+ WidenUpper(radians) * scale
+ );
}
public static TVectorDouble RoundDouble(TVectorDouble vector, MidpointRounding mode)
@@ -2511,6 +2674,40 @@ private static TVectorDouble CosSingleSmall(TVectorDouble x)
return TVectorDouble.MultiplyAddEstimate(CosSinglePoly(x), x4, s);
}
+ [MethodImpl(MethodImplOptions.AggressiveInlining)]
+ private static TVectorDouble FusedMultiplyAdd(TVectorDouble left, TVectorDouble right, TVectorDouble addend)
+ where TVectorDouble : unmanaged, ISimdVector
+ {
+ Unsafe.SkipInit(out TVectorDouble result);
+
+ if (typeof(TVectorDouble) == typeof(Vector))
+ {
+ result = (TVectorDouble)(object)Vector.FusedMultiplyAdd((Vector)(object)left, (Vector)(object)right, (Vector)(object)addend);
+ }
+ else if (typeof(TVectorDouble) == typeof(Vector64))
+ {
+ result = (TVectorDouble)(object)Vector64.FusedMultiplyAdd((Vector64)(object)left, (Vector64)(object)right, (Vector64)(object)addend);
+ }
+ else if (typeof(TVectorDouble) == typeof(Vector128))
+ {
+ result = (TVectorDouble)(object)Vector128.FusedMultiplyAdd((Vector128)(object)left, (Vector128)(object)right, (Vector128)(object)addend);
+ }
+ else if (typeof(TVectorDouble) == typeof(Vector256))
+ {
+ result = (TVectorDouble)(object)Vector256.FusedMultiplyAdd((Vector256)(object)left, (Vector256)(object)right, (Vector256)(object)addend);
+ }
+ else if (typeof(TVectorDouble) == typeof(Vector512))
+ {
+ result = (TVectorDouble)(object)Vector512.FusedMultiplyAdd((Vector512)(object)left, (Vector512)(object)right, (Vector512