Skip to content

Added SIMD path for Vector3, Quaternion and Matrix4x4 operations #99547

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 7 commits into from
Closed

Added SIMD path for Vector3, Quaternion and Matrix4x4 operations #99547

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
7 commits
Select commit Hold shift + click to select a range
38bb586
Added SIMD paths for:
martenf Mar 11, 2024
0e14da8
Fixed typo
martenf Mar 11, 2024
cfa24d9
fixed missing conversion back to Vector3
martenf Mar 11, 2024
8d71e34
Implemented feedback
martenf Apr 8, 2024
27bd42c
fixed non static field
martenf Apr 8, 2024
1f35b22
removed whitespace
martenf Apr 8, 2024
00199d2
Use Vector4.Transform
martenf Apr 8, 2024
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
29 changes: 4 additions & 25 deletions src/libraries/System.Private.CoreLib/src/System/Numerics/Matrix4x4.Impl.cs
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -166,31 +166,10 @@ public Vector3 Translation
public static Impl operator *(in Impl left, in Impl right)
{
Impl result;

// result.X = Transform(left.X, in right);
result.X = right.X * left.X.X;
result.X += right.Y * left.X.Y;
result.X += right.Z * left.X.Z;
result.X += right.W * left.X.W;

// result.Y = Transform(left.Y, in right);
result.Y = right.X * left.Y.X;
result.Y += right.Y * left.Y.Y;
result.Y += right.Z * left.Y.Z;
result.Y += right.W * left.Y.W;

// result.Z = Transform(left.Z, in right);
result.Z = right.X * left.Z.X;
result.Z += right.Y * left.Z.Y;
result.Z += right.Z * left.Z.Z;
result.Z += right.W * left.Z.W;

// result.W = Transform(left.W, in right);
result.W = right.X * left.W.X;
result.W += right.Y * left.W.Y;
result.W += right.Z * left.W.Z;
result.W += right.W * left.W.W;

result.X = Vector4.Transform(left.X, in right);
result.Y = Vector4.Transform(left.Y, in right);
result.Z = Vector4.Transform(left.Z, in right);
result.W = Vector4.Transform(left.W, in right);
return result;
}

Expand Down
67 changes: 2 additions & 65 deletions src/libraries/System.Private.CoreLib/src/System/Numerics/Quaternion.cs
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -118,40 +118,7 @@ public readonly bool IsIdentity
/// <remarks>The <see cref="op_Division" /> method defines the division operation for <see cref="Quaternion" /> objects.</remarks>
public static Quaternion operator /(Quaternion value1, Quaternion value2)
{
Quaternion ans;

float q1x = value1.X;
float q1y = value1.Y;
float q1z = value1.Z;
float q1w = value1.W;

//-------------------------------------
// Inverse part.
float ls = value2.X * value2.X + value2.Y * value2.Y +
value2.Z * value2.Z + value2.W * value2.W;
float invNorm = 1.0f / ls;

float q2x = -value2.X * invNorm;
float q2y = -value2.Y * invNorm;
float q2z = -value2.Z * invNorm;
float q2w = value2.W * invNorm;

//-------------------------------------
// Multiply part.

// cross(av, bv)
float cx = q1y * q2z - q1z * q2y;
float cy = q1z * q2x - q1x * q2z;
float cz = q1x * q2y - q1y * q2x;

float dot = q1x * q2x + q1y * q2y + q1z * q2z;

ans.X = q1x * q2w + q2x * q1w + cx;
ans.Y = q1y * q2w + q2y * q1w + cy;
ans.Z = q1z * q2w + q2z * q1w + cz;
ans.W = q1w * q2w - dot;

return ans;
return value1 * Inverse(value2);
}

