[Edit: 4th time's the charm, at last something sensible]
I got to this backwards: I started with another answer showing the filter-based approach, and used that to generate all the valid combinations for a series of values of $n$, and looked the sequence up on the online integer sequence database. The number of combinations is $n^2(n^2+3)$, which looks unlikely (why the 3?). It is also $t(t(n))+t(t(n-1))$ where $t(a)$ is the triangular number of $a$, $t(a)=a(a+1)/2$. Having got this, we need to know why.
The first term is simpler - the pairs of pairs where the $i \ge j$ and $tid(i,j) \ge tid(k,l)$, where $tid(a,b)$ is the triangular index of $a,b$. That is fulfilled by a function like this:
def ascendings(n):
idx = 0
for i in range(1,n+1):
for j in range(1,i+1):
for k in range(1,i):
for l in range(1,k+1):
idx = idx + 1
print(i,j,k,l)
k=i
for l in range(1,j+1):
idx = idx + 1
print(i,j,k,l)
return idx
where the second $l$ loop is because we cannot nest the $l$ loop entirely inside the $k$ loop without an if/skip to check the triangular indices.
The second term $t(t(n-1))$ is where the first pair is ascending, and the second pair is descending (note, no ==, as they are handled above).
def mixcendings(n):
idx = 0
for j in range(2,n+1):
for i in range(1,j):
for k in range(1,j):
for l in range(1,k):
print(i,j,k,l)
idx = idx + 1
k=j
for l in range(1,i+1):
print(i,j,k,l)
idx = idx + 1
return idx
The combination of both of these gives the complete set, so putting both loops together gives us the complete set of indices.
A significant problem is that these patterns are difficult to calculate for an arbitrary i,j,k,l. So I would suggest a map that yields the index given i,j,k,l. Frankly, if you are doing this at all, you might as well use the generate+filter approach, because you only need to do that once for a given $n$. The plus side of the above method is that you at least have a predictable loop structure.
In python we can write the following iterator to give us the idx and i,j,k,l values for each different scenario:
def iterate_quad(n):
idx = 0
for i in range(1,n+1):
for j in range(1,i+1):
for k in range(1,i):
for l in range(1,k+1):
idx = idx + 1
yield (idx,i,j,k,l)
#print(i,j,k,l)
k=i
for l in range(1,j+1):
idx = idx + 1
yield (idx,i,j,k,l)
for i in range(2,n+1):
for j in range(1,i):
for k in range(1,i):
for l in range(1,k):
idx = idx + 1
yield (idx,i,j,k,l)
k=i
for l in range(1,j+1):
idx = idx + 1
yield (idx,i,j,k,l)
In fortran we would just have to run the loop and store the values. We can use a simple index to store the i,j,k,l combination as a single value ($i n^3 + j n^2 + k n + l$), and store these values in an array whose index is the same as the index above. We can then iterate over this array and retrieve the i,j,k,l from the values. To get the idx for arbitrary i,j,k,l would require a reverse map and a filter to handle the symmetry, although we could probably construct a function from the above structure. The idx array generation function in fortran would be:
integer function squareindex(i,j,k,l,n)
integer,intent(in)::i,j,k,l,n
squareindex = (((i-1)*n + (j-1))*n + (k-1))*n + l
end function
integer function generate_order_array(n,arr)
integer,intent(in)::n,arr(*)
integer::total,idx,i,j,k,l
total = n**2 * (n**2 + 3)
reshape(arr,total)
idx = 0
do i=1,n
do j=1,i
do k=1,i-1
do l=1,k
idx = idx+1
arr(idx) = squareindex(i,j,k,l,n)
end do
end do
k=i
do l=1,j
idx = idx+1
arr(idx) = squareindex(i,j,k,l,n)
end do
end do
end do
do i=2,n
do j=1,i-1
do k=1,i-1
do l=1,j
idx = idx+1
arr(idx) = squareindex(i,j,k,l,n)
end do
end do
k=i
do l=1,j
idx = idx+1
arr(idx) = squareindex(i,j,k,l,n)
end do
end do
end do
generate_order_array = idx
end function
And then loop over it thus:
maxidx = generate_order_array(n,arr)
do idx=1,maxidx
i = idx/(n**3) + 1
t_idx = idx - (i-1)*n**3
j = t_idx/(n**2) + 1
t_idx = t_idx - (j-1)*n**2
k = t_idx/n + 1
t_idx = t_idx - (k-1)*n
l = t_idx
! now have i,j,k,l, so do stuff
! ...
end do
