floor(fsolve(1/(n+1)! = 0.0001, n)) + 1;

You can't use Dirac like that. Dirac is not a function; it is a distribution and must be handled with great care, and only if the user knows exactly its meaning. In Maple the use of the Dirac "function" is restricted to some contexts, e.g. as in the input and output of dsolve and in integral transforms (Fourier, Laplace).

The answer is no.
You will have to use the theory of Parametric LP.

The modp1 arithmetic is very fast in maple (but especially when the modulus is a hardware integer).
The usual mod operations use modp1, so they are fast too.
Using modp1 directly, the gain is modest:

restart;
m := 8009;
N := ceil(m/2) + 2;
read "nn.txt":
st := time[real]():
r := modpol(a*a, f_t, t, 2^N):
t1 := time[real]() - st;

N2:=2^N:
A:=modp1(ConvertIn(a,t),N2):
F:=modp1(ConvertIn(f_t,t),N2):

st := time[real]():
R:=modp1(Rem(Multiply(A,A),F), N2):
t2 := time[real]() - st;

rr:=modp1(ConvertOut(R,t),N2):
r-rr; # 0 
t2/t1;

st := time[real]():
modp1(Rem(Power(A,2),F), N2):
t3 := time[real]() - st;
t3/t1;

So, probably in your case that's all Maple can do.

For a general approach, you can use the package QuadraticInt (for working in the ring Z[sqrt(d)]); see the Application Center.

Your Bst is the parametric representation of a surface.
But it is in the plane z = 0. So you obviously have a single "contour", but this "contour" is actually a 2-dimensional set in R^2 (it has interior points).

isolve is not able to manage inequalities in this case. Select the desired solutions later:

isolve(29 = x^2 + y^2):
select(u -> (eval(x,u)>=1 and eval(y,u)>=1), [%]);

It does not work because eval needs subexpressions (i.e. operands of some level).
[E.g.  a*b  is not subexpression in a*b*c;  for such substitution there is algsubs, but in Physics almost everything is redefined, so this probably fails too.]

But here you can simply use:

eval((2), Dagger(A[i]) = -A[i]);

If you just want symbolic computations, you may use the Physics package:
with(Physics):
Setup(noncommutativeprefix = {A,B});

Now, the adjoint (Hermitian conjugate) ^* will act (almost) as a generic involution (displayed: Dagger).

(A^(-1) * B)^*; # ==>  B^*  * (A^*)^(-1)

Download mimiz.mw

A faster version.

F := proc(N::posint)
  local V:=Vector(4), p:=2;
  while (p:=nextprime(p))<N do V[(irem(p,8)+1)/2]++ od;
  [seq(V)]
end proc:

A := Matrix(2, 3, [[1, 2, 3], [3, 1, 2]]); B := Matrix(2, 3, [[1, x, 3], [y, 1, z]]);
solve(Equate(A,B), {x,y,z});

Note that solve needs a list or set for its first argument, so solve(A=B, {x,y,z}); will not work.
The second argument is optional, so, solve(Equate(A,B)); is enough.

simplify ... assuming ... should be enough.
Unfortunately an extra step is needed, for example:

arctan(cos(beta)/sin(beta));
simplify(convert(%, tan)) assuming beta>0, beta<Pi/2;

   Pi/2 - beta