/// <summary>Returns a value that indicates whether two quaternions are equal.</summary>
Expand Down Expand Up @@ -192,7 +159,6 @@ public readonly bool IsIdentity
{
var left = value1.AsVector128();
var right = value2.AsVector128();

var result = right * left.GetElementUnsafe(3);
result += (Vector128.Shuffle(right, Vector128.Create(3, 2, 1, 0)) * left.GetElementUnsafe(0)) * Vector128.Create(+1.0f, -1.0f, +1.0f, -1.0f);
result += (Vector128.Shuffle(right, Vector128.Create(2, 3, 0, 1)) * left.GetElementUnsafe(1)) * Vector128.Create(+1.0f, +1.0f, -1.0f, -1.0f);
Expand Down Expand Up @@ -289,36 +255,7 @@ public static Quaternion Add(Quaternion value1, Quaternion value2)
/// <param name="value1">The first quaternion rotation in the series.</param>
/// <param name="value2">The second quaternion rotation in the series.</param>
/// <returns>A new quaternion representing the concatenation of the <paramref name="value1" /> rotation followed by the <paramref name="value2" /> rotation.</returns>
public static Quaternion Concatenate(Quaternion value1, Quaternion value2)
{
Quaternion ans;

// Concatenate rotation is actually q2 * q1 instead of q1 * q2.
// So that's why value2 goes q1 and value1 goes q2.
float q1x = value2.X;
float q1y = value2.Y;
float q1z = value2.Z;
float q1w = value2.W;

float q2x = value1.X;
float q2y = value1.Y;
float q2z = value1.Z;
float q2w = value1.W;

// cross(av, bv)
float cx = q1y * q2z - q1z * q2y;
float cy = q1z * q2x - q1x * q2z;
float cz = q1x * q2y - q1y * q2x;

float dot = q1x * q2x + q1y * q2y + q1z * q2z;

ans.X = q1x * q2w + q2x * q1w + cx;
ans.Y = q1y * q2w + q2y * q1w + cy;
ans.Z = q1z * q2w + q2z * q1w + cz;
ans.W = q1w * q2w - dot;

return ans;
}
public static Quaternion Concatenate(Quaternion value1, Quaternion value2) => value2 * value1;

/// <summary>Returns the conjugate of a specified quaternion.</summary>
/// <param name="value">The quaternion.</param>
Expand Down
39 changes: 34 additions & 5 deletions src/libraries/System.Private.CoreLib/src/System/Numerics/Vector3.cs
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -302,11 +302,23 @@ public static Vector3 Clamp(Vector3 value1, Vector3 min, Vector3 max)
[MethodImpl(MethodImplOptions.AggressiveInlining)]
public static Vector3 Cross(Vector3 vector1, Vector3 vector2)
{
return new Vector3(
(vector1.Y * vector2.Z) - (vector1.Z * vector2.Y),
(vector1.Z * vector2.X) - (vector1.X * vector2.Z),
(vector1.X * vector2.Y) - (vector1.Y * vector2.X)
);
if (Vector128.IsHardwareAccelerated)
{
var v1 = vector1.AsVector128();
var v2 = vector2.AsVector128();
return ((Vector128.Shuffle(v1, Vector128.Create(1, 2, 0, 3)) *
Vector128.Shuffle(v2, Vector128.Create(2, 0, 1, 3))) -
(Vector128.Shuffle(v1, Vector128.Create(2, 0, 1, 3)) *
Vector128.Shuffle(v2, Vector128.Create(1, 2, 0, 3)))).AsVector3();
}
else
{
return new Vector3(
(vector1.Y * vector2.Z) - (vector1.Z * vector2.Y),
(vector1.Z * vector2.X) - (vector1.X * vector2.Z),
(vector1.X * vector2.Y) - (vector1.Y * vector2.X)
);
}
}

/// <summary>Computes the Euclidean distance between the two given points.</summary>
Expand Down Expand Up @@ -516,6 +528,23 @@ public static Vector3 Transform(Vector3 position, Matrix4x4 matrix)
[MethodImpl(MethodImplOptions.AggressiveInlining)]
public static Vector3 Transform(Vector3 value, Quaternion rotation)
{
if (Vector128.IsHardwareAccelerated)
{
var vVector = value.AsVector128();
var rVector = Unsafe.BitCast<Quaternion, Vector128<float>>(rotation);

// v + 2 * (q x (q.W * v + q x v)
return (vVector + Vector128.Create(2f) * Cross(rVector, Vector128.Shuffle(rVector, Vector128.Create(3, 3, 3, 3)) * vVector + Cross(rVector, vVector))).AsVector3();

static Vector128<float> Cross(Vector128<float> v1, Vector128<float> v2)
{
return (Vector128.Shuffle(v1, Vector128.Create(1, 2, 0, 3)) *
Vector128.Shuffle(v2, Vector128.Create(2, 0, 1, 3))) -
(Vector128.Shuffle(v1, Vector128.Create(2, 0, 1, 3)) *
Vector128.Shuffle(v2, Vector128.Create(1, 2, 0, 3)));
}
}
Comment on lines +531 to +546
Copy link
Member

Choose a reason for hiding this comment

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

I'm not immediately following where this implementation is from...

The standard naive algorithm is typically broken down to something like:

return ((Quaternion.Conjugate(rotation) * new Quaternion(value, 0.0f)) * rotation).AsVector128().AsVector3();

(noting that Quaternion.operator *(Quaternion, Quaterion) is not currently marked as aggressively inline)

Copy link
Contributor Author

Choose a reason for hiding this comment

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

I got it from a private C library I wrote a decade ago which quotes a dead link.
The closest I could find on the internet after a quick search is this forum post which gets 95% of the way there and just misses the last 'simplification for readability' step.

Copy link
Member

Choose a reason for hiding this comment

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

Can we update this to follow the "standard" algorithm I gave above.

I expect we'll get better support for things like constant folding, more consistent performance across a range of hardware, and we can accurately quote the root implementation to a known source (which avoids any potential licensing issues): https://github.com/microsoft/DirectXMath/blob/main/Inc/DirectXMathVector.inl#L10307-L10317

IIRC, glm uses a similar algorithm for it's implementation

Copy link
Contributor Author

Choose a reason for hiding this comment

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

The naive implementation is even slower then the dotnet8 implementation, even with the dotnet9 improvements to quaternion multiplication. It will also always be slower then the reduced version because it contains significantly more math operation.

I also don't see any potential licensing problem. I wrote the c# code. Btw it looks nothing like c code I based it on which I also wrote. Both version are based on math that cannot be copyrighted and that has been floating around the internet for at least decade at this point.

To reduce any variation between different archs the implementation could be changed to

return (vVector + 2 * Cross(rVector, Vector128.Create(rotation.W) * vVector + Cross(rVector, vVector))).AsVector3();

which would give the JIT a bit more leeway, but I doubt it makes any differens.

Copy link
Member

Choose a reason for hiding this comment

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

The naive implementation is even slower then the dotnet8 implementation, even with the dotnet9 improvements to quaternion multiplication

This is in part because quaternion multiplication is not being inlined currently. But, also because it's doing some "unnecessary" work, as you indicated. But that can be accounted for as well.

I also don't see any potential licensing problem. I wrote the c# code. Btw it looks nothing like c code I based it on which I also wrote. Both version are based on math that cannot be copyrighted and that has been floating around the internet for at least decade at this point.

In general we require correct attribution to exist and based on the statement given above, the code (even if you authored it) was itself based on some internet blog/forum/site/etc post that can no longer be referenced (due to a dead link). Such a post may have itself may have been based on something else which overall makes it very difficult for us to use "safely" as the chain of citations is broken.

The actual requirements for attribution, whether or not something can be copyrighted, and the like can get quite complex. When there isn't a clear chain anymore, then it requires additional effort on our part, potentially even requiring checks with legal, so an actual lawyer can make the determination on whether or not it is safe to do.

.NET is a huge open source project depending on by millions. We must do the due diligence and ensure that all the boxes are checked.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

What makes the result of Q*vQ any more correct that any of the other ways you could use to rotate a vector by the rotation saved in a quaternion? Rodrigues' formula is actually older than Hamilton's, but it's not like any of this is a universal law written in the stars. This entire conversation feels pointless so I will stop wasting everyone's time, close this pull request and just use a private library.

Copy link
Member

Choose a reason for hiding this comment

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

From a mathematical perspective they all produce the "same" result. The consideration is that they are not the same from a programmatic perspective due to additional considerations.

The intent was to push this PR in a direction where we could still improve the performance but without changing certain invariants that are extremely important to maintain. That requires more in depth consideration and documentation to show how those requirements are being met.

The intent was not to waste time or get into an argument over what is right vs wrong. The contribution was appreciated and was otherwise nearly in the right shape to be taken.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

I know and appreciate that you weren't trying to waste any time, but I just don't have the spare time to have long debates about computational vs mathematical correctness right now, so I think it's just easier if I leave it at that so that someone else can make a pull request if they run into the same bottlenecks.

Copy link
Member

@tannergooding tannergooding Jun 14, 2024 •
edited
Loading

Choose a reason for hiding this comment

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

👍, are you fine with me cherry-picking your commits into a PR so I can fix the last couple bits of feedback and get merged? That way you can still get credit for the overall contribution

Copy link
Contributor Author

Choose a reason for hiding this comment

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

Sure


float x2 = rotation.X + rotation.X;
float y2 = rotation.Y + rotation.Y;
float z2 = rotation.Z + rotation.Z;
Expand Down
Loading