diff --git a/src/libraries/System.Numerics.Tensors/src/System/Numerics/Tensors/netcore/TensorPrimitives.Exp.cs b/src/libraries/System.Numerics.Tensors/src/System/Numerics/Tensors/netcore/TensorPrimitives.Exp.cs index 1147ec5bea5f7..c89ad27074270 100644 --- a/src/libraries/System.Numerics.Tensors/src/System/Numerics/Tensors/netcore/TensorPrimitives.Exp.cs +++ b/src/libraries/System.Numerics.Tensors/src/System/Numerics/Tensors/netcore/TensorPrimitives.Exp.cs @@ -157,9 +157,6 @@ public static Vector512 Invoke(Vector512 x) private const ulong V_ARG_MAX = 0x40862000_00000000; private const ulong V_DP64_BIAS = 1023; - private const double V_EXPF_MIN = -709.782712893384; - private const double V_EXPF_MAX = +709.782712893384; - private const double V_EXPF_HUGE = 6755399441055744; private const double V_TBL_LN2 = 1.4426950408889634; @@ -183,155 +180,145 @@ public static Vector512 Invoke(Vector512 x) public static Vector128 Invoke(Vector128 x) { - // x * (64.0 / ln(2)) - Vector128 z = x * Vector128.Create(V_TBL_LN2); - - Vector128 dn = z + Vector128.Create(V_EXPF_HUGE); + // Check if -709 < vx < 709 + if (Vector128.LessThanOrEqualAll(Vector128.Abs(x).AsUInt64(), Vector128.Create(V_ARG_MAX))) + { + // x * (64.0 / ln(2)) + Vector128 z = x * Vector128.Create(V_TBL_LN2); - // n = (int)z - Vector128 n = dn.AsUInt64(); + Vector128 dn = z + Vector128.Create(V_EXPF_HUGE); - // dn = (double)n - dn -= Vector128.Create(V_EXPF_HUGE); + // n = (int)z + Vector128 n = dn.AsUInt64(); - // r = x - (dn * (ln(2) / 64)) - // where ln(2) / 64 is split into Head and Tail values - Vector128 r = x - (dn * Vector128.Create(V_LN2_HEAD)) - (dn * Vector128.Create(V_LN2_TAIL)); + // dn = (double)n + dn -= Vector128.Create(V_EXPF_HUGE); - Vector128 r2 = r * r; - Vector128 r4 = r2 * r2; - Vector128 r8 = r4 * r4; + // r = x - (dn * (ln(2) / 64)) + // where ln(2) / 64 is split into Head and Tail values + Vector128 r = x - (dn * Vector128.Create(V_LN2_HEAD)) - (dn * Vector128.Create(V_LN2_TAIL)); - // Compute polynomial - Vector128 poly = ((Vector128.Create(C12) * r + Vector128.Create(C11)) * r2 + - Vector128.Create(C10) * r + Vector128.Create(C9)) * r8 + - ((Vector128.Create(C8) * r + Vector128.Create(C7)) * r2 + - (Vector128.Create(C6) * r + Vector128.Create(C5))) * r4 + - ((Vector128.Create(C4) * r + Vector128.Create(C3)) * r2 + (r + Vector128.One)); + Vector128 r2 = r * r; + Vector128 r4 = r2 * r2; + Vector128 r8 = r4 * r4; - // m = (n - j) / 64 - // result = polynomial * 2^m - Vector128 ret = poly * ((n + Vector128.Create(V_DP64_BIAS)) << 52).AsDouble(); + // Compute polynomial + Vector128 poly = ((Vector128.Create(C12) * r + Vector128.Create(C11)) * r2 + + Vector128.Create(C10) * r + Vector128.Create(C9)) * r8 + + ((Vector128.Create(C8) * r + Vector128.Create(C7)) * r2 + + (Vector128.Create(C6) * r + Vector128.Create(C5))) * r4 + + ((Vector128.Create(C4) * r + Vector128.Create(C3)) * r2 + (r + Vector128.One)); - // Check if -709 < vx < 709 - if (Vector128.GreaterThanAny(Vector128.Abs(x).AsUInt64(), Vector128.Create(V_ARG_MAX))) + // m = (n - j) / 64 + // result = polynomial * 2^m + return poly * ((n + Vector128.Create(V_DP64_BIAS)) << 52).AsDouble(); + } + else { - // (x > V_EXPF_MAX) ? double.PositiveInfinity : x - Vector128 infinityMask = Vector128.GreaterThan(x, Vector128.Create(V_EXPF_MAX)); - - ret = Vector128.ConditionalSelect( - infinityMask, - Vector128.Create(double.PositiveInfinity), - ret - ); + return ScalarFallback(x); - // (x < V_EXPF_MIN) ? 0 : x - ret = Vector128.AndNot(ret, Vector128.LessThan(x, Vector128.Create(V_EXPF_MIN))); + static Vector128 ScalarFallback(Vector128 x) => + Vector128.Create(Math.Exp(x.GetElement(0)), + Math.Exp(x.GetElement(1))); } - - return ret; } public static Vector256 Invoke(Vector256 x) { - // x * (64.0 / ln(2)) - Vector256 z = x * Vector256.Create(V_TBL_LN2); - - Vector256 dn = z + Vector256.Create(V_EXPF_HUGE); + // Check if -709 < vx < 709 + if (Vector256.LessThanOrEqualAll(Vector256.Abs(x).AsUInt64(), Vector256.Create(V_ARG_MAX))) + { + // x * (64.0 / ln(2)) + Vector256 z = x * Vector256.Create(V_TBL_LN2); - // n = (int)z - Vector256 n = dn.AsUInt64(); + Vector256 dn = z + Vector256.Create(V_EXPF_HUGE); - // dn = (double)n - dn -= Vector256.Create(V_EXPF_HUGE); + // n = (int)z + Vector256 n = dn.AsUInt64(); - // r = x - (dn * (ln(2) / 64)) - // where ln(2) / 64 is split into Head and Tail values - Vector256 r = x - (dn * Vector256.Create(V_LN2_HEAD)) - (dn * Vector256.Create(V_LN2_TAIL)); + // dn = (double)n + dn -= Vector256.Create(V_EXPF_HUGE); - Vector256 r2 = r * r; - Vector256 r4 = r2 * r2; - Vector256 r8 = r4 * r4; + // r = x - (dn * (ln(2) / 64)) + // where ln(2) / 64 is split into Head and Tail values + Vector256 r = x - (dn * Vector256.Create(V_LN2_HEAD)) - (dn * Vector256.Create(V_LN2_TAIL)); - // Compute polynomial - Vector256 poly = ((Vector256.Create(C12) * r + Vector256.Create(C11)) * r2 + - Vector256.Create(C10) * r + Vector256.Create(C9)) * r8 + - ((Vector256.Create(C8) * r + Vector256.Create(C7)) * r2 + - (Vector256.Create(C6) * r + Vector256.Create(C5))) * r4 + - ((Vector256.Create(C4) * r + Vector256.Create(C3)) * r2 + (r + Vector256.One)); + Vector256 r2 = r * r; + Vector256 r4 = r2 * r2; + Vector256 r8 = r4 * r4; - // m = (n - j) / 64 - // result = polynomial * 2^m - Vector256 ret = poly * ((n + Vector256.Create(V_DP64_BIAS)) << 52).AsDouble(); + // Compute polynomial + Vector256 poly = ((Vector256.Create(C12) * r + Vector256.Create(C11)) * r2 + + Vector256.Create(C10) * r + Vector256.Create(C9)) * r8 + + ((Vector256.Create(C8) * r + Vector256.Create(C7)) * r2 + + (Vector256.Create(C6) * r + Vector256.Create(C5))) * r4 + + ((Vector256.Create(C4) * r + Vector256.Create(C3)) * r2 + (r + Vector256.One)); - // Check if -709 < vx < 709 - if (Vector256.GreaterThanAny(Vector256.Abs(x).AsUInt64(), Vector256.Create(V_ARG_MAX))) + // m = (n - j) / 64 + // result = polynomial * 2^m + return poly * ((n + Vector256.Create(V_DP64_BIAS)) << 52).AsDouble(); + } + else { - // (x > V_EXPF_MAX) ? double.PositiveInfinity : x - Vector256 infinityMask = Vector256.GreaterThan(x, Vector256.Create(V_EXPF_MAX)); + return ScalarFallback(x); - ret = Vector256.ConditionalSelect( - infinityMask, - Vector256.Create(double.PositiveInfinity), - ret - ); - - // (x < V_EXPF_MIN) ? 0 : x - ret = Vector256.AndNot(ret, Vector256.LessThan(x, Vector256.Create(V_EXPF_MIN))); + static Vector256 ScalarFallback(Vector256 x) => + Vector256.Create(Math.Exp(x.GetElement(0)), + Math.Exp(x.GetElement(1)), + Math.Exp(x.GetElement(2)), + Math.Exp(x.GetElement(3))); } - - return ret; } public static Vector512 Invoke(Vector512 x) { - // x * (64.0 / ln(2)) - Vector512 z = x * Vector512.Create(V_TBL_LN2); - - Vector512 dn = z + Vector512.Create(V_EXPF_HUGE); + // Check if -709 < vx < 709 + if (Vector512.LessThanOrEqualAll(Vector512.Abs(x).AsUInt64(), Vector512.Create(V_ARG_MAX))) + { + // x * (64.0 / ln(2)) + Vector512 z = x * Vector512.Create(V_TBL_LN2); - // n = (int)z - Vector512 n = dn.AsUInt64(); + Vector512 dn = z + Vector512.Create(V_EXPF_HUGE); - // dn = (double)n - dn -= Vector512.Create(V_EXPF_HUGE); + // n = (int)z + Vector512 n = dn.AsUInt64(); - // r = x - (dn * (ln(2) / 64)) - // where ln(2) / 64 is split into Head and Tail values - Vector512 r = x - (dn * Vector512.Create(V_LN2_HEAD)) - (dn * Vector512.Create(V_LN2_TAIL)); + // dn = (double)n + dn -= Vector512.Create(V_EXPF_HUGE); - Vector512 r2 = r * r; - Vector512 r4 = r2 * r2; - Vector512 r8 = r4 * r4; + // r = x - (dn * (ln(2) / 64)) + // where ln(2) / 64 is split into Head and Tail values + Vector512 r = x - (dn * Vector512.Create(V_LN2_HEAD)) - (dn * Vector512.Create(V_LN2_TAIL)); - // Compute polynomial - Vector512 poly = ((Vector512.Create(C12) * r + Vector512.Create(C11)) * r2 + - Vector512.Create(C10) * r + Vector512.Create(C9)) * r8 + - ((Vector512.Create(C8) * r + Vector512.Create(C7)) * r2 + - (Vector512.Create(C6) * r + Vector512.Create(C5))) * r4 + - ((Vector512.Create(C4) * r + Vector512.Create(C3)) * r2 + (r + Vector512.One)); + Vector512 r2 = r * r; + Vector512 r4 = r2 * r2; + Vector512 r8 = r4 * r4; - // m = (n - j) / 64 - // result = polynomial * 2^m - Vector512 ret = poly * ((n + Vector512.Create(V_DP64_BIAS)) << 52).AsDouble(); + // Compute polynomial + Vector512 poly = ((Vector512.Create(C12) * r + Vector512.Create(C11)) * r2 + + Vector512.Create(C10) * r + Vector512.Create(C9)) * r8 + + ((Vector512.Create(C8) * r + Vector512.Create(C7)) * r2 + + (Vector512.Create(C6) * r + Vector512.Create(C5))) * r4 + + ((Vector512.Create(C4) * r + Vector512.Create(C3)) * r2 + (r + Vector512.One)); - // Check if -709 < vx < 709 - if (Vector512.GreaterThanAny(Vector512.Abs(x).AsUInt64(), Vector512.Create(V_ARG_MAX))) + // m = (n - j) / 64 + // result = polynomial * 2^m + return poly * ((n + Vector512.Create(V_DP64_BIAS)) << 52).AsDouble(); + } + else { - // (x > V_EXPF_MAX) ? double.PositiveInfinity : x - Vector512 infinityMask = Vector512.GreaterThan(x, Vector512.Create(V_EXPF_MAX)); - - ret = Vector512.ConditionalSelect( - infinityMask, - Vector512.Create(double.PositiveInfinity), - ret - ); - - // (x < V_EXPF_MIN) ? 0 : x - ret = Vector512.AndNot(ret, Vector512.LessThan(x, Vector512.Create(V_EXPF_MIN))); + return ScalarFallback(x); + + static Vector512 ScalarFallback(Vector512 x) => + Vector512.Create(Math.Exp(x.GetElement(0)), + Math.Exp(x.GetElement(1)), + Math.Exp(x.GetElement(2)), + Math.Exp(x.GetElement(3)), + Math.Exp(x.GetElement(4)), + Math.Exp(x.GetElement(5)), + Math.Exp(x.GetElement(6)), + Math.Exp(x.GetElement(7))); } - - return ret; } } diff --git a/src/libraries/System.Numerics.Tensors/tests/TensorPrimitivesTests.cs b/src/libraries/System.Numerics.Tensors/tests/TensorPrimitivesTests.cs index 278ab39938b41..920ed00af7e90 100644 --- a/src/libraries/System.Numerics.Tensors/tests/TensorPrimitivesTests.cs +++ b/src/libraries/System.Numerics.Tensors/tests/TensorPrimitivesTests.cs @@ -163,9 +163,17 @@ protected T NextRandom(T avoid) /// the value is stored into a random position in , and the original /// value is subsequently restored. /// - protected void RunForEachSpecialValue(Action action, BoundedMemory x) + protected void RunForEachSpecialValue(Action action, BoundedMemory x) => + RunForEachSpecialValue(action, x, GetSpecialValues()); + + /// + /// Runs the specified action for each special value. Before the action is invoked, + /// the value is stored into a random position in , and the original + /// value is subsequently restored. + /// + protected void RunForEachSpecialValue(Action action, BoundedMemory x, IEnumerable specialValues) { - Assert.All(GetSpecialValues(), value => + Assert.All(specialValues, value => { int pos = Random.Next(x.Length); T orig = x[pos]; @@ -1021,6 +1029,17 @@ public void Exp_SpecialValues() using BoundedMemory x = CreateAndFillTensor(tensorLength); using BoundedMemory destination = CreateTensor(tensorLength); + T[] additionalSpecialValues = + [ + typeof(T) == typeof(float) ? (T)(object)-709.7f : + typeof(T) == typeof(double) ? (T)(object)-709.7 : + default, + + typeof(T) == typeof(float) ? (T)(object)709.7f : + typeof(T) == typeof(double) ? (T)(object)709.7 : + default, + ]; + RunForEachSpecialValue(() => { Exp(x, destination); @@ -1028,7 +1047,7 @@ public void Exp_SpecialValues() { AssertEqualTolerance(Exp(x[i]), destination[i]); } - }, x); + }, x, GetSpecialValues().Concat(additionalSpecialValues)); }); } 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 bc0d1b7c82f84..fbd5af0fa1b78 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 @@ -391,9 +391,6 @@ public static TVectorDouble ExpDouble(TVectorDoubl const ulong V_ARG_MAX = 0x40862000_00000000; const ulong V_DP64_BIAS = 1023; - const double V_EXPF_MIN = -709.782712893384; - const double V_EXPF_MAX = +709.782712893384; - const double V_EXPF_HUGE = 6755399441055744; const double V_TBL_LN2 = 1.4426950408889634; @@ -411,66 +408,69 @@ public static TVectorDouble ExpDouble(TVectorDoubl const double C11 = 2.7632293298250954E-07; const double C12 = 2.499430431958571E-08; - // x * (64.0 / ln(2)) - TVectorDouble dn = TVectorDouble.MultiplyAddEstimate(x, TVectorDouble.Create(V_TBL_LN2), TVectorDouble.Create(V_EXPF_HUGE)); + // Check if -709 < vx < 709 + if (TVectorUInt64.LessThanOrEqualAll(Unsafe.BitCast(TVectorDouble.Abs(x)), TVectorUInt64.Create(V_ARG_MAX))) + { + // x * (64.0 / ln(2)) + TVectorDouble dn = TVectorDouble.MultiplyAddEstimate(x, TVectorDouble.Create(V_TBL_LN2), TVectorDouble.Create(V_EXPF_HUGE)); - // n = (int)z - TVectorUInt64 n = Unsafe.BitCast(dn); + // n = (int)z + TVectorUInt64 n = Unsafe.BitCast(dn); - // dn = (double)n - dn -= TVectorDouble.Create(V_EXPF_HUGE); + // dn = (double)n + dn -= TVectorDouble.Create(V_EXPF_HUGE); - // r = x - (dn * (ln(2) / 64)) - // where ln(2) / 64 is split into Head and Tail values - TVectorDouble r = TVectorDouble.MultiplyAddEstimate(dn, TVectorDouble.Create(-V_LN2_HEAD), x); - r = TVectorDouble.MultiplyAddEstimate(dn, TVectorDouble.Create(-V_LN2_TAIL), r); + // r = x - (dn * (ln(2) / 64)) + // where ln(2) / 64 is split into Head and Tail values + TVectorDouble r = TVectorDouble.MultiplyAddEstimate(dn, TVectorDouble.Create(-V_LN2_HEAD), x); + r = TVectorDouble.MultiplyAddEstimate(dn, TVectorDouble.Create(-V_LN2_TAIL), r); - TVectorDouble r2 = r * r; - TVectorDouble r4 = r2 * r2; - TVectorDouble r8 = r4 * r4; + TVectorDouble r2 = r * r; + TVectorDouble r4 = r2 * r2; + TVectorDouble r8 = r4 * r4; - // Compute polynomial - TVectorDouble poly = TVectorDouble.MultiplyAddEstimate( - TVectorDouble.MultiplyAddEstimate( - TVectorDouble.MultiplyAddEstimate(TVectorDouble.Create(C12), r, TVectorDouble.Create(C11)), - r2, - TVectorDouble.MultiplyAddEstimate(TVectorDouble.Create(C10), r, TVectorDouble.Create(C09))), - r8, - TVectorDouble.MultiplyAddEstimate( + // Compute polynomial + TVectorDouble poly = TVectorDouble.MultiplyAddEstimate( TVectorDouble.MultiplyAddEstimate( - TVectorDouble.MultiplyAddEstimate(TVectorDouble.Create(C08), r, TVectorDouble.Create(C07)), + TVectorDouble.MultiplyAddEstimate(TVectorDouble.Create(C12), r, TVectorDouble.Create(C11)), r2, - TVectorDouble.MultiplyAddEstimate(TVectorDouble.Create(C06), r, TVectorDouble.Create(C05))), - r4, + TVectorDouble.MultiplyAddEstimate(TVectorDouble.Create(C10), r, TVectorDouble.Create(C09))), + r8, TVectorDouble.MultiplyAddEstimate( - TVectorDouble.MultiplyAddEstimate(TVectorDouble.Create(C04), r, TVectorDouble.Create(C03)), - r2, - r + TVectorDouble.One + TVectorDouble.MultiplyAddEstimate( + TVectorDouble.MultiplyAddEstimate(TVectorDouble.Create(C08), r, TVectorDouble.Create(C07)), + r2, + TVectorDouble.MultiplyAddEstimate(TVectorDouble.Create(C06), r, TVectorDouble.Create(C05))), + r4, + TVectorDouble.MultiplyAddEstimate( + TVectorDouble.MultiplyAddEstimate(TVectorDouble.Create(C04), r, TVectorDouble.Create(C03)), + r2, + r + TVectorDouble.One + ) ) - ) - ); - - // m = (n - j) / 64 - // result = polynomial * 2^m - TVectorDouble result = poly * Unsafe.BitCast((n + TVectorUInt64.Create(V_DP64_BIAS)) << 52); + ); - // Check if -709 < vx < 709 - if (TVectorUInt64.GreaterThanAny(Unsafe.BitCast(TVectorDouble.Abs(x)), TVectorUInt64.Create(V_ARG_MAX))) + // m = (n - j) / 64 + // result = polynomial * 2^m + return poly * Unsafe.BitCast((n + TVectorUInt64.Create(V_DP64_BIAS)) << 52); + } + else { - // (x > V_EXPF_MAX) ? double.PositiveInfinity : x - TVectorDouble infinityMask = TVectorDouble.GreaterThan(x, TVectorDouble.Create(V_EXPF_MAX)); + return ScalarFallback(x); - result = TVectorDouble.ConditionalSelect( - infinityMask, - TVectorDouble.Create(double.PositiveInfinity), - result - ); + static TVectorDouble ScalarFallback(TVectorDouble x) + { + TVectorDouble expResult = TVectorDouble.Zero; - // (x < V_EXPF_MIN) ? 0 : x - result = TVectorDouble.AndNot(result, TVectorDouble.LessThan(x, TVectorDouble.Create(V_EXPF_MIN))); - } + for (int i = 0; i < TVectorDouble.Count; i++) + { + double expScalar = double.Exp(x[i]); + expResult = expResult.WithElement(i, expScalar); + } - return result; + return expResult; + } + } } public static TVectorSingle ExpSingle(TVectorSingle x)