pairwise distribution function for 1D waves

Is there a built-in function to calculate the pairwise distribution function for points in a pair of 1D waves (especially between the same wave, i.e. an auto-pairwise distribution function)? I came up with a simple program to do so (see below), but it runs slowly (with waves that are ~10,000 points) and I expect there already is a function that does this but wasn't able to find it in the manual.

function test_distribution_function()
    wave test
    make/o/n=(numpnts(test)^2), test_DF
    variable i, j
    for(i=0;i<numpnts(test);i+=1)
        for(j=0;j<numpnts(test);j+=1)
            test_df[i*numpnts(test)+j] = test[i] - test[j]
        endfor
    endfor
   
end


Note that, despite my unfortunate naming choice, "test_DF" is not the actual distribution function. A histogram of test_DF gives the correct distribution function.
I saw the discussion in the other thread about multiplying two 1D waves to give a matrix. Here you're asking something similar (except the operation is different) and you want a 1D wave at the end. Since the matrixops are faster, this could be worth exploring to improve the speed of this procedure. Sorry I can't help more, the procedure you have is the way I'd write it.
MatrixOp is definitly worth exploring as that give you the best solutions.

Without it I'd write something like
Function test_distribution_function()

    Make/O/N=128 test = p
    variable numRows = numpnts(test)

    Make/O/N=(numRows, numRows) test_DF = test[p] - test[q]
    Redimension/N=(numRows * numRows) test_DF
 end