$$\begin{align} px+qy+rz &=(p^2+q^2+r^2)\cos\theta \\[4pt] ax+by+cz &=0 \\[4pt] \sqrt{x^2+y^2+z^2} &=\sqrt{p^2+q^2+r^2} \end{align}$$ $$ap+bq+cr=0$$
How would one express the solution $(x,y,z)$ in terms of real constants $a,b,c,p,q,r$ and $\theta$?
What I have so far: If I choose z to be a free parameter, then $$ \begin{bmatrix} p & q \\ a & b \end{bmatrix} \begin{bmatrix} x \\ y \end{bmatrix}=\begin{bmatrix} (p^2+q^2+r^2)cos(\theta)-rz \\ -cz \end{bmatrix} $$ and $$ \begin{bmatrix} x \\ y \end{bmatrix}=\begin{bmatrix} b((p^2+q^2+r^2)cos(\theta)-rz) +cqz \\ -a((p^2+q^2+r^2)cos(\theta)-rz)-cpz \end{bmatrix}1/(bp-aq) $$ so $$\begin{bmatrix} x \\ y \\ z\end{bmatrix}=\begin{bmatrix} b((p^2+q^2+r^2)cos(\theta)-rz) +cqz \\ -a((p^2+q^2+r^2)cos(\theta)-rz)-cpz \\ z(bp-aq)\end{bmatrix}1/(bp-aq)$$ Then, the rest is a straightforward polynomial in z $\sqrt{x^2+y^2+z^2}=\sqrt{p^2+q^2+r^2}=\sqrt{(b((p^2+q^2+r^2)cos(\theta)-rz) +cqz)^2/(bp-aq)^2+(-a((p^2+q^2+r^2)cos(\theta)-rz) -cpz)^2/(bp-aq)^2+z^2}$
When the algebra becomes messy, my method is to repeatedly define groups of constants as new constants. In your problem, you only care about $x,y,z$: you don't need to worry about the interactions of the remaining constants until the end where you express the final solution in terms of them. So for example, you could define common constant expressions in your final equation as such $$d=\sqrt{(b(d^2 cos(\theta)-rz) +cqz)^2/e +(-a(d^2 cos(\theta)-rz) -cpz)^2/e+z^2} \\ d = ||(p,q,r)|| = \sqrt{p^2+q^2+r^2}, \quad e = (bp-aq)^2 $$
But this is still ugly, so if you square both sides, expand the $b()$ and $a()$ then define some more, you get this relatively cleaner equation.
$$f_1 = (f_2 z + f_3)^2 +(f_4 z + f_5)^2 + e z^2 \\ f_1 = d^2e, \quad f_2 = cq - br, \quad f_3 = bd^2 cos(\theta), \quad f_4 = a r -cp, \quad f_5 = -a d^2 cos(\theta)$$
Hopefully the $f$ equation is now clean enough for you to solve for $z$ in terms of $f$'s and $e$. Then you could substitute for $z$ in your $x,y$ equations, then substitute the $f$'s, then $d,e$. I agree that the final solution will still be ugly, but maybe it won't be as hard for you to get there now.
Alternative thoughts in case they lead to a simpler solution approach:
Interpreting the algebraic system of equations as vector constraints, you want to solve for $\vec{X}$ given $\vec{P}, \vec{A}, \theta$ in the following. $$\vec{P} \cdot \vec{X} = ||P||^2 \cos \theta \\ \vec{A} \cdot \vec{X} = 0 \\ ||X|| = ||P||$$
Geometrically, you want to find a vector with same magnitude as $\vec{P}$ but forming an angle of $\theta$ with it and is also perpendicular to $\vec{A}$. So in line with your comment, you would want to find something like, $$\vec{X} = M \vec{P}$$
where $M$ is a rotation matrix that rotates by $\theta$ in a rotational plane containing $\vec{A}$ such that the resulting vector is in the plane normal to $\vec{A}$.