How do you compute the primes of bad reduction? Suppose that I am given a subscheme $Y$ of $\mathbf{P}^n_{\mathbf{Z}}$, flat over $\operatorname{Spec}\mathbf{Z}$ and with smooth generic fiber $Y_{\mathbf{Q}}$, defined by the vanishing of some homogeneous polynomials
$$
F_1, \ldots, F_k \in \mathbf{Z}[X_0,\ldots,X_n].
$$
How does one determine the set $S$ consisting of primes $p$ such that the fiber $Y_{\mathbf{F}_p}$ is non-smooth? 
The way I'm going about the problem seems a bit crude to me, and has some problems. In particular, I can't prove that my algorithm solves the problem.
For simplicity suppose that $k=1$, and write $F=F_1$. In a nutshell, I compute the derivatives $f_i := \frac{\partial F}{\partial X_i}$, and determine an integer $N$ such that $N$ is contained in the ideal $I = (F, f_1, \ldots, f_n) \subset \mathbf{Z}[X_0,\ldots,X_n]$  by repeatedly taking resultants of polynomials in $I$. Then at least for all $p$ not dividing $N$, the fiber $Y_{\mathbf{F}_p}$ is smooth, since for such $p$ the ideal $I_p \subset \mathbf{F}_p[X_0,\ldots,X_n]$ generated by the reduced polynomials $\widetilde{F}$, $\widetilde{f}_1$,$\ldots$,$\widetilde{f}_n$ contains $1$.
One problem is that it's quite expensive computationally. I'm using resultants to eliminate the variables $X_i$ one by one, say in the order $i=0,1,\ldots,n$. Then every time I've eliminated one of the $X_i$s, I get a lot of polynomials that don't contain $X_0$ up to $X_i$, and I take them all into the next elimination round. (You have to be careful when the resultant comes out $0$, and perhaps there are some other subtleties that I'm forgetting, but this is basically it.) 
Another problem with this is that the final outcome seems to depend on the order in which I take the $i$s. (This may have something to do with singularities being located at infinity with respect to one of the $X_i$s [edit: on second thought, this doesn't make any sense; see comment by François], but I really can't see the geometry of what I'm doing clearly enough to be confident about this.) And I don't see why my method necessarily gives me the whole $S$ - instead of some upper bound on it - even when all orderings of the $X_i$ are taken into account. 
My questions: 


*

*Does my method give the whole set $S$? 

*If yes, why? If no, what method does work?

*Regardless of the answers to 1. and 2., is there a faster or more natural way of finding the primes of bad reduction?

 A: Assume for simplicity that $Y$ has pure relative dimension $d$. By considering the standard affine cover of $\mathbf{P}^n_{\mathbf{Z}}$, one easily reduces to the case where $Y=V(F_1,\ldots,F_k)$ is a closed subscheme of $\mathbf{A}^n_{\mathbf{Z}}$. Then the special fiber $Y_p = V(F_{1,p},\ldots,F_{k,p}) \subset \mathbf{A}^n_{\mathbf{F}_p}$ is smooth if and only if for every $x \in Y_p(\overline{\mathbf{F}}_p)$, the Jacobian matrix $(\frac{\partial F_{i,p}}{\partial X_j}(x))$ has rank $n-d$. Note that $Y_p$ is $d$-dimensional by the flatness assumption, so all geometric tangent spaces have dimension $\geq  d$, which implies that the rank of the Jacobian matrix is everywhere $\leq n-d$. Thus $Y_p$ is smooth if and only if the ideal generated by $F_{1,p},\ldots,F_{k,p}$ together with all $(n-d)$-minors of the Jacobian matrix is equal to $\mathbf{F}_p[X_1,\ldots,X_n]$.
Now consider the ideal $I$ of $R=\mathbf{Z}[X_1,\ldots,X_n]$ generated by $F_1,\ldots,F_k$ together with all $(n-d)$-minors of the Jacobian matrix. Since the generic fiber of $Y$ is smooth, we have $I \cap \mathbf{Z}=N\mathbf{Z}$ for some integer $N \geq 1$, and the prime factors of $N$ are precisely the primes of bad reduction of $Y$.
Proof. If $p$ doesn't divide $N$ then $(I,p)$ contains $N$ and $p$, so $(I,p)=R$ and $Y_p$ is smooth. If $p$ divides $N$, write $N=p^k M$ with $p$ not dividing $M$. If we had $(I,p)=R$, then we would also have $(I,p^k)=R$. Then $M \in (MI,p^k M)=(MI,N) \subset I$, a contradiction. Thus $Y_p$ is not smooth.
You can compute the integer $N$ using Gröbner bases for polynomials over $\mathbf{Z}$, see e.g. http://magma.maths.usyd.edu.au/magma/handbook/text/1112#12189 for the Magma implementation. I'm not expert in Gröbner bases, so I would appreciate if someone could confirm whether it is always possible to find $I \cap \mathbf{Z}$ when a set of generators of $I$ is given, and whether the Magma implementation works in all cases.
