Let me first state that functional determinants are not the same as determinants. So point A is wrong. About B, there is no single best way to do it, it is generally a hard task. Let me know elaborate.
So a functional determinant is the determinant of a differential operators, that is the determinant of of a linear functional in a function space, that means we are very likely dealing with an infinite dimensional vector space (our space of functions).
One way to compute a determinant is indeed by taking a product of eigenvalues, so you would expect
$$\det(\mathcal{D}) = \prod_\lambda \lambda $$
where $\lambda$'s are the eigenvalues of the differential operator $\mathcal{D}$. However the situation is generally not so simple since, we no longer have necessarily discrete eigenvalues and moreover zero-modes would imply the determinant is zero.
So we are looking for functions $\phi(x)$ such that
$$\mathcal{D}\phi(x) = \lambda \phi(x)$$
if you take $\mathcal{D} = \partial^2 + m^2$ you see why $\pm im$ are not the eigenvalues.
A simple physicist argument for the case of the Klein-Gordon operator is to consider plain waves. Assume $phi_k(x) = e^{-i k\cdot x}$, where $k=(k^0,k^1,k^2,k^3)$ is for now some four-vector. Trying this solutions one gets
$$\mathcal{D}\phi(x) = (-k^2+ m^2)e^{-ik\cdot x}$$
you can then identify the eigenvalues to be $-k^2+m^2$ (notice, we are not looking for solutions to the Klein-Gordon equation, but fur eigenvalues, therefore I don't need to impose any on-shell condition on $k$).
At this point you would be tempted to write
$$\det{\mathcal{D}} = \prod_{k}(-k^2 + m^2) $$
here you can already spot there are several issues. There will be some $k$'s which make that vanish and therefore the whole thing vanish and even if extracted the product doesn't seem to converge at all.
Each of this problems have to be solved with care. One can for example either extract the zero modes when they are discrete and exchange them for collective coordinates or perturb the operator to eliminate the zero-modes and later take the limit. For the divergence behavior some sort of regularization is needed. In applications one usually regularizes by taking ratios of functional determinants, which can have a more physical interpretation since they become finite. E.g.
$$ \frac{\det\mathcal{D}\phantom{m}}{\det\mathcal{D}_{free}}$$
Alternatively you could still obtained a regularized expression for the logarithm of the determinant for free Klein-Gordon op. by dimensional regularization, as shown in Peskin & Schröder p374, where a Wick rotation is performed and then the integral solved for $d$-dimensional Euclidean space. That should make clear what the computation more or less should look like.
To answer your point B, I can only say generally you want to first identify what sort of spectrum you are dealing with, whether it has a discrete or a continuous part or both. And then you can attempt to sum/multiply them in some smart way. There are methods to compute functional determinants that rely only on zero-eigenvalues such as the Gelfand-Yaglom method, or methods were it is enough to know that the spectrum will be continuous so that one can integrate over them (for the case were one is after the $\log$ ). All in all it is a delicate matter and as far as I know there is no general recipe.
