The Cross-Product

The cross-product is a strange beast.

Do it that way

It's usually taught to you in 6^th grade where you are told its 2D definition first:

$$\boldsymbol{u} = ( u_x, u_y, 0 ) ~~~~~~~~~~ \boldsymbol{v} = ( v_x, v_y, 0 ) ~~~~~~~~~~ \boldsymbol{w}= \boldsymbol{u} \times \boldsymbol{v} = (0, 0, u_x v_y - u_y v_x )$$

You are also told that $\Vert\boldsymbol{w}\Vert = \Vert\boldsymbol{u}\Vert \Vert\boldsymbol{v}\Vert \sin(\alpha)$ and that the resulting vector $\boldsymbol{w}$ is orthogonal to its base constituents $\boldsymbol{u}$ and $\boldsymbol{v}$.

Then everyone in the classroom starts doing some strange gestures with their fingers to understand which way the vector should point to, something about the "right hand rule", clockwise rotation, counter-clockwise rotation, etc.

And then you start using it everywhere because as soon as you need something orthogonal to another, the cross-product immediately pops into your mind.

But why?

Why introduce such a tool? Where does it come from? Is it really necessary? Isn't the $u_x v_y - u_y v_x$ part a bit arbitrary?

The history section of the Wikipedia page says the dot and cross product were apparently introduced by Lagrange in order to study the tetrahedron in three dimensions. Obviously, you start imagining the volume of a tetrahedron, some parallelograms come to mind, the wedge product, areas, volumes, tensors, group theory and whatnots, etc. (By the way, I encourage you to read the very interesting article ¹ by Nathan Reed about the Grassman Algebra that makes a nice echo to what I am writing here) (and he published it pretty much the same day I finished writing this )

But you're no closer to an idea than when you started... I won't pretend I will provide a brilliant explanation of "why the cross product" either, just adding oil to the thought process here.

Factorizing the Tools

As a programmer, there is nothing I hate more than duplicate code... And I claim that the cross product is nothing more than a duplicate tool: the essential original tool being the dot product!

Why bother learning a new notation? Only because it's found in many places? Why not calling it the "orthogonal dot product" at the very least? Because it would seem that most people forgot the cross product is really the dot product but rewrote in another way. Let me explain.

When we write:

$$u_x v_y - u_y v_x$$

Aren't we actually doing the dot product?

$$u_x (v_y) + u_y (-v_x) = \boldsymbol{u} \cdot \boldsymbol{n}$$

With $\boldsymbol{n}=( v_y, -v_x, 0 )$

I chose the symbol $\boldsymbol{n}$ on purpose, to show that this vector is actually the normal vector to our original $\boldsymbol{v}$ vector:

We know that $\boldsymbol{u} \cdot \boldsymbol{n} = \Vert\boldsymbol{u}\Vert \Vert\boldsymbol{n}\Vert \cos(\alpha')$ and, from the figure we can see that $\alpha' = \alpha-\frac{\pi}{2}$ and thus, indeed, $\cos(\alpha') = \cos(\alpha-\frac{\pi}{2}) = \sin(\alpha)$ and we find back the original expression:

$$\Vert\boldsymbol{w}\Vert = u_x v_y - u_y v_x = \boldsymbol{u} \cdot \boldsymbol{n} = \Vert\boldsymbol{u}\Vert \Vert\boldsymbol{n}\Vert \cos(\alpha') = \Vert\boldsymbol{u}\Vert \Vert\boldsymbol{v}\Vert \sin(\alpha)$$

So what is it then?

We just saw that the cross-product is nothing more than a mere dot product with a vector $\boldsymbol{n}$ normal to the original vector $\boldsymbol{v}$. As the dot product measures the alignment of 2 vectors (i.e. the alignment being maximum when the dot product is 1 in absolute value), and since we're measuring the alignment of a vector $\boldsymbol{n}$ that is orthogonal to $\boldsymbol{v}$, then we're actually measuring the orthogonality of the two vectors $\boldsymbol{u}$ and $\boldsymbol{v}$.

Claim: The cross product is just a twist of the dot product to measure "how much two vectors are orthogonal from each other".

Three Dimensions

In 3 dimensions, the arbitrary nature of the operations you are told to remember is at its maximum!

Indeed, if you have:

$$\boldsymbol{u} = ( u_x, u_y, u_z ) ~~~~~~~~~~ \boldsymbol{v} = ( v_x, v_y, v_z )$$

Then, get ready... The cross product $\boldsymbol{u} \times \boldsymbol{v}$ is given by:

$$~~~~~~~~~~ \boldsymbol{w}= \boldsymbol{u} \times \boldsymbol{v} =\begin{cases} u_y v_z - u_z v_y\\\\ u_z v_x - u_x v_z\\\\ u_x v_y - u_y v_x\\\\ \end{cases}$$

You easily see the patterns, the rotations and permutations that are necessary to retrieve the result. You memorize it and you retrieve it each time you need it.

But one more time, I claim these are just 3 little dot products in disguise:

$$~~~~~~~~~~ \boldsymbol{w}= \boldsymbol{u} \times \boldsymbol{v} =\begin{cases} u_y (v_z) + u_z (-v_y) = u_y n_y + u_z n_z = \boldsymbol{u_{yz}} \cdot \boldsymbol{n_{yz}}\\\\ u_z (v_x) + u_x (-v_z) = u_z n_z + u_x n_x = \boldsymbol{u_{zx}} \cdot \boldsymbol{n_{zx}}\\\\ u_x (v_y) + u_y (-v_x) = u_x n_x + u_y n_y = \boldsymbol{u_{xy}} \cdot \boldsymbol{n_{xy}}\\\\ \end{cases}$$

With $\boldsymbol{u_{yz}} = \begin{cases}u_y\\\\u_z\end{cases} ~~~~~$ $\boldsymbol{u_{zx}} = \begin{cases}u_z\\\\u_x\end{cases} ~~~~~$ $\boldsymbol{u_{xy}} = \begin{cases}u_x\\\\u_y\end{cases}$

And $\boldsymbol{n_{yz}} = \begin{cases}v_z\\\\-v_y\end{cases} ~~~~~$ $\boldsymbol{n_{zx}} = \begin{cases}v_x\\\\-v_z\end{cases} ~~~~~$ $\boldsymbol{n_{xy}} = \begin{cases}v_y\\\\-v_x\end{cases}$

Each little dot product is actually measuring how much 2D sub-vectors are orthogonal from each other along each axis, each in their own 2D sub-space XY, YZ, ZX (once again, Nathan Reed uses better words than I am, in his article about Grassman algebra ¹).

• The X coordinate equals how much the 2D vectors are orthogonal from each other in the complementary 2D sub-space YZ
• The Y coordinate equals how much the 2D vectors are orthogonal from each other in the complementary 2D sub-space ZX
• The Z coordinate equals how much the 2D vectors are orthogonal from each other in the complementary 2D sub-space XY

Orthogonality

Building a new 3D vector this way is actually yielding a vector that is orthogonal to the other two base vectors (unless the 2 base vectors are colinear). Obviously, there are only 6 possible permutations of different triplets of vector components for arbitrary vectors $\boldsymbol{u}$, $\boldsymbol{v}$ and $\boldsymbol{w}$:

$~~~~u_x v_y w_z$

$~~~~u_y v_z w_x$

$~~~~u_z v_x w_y$

$~~~~u_x v_z w_y$

$~~~~u_y v_x w_z$

$~~~~u_z v_y w_x$

And they all find their expression as positive and negative quantities when writing the dot product of $\boldsymbol{w}$ with the resulting cross product $\boldsymbol{u} \times \boldsymbol{v}$:

$\boldsymbol{w} \cdot (\boldsymbol{u} \times \boldsymbol{v}) = u_y v_z w_x - u_z v_y w_z + u_z v_x w_y - u_x v_z w_y + u_x v_y w_z - u_y v_x w_z$

Whenever $\boldsymbol{w}$ turns out to be $\boldsymbol{w} = \boldsymbol{u}$ or $\boldsymbol{w} = \boldsymbol{v}$, the quantities find balancing positive and negative quantities that eventually vanish, thus proving that the vector resulting from the cross product construction is indeed orthogonal to its base construction vectors $\boldsymbol{u}$ and $\boldsymbol{v}$.

Higher Dimensions

So we gathered interesting information so far:

1. The cross product is actually a collection of $N$ dot products, $N$ being the amount of dimensions of our space
2. The $N$ dot products from point 1. are performed between vectors of dimension $N-1$, each of these vectors $\boldsymbol{u_i}$, $\boldsymbol{v_i}$ and $\boldsymbol{n_i}$ being expressed in the sub-space $i \in [1,N]$ excluding the dimension $i$ itself (e.g. in 3D, $\boldsymbol{n_1}=\boldsymbol{n_{yz}}$ is a combination of components from axes $Y$ and $Z$, excluding axis $X$)
3. Each dot product is performed between the original vector $\boldsymbol{u}$ and a "normal vector" $\boldsymbol{n}$ that is derived from the original vector $\boldsymbol{v}$ by permutations and optional sign changes of the components of $\boldsymbol{v}$ so that $\boldsymbol{v} \cdot \boldsymbol{n} = 0$
4. The combination of permutations and sign changes from point 3. stays the same for all $N$ sub-space normal vectors
5. Dotting the vector $\boldsymbol{w} = \boldsymbol{u} \times \boldsymbol{v}$ with either one of the base vectors $\boldsymbol{u}$ and $\boldsymbol{v}$ must yield 0, guaranteeing the orthogonality of the result
6. The algebra must also be anti-commutative so that $\boldsymbol{u} \times \boldsymbol{v} = -\left(\boldsymbol{v} \times \boldsymbol{u}\right)$ (but is this a consequence or a condition? Couldn't we accept a commutative cross-product?)

Assuming these rules apply for any dimension, let's try and build a 4D cross-product!

Note

For dimensions > 3, we will drop the scalar notation of the individual components of vector $\boldsymbol{u_{xyz}} = \left( u_x, u_y, u_z \right)$ in favor of the notation $u_0$, $u_1$, $u_2 \dots u_n$

And the vector notation $\boldsymbol{u_{xyz}}$ will be rewritten $\boldsymbol{u_{012}}$

Odd Dimensions

From point 3., we know that $\boldsymbol{n} \cdot \boldsymbol{v} = 0$.

In order for this to be possible, we use point 2. to notice that we must have sub-vectors of even dimension for terms coupling and canceling to occur!

Indeed, for example in 4 dimensions, we would build our cross product as:

$$\boldsymbol{w} = \begin{cases} \boldsymbol{u_{123}} \cdot \boldsymbol{n_{123}}\\\\ \boldsymbol{u_{230}} \cdot \boldsymbol{n_{230}}\\\\ \boldsymbol{u_{301}} \cdot \boldsymbol{n_{301}}\\\\ \boldsymbol{u_{012}} \cdot \boldsymbol{n_{012}} \end{cases}$$

And for any of these 4 cases, we could write something like this expression:

$$\boldsymbol{v_{012}} \cdot \boldsymbol{n_{012}} = 0$$

But whatever the way we may construct the vector $\boldsymbol{n_{012}}$, we will find no combination so that:

$$v_0 \cdot ? + v_1 \cdot ? + v_2 \cdot ? = 0$$

Simply because we need an even amount of terms for that to be possible (i.e. so that an even amount of + and - terms cancel each other).

Note

This leads to the very interesting remark that cross products only exist for odd dimensions 3, 5, 7 and so on. (you may object the 2D case but we're actually doing a 3D cross product with planar 2D vectors really, the true 2D cross-product doesn't exist) (interestingly, the 1D cross product is simply the product of two scalars )

And indeed, for odd dimensions you get sub-space vectors of even dimensions and you can then have terms coupling that cancel each other.

Permutations

For example, for 5 dimensions you can write the cross product as:

$$\boldsymbol{w} = \begin{cases} \boldsymbol{u_{1234}} \cdot \boldsymbol{n_{1234}}\\\\ \boldsymbol{u_{2340}} \cdot \boldsymbol{n_{2340}}\\\\ \boldsymbol{u_{3401}} \cdot \boldsymbol{n_{3401}}\\\\ \boldsymbol{u_{4012}} \cdot \boldsymbol{n_{4012}}\\\\ \boldsymbol{u_{0123}} \cdot \boldsymbol{n_{0123}} \end{cases}$$

Selecting the simple case $\boldsymbol{u_{0123}} \cdot \boldsymbol{n_{0123}}$, we first notice from point 2. and 3. that $\boldsymbol{n_{0123}}$ needs to be a derangement of the indices $0123$ to ensure that no component index from $\boldsymbol{u}$ gets multiplied by the same component index from $\boldsymbol{n} = derangement\left(\boldsymbol{v}\right)$ (i.e. $u_0$ can't get multiplied by $v_0$, $u_1$ can't get multiplied by $v_1$ and so on).

The amount of possible derangements in 4 dimensions, noted $!4$, is 9:

$$\begin{align} (1032) (1302) (1203) \\\\ (2031) (2301) (2310) \\\\ (3012) (3201) (3210) \\\\ \end{align}$$

We can further reduce the amount of possibilities by remembering that $\boldsymbol{v_{0123}} \cdot \boldsymbol{n_{0123}} = 0$, implying the construction of the normal vector must satisfy:

$$\boldsymbol{v_{0123}} \cdot \boldsymbol{n_{0123}} = v_0 \cdot ? + v_1 \cdot ? + v_2 \cdot ? + v_3 \cdot ? = 0$$

And this can only happen in the following "mirror" cases:

$$\boldsymbol{n_{0123}} = \begin{cases} \pm v_1 \\\\ \mp v_0 \\\\ \pm v_3 \\\\ \mp v_2 \end{cases} \quad \text{or} \quad \boldsymbol{n_{0123}} = \begin{cases} \pm v_2 \\\\ \pm v_3 \\\\ \mp v_0 \\\\ \mp v_1 \end{cases} \quad \text{or} \quad \boldsymbol{n_{0123}} = \begin{cases} \pm v_3 \\\\ \pm v_2 \\\\ \mp v_1 \\\\ \mp v_0 \end{cases}$$

This is because if you decide to couple $v_0$ with $v_2$, then you necessarily need to couple $v_2$ with $v_0$ in mirror, with a negative sign to make both terms cancel each other.
Moreover, coupling $v_0$ with $v_2$ also forces the mirror coupling of $v_1$ with $v_3$...

As for the signs coupling, with 4D vectors you only have the following mirror choices, in the case of permutation $\boldsymbol{n_{0123}} = \boldsymbol{v_{1032}}$:

$$+-+- \\\\ -+-+ \\\\ +--+ \\\\ -++-$$

It's now quite easy to notice that for any odd dimension $N=2k+1$, the number of possible combinations for constructing a normal vector is equal to $2^k \cdot (2k-1)!!$.

→ The term $2^k$ is the amount of binary sign permutations (so $2^x$) for each pair (so $\frac{N-1}{2} = \frac{2k}{2} = k$ as exponent)

→ The term $(2k-1)!!$ is due to the number of possible swaps of even sub-space dimensions (i.e. $2k$) without including your own dimension (i.e. in 3D, $x$ can use $y$ and $z$ but not itself, so $-1$).
Notice that it involves the double factorial $n!! = n(n-2)(n-4) \cdots 5\cdot3\cdot1$ because we're dealing with recursion here: if it takes only 3 valid derangements in 5D, then in 7D it will take 5 times the amount of valid derangements so $5 \times 3$, and in 9D it will take $7 \times 5 \times 3$ valid derangements, etc.

Finding the Right One

So all in all, for 5D vectors, we saw that we can choose among $(2\cdot2-1) \cdot 2^2 = 12$ possible ways to rewrite any sub-space 4D vector like $\boldsymbol{v_{0123}}$.

All we need now is to enforce point 5. to guarantee the orthogonality of the constructed vector.

So we write the final required expression:

$$\boldsymbol{u} \cdot (\boldsymbol{u} \times \boldsymbol{v}) = \boldsymbol{v} \cdot (\boldsymbol{u} \times \boldsymbol{v}) = \boldsymbol{0}$$

This expression can be developped into a pretty nightmarish result, for example with the 5D vectors $\boldsymbol{u_{01234}}$ and $\boldsymbol{v_{01234}}$ and using the combination $\boldsymbol{n_{0123}} = (v_1,-v_0,v_3,-v_2)$ for each sub-space vector we get:

$$\begin{align} \boldsymbol{u} \cdot (\boldsymbol{u} \times \boldsymbol{v}) &= u_0.u_1.v_2 - u_0.u_2.v_1 + u_0.u_3.v_2 + u_0.u_2.v_3 \\\\ &+ u_1.u_2.v_3 - u_1.u_3.v_2 + u_1.u_4.v_0 + u_1.u_0.v_4 \\\\ &+ u_2.u_3.v_4 - u_2.u_4.v_3 + u_2.u_0.v_1 + u_2.u_1.v_0 \\\\ &+ u_3.u_4.v_0 - u_3.u_0.v_4 + u_3.u_1.v_2 + u_3.u_2.v_1 \\\\ &+ u_4.u_0.v_1 - u_4.u_1.v_0 + u_4.u_2.v_3 + u_4.u_3.v_2 \\\\ \end{align}$$

Focusing on the single first term $u_0.u_1.v_2$ we easily see that it cannot be found again anywhere else in the expression, either as a positive or negative quantity, and is thus not cancelled.
The combination $\boldsymbol{n_{0123}} = (v_1,-v_0,v_3,-v_2)$ is thus not a good candidate for creating the 5D cross product operator.

I wrote a simple program that tries all possible combinations and accumulates $uuv$ triplets as positive and negative quantities. If all the resulting factors for each triplet are 0 then the combination is marked as valid for the cross product.

Code for finding a 5D Cross Product (C#)

void    TestCross5D() {
    int[][] derangements = new int[3][] {
        new int[4] { 1, 0, 3, 2 },      // Couple {x,y} and {z,w}
        new int[4] { 2, 3, 0, 1 },      // Couple {x,z} and {y,w}
        new int[4] { 3, 2, 1, 0 },      // Couple {x,w} and {y,z}
    };
    int[][] couples = new int[3][] {
        new int[4] { 1, -1, 2, -2 },    // Couple {x,y} and {z,w}
        new int[4] { 1, 2, -1, -2 },    // Couple {x,z} and {y,w}
        new int[4] { 1, 2, -2, -1 },    // Couple {x,w} and {y,z}
    };

    int[,]  permutations = new int[5,4] {
        { 1, 2, 3, 4 },                 // Row x uses yzwt
        { 2, 3, 4, 0 },                 // Row y uses zwtx
        { 3, 4, 0, 1 },                 // Row z uses wtxy
        { 4, 0, 1, 2 },                 // Row w uses txyz
        { 0, 1, 2, 3 },                 // Row t uses xyzw
    };

    int[]   signs = new int[4];

    List< Tuple<int,int> >  solutionsU = new List<Tuple<int, int>>();
    List< Tuple<int,int> >  solutionsV = new List<Tuple<int, int>>();
    List< Tuple<int,int> >  solutions = new List<Tuple<int, int>>();

    // Try all possible derangements
    for ( int i=0; i < 3; i++ ) {
        int[]   derangement = derangements[i];
        int[]   couple = couples[i];

        // Try all possible sign combinations
        for ( int s=0; s < 4; s++ ) {
            int s0 = (s & 1) != 0 ? 1 : -1;
            int s1 = (s & 2) != 0 ? 1 : -1;
            signs[0] = (couple[0] < 0 ? -1 : 1) * (Math.Abs(couple[0]) == 1 ? s0 : s1);
            signs[1] = (couple[1] < 0 ? -1 : 1) * (Math.Abs(couple[1]) == 1 ? s0 : s1);
            signs[2] = (couple[2] < 0 ? -1 : 1) * (Math.Abs(couple[2]) == 1 ? s0 : s1);
            signs[3] = (couple[3] < 0 ? -1 : 1) * (Math.Abs(couple[3]) == 1 ? s0 : s1);

            // Estimate each line of w.(u x v)
            Array.Clear( m_sumTripletsU, 0, 5*5*5 );
            Array.Clear( m_sumTripletsV, 0, 5*5*5 );
            for ( int d=0; d < 5; d++ ) {
                int index_w_0 = d;
                int index_u_0 = permutations[d,0];
                int index_v_0 = permutations[d,derangement[0]];
                int sign_0 = signs[0];
                AddTripletU( index_w_0, index_u_0, index_v_0, sign_0 );
                AddTripletV( index_w_0, index_u_0, index_v_0, sign_0 );

                int index_w_1 = d;
                int index_u_1 = permutations[d,1];
                int index_v_1 = permutations[d,derangement[1]];
                int sign_1 = signs[1];
                AddTripletU( index_w_1, index_u_1, index_v_1, sign_1 );
                AddTripletV( index_w_1, index_u_1, index_v_1, sign_1 );

                int index_w_2 = d;
                int index_u_2 = permutations[d,2];
                int index_v_2 = permutations[d,derangement[2]];
                int sign_2 = signs[2];
                AddTripletU( index_w_2, index_u_2, index_v_2, sign_2 );
                AddTripletV( index_w_2, index_u_2, index_v_2, sign_2 );

                int index_w_3 = d;
                int index_u_3 = permutations[d,3];
                int index_v_3 = permutations[d,derangement[3]];
                int sign_3 = signs[3];
                AddTripletU( index_w_3, index_u_3, index_v_3, sign_3 );
                AddTripletV( index_w_3, index_u_3, index_v_3, sign_3 );
            }

            // Check if the sum only contains zeroes
            bool    allZeroesU = true;
            bool    allZeroesV = true;
            for ( int j=0; j < 5*5*5; j++ ) {
                int z = j;
                int x = z / (5*5);
                z -= 5*5 * x;
                int y = z / 5;
                z -= 5* y;

                int sumU = m_sumTripletsU[x,y,z];
                int sumV = m_sumTripletsV[x,y,z];
                if ( sumU != 0 ) {
                    allZeroesU = false;
                }
                if ( sumV != 0 ) {
                    allZeroesV = false;
                }
            }
            if ( allZeroesU && allZeroesV ) {
                solutions.Add( new Tuple<int,int>( i, s ) );
            } else if ( allZeroesU ) {
                solutionsU.Add( new Tuple<int,int>( i, s ) );
            } else if ( allZeroesV ) {
                solutionsV.Add( new Tuple<int,int>( i, s ) );
            }
        }
    }

    if ( solutions.Count == 0 )
        throw new Exception( "No solution!" );
}

int[,,] m_sumTripletsU = new int[5,5,5];
void    AddTripletU( int _u0, int _u1, int _v, int _sign ) {
    if ( _u0 == _u1 || _u1 == _v ) throw new Exception(  "Can't have identical indices!" );
    m_sumTripletsU[_u0, _u1, _v] += _sign;
    m_sumTripletsU[_u1, _u0, _v] += _sign;
}
int[,,] m_sumTripletsV = new int[5,5,5];
void    AddTripletV( int _v0, int _u, int _v1, int _sign ) {
    if ( _v0 == _u || _u == _v1 ) throw new Exception(  "Can't have identical indices!" );
    m_sumTripletsV[_v0, _u, _v1] += _sign;
    m_sumTripletsV[_v1, _u, _v0] += _sign;
}

Unfortunately, this code returns no valid solution so, unless I made a mistake, it seems that there is no 4D combination of component and sign swaps that creates an algebra that induces the cross-product operator.

Let's try 7 dimensions then

According to wikipedia, the next space where a cross product is available is in 7D so I'm expecting to find solutions with this new code that checks the $5!! \cdot 2^3 = 120$ possible indices permutations for 6D sub-space vectors:

Code for finding a 7D Cross Product (C#)

void    TestCross7D() {
    // Manual creation of the 5*3 possible valid derangements
    const int   x=0, y=1, z=2, w=3, s=4, t=5;
    int[][] derangements = new int[5*3][] {
        new int[6] { y, x, w, z, t, s },
        new int[6] { y, x, s, t, z, w },
        new int[6] { y, x, t, s, w, z },

        new int[6] { z, w, x, y, t, s },
        new int[6] { z, s, x, t, y, w },
        new int[6] { z, t, x, s, w, y },    // SOLUTION #1!

        new int[6] { w, z, y, x, t, s },
        new int[6] { w, s, t, x, y, z },
        new int[6] { w, t, s, x, z, y },

        new int[6] { s, z, y, t, x, w },    // SOLUTION #2!
        new int[6] { s, w, t, y, x, z },
        new int[6] { s, t, w, z, x, y },

        new int[6] { t, z, y, s, w, x },
        new int[6] { t, w, s, y, z, x },
        new int[6] { t, s, w, z, y, x },
    };
    int[][] coupless = new int[5*3][];
    for ( int i=0; i < derangements.Length; i++ ) {     // Automate couples creation to avoid mistakes...
        int[]   couples = new int[6];
        coupless[i] = couples;
        int coupleIndex = 1;
        for ( int j=0; j < 6; j++ ) {
            if ( couples[j] == 0 ) {
                couples[j] = coupleIndex;
                couples[derangements[i][j]] = -coupleIndex;
                coupleIndex++;
            }
        }
    }

    int[,]  permutations = new int[7,6] {
        { 1, 2, 3, 4, 5, 6 },           // Row i uses jklmno
        { 2, 3, 4, 5, 6, 0 },           // Row j uses klmnoi
        { 3, 4, 5, 6, 0, 1 },           // Row k uses lmnoij
        { 4, 5, 6, 0, 1, 2 },           // Row l uses mnoijk
        { 5, 6, 0, 1, 2, 3 },           // Row m uses noijkl
        { 6, 0, 1, 2, 3, 4 },           // Row n uses oijklm
        { 0, 1, 2, 3, 4, 5 },           // Row o uses ijklmn
    };

    int[]   signs = new int[6];

    List< Tuple<int,int> >  solutionsU = new List<Tuple<int, int>>();
    List< Tuple<int,int> >  solutionsV = new List<Tuple<int, int>>();
    List< Tuple<int,int> >  solutions = new List<Tuple<int, int>>();

    // Try all possible derangements
    for ( int i=0; i < derangements.Length;i++ ) {
        int[]   derangement = derangements[i];
        int[]   couples = coupless[i];

        // Try all possible sign combinations
        for ( int sign=0; sign < 8; sign++ ) {
            int[]   signBits = new int[] {
                (sign & 1) != 0 ? -1 : 1,
                (sign & 2) != 0 ? -1 : 1,
                (sign & 4) != 0 ? -1 : 1
            };
            signs[0] = (couples[0] < 0 ? -1 : 1) * signBits[Math.Abs(couples[0]) - 1];
            signs[1] = (couples[1] < 0 ? -1 : 1) * signBits[Math.Abs(couples[1]) - 1];
            signs[2] = (couples[2] < 0 ? -1 : 1) * signBits[Math.Abs(couples[2]) - 1];
            signs[3] = (couples[3] < 0 ? -1 : 1) * signBits[Math.Abs(couples[3]) - 1];
            signs[4] = (couples[4] < 0 ? -1 : 1) * signBits[Math.Abs(couples[4]) - 1];
            signs[5] = (couples[5] < 0 ? -1 : 1) * signBits[Math.Abs(couples[5]) - 1];

            // Estimate each line of w.(u x v)
            Array.Clear( m_sumTripletsU, 0, 7*7*7 );
            Array.Clear( m_sumTripletsV, 0, 7*7*7 );
            for ( int d=0; d < 7; d++ ) {
                int index_w_0 = d;
                int index_u_0 = permutations[d,0];
                int index_v_0 = permutations[d,derangement[0]];
                int sign_0 = signs[0];
                AddTripletU( index_w_0, index_u_0, index_v_0, sign_0 );
                AddTripletV( index_w_0, index_u_0, index_v_0, sign_0 );

                int index_w_1 = d;
                int index_u_1 = permutations[d,1];
                int index_v_1 = permutations[d,derangement[1]];
                int sign_1 = signs[1];
                AddTripletU( index_w_1, index_u_1, index_v_1, sign_1 );
                AddTripletV( index_w_1, index_u_1, index_v_1, sign_1 );

                int index_w_2 = d;
                int index_u_2 = permutations[d,2];
                int index_v_2 = permutations[d,derangement[2]];
                int sign_2 = signs[2];
                AddTripletU( index_w_2, index_u_2, index_v_2, sign_2 );
                AddTripletV( index_w_2, index_u_2, index_v_2, sign_2 );

                int index_w_3 = d;
                int index_u_3 = permutations[d,3];
                int index_v_3 = permutations[d,derangement[3]];
                int sign_3 = signs[3];
                AddTripletU( index_w_3, index_u_3, index_v_3, sign_3 );
                AddTripletV( index_w_3, index_u_3, index_v_3, sign_3 );

                int index_w_4 = d;
                int index_u_4 = permutations[d,4];
                int index_v_4 = permutations[d,derangement[4]];
                int sign_4 = signs[4];
                AddTripletU( index_w_4, index_u_4, index_v_4, sign_4 );
                AddTripletV( index_w_4, index_u_4, index_v_4, sign_4 );

                int index_w_5 = d;
                int index_u_5 = permutations[d,5];
                int index_v_5 = permutations[d,derangement[5]];
                int sign_5 = signs[5];
                AddTripletU( index_w_5, index_u_5, index_v_5, sign_5 );
                AddTripletV( index_w_5, index_u_5, index_v_5, sign_5 );
            }

            // Check if the sum only contains zeroes
            bool    allZeroesU = true;
            bool    allZeroesV = true;
            for ( int j=0; j < 7*7*7; j++ ) {
                int Z = j;
                int X = Z / (7*7);
                Z -= 7*7 * X;
                int Y = Z / 7;
                Z -= 7 * Y;

                int sumU = m_sumTripletsU[X,Y,Z];
                int sumV = m_sumTripletsV[X,Y,Z];
                if ( sumU != 0 ) {
                    allZeroesU = false;
                }
                if ( sumV != 0 ) {
                    allZeroesV = false;
                }
            }
            if ( allZeroesU && allZeroesV ) {
                solutions.Add( new Tuple<int,int>( i, sign ) );
            } else if ( allZeroesU ) {
                solutionsU.Add( new Tuple<int,int>( i, sign ) );
            } else if ( allZeroesV ) {
                solutionsV.Add( new Tuple<int,int>( i, sign ) );
            }
        }
    }

    if ( solutions.Count == 0 )
        throw new Exception( "No solution!" );
}

int[,,] m_sumTripletsU = new int[7,7,7];
void    AddTripletU( int _u0, int _u1, int _v, int _sign ) {
    if ( _u0 == _u1 || _u1 == _v ) throw new Exception(  "Can't have identical indices!" );
    m_sumTripletsU[_u0, _u1, _v] += _sign;
    m_sumTripletsU[_u1, _u0, _v] += _sign;
}
int[,,] m_sumTripletsV = new int[7,7,7];
void    AddTripletV( int _v0, int _u, int _v1, int _sign ) {
    if ( _v0 == _u || _u == _v1 ) throw new Exception(  "Can't have identical indices!" );
    m_sumTripletsV[_v0, _u, _v1] += _sign;
    m_sumTripletsV[_v1, _u, _v0] += _sign;
}

We indeed find 2 solutions (actually 4 but I ignored the signs mirroring) with the following 2 derangements:

• $\boldsymbol{n_{012345}} = ( +v_2, +v_5, -v_0, +v_4, -v_3, -v_1 )$

• $\boldsymbol{n_{012345}} = ( +v_4, +v_2, -v_1, +v_5, -v_0, -v_3 )$

Fully expanding the 1^st solution we get:

$$\begin{align} \boldsymbol{u} \cdot (\boldsymbol{u} \times \boldsymbol{v}) &= u_0.u_1.v_3 + u_0.u_2.v_6 - u_0.u_3.v_1 + u_0.u_4.v_5 - u_0.u_5.v_4 - u_0.u_6.v_2 \\\\ &+ u_1.u_2.v_4 + u_1.u_3.v_0 - u_1.u_4.v_2 + u_1.u_5.v_6 - u_1.u_6.v_5 - u_1.u_0.v_3 \\\\ &+ u_2.u_3.v_5 + u_2.u_4.v_1 - u_2.u_5.v_3 + u_2.u_6.v_0 - u_2.u_0.v_6 - u_2.u_1.v_4 \\\\ &+ u_3.u_4.v_6 + u_3.u_5.v_2 - u_3.u_6.v_4 + u_3.u_0.v_1 - u_3.u_1.v_0 - u_3.u_2.v_5 \\\\ &+ u_4.u_5.v_0 + u_4.u_6.v_3 - u_4.u_0.v_5 + u_4.u_1.v_2 - u_4.u_2.v_1 - u_4.u_3.v_6 \\\\ &+ u_5.u_6.v_1 + u_5.u_0.v_4 - u_5.u_1.v_6 + u_5.u_2.v_3 - u_5.u_3.v_2 - u_5.u_4.v_0 \\\\ &+ u_6.u_0.v_2 + u_6.u_1.v_5 - u_6.u_2.v_0 + u_6.u_3.v_4 - u_6.u_4.v_3 - u_6.u_5.v_1 \\\\ \end{align}$$

We can verify the terms indeed cancel each other...

Only special odd dimensions

We can generalize the previous programs for any odd dimension:

Generic Code for Finding a N-Dimensional Cross Product (C#)

void    TestCrossGeneric() {

    List< Tuple< int, int > >[] allSolutions = new List<Tuple<int, int>>[(21-3) / 2];
    int[]                       totalValidDerangementsCounts = new int[(21-3) / 2];

    // Try all possible odd dimensions
    for ( int N=3; N < 21; N+=2 ) {

        // Recursively generate all possible derangements
        int[][] derangements = null;
        int[][] coupless = null;
        int[]   originalOrder = new int[N-1];
        for ( int i=0; i < N-1; i++ )
            originalOrder[i] = i;
        RecurseGenerateDerangements( originalOrder, out derangements, out coupless );

        totalValidDerangementsCounts[(N-3) / 2] = derangements.Length;

        int[,]  permutations = new int[N,N-1];
        for ( int d=0; d < N; d++ ) {
            for ( int x=0; x < N-1; x++ ) {
                permutations[d,x] = (d+1 + x) % N;
            }
        }

        int     signBitsCount = (N-1) / 2;
        int     totalSignsCombinations = 1 << signBitsCount;
        int[]   signs = new int[N-1];

        List< Tuple<int,int> >  solutionsU = new List<Tuple<int, int>>();
        List< Tuple<int,int> >  solutionsV = new List<Tuple<int, int>>();
        List< Tuple<int,int> >  solutions = new List<Tuple<int, int>>();

        // Try all possible derangements
        for ( int derangementIndex=0; derangementIndex < derangements.Length; derangementIndex++ ) {
            int[]   derangement = derangements[derangementIndex];
            int[]   couples = coupless[derangementIndex];

            // Try all possible sign combinations
            for ( int signCombination=0; signCombination < totalSignsCombinations; signCombination++ ) {

                // Build signs for each component
                for ( int componentIndex=0; componentIndex < N-1; componentIndex++ ) {
                    int bitIndex = Math.Abs( couples[componentIndex] ) - 1;
                    int signValue = (signCombination & (1 << bitIndex)) != 0 ? -1 : 1;
                        signValue *= couples[componentIndex] > 0 ? 1 : -1;
                    signs[componentIndex] = signValue;
                }

                // Estimate each line of w.(u x v)
                int[,,] sumTripletsU = new int[N,N,N];
                int[,,] sumTripletsV = new int[N,N,N];
                for ( int d=0; d < N; d++ ) {
                    for ( int x=0; x < N-1; x++ ) {
                        int index_w = d;
                        int index_u = permutations[d,x];
                        int index_v = permutations[d,derangement[x]];
                        int sign = signs[x];
                        AddTripletU( sumTripletsU, index_w, index_u, index_v, sign );
                        AddTripletV( sumTripletsV, index_w, index_u, index_v, sign );
                    }
                }

                // Check if the sum only contains zeroes
                bool    allZeroesU = true;
                bool    allZeroesV = true;
                int     totalTriplets = N*N*N;
                for ( int j=0; j < totalTriplets; j++ ) {
                    int Z = j;
                    int X = Z / (N*N);
                        Z -= N*N * X;
                    int Y = Z / N;
                        Z -= N * Y;

                    int sumU = sumTripletsU[X,Y,Z];
                    int sumV = sumTripletsV[X,Y,Z];
                    if ( sumU != 0 ) {
                        allZeroesU = false;
                    }
                    if ( sumV != 0 ) {
                        allZeroesV = false;
                    }
                }
                if ( allZeroesU && allZeroesV ) {
                    solutions.Add( new Tuple<int,int>( derangementIndex, signCombination ) );
                } else if ( allZeroesU ) {
                    solutionsU.Add( new Tuple<int,int>( derangementIndex, signCombination ) );
                } else if ( allZeroesV ) {
                    solutionsV.Add( new Tuple<int,int>( derangementIndex, signCombination ) );
                }
            }
        }

        allSolutions[(N-3) / 2] = solutions;
    }
}

void    AddTripletU( int[,,] _sumTriplets, int _u0, int _u1, int _v, int _sign ) {
    if ( _u0 == _u1 || _u1 == _v ) throw new Exception(  "Can't have identical indices!" );
    _sumTriplets[_u0, _u1, _v] += _sign;
    _sumTriplets[_u1, _u0, _v] += _sign;
}
void    AddTripletV( int[,,] _sumTriplets, int _v0, int _u, int _v1, int _sign ) {
    if ( _v0 == _u || _u == _v1 ) throw new Exception(  "Can't have identical indices!" );
    _sumTriplets[_v0, _u, _v1] += _sign;
    _sumTriplets[_v1, _u, _v0] += _sign;
}

void    RecurseGenerateDerangements( int[] _sequence, out int[][] _derangements, out int[][] _coupless ) {
    if ( _sequence.Length == 2 ) {
        // The only sensible choice...
        _derangements = new int[][] { new int[2] { _sequence[1], _sequence[0] } };
        _coupless = new int[][] { new int[2] { 1, -1 } };
        return;
    }

    int             N = _sequence.Length;
    List< int[] >   derangements = new List<int[]>();
    List< int[] >   coupless = new List<int[]>();
    for ( int i=1; i < N; i++ ) {

        // Build the original sub-sequence containing the remaining terms (i.e. all terms except the terms at index 0 and i)
        int[]   originalSubSequence = new int[N-2];
        for ( int j=1; j <= N-2; j++ ) {
            originalSubSequence[j-1] = j < i ? _sequence[j] : _sequence[j+1];   // Skip i^th term
        }
        int[][] subSequences, subCoupless;
        RecurseGenerateDerangements( originalSubSequence, out subSequences, out subCoupless );

        // Insert the sub-sequences back
        for ( int subSequenceIndex=0; subSequenceIndex < subSequences.Length; subSequenceIndex++ ) {
            int[]   subSequence = subSequences[subSequenceIndex];
            int[]   subCouples = subCoupless[subSequenceIndex];

            int[]   derangement = new int[N];
            int[]   couples = new int[N];

            // Initial derangement is always between the first term of the sequence and any other term
            derangement[0] = _sequence[i];
            derangement[i] = _sequence[0];
            couples[0] = 1;     // First, positive couple
            couples[i] = -1;    // First, negative couple

            // Insert remaining terms
            for ( int j=1; j <= N-2; j++ ) {
                derangement[j < i ? j : j+1] = subSequence[j-1];
                couples[j < i ? j : j+1] = subCouples[j-1] < 0 ? subCouples[j-1] - 1 : 1 + subCouples[j-1];
            }

            // Store new combination
            derangements.Add( derangement );
            coupless.Add( couples );
        }
    }

    _derangements = derangements.ToArray();
    _coupless = coupless.ToArray();
}

It's not optimized at all and you start feeling it when reaching 13 dimensions because the amount of valid derangements grows like crazy...

Dimension	Derangements	Solutions
3	1	1
5	3	0
7	15	2
9	105	0
11	945	0
13	10,395	8
15	135,135	16
17	2,027,025	0
19	34,459,425	Computer Died

I was kind of expecting solutions in 15D so that would match the algebra of quaternions, octonions, sedenions, etc(?) and indeed there are but, strangely enough I wasn't prepared for finding 8 solutions in 13D...

Solutions

I give the solutions below in case you're interested. It's in the output form of my search program as a 2-tuple (derangement combination index, sign bits combination):

Solutions for 13D:

(2307, 0)
(2307, 26)

(2405, 0)
(2405, 26) ←

(7989, 0)
(7989, 26) ←

(8089, 0)
(8089, 26) ← Always the same bit combinations 26!

Interestingly enough, we notice that $2307 + 8089 = 10,396 = \#derangements+1$ and that $2405 + 7989 = 10,394 = \#derangements-1$ so these derangement combinations must be mirrors of each other...

Solutions for 15D:

(25467, 0)
(25467, 8) ←
(25467, 50) ←
(25467, 58) ←

(26405, 0)
(26405, 8) ←
(26405, 50) ←
(26405, 58) ←

(108729, 0)
(108729, 8) ←
(108729, 50) ←
(108729, 58) ←

(109669, 0)
(109669, 8) ←
(109669, 50) ←
(109669, 58) ← Always the same bit combinations!

Same as in 13D: we notice that $25467 + 109669 = 135,136 = \#derangements+1$ and that $26405 + 108729 = 135,134 = \#derangements-1$ and these derangement combinations must also be mirrors of each other...

Remarks

It would seem the amount of solutions for each valid dimension is somewhat "doubling" each time, and if we look closely, we notice that there really are only 4 solutions both in 13D and in 15D corresponding to mirror derangements, so technically we only have 2 distinct derangement solutions.
The rest are just sign variations working for each one of these 2 solutions.

So in 13D and in 15D there should really be only 2 solutions, but 13D offers 2 sign variations on each one, while 15D offers 4 sign variations.

Going back to 7D, looking back at the exact same derangement mirroring we notice there is only 1 true solution and 1 sign variation.

So maybe the actual progression is something expanding on this scheme:

Dimension	Derangements	Actual Solutions	Mirror Solutions	Sign Variations	Total Solutions = Mirror Solutions * Sign Variations
3	1	1	1	1	1
7	15	1	2	1	2
13	10,395	2	4	2	8
15	135,135	2	4	4	16

But in which way? Will the next valid dimension have 4 actual solutions and 4 sign variations to start with? Or maybe there will still be 2 actual solutions but 8 sign variations? Maybe both grow in parallel in a strange fashion? Argh. I'm dying to understand the logic behind all this but enough! Time to get back to our good old reality: 4 dimensions!

Anti-commutativity

As a final note, for each of the solutions found in all the dimensions where a cross-product-generating algebra exists, it would seem that the anti-commutativity property of the cross product (i.e. $\boldsymbol{u} \times \boldsymbol{v} = -\boldsymbol{v} \times \boldsymbol{u}$) is always a direct consequence of the algebra.

I couldn't find a dimension (up until 15D anyway) with a commutative solution but, unless we verify for all possible cases or deal with a terrifying amount of hair-pulling algebra, there's no way to really be sure (and that's the kind of research subject I know people spent a good part of their life looking into so I won't linger there, I already spent much too much time on this!)...

References

---

1. Reed, N. "Normals and the Inverse Transpose, Part 1: Grassmann Algebra" ↩↩