@Joe Riel  Of course I agree the help page is confusing. At one point, it says that all entries of an rtable are "scanned" (as opposed to all opererands), and later that for sparse rtables only the nonzero ones are "scanned", so they can't both be correct.

But I expect (1) any mathematical operation on an rtable to give the same result, whether sparse of not, and (2) some algorithmic efficiencies  may be be applied if it is sparse.

So I think the first point (and the first usage of scan in terms of all entries) means I expect to see seq return the number of elements of the rectangular size. In the case of mul (which has the same strange wording about "scan"), mul(V) returns 0, (as I expect) - it involves the same result as if non sparse, and must involve more that just the stored entries.

I took the second usage of scan (that for sparse rtables only the non-zero entries are scanned), to be associated with the second point, that somehow seq uses a more efficient algorithm, and can be faster because it doesn't attempt to access the non-zero entries. 

But perhaps this is just me being optimistic that there is some resonable meaning behind the second usage (that @vv pointed to as strange)...

@eslamelidy You keep changing the worksheet so each seems unrelated to the previous ones, so suggest you start a new post. The error meesage is telling you there is a singularity, so you can only solve up to just before that point, so you need to look at the equations and figure out why - that is a math question and not a Maple question.

@vv Thanks for the improvement. I did expect there to be duplicates, but mainly because the role of the "permutation matrix" 1s and the "extra 1s" could be reversed. Removing duplicates in this case does produce 40320 = 8*7!.

There are other things that could be checked. For example, if there is a 1 on the diagonal, then there must be an extra 1 in that column and an extra 1 in that row ("guard 1s") to prevent the diagonal element being exposed as eigenvalue (which would make the polynomial reducible). In fact for the polynomial given, since the trace must be 1, and there are only two extra 1s, I could devise a more efficient scheme, though it wouldn't be general.

@dharr Thanks for the update. In the meantime, I figured out how to test this case exhaustively in about an hour, assuming I haven't made a logic error. I'll post it as an answer after I tidy it up and rerun it.

@Magma It is interesting, but perhaps not surprising that the coefficients make a big difference. (I have a hard time thinking about what carries over in the field case from regular matrices.) In your second case the coefficient of x^7 is zero, so the trace is zero, which implies an even number of 1's on the diagonal. So that is quite different from the coefficient of x^6 being zero, which is related to sums of products taken two at a time.

I agree with @vv that using theoretical ideas to limit things might make it tractable to search various subclasses, If you are interested in the general problem, why not start with smaller matrices than 8x8, in which you can do exhaustive tests.

A lot also depends on whether you are interested in the number of solutions, which is in general a hard problem, or just finding examples.

@radaar My code was not Metropolis and is not much use without the whole project. But basically I used

Statistics:-Sample(Statistics:-RandomVariable(Normal(0,sdev)),nfreqs) in a loop (I wanted Normal distribution)

Note that if you rerun a worksheet you get the same numbers, but you can use randomize() to get different numbers every time you run the worksheet. The actual sequence is very long.

@Magma You say "and also the determinant of A should be non-zero" but I'm not clear whether than means it follows automatically or it has to be another test.

I would have thought if you have an irreducible polynomial of degree n which doesn't have zero as a root (such as the one you gave), then the determinant is automatically going to be nonzero? So you only have to test if the characteristic poiynomial is equal to the required polynomial, and then you are done (or could do the determinant test only when the polynomials match to improve efficiency).

@radaar I recently used the Statistics method for a Monte Carlo method of estimating parameter errors in nonlinear regression and it seemed to work well. 

But I don't understand - your timing tests suggested it was faster than rand() and now you are saying rand is faster? Any of them can be used for Metropolis, surely.

Why don't you think the Statistics package method is random? This is what I would use.

Nice worksheet!

To simplify the log of the product to a sum you need to tell Maple N::posint (not just integer). The following works:

LL := expand(log(L)) assuming positive;

I'm assuming in the first part you didn't want to use the Statistics package, but to avoid typing errors for the PDF, I might have done

Statistics:-PDF(Normal(mu,sigma),x);

@Magma A final suggestion might be to look into the literature on totally positive matrices (all minors nonnegative) or strictly totally positive matrices (all minors positive), where there are improvements on how many determinants you need to check - for example you only need check n^2 minors for totally positive matrices, not all of them. These originate in relationships between various determinants, which might help for the MDS case though it seems the MDS case will likely be more difficult. On the other hand, all strictly totally positive matrices are MDS, so if you just want to find some MDS matrices, this might be helpful.

A good starting point is J. Pen~a, SeMA Journal, 62 (2013) 61. 

http://dx.doi.org/10.1007/s40324-013-0008-z
 

@Magma My point was mainly that doing the inverse is slower than the adjoint, and you can keep everything in integer arithmetic. And you're right that mul is probably not the best way to check all zeroes. But thinking about it, adjoint has to do the same determinants anyway, so I think the inverse idea really is not any better. I tried to do it recursively (often efficient in Maple), but it took forever.

But note that Determinant has the possibility of using many different algorithms, so you could optimize the routine to some extent. And see ?LinearAlgebra,Generic,Determinant for working in different fields. 

@Magma If Determinant(A) is 0, then it's not MDS and you are done, but otherwise just calculate the Adjoint rather than the inverse. Then for an integer matrix the Adjoint is another integer matrix. Multiplying all the entries and checking if the answer is zero is also integer-only operations. So for 16x16 this is certainly feasible. Surprisingly, mul(Adjoint(A)) is much faster than working out the Adjoint and then using mul in a seperate step [edit: no, the display just slows things down], and takes about 0.02 s on my laptop. 

 Then you can work down to the smaller matrices.

Adjoint.mw

@Thomas Richard that's embarrasing; its years since I made that mistake in my own worksheets. Perhaps pi should be prettyprinted in grey and Pi in blue, since it is a common error.

@vv taylor can do this, but I agree it shouldn't. I guess it just passes it to asympt.

Download asympt.mw