Your question defines a set of matrices you want to sample from. However as pointed out in the comments, the notion of sampling depends on some distribution or probability measure over the set.

As a simple example, consider the $\delta$ distribution around some matrix $M_0$ that satisfies your constraints - this is a perfectly valid distribution and it answers your requirements because 'sampling' from it will always give you the matrix $M_0$, but this is probably not what you want.

Now having that in mind we can try to come up with of with some distribution that might be reasonable for your purpose (with the caveat that without knowing the purpose of this exercise it is hard to tell what is 'reasonable'): note that your specific set of constraints defines matrices that are close to unitary, since if $M$ is unitary ($M^TM=I$) then both $f(M)=1$ and $g(M)=0$. So, let

$$M = U(I + \epsilon X)$$

Where $U$ is unitary, $X$ is a random matrix (with some distribution) and $\epsilon$ is a small real number, therefore

$$ \det(M)=1 + \epsilon \text{tr}(X) + O(\epsilon^2) $$

$$M^TM - I = \epsilon(X+X^T) + O(\epsilon^2) $$

For a given $X$ this defines a possible range for $\epsilon$, since you require that

$$ -\epsilon_1 < \epsilon \text{tr}(X) < \epsilon_2 $$ $$ -\epsilon_1 < |\epsilon| \cdot \|\frac{X+X^T}{2}\|_{\infty} < \epsilon_2 $$

(For the second constraints only the upper bound is actually relevant since the norm can only be positive). So, generate a random matrix $X$ form some distribution, then generate a random $\epsilon$ from the range defined above (or alternatively you can use different distributions for the diagonal and off diagonal elements of $X$ in a way that will automatically satisfy the constraints): this will give you a matrix $M$ that satisfies the constraints with high probability (not surely because of the $O(\epsilon^2)$ terms).

A random unitary matrix $U$ can be generated in many ways, for example using Cayley transform: $U = (I - A)(I + A)^{-1}$ is unitary if $A$ is anti-symmetric. Therefore generate another random matrix $X$ and use its anti-symmetric part $X-X^T$ to get $U$.

You can generalize the above procedure to other functions $f$ and $g$ by expanding them in a similar way with a small parameter $\epsilon$ around $M_0$ that satisfies $f(M_0)-1=g(M_0)=0$ (because you assume $\epsilon_{1,2}$ are small). However if this condition defines a set of matrices that you want to sample from, you will have to find some parametrization that will allow you to do that (the equivalent of generating the unitary matrices in your case). There is probably no universal recipe for that so it will have to be worked out case by case.