This doesn't answer your question, but it might gives some intuition as to why this is true. Given the fact that there exists a solution to the equation $x^2-Py^2=-1$, one sees that the matrix $$\left[\begin{matrix}x & Py \\ y & x\end{matrix}\right]$$ fixes $\pm\sqrt{P}$. The conjugacy class of a matrix in $GL_2(\mathbb{Z})$ (which is not reducible or finite order) is determined by the closed geodesic on the modular orbifold that it represents. This in turn is determined by a sequence of triangles which the geodesic crosses in the Farey graph. ![alt text][1] These triangles come in bunches sharing a common vertex, where the number in each bunch corresponds to coefficients of the continued fraction expansion. The matrix is conjugate to $$\pm \left[\begin{array}{cc}1 & a_1 \\ 0 & 1\end{array}\right] \left[\begin{array}{cc}1 & 0 \\ 1 & a_2\end{array}\right] \cdots \left[\begin{array}{cc}1 & a_{2n} \\ 0 & 1\end{array}\right]$$ if the determinant is 1, and to $$\pm \left[\begin{array}{cc}1 & a_1 \\ 0 & 1\end{array}\right] \left[\begin{array}{cc}1 & 0 \\ 1 & a_2\end{array}\right] \cdots \left[\begin{array}{cc}1 & 0 \\ a_{2n-1} & 1\end{array}\right] \left[\begin{array}{cc}0 & 1 \\ 1 & 0\end{array}\right] $$ if the determinant is $-1$. The number of such factors corresponds precisely to period the continued fraction expansion of fixed points of the matrix, since the closed geodesic is asymptotic in $\mathbb{H}^2$ to the geodesic connecting $\infty$ to $\sqrt{P}$. The number is even if and only if the matrix is orientation preserving, which is if and only if the determinant is 1. So the determinant is $1$ if and only if the continued fraction has even period, and the determinant is $-1$ if and only if the continued fraction has odd period. [1]: http://math.berkeley.edu/~ianagol/farey.jpg