From the documentation it looks like the * operator will perform a cross product and ;* will perform a dot product. I want to use a cross product to determine the vector normal to a pair of known vectors like this:
p1 = [1, 2, 4];
p2 = [2, 1, 4];
p3 = [2, 2, 4];
v1 = p3 - p1;
v2 = p2 - p1;
xProd = v1*v2;
3 points define a plane, two vectors taken from those points do the same, the cross-product should give the normal vector of that plane. The problem seems to be that the dimensions of the two matrices being multiplied (v1 and v2) do not agree -- nRows(v1) should equal nCols(v2). I can implement this manually with a simplified formula for vector cross products but does anyone know how to format this differently to work with matrix multiplication?
The inferior performance doesn't surprise... thanks for inserting the sign adjustment. I omitted the "sign checkerboard" inadvertently.
IMO this is one of those instances where "just because you can, doesn't mean you should" applies... {cool, fun, short} doesn't outweigh {understandable, maintainable, speedy}.
@brady_brady Cool! It needs a tweak if you are going to use it:
i = 1;
j = -1;
J( 3, 1, (j *= -1) * Det( (v1 || v2)[Remove( [1, 2, 3], i++ ), 0] ) );
</CODE><A href="https://mathinsight.org/cross_product_formula" target="_blank" rel="nofollow noopener noreferrer">https://mathinsight.org/cross_product_formula</A> shows how the middle element needs to be negative. Above, j will alternate sign (1, -1, 1), and i increases (1, 2, 3). The J(...) function has an implicit loop that evaluates the 3rd argument for each element.</P>
<P> </P>
<P>The straight forward corrected formula from <A href="https://community.jmp.com/t5/user/viewprofilepage/user-id/11962" target="_blank">@klk</A> is about 3 times faster and about N times easier to understand.</P>
<P> </P>
<P>Test code.</P>
<DIV class="lia-spoiler-container"><A class="lia-spoiler-link" href="#" rel="nofollow noopener noreferrer" target="_blank">View more...</A><NOSCRIPT><DIV class="lia-spoiler-content"><BR /><PRE><CODE class=" language-jsl">v1 = [0, 0, 1];
v2 = [0, 1, 0];
start = HP Time();
For( t = 1, t < 1e6, t += 1,
i = 1;
j = -1;
J( 3, 1, (j *= -1) * Det( (v1 || v2)[Remove( [1, 2, 3], i++ ), 0] ) );
);
stop = HP Time();
Show( (stop - start) / 1e6 );//9.3s
start = HP Time();
For( t = 1, t < 1e6, t += 1,
Matrix( {v1[2] * v2[3] - v1[3] * v2[2], v1[3] * v2[1] - v1[1] * v2[3], v1[1] * v2[2] - v1[2] * v2[1]} )
);
stop = HP Time();
Show( (stop - start) / 1e6 );//3.2s
For( q = 1, q < 10000, q += 1,
v1 = J( 3, 1, Random Uniform( -1, 1 ) );
v2 = J( 3, 1, Random Uniform( -1, 1 ) );
If(
All(
Round( Matrix( {v1[2] * v2[3] - v1[3] * v2[2], v1[3] * v2[1] - v1[1] * v2[3], v1[1] * v2[2] - v1[2] * v2[1]} ), 10 ) //
== //
(i = 1; j = -1 ; Round( J( 3, 1, (j *= -1) * Det( (v1 || v2)[Remove( [1, 2, 3], i++ ), 0] ) ), 10 ) ; ) //
) == 0//
, //
Throw( Char( q ) )
);
);Curiously, rounding to 11 places occasionally fails to match results. Probably nearly collinear vectors.
Well, the revived zombie caught me. Here is a way to do it via determinants as @Craige_Hales mentioned. It is short, but I don't know... it isn't very simple.
crossP = Function( {x, y, i = 1},
Return( J( 3, 1, Det( (x || y)[Remove( [1, 2, 3], i++ ), 0] ) ) )
);
i.e., using the originally-supplied data:
Sorry for reviving this zombie thread, but there is a sign error in the second element here. I believe it should be
Matrix(
{v1[2] * v2[3] - v1[3] * v2[2], v1[3] * v2[1] - v1[1] * v2[3], v1[1] * v2[2] - v1[2] * v2[1]}<BR />)Just want to leave this here for anyone else who comes along for a quick copy-paste.
Thanks cwillden, this is essentially what I have implemented. It does work for me since I am only concerned with 3D vectors, but I am supprised this isn't a built in function. I'm going to leave this open for a bit to see if anyone else knows of such a way.
I hope someone comes up with a clever way to do this. The det() function is closely related, but I don't understand how to use det(3x3 matrix) to get back more than a scalar value. I think you could use det(2x2 sub-matrix) three times, but I think the straight-forward xyzzy approach is simpler and just as fast.
I use it for the surface normals in 3D scenes, there is some JSL in https://community.jmp.com/t5/Uncharted/Custom-Visualization/ba-p/191559 .
If you're only ever dealing 3x1 (or 1x3) vectors, then I would think computational efficiency or elegance is not too important. Would this be sufficient?
cross_prod = function({v1, v2},
matrix({v1[2]*v2[3] - v1[3]*v2[2], v1[1]*v2[3] - v1[3]*v2[1], v1[1]*v2[2] - v1[2]*v2[1]});
);
cross_prod(v1, v2);
