@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).

@Carl Love 

OK, but we cannot assume that x-1 and x-3 are the only OP's polynomials.