$\newcommand{\la}{\lambda}\newcommand{\R}{\mathbb{R}}$Let $I:=(0,1)$ and $Q=I^n$. Let $\la$ denote the Lebesgue measure over $\R^n$. Let $P_j$ denote the partial derivative of $P$ wrt to its $j$th argument. Because the polynomial $P=P(x_1,\dots,x_n)$ is non-constant, at least one of the $P_j$'s is a nonzero polynomial. Without loss of generality, $P_1$ is a nonzero polynomial.
It follows by the main theorem of algebra and the Tonelli theorem that the Lebesgue measure $\la(N)$ of the (closed) set \begin{equation*} N:=\{(x_1,\dots,x_n)\in \R^n\colon P_1(x_1,\dots,x_n)=0\} \tag{10}\label{10} \end{equation*} of the zeros of the nonzero polynomial $P_1$ is $0$.
Consider now the transformation
$$\R^n\ni(x_1,\dots,x_n)\mapsto g(x_1,\dots,x_n):=(P(x_1,\dots,x_n),x_2,\dots,x_n)\in\R^n, $$
the polynomials
\begin{equation*}
p_{(y_1,\dots,y_n)}(x_1):=P(x_1,y_2,\dots,y_n)-y_1,
\end{equation*}
the sets
\begin{equation*}
Z_{(y_1,\dots,y_n)}:=\{x_1\in \R\colon(x_1,y_2,\dots,y_n)\notin N, p_{(y_1,\dots,y_n)}(x_1)=0\},
\end{equation*}
for $(y_1,\dots,y_n)\in\R^n$,
the set
\begin{equation*}
Y:=\{(y_1,\dots,y_n)\in\R^n\colon p_{(y_1,\dots,y_n)} \text{ has no multiple roots}\},
\end{equation*}
and the sets
\begin{equation*}
A_m:=\{(y_1,\dots,y_n)\in Y\colon\text{card}\,Z_{(y_1,\dots,y_n)}=m\}
\end{equation*}
for $m=0,\dots,d$, where $\text{card}$ denotes the cardinality and $d$ is the degree of $P$ wrt $x_1$. Note that $(A_0,\dots,A_d)$ is a partition of $Y$, and the sets $A_1,\dots,A_m$ are open.
By the main theorem of algebra and the Tonelli theorem, the open set $Y$ is of full Lebesgue measure.
By the implicit function theorem, for each $m\in\{0,\dots,d\}$, each $(y_1,\dots,y_n)\in A_m$, and each integer $k\in[1,m]$, there is a real-analytic function $X_{m,k}$ defined on a neighborhood $U_{(y_1,\dots,y_n)}$ of the point $(y_1,\dots,y_n)\in A_m$ such that \begin{equation*} g(X_{m,k}(z_1,\dots,z_n),z_2,\dots,z_n)=(z_1,\dots,z_n) \end{equation*} for all $(z_1,\dots,z_n)\in U_{(y_1,\dots,y_n)}$, and the value of the Jacobian determinant of the map \begin{equation*} U_{(y_1,\dots,y_n)}\ni(z_1,\dots,z_n)\mapsto(X_{m,k}(z_1,\dots,z_n),z_2,\dots,z_n) \end{equation*} at the point $(z_1,\dots,z_n)$ is \begin{equation*} \partial_{z_1}X_{m,k}(z_1,\dots,z_n)=\frac1{P_1(X_{m,k}(z_1,\dots,z_n),z_2,\dots,z_n)}. \end{equation*}
So, by the change of variables in multiple integrals, the density (say $f$) of the pushforward $g_\#(\la\downharpoonright Q)$ under $g$ of the measure $\la\downharpoonright Q$ wrt to $\la$ is given by the formula \begin{equation*} f(y_1,\dots,y_n)=\sum_{m=1}^d \sum_{k=1}^m\frac{1((y_1,\dots,y_n)\in B_{m,k})}{|P_1(X_{m,k}(y_1,\dots,y_n),y_2,\dots,y_n)|} \tag{20}\label{20} \end{equation*} for $(y_1,\dots,y_n)\in\R^n\setminus g(N)$, where \begin{equation*} B_{m,k}:=\{(y_1,\dots,y_n)\in A_m\colon(X_{m,k}(y_1,\dots,y_n),y_2,\dots,y_n)\in Q\}. \end{equation*} Note that $\la(g(N))=0$, since $\la(N)=0$ and $g$ is Lipschitz. So, we can let $f:=0$ on $g(N)$.
Now, the density (say $f_1$) of the pushforward $P_\#(\la\downharpoonright Q)$ of the measure $\la\downharpoonright Q$ under $P$ wrt to $\la$ is given by the formula \begin{equation*} f_1(y_1)=\int_{\R^{n-1}}dy_2\,\cdots\, dy_n\,f(y_1,\dots,y_n) \tag{30}\label{30} \end{equation*} for $y_1\in\R$.
More detailed information can be obtained using cylindrical algebraic decomposition and the Tarski–Seidenberg theorem.
Indeed, using cylindrical algebraic decomposition, one can partition $\R^n$ into finitely many connected semialgebraic sets $(C_i)_{i\in I}$ called cells, on which the polynomial $P_1$ has constant sign, either $+$, $-$ or $0$, such that for any $i$ and $j$ in $I$ one has either $\pi(C_i)=\pi(C_j)$ or $\pi(C_i)\cap\pi(C_j)=\emptyset$, where $\pi$ is the projection of $\R^n$ onto $\R$ consisting in removing the last $n-1$ coordinates.
Moreover, by the Tarski–Seidenberg theorem, $\pi(C_i)$ is a semialgebraic set for each $i\in I$.
Thus, we have a partition $(I_k)_{k\in K}$ of the finite set $I$ and a partition $(D_k)_{k\in K}$ of $\R$ into semialgebraic sets $D_k$ such that for each $k\in K$ \begin{equation*} D_k=\pi(C_i)\text{ for }i\in I_k. \end{equation*}
For each $k\in K$, the semialgebraic subset $(0,1)\cap D_k$ of $\R$ is the finite disjoint union of intervals $D_{k,l}$: \begin{equation*} (0,1)\cap D_k=\;\cdot \hspace{-10pt}\bigcup_{l=1}^{L_k} D_{k,l}. \end{equation*} Introducing now the connected semialgebraic sets \begin{equation*} C_{k,i,l}:=C_i\cap\pi^{-1}(D_{k,l})\text{ for }i\in I_k \end{equation*} $k\in K$, and $l\in[L_k]$, we get the finite partition $(C_{k,i,l}\colon k\in K, i\in I_k, l\in[L_k])$ of $(0,1)\times\R^{n-1}$ into connected semialgebraic sets and the finite partition $(D_{k,l}\colon k\in K, l\in[L_k])$ of the interval $(0,1)$ into intervals $D_{k,l}$ such that \begin{equation*} D_{k,l}=\pi(C_{k,i,l}) \end{equation*} for all $k\in K, i\in I_k, l\in[L_k]$. In each of the cells $C_{k,i,l}$, the sign (say $s_{k,i,l}$) is constant. Also, the union of all $C_{k,i,l}$'s with $s_{k,i,l}=0$ is contained in the set $N$ defined by \eqref{10}, and $\la(N)=0$, as was established previously.
So (cf. \eqref{20}),
\begin{equation*}
f(y_1,\dots,y_n)=\sum_{k=1}^K \sum_{i\in I_k}\sum_{l=1}^{L_k}
\frac{1((y_1,\dots,y_n)\in C_{k,i,l})}{|P_1(X_{k,i,l}(y_1,\dots,y_n),y_2,\dots,y_n)|} \tag{40}\label{40}
\end{equation*}
for certain real-analytic functions $X_{k,i,l}$ defined on the corresponding cells $C_{k,i,l}$
$(y_1,\dots,y_n)\in(0,1)\times\R^{n-1}\setminus g(N)$. Again, we can let $f:=0$ on $(0,1)\times\R^{n-1}\cap g(N)$.
It follows from \eqref{40} that $f$ can explode near the boundaries of cells $C_{k,i,l}$ only polynomially. So, by \eqref{30}, $f_1$ can explode near the endpoints of the intervals $D_{k,l}$ only polynomially.