Here's a solution using lattice reduction:
1) Find degree $d$ polynomials $p_i(x)$ such that $p_i(x_j) = |M(x_i)| \delta_{i j}$.
2) Let $c_i$ be the coefficient of $x^d$ in $p_i(x)$, and $c$ the $d+1$ long column vector whose coordinates are $c_i$.
3) Find a matrix $U \in SL_{d+1}(\mathbb{Z})$ such that $U c = e$, where $e$ is the $d+1$ long column vector with a 1 in the first coordinate and zeros elsewhere. Added later:
not quite. Let $D$ be a common denominator of all the elements of $c$, and form a $d+1 \times d+1$ matrix, $A$, whose first column is $cD$ and the lower $d \times d$ block is the identity matrix (with the rest of the top row 0). Find the Hermite normal form: $U \in SL_{d+1}(\mathbb{Z})$, $U A = H$. The first column of $H$ will be zeros below the first entry, which will be a positive integer $r$. In order for there to be a solution it is necessary that $r | D$. Form a new matrix $U'$ by multiplying the top row of $U$ by $D/r$.
4) The answer (see below) is a vector in the $\mathbb{Z}$-lattice generated by the bottom $d$ rows of $U'$ which is close to the top row.
Namely, form the matrix $W$ by adjoining a $d+1 \times d+1$ identity matrix to the right of $c$. Since only the coefficient of $x^d$ matters in the answer, we can see that an answer will be given by some integer linear combination of the rows of $W$ which has $\pm 1$ in the last $d$ coordinates. The squared Euclidean length of that vector will be $d+1$, which is quite short. There are a number of algorithms for finding a closest vector (in theory for general lattices it's a hard problem, but in practice in a lattice like this it's not too hard). For a nice account of how to do it look in http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.81.8089 starting around page 14.
The idea here is to prepend a column to $U'$ which is a unit vector with 1 in the first position and 0's everywhere else. Now multiply the rest of the matrix (all columns but the first) by a large scaling factor, $s$. This will make sure that the first row will show up in the linear combination forming the shortest basis vector.
Lattice reduction will supply a short vector in the lattice which we know that our answer is. We then read off the coefficients of the $p_i$ in the last $d+1$ coordinates.
I've programmed this, and tested it on random polynomials of degree 20, and it successfully finds the $\pm 1$ combination leading to a monic polynomial.
In Tony Scholl's example of three different polynomials having the same data, the lattice generated has a lot of short vectors, so in that case one needs to enumerate short vectors to pick out the answers.
$\{0,1,2\}$
. You will run in to this same problem whenever you try to find a quadratic with $f(x_i)=\pm(x_i-x_j)(x_i-x_k)$ - there are 3 monic polynomials satisfying this permuted by $S_3$. $\endgroup$$\mathbb{Q}[X]/<(x-a_1) \cdots (x-a_d)>$
as vectors and add more coordinates for each to keep the coefficients under control. $\endgroup$