Maxima function to calculate the second additive compound of a square symbolic matrix

In this post, I present a function to calculate the second additive compound of a square matrix using the symbolic algebra program Maxima.

The second additive compound of a square matrix has applications, among other places, in Michael Li and James Muldowney’s autonomous convergence theorem, and in identifying Hopf bifurcations in reaction networks. A fair amount of my research has required analysing the second additive compound of symbolic matrices, and as it’s a fairly niche application I haven't been able to find any existing implementations. For this reason, I hope that the function below will prove useful to other researchers working on similar topics!

Maxima’s function declaration is rather strange (to my mind), which goes some way to explaining the weird indentation. It works as it appears below, so I thought it best not to move things around too much (although it also works when formatted on one line, as implemented in CoNtRol). The $ at the end of the function declaration tells Maxima to suppress output for that line, which I didn’t know was possible until discovering it while writing code for CoNtRol.

The following code is hereby released into the public domain.

second_additive_compound(testm) := block([testm:testm, i:1, j:1],
	matsize:matrix_size(testm)[2],
	matsize2:binomial(matsize,2),
	testm2:zeromatrix(matsize2,matsize2),
	for i1:1 while i1 <= matsize do
	    for i2:i1+1 while i2 <= matsize do
	    (
		for j1:1 while j1 <= matsize do
		    for j2:j1+1 while j2 <= matsize do
		    (
		        if i1=j1 then
		        (
		            if i2=j2 then testm2[i][j]:testm[i1][j1]+testm[i2][j2]
		            else testm2[i][j]:testm[i2][j2]
		        ) elseif i1=j2 then
		        (
		            testm2[i][j]:-testm[i2][j1]
		        ) elseif i2=j1 then
		        (
		            testm2[i][j]:-testm[i1][j2]
		        ) elseif i2=j2 then
		        (
		            testm2[i][j]:testm[i1][j1]
		        ),
		        j:j+1
		    ),
		i:i+1,j:1
	    )
	,
	return(testm2))$

To test it:

a:matrix([a11,a12,a13],[a21,a22,a23],[a31,a32,a33]);
second_additive_compound(a);