$\def\diag{\operatorname{diag}}$
The suggestion of Qiaochu Yuan is just what you need. You can do this using Smith normal form, and you are not likely to be able to do this without using Smith normal form. I assume you are only interested in the case where the determinant is non zero. In that case, the theorem regarding Smith normal form is that there exist integer matrices $U$ and $V$ with determinant equal to $\pm 1$, and a diagonal integer matrix $D = \operatorname{diag}(d_1, d_2, \dots, d_n)$ such that $A = U D V$. (Moreover, the diagonal entries satisfy $d_j$ divides $d_{j+1}$ for all $j$.) Thus $\det(A) = d_1 d_2 \cdots d_n$. Now you can get your result by factoring the diagonal matrix $D$. You can factor $D$ into (mutually commuting) diagonal matrices, each with only one diagonal entry different from 1, and with that special entry equal to some prime factor of $\det(A)$. Then you can regroup these factors to obtain the desired matrix factors with determinants equal to your chosen factors $b_j$. For example, with $D = \diag(1, 4, 12)$, we have $\det(D) = 48$ and $D$ factors as
$$D = \diag(1,2,1)\diag(1,2,1) \diag(1,1,2) \diag(1,1,2) \diag(1,1,3).$$
If your desired factors are $b_1 = 8$ and $b_2 = 6$, then you can re-group the factors above as
$$D = \diag(1, 4, 2) \diag(1, 1, 6).$$
So finally, you would take $B_1 = U \diag(1,4,2)$ and $B_2 = \diag(1,1,6)V$.
A small amount of additional fiddling might be necessary to take care of unwanted signs $\pm 1$, ask me if you are concerned about that.
If you need to implement this method computationally, you can do that, because there is an algorithm for computing Smith normal form, i.e. for computing $U, V, D$.