@Carl Love 

It works simply because the inner evalf has Digits=Digits0=10 >= 5.

But a sigfig is not simple. evalf[n] for a mathematical function gives in general n significant correct digits, but for an expession it's difficult because we must estimate the needed precision (to compensate e.g. the catastrophic cancellations).
Mathematica tries to do this, but fails sometimes. Maple prefers to let the user choose the precision for the desired accuracy. However for some commands (e.g. Int) there is an "epsilon" for accuracy.

f := x -> sqrt(x^2+1)-x + 1/3:
g := x -> 1/(sqrt(x^2+1)+x):
a:=10^9/3: # exact
evalf[15](f(a)-g(a)); # exact = 1/3;
                       0.333332998500000
a:=10^15/3: # exact
evalf[15](f(a)-g(a)); # exact = 1/3;
                    -1.50000000000000*10^(-15)  

@Glowing 

restart
p1 := 1007.0:
p2 := 1014.0:
evalf(evalf(p2 - p1), 2);   
# The outer evalf passes Digits=2 to the inner evalf
# For Digits=2 we have p1 - p2 = 1e3 - 1e3 = 0.
# So, the outer evalf does evalf[2](0.) = 0.  

                               0.
evalf(evalf(p2 - p1, 4), 2);
# Now the inner evalf sets Digits from 2 to 4 ==> 7.0
                               7.

Note that

evalf[2](1014.0-1007.0);
           7.

due to automatic simplification, which is done with Digits=10.

@Kitonum 

We cannot compare evalf[n] with Mathematica's N because N attempts to give the result with n digits precision
while evalf[n] uses n digits precision.
I think that it would be useful to explain how Maple obtains the following answers (evaluation rules, option remember):

restart;
m:=206:  x:=210.0:
evalf[2](`+`(x,-m)) , (evalf[2]@`+`)(x,-m), (evalf[2]@`+`)('x',-m);
                       0., 4.0, x - 210.
 

restart;
m:=206:  n:=210:  x:=210.0:
evalf[2](x-m), evalf[2](n-m);  
                             0., 4.

restart;
m:=206:  n:=210: z:=0.0:
evalf[2](n-m+z), evalf[2](2*z+n-m);  
                             4., 0.
restart;
m:=206:  n:=210: z:=0.0:
evalf[2](n-m+z), evalf[2](z+n-m);  
                             4., 4.
restart;
m:=206:  n:=210: z:=0.0:
evalf[2](z+n-m), evalf[2](n-m+z);  
                             0., 0.
restart;
m:=206:  n:=210: z:=0.0:
evalf[2](2*z+n-m), evalf[2](n-m+z);
                             0., 4.

@Carl Love 

1.

Yes, it's curious that RowReduce is faster, because for Determinant, the default method is REF.
BTW, replacing in MDS Determinant by RowReduce, in Maple 2019 I get [each after restart]:

CodeTools:-Usage(MDS2(A80,8)):
memory used=1.02GiB, alloc change=32.00MiB, cpu time=15.96s, real time=15.96s, gc time=265.20ms
CodeTools:-Usage(IsMDS(A80,8)):
memory used=2.65GiB, alloc change=60.50MiB, cpu time=15.27s, real time=15.20s, gc time=795.61ms
(In Maple 2018 the results are the same.)

2.

n <= m <= n^2 - n +1

both limits being attained.

@Carl Love 

But of course, if A has integer[], Mod is not necessary, e.g.
A:= LA:-RandomMatrix(64$2, generator= rand(0..1), datatype=integer[kernelopts(wordsize)/8]):

Merry Christmas !

@Magma 

LinearAlgebra:-Modular:-MatrixPower(2, A, k)

@Magma

It also appears in the definition of B[i,j].
I am curious if it takes <1 min. I do not have a 64 MDS matrix to check.

@Magma 

It seems that you have a 32 bit Maple. Please replace integer[8] with integer[4] (twice)

@Magma 

I think you have copied the code while I was editing. Please copy again.

@Carl Love 

Yes. Actually division ring is enough (commutativity is automatic) by Wedderburn's theorem.

@Magma 

Thanks, everything is clear now. So, A is not an arbitrary binary matrix, it is probably derived from a matrix over GF.
I wonder whether this approach is faster than computing the determinants over GF(2^r).
I'll try to find an enhanced version soon.

@Magma

The MDS definition is not in terms of GF. Actually GF is not mentioned at all in the first paper.
That is why I have asked. So, do you know the answers (mainly the first one)?

@Magma 

I did not know about this technique. So, it seems that each r x r block of B corresponds to an element of GF(2^r).

Is it known an embedding of GF(2^r) into the M_r(Z2) ? Or, here it is used another definition for MDS?

@Jjjones98 

I don't see how such a solution could be useful. It is inherently huge (almost each coefficient being a symbolic expression).

Not to mention that for some values of the parameters the system will be inconsistent or undetermined, i.e. the obtained solution is generic only. A complete solution (with all possibilities) would be much much longer.
 

@CyberRob 

f is already in the form you describe. Or, use
collect(f, [k12, k21, Ve], distributed);
It would be better to present mathematically the entire problem (preferably with generic notations: x,y,...  the unknowns,  a,b,... the parameters).