Skip to content

Vectorize TensorPrimitives.Acosh/sinh/tanh in terms of Log #98302

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Closed
wants to merge 1 commit into from
Closed

Vectorize TensorPrimitives.Acosh/sinh/tanh in terms of Log #98302

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
154 changes: 142 additions & 12 deletions ...s/System.Numerics.Tensors/src/System/Numerics/Tensors/netcore/TensorPrimitives.netcore.cs
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -15868,11 +15868,54 @@ public static Vector512<T> Invoke(Vector512<T> x)
internal readonly struct AcoshOperator<T> : IUnaryOperator<T, T>
where T : IHyperbolicFunctions<T>
{
public static bool Vectorizable => false; // TODO: Vectorize
public static bool Vectorizable => typeof(T) == typeof(float) || typeof(T) == typeof(double);

public static T Invoke(T x) => T.Acosh(x);
public static Vector128<T> Invoke(Vector128<T> x) => throw new NotSupportedException();
public static Vector256<T> Invoke(Vector256<T> x) => throw new NotSupportedException();
public static Vector512<T> Invoke(Vector512<T> x) => throw new NotSupportedException();

public static Vector128<T> Invoke(Vector128<T> x)
{
if (typeof(T) == typeof(float))
{
Vector128<float> f = x.AsSingle();
return LogOperator<float>.Invoke(f + Vector128.Sqrt(MultiplyAddEstimateOperator<float>.Invoke(f, f, Vector128.Create(-1f)))).As<float, T>();
}
else
{
Debug.Assert(typeof(T) == typeof(double));
Vector128<double> d = x.AsDouble();
return LogOperator<double>.Invoke(d + Vector128.Sqrt(MultiplyAddEstimateOperator<double>.Invoke(d, d, Vector128.Create(-1d)))).As<double, T>();
}
}

public static Vector256<T> Invoke(Vector256<T> x)
{
if (typeof(T) == typeof(float))
{
Vector256<float> f = x.AsSingle();
return LogOperator<float>.Invoke(f + Vector256.Sqrt(MultiplyAddEstimateOperator<float>.Invoke(f, f, Vector256.Create(-1f)))).As<float, T>();
}
else
{
Debug.Assert(typeof(T) == typeof(double));
Vector256<double> d = x.AsDouble();
return LogOperator<double>.Invoke(d + Vector256.Sqrt(MultiplyAddEstimateOperator<double>.Invoke(d, d, Vector256.Create(-1d)))).As<double, T>();
}
}

public static Vector512<T> Invoke(Vector512<T> x)
{
if (typeof(T) == typeof(float))
{
Vector512<float> f = x.AsSingle();
return LogOperator<float>.Invoke(f + Vector512.Sqrt(MultiplyAddEstimateOperator<float>.Invoke(f, f, Vector512.Create(-1f)))).As<float, T>();
}
else
{
Debug.Assert(typeof(T) == typeof(double));
Vector512<double> d = x.AsDouble();
return LogOperator<double>.Invoke(d + Vector512.Sqrt(MultiplyAddEstimateOperator<double>.Invoke(d, d, Vector512.Create(-1d)))).As<double, T>();
}
}
}

/// <summary>T.AcosPi(x)</summary>
Expand Down Expand Up @@ -15901,11 +15944,55 @@ public static Vector512<T> Invoke(Vector512<T> x)
internal readonly struct AsinhOperator<T> : IUnaryOperator<T, T>
where T : IHyperbolicFunctions<T>
{
public static bool Vectorizable => false; // TODO: Vectorize

public static bool Vectorizable => typeof(T) == typeof(float) || typeof(T) == typeof(double);

public static T Invoke(T x) => T.Asinh(x);
public static Vector128<T> Invoke(Vector128<T> x) => throw new NotSupportedException();
public static Vector256<T> Invoke(Vector256<T> x) => throw new NotSupportedException();
public static Vector512<T> Invoke(Vector512<T> x) => throw new NotSupportedException();

public static Vector128<T> Invoke(Vector128<T> x)
{
if (typeof(T) == typeof(float))
{
Vector128<float> f = x.AsSingle();
return LogOperator<float>.Invoke(f + Vector128.Sqrt(MultiplyAddEstimateOperator<float>.Invoke(f, f, Vector128<float>.One))).As<float, T>();
}
else
{
Debug.Assert(typeof(T) == typeof(double));
Vector128<double> d = x.AsDouble();
return LogOperator<double>.Invoke(d + Vector128.Sqrt(MultiplyAddEstimateOperator<double>.Invoke(d, d, Vector128<double>.One))).As<double, T>();
}
}

public static Vector256<T> Invoke(Vector256<T> x)
{
if (typeof(T) == typeof(float))
{
Vector256<float> f = x.AsSingle();
return LogOperator<float>.Invoke(f + Vector256.Sqrt(MultiplyAddEstimateOperator<float>.Invoke(f, f, Vector256<float>.One))).As<float, T>();
}
else
{
Debug.Assert(typeof(T) == typeof(double));
Vector256<double> d = x.AsDouble();
return LogOperator<double>.Invoke(d + Vector256.Sqrt(MultiplyAddEstimateOperator<double>.Invoke(d, d, Vector256<double>.One))).As<double, T>();
}
}

public static Vector512<T> Invoke(Vector512<T> x)
{
if (typeof(T) == typeof(float))
{
Vector512<float> f = x.AsSingle();
return LogOperator<float>.Invoke(f + Vector512.Sqrt(MultiplyAddEstimateOperator<float>.Invoke(f, f, Vector512<float>.One))).As<float, T>();
}
else
{
Debug.Assert(typeof(T) == typeof(double));
Vector512<double> d = x.AsDouble();
return LogOperator<double>.Invoke(d + Vector512.Sqrt(MultiplyAddEstimateOperator<double>.Invoke(d, d, Vector512<double>.One))).As<double, T>();
}
}
}

/// <summary>T.AsinPi(x)</summary>
Expand Down Expand Up @@ -15934,11 +16021,54 @@ public static Vector512<T> Invoke(Vector512<T> x)
internal readonly struct AtanhOperator<T> : IUnaryOperator<T, T>
where T : IHyperbolicFunctions<T>
{
public static bool Vectorizable => false; // TODO: Vectorize
public static bool Vectorizable => typeof(T) == typeof(float) || typeof(T) == typeof(double);

public static T Invoke(T x) => T.Atanh(x);
public static Vector128<T> Invoke(Vector128<T> x) => throw new NotSupportedException();
public static Vector256<T> Invoke(Vector256<T> x) => throw new NotSupportedException();
public static Vector512<T> Invoke(Vector512<T> x) => throw new NotSupportedException();

public static Vector128<T> Invoke(Vector128<T> x)
{
if (typeof(T) == typeof(float))
{
Vector128<float> f = x.AsSingle();
return (Vector128.Create(0.5f) * LogOperator<float>.Invoke((Vector128<float>.One + f) / (Vector128<float>.One - f))).As<float, T>();
Copy link
Contributor

@lilinus lilinus Feb 12, 2024 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

A more numerically stable formula for small f would be to use logp1
Vector128.Create(0.5f) * LogP1Operator<float>.Invoke(Vector128<float>.Create(2f) * f / (Vector128<float>.One - f))

Although I think logp1 is just implemented as log(x+1) at the moment...

}
else
{
Debug.Assert(typeof(T) == typeof(double));
Vector128<double> d = x.AsDouble();
return (Vector128.Create(0.5d) * LogOperator<double>.Invoke((Vector128<double>.One + d) / (Vector128<double>.One - d))).As<double, T>();
}
}

public static Vector256<T> Invoke(Vector256<T> x)
{
if (typeof(T) == typeof(float))
{
Vector256<float> f = x.AsSingle();
return (Vector256.Create(0.5f) * LogOperator<float>.Invoke((Vector256<float>.One + f) / (Vector256<float>.One - f))).As<float, T>();
}
else
{
Debug.Assert(typeof(T) == typeof(double));
Vector256<double> d = x.AsDouble();
return (Vector256.Create(0.5d) * LogOperator<double>.Invoke((Vector256<double>.One + d) / (Vector256<double>.One - d))).As<double, T>();
}
}

public static Vector512<T> Invoke(Vector512<T> x)
{
if (typeof(T) == typeof(float))
{
Vector512<float> f = x.AsSingle();
return (Vector512.Create(0.5f) * LogOperator<float>.Invoke((Vector512<float>.One + f) / (Vector512<float>.One - f))).As<float, T>();
}
else
{
Debug.Assert(typeof(T) == typeof(double));
Vector512<double> d = x.AsDouble();
return (Vector512.Create(0.5d) * LogOperator<double>.Invoke((Vector512<double>.One + d) / (Vector512<double>.One - d))).As<double, T>();
}
}
}

/// <summary>T.AtanPi(x)</summary>
Expand Down
2 changes: 1 addition & 1 deletion src/libraries/System.Numerics.Tensors/tests/TensorPrimitives.Generic.cs
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -365,7 +365,7 @@ public static IEnumerable<object[]> SpanDestinationFunctionsToTest()
yield return new object[] { new SpanDestinationDelegate(TensorPrimitives.Acosh), new Func<T, T>(T.Acosh) };
yield return new object[] { new SpanDestinationDelegate(TensorPrimitives.AcosPi), new Func<T, T>(T.AcosPi) };
yield return new object[] { new SpanDestinationDelegate(TensorPrimitives.Acos), new Func<T, T>(T.Acos) };
yield return new object[] { new SpanDestinationDelegate(TensorPrimitives.Asinh), new Func<T, T>(T.Asinh) };
yield return new object[] { new SpanDestinationDelegate(TensorPrimitives.Asinh), new Func<T, T>(T.Asinh), T.CreateTruncating(0.001) };
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This is quite a lot of variance being allowed, we're basically saying we can't get within 3 digits of accuracy for some inputs.

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

What do you suggest?

(Note that the default being used everywhere right now is 0.0001.)

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

What I'd like to understand is if this a variance being introduced more generally, or just for the special values we test like MaxValue/MinValue. If it's all or most values, then I don't think we can take the implementation. If its only for these extreme boundary values, then I think we're perfectly fine.

For a TL;DR. I think the upper boundary for a normal input is:

  • float - 2^17
  • double - 2^46

Beyond that we've already started with enough precision loss that its not going to be super meaningful no matter what we do.

For float we can exactly represent any integer up to 2^24. We can only track any fractional portion up to 2^23 (where the distance between values is 0.5). At 2^8 (256) the distance between values is 0.000030517578125, which is just enough to represent 5 significant digits of PI. For double, the same distance between values occurs at 2^37.

The ability to represent 5 digits for PI (3.1415) is a reasonably boundary at which these types of trigonometric functions can be expected to be most accurate, regardless of underlying implementation. On the more extreme end, we can represent 3 digits for PI (3.14) at 2^17 (float) or 2^46 (double). With any fewer digits we've lost enough precision that results are less meaningful, so I think this is a good maximum range for expected inputs to be accurate.

For typical inputs, we really want to expect 6-9 digits of accuracy for float and 15-17 digits of accuracy for double. With an allowance for this accuracy to taper off as the input grows in magnitude.

yield return new object[] { new SpanDestinationDelegate(TensorPrimitives.AsinPi), new Func<T, T>(T.AsinPi) };
yield return new object[] { new SpanDestinationDelegate(TensorPrimitives.Asin), new Func<T, T>(T.Asin) };
yield return new object[] { new SpanDestinationDelegate(TensorPrimitives.Atanh), new Func<T, T>(T.Atanh) };
Expand Down