Unexpected result with `Vector4.Exp`

[kristiker]

Does anyone know why vector exponential

Vector4.Exp(Vector4.Zero);

gives <1.0000001, 1.0000001, 1.0000001, 1.0000001>

while

MathF.Exp(0);

gives 1.0f

Where is the error coming from?

Check the results for youself here:
https://sharplab.io/#v2:C4LgTgrgdgNAJiA1AHwAICYCMBYAUKgBgAJVMA6AOQgFsBTMASwGMBnAbj0JPICVpgGdMgEkowRlBbN2eTpgCcACgCyAQ2AALAGJkAogA8ADooIBKUx3wKV6jXqMnzluUoBqtJsAD2YACz3jd08ffwAtei8nF0Ug7zBMdAAOAJiPOLCIsgBBFlifBMTFc2zctJCiizwgA===

[EgorBo]

cc @tannergooding

[dotnet-policy-service]

Tagging subscribers to this area: @dotnet/area-system-numerics
See info in area-owners.md if you want to be subscribed.

[huoyaoyuan]

They are using different math implementations. Vector4.Exp ends up calling here:

runtime/src/libraries/System.Private.CoreLib/src/System/Runtime/Intrinsics/VectorMath.cs

Lines 475 to 590 in cff2063

	public static TVectorSingle ExpSingle<TVectorSingle, TVectorUInt32, TVectorDouble, TVectorUInt64>(TVectorSingle x)
	where TVectorSingle : unmanaged, ISimdVector<TVectorSingle, float>
	where TVectorUInt32 : unmanaged, ISimdVector<TVectorUInt32, uint>
	where TVectorDouble : unmanaged, ISimdVector<TVectorDouble, double>
	where TVectorUInt64 : unmanaged, ISimdVector<TVectorUInt64, ulong>
	{
	// This code is based on `vrs4_expf` from amd/aocl-libm-ose
	// Copyright (C) 2019-2022 Advanced Micro Devices, Inc. All rights reserved.
	//
	// Licensed under the BSD 3-Clause "New" or "Revised" License
	// See THIRD-PARTY-NOTICES.TXT for the full license text
	// Implementation Notes:
	// 1. Argument Reduction:
	// e^x = 2^(x/ln2) --- (1)
	//
	// Let x/ln(2) = z --- (2)
	//
	// Let z = n + r , where n is an integer --- (3)
	// |r| <= 1/2
	//
	// From (1), (2) and (3),
	// e^x = 2^z
	// = 2^(N+r)
	// = (2^N)*(2^r) --- (4)
	//
	// 2. Polynomial Evaluation
	// From (4),
	// r = z - N
	// 2^r = C1 + C2*r + C3*r^2 + C4*r^3 + C5 *r^4 + C6*r^5
	//
	// 4. Reconstruction
	// Thus,
	// e^x = (2^N) * (2^r)
	const uint V_ARG_MAX = 0x42AE0000;
	const float V_EXPF_MIN = -103.97208f;
	const float V_EXPF_MAX = +88.72284f;
	const double V_EXPF_HUGE = 6755399441055744;
	const double V_TBL_LN2 = 1.4426950408889634;
	const double C1 = 1.0000000754895704;
	const double C2 = 0.6931472254087585;
	const double C3 = 0.2402210737432219;
	const double C4 = 0.05550297297702539;
	const double C5 = 0.009676036358193323;
	const double C6 = 0.001341000536524434;
	TVectorSingle result;
	if (TVectorSingle.ElementCount == TVectorDouble.ElementCount)
	{
	result = Narrow<TVectorDouble, TVectorSingle>(
	CoreImpl(Widen<TVectorSingle, TVectorDouble>(x))
	);
	}
	else
	{
	result = Narrow<TVectorDouble, TVectorSingle>(
	CoreImpl(WidenLower<TVectorSingle, TVectorDouble>(x)),
	CoreImpl(WidenUpper<TVectorSingle, TVectorDouble>(x))
	);
	}
	// Check if -103 < |x| < 88
	if (TVectorUInt32.GreaterThanAny(Unsafe.BitCast<TVectorSingle, TVectorUInt32>(TVectorSingle.Abs(x)), TVectorUInt32.Create(V_ARG_MAX)))
	{
	// (x > V_EXPF_MAX) ? float.PositiveInfinity : x
	TVectorSingle infinityMask = TVectorSingle.GreaterThan(x, TVectorSingle.Create(V_EXPF_MAX));
	result = TVectorSingle.ConditionalSelect(
	infinityMask,
	TVectorSingle.Create(float.PositiveInfinity),
	result
	);
	// (x < V_EXPF_MIN) ? 0 : x
	result = TVectorSingle.AndNot(result, TVectorSingle.LessThan(x, TVectorSingle.Create(V_EXPF_MIN)));
	}
	return result;
	[MethodImpl(MethodImplOptions.AggressiveInlining)]
	static TVectorDouble CoreImpl(TVectorDouble x)
	{
	// x * (64.0 / ln(2))
	TVectorDouble z = x * TVectorDouble.Create(V_TBL_LN2);
	TVectorDouble v_expf_huge = TVectorDouble.Create(V_EXPF_HUGE);
	TVectorDouble dn = z + v_expf_huge;
	// n = (int)z
	TVectorUInt64 n = Unsafe.BitCast<TVectorDouble, TVectorUInt64>(dn);
	// r = z - n
	TVectorDouble r = z - (dn - v_expf_huge);
	TVectorDouble r2 = r * r;
	TVectorDouble r4 = r2 * r2;
	TVectorDouble poly = TVectorDouble.MultiplyAddEstimate(
	r4,
	TVectorDouble.MultiplyAddEstimate(TVectorDouble.Create(C6), r, TVectorDouble.Create(C5)),
	TVectorDouble.MultiplyAddEstimate(
	r2,
	TVectorDouble.MultiplyAddEstimate(TVectorDouble.Create(C4), r, TVectorDouble.Create(C3)),
	TVectorDouble.MultiplyAddEstimate(TVectorDouble.Create(C2), r, TVectorDouble.Create(C1))
	)
	);
	// result = poly + (n << 52)
	return Unsafe.BitCast<TVectorUInt64, TVectorDouble>(Unsafe.BitCast<TVectorDouble, TVectorUInt64>(poly) + (n << 52));
	}
	}

Math.Exp/MathF.Exp call the platform CRT math functions. The actual result can vary among platforms.

[kristiker]

It's really weird that the implementation does not optimize for key points such as 0 and 1.

Fyi the error at 0 is +0.00000011920928955078125, at 1 it is doubled +0.0000002384185791015625

Edit: Also double clicking in my IDE leads me to VectorMath.ExpDouble, but I'm not sure if that is accurate.

[tannergooding]

It's really weird that the implementation does not optimize for key points such as 0 and 1.

It's really not. That is incredibly atypical for an algorithm to do, because it's not typically necessary.

The general issue here is that what scalar logic can do and what vectorized logic can do is a bit different. The scalar algorithm ends up being driven by a lookup table: https://github.com/amd/aocl-libm-ose/blob/master/src/optimized/expf.c. This works because there is a single value and so the cost is kept down while still being accurate.

Vectorization cannot do that trivially, however. It would require a "gather" operation which then negates most of the benefits of vectorization in typical scenarios: https://github.com/amd/aocl-libm-ose/blob/master/src/optimized/vec/vrs4_expf.c. It uses an alternative algorithm instead, which is more akin to a Taylor Series (which is also how most math algorithms are implemented). The number of polynomial steps it's doing isn't sufficient to preserve the accuracy and so it ends up with this off by one and returns bitwise 0x3F800001 when it should've returned 0x3F800000.

This off by one is typically not an issue for vectorization where perf is typically of more importance than accuracy. However, it is still undesirable if it could otherwise be trivially handled. There likely needs to be a single extra step in the polynomial handling to handle such edges, but it will require some non-trivial digging to find the right adjustment to handle all cases.