I recommend you get Binary Quadratic Forms by Duncan A. Buell.
Theorem 4.23 on page 74 is, when $p$ is an odd prime with $(\Delta | p) = 1,$ and $b$ is any integral solution to $b^2 \equiv \Delta \pmod {4p},$ the $p$ is primitively represented by the binary quadratic forms belonging to the classes of $\langle p, b, \ast \rangle$ and its opposite $\langle p, -b, \ast \rangle$
In you case, the discriminant $\Delta$ of $x^2 - 2 y^2$ is 8. There are no imprimitive forms of this discriminant, so we may ignore the word "primitively." On page 30, we find that there is only one equivalence class of forms of this discriminant. As a result, any prime $p \equiv \pm 1 \pmod 8$ is represented by $x^2 - 2 y^2.$
How to find a representation? Using TONELLI-SHANKS solve $b^2 \equiv 8 \pmod p.$ If $b$ is odd, replace it by $p-b$ to get an even value, so now $b^2 \equiv 8 \pmod {4p}.$ That is, $b^2 = 8 + 4 p t,$ Which is to say $$ b^2 - 4 pt = 8.$$ So we have an indefinite form, almost certainly not "reduced,"
$$\langle p, b, t \rangle.$$
Pages 21-26 in Buell tell you everything necessary about indefinite reduced forms, with some positive discriminant $D$ not a square, and some coefficients $\langle a, b, c \rangle$ with $b^2 - 4 a c = D,$ the form is reduced if
$$ 0 < b < \sqrt D \; \; \mbox{and} \; \; \sqrt D - b < 2 |a| < \sqrt D + b $$
The reduced form of discriminant 8 are $\langle 1,2,-1 \rangle$ and $\langle -1,2,1 \rangle,$ these being equivalent.
Finally, on page 22, we have Proposition 3.3, any indefinite form is equivalent to a reduced form of the same discriminant. Furthermore, in the body of the proof, Buell shows exactly how to get from non-reduced $\langle p, b, t \rangle$ to one of the reduced forms mentioned, also resulting in a 2 by 2 matrix of determinant 1, let us call it $R,$ which provides the (invertible) change of variables. Inverting $R$ tels us how to get from the reduced form to $\langle p, b, t \rangle.$ A little extra care lets us find a solution to $v^2 - 2 w^2 = p.$
The fact that $x^2 - 2 y^2$ integrally represents all primes $p \equiv \pm 1 \pmod 8$ is due to people such as Lagrange (Franz would know). Lagrange would have seen this material in the way I describe, if you can find it see Introduction to the Theory of Numbers by Leonard Eugene Dickson.
I have software I wrote that takes an indefinite form and reduces it, then shows the entire "chain" or "cycle" of reduced equivalent forms. Email me if you would like a copy (it is in C++). Also, your name seems very familiar, perhaps you are Marvin Greenberg's friend who wrote to me once or twice? Not sure.
A note about speed. This all works very rapidly. It is essentially finding the repeated part of a continued fraction for a quadratic irrational. Indeed, the special thing about reduced forms is that the continued fraction for the largest root is purely periodic.
EDIT, 25 January 2012: I finally got the reduction step cleaned up with a nice printout. I find a pretty example, $p = 159287$ and $45991^2 \equiv 8 \pmod {159287}.$ Then
$159287 - 45991 = 113296$ is even, and the constructed indefinite form of discriminant $D = 8$ is
$$ \langle 159287, \; 113296, \; 20146 \rangle. $$
The reduction step is this: we find the (unique) integral $\delta$ such that
$$ \sqrt D - 2 |c| < -b + 2 c \delta < \sqrt D. $$
Then the next form is
$$ \langle c, \; -b + 2 c \delta , \; a - b \delta + c \delta^2 \rangle. $$
The change of variables matrix is just
$$ \left( \begin{array}{cc} 0 & -1 \\ 1 & \delta \end{array} \right) $$
As we follow through several steps, we just keep a running account of the matrix. At the end, we have the matrix, written on the right, that takes the original form to the resulting reduced form. The inverse takes the reduced form back. Here is the output:
jagy@phobeusjunior:~/old drive/home/jagy/Cplusplus$ ./indefCycle
Input three coefficients a b c for indef f(x,y)= a x^2 + b x y + c y^2
159287 113296 20146
0 form 159287 113296 20146 delta 2
1 form 20146 -32712 13279 delta -2
2 form 13279 -20404 7838 delta -2
3 form 7838 -10948 3823 delta -2
4 form 3823 -4344 1234 delta -2
5 form 1234 -592 71 delta -5
6 form 71 -118 49 delta -2
7 form 49 -78 31 delta -2
8 form 31 -46 17 delta -2
9 form 17 -22 7 delta -2
10 form 7 -6 1 delta -2
11 form 1 2 -1
85 -101
-239 284
To Return
284 101
239 85
0 form 1 2 -1 delta -2
1 form -1 2 1 delta 2
2 form 1 2 -1
minimum was 1rep 1 0 disc 8 dSqrt 2.8284271247 M_Ratio 8
Automorph, written on right of Gram matrix:
-1 -2
-2 -5
Trace: -6 gcd(a21, a22 - a11, a12) : 2
=========================================
jagy@phobeusjunior:~/old drive/home/jagy/Cplusplus$
The matrix of the form $ \langle 1, \; 2, \; -1 \rangle $ is
$$ \left( \begin{array}{cc} 1 & 1 \\ 1 & -1 \end{array} \right) $$
and the part about returning is
$$ \left( \begin{array}{cc} 284 & 239 \\ 101 & 85 \end{array} \right) \cdot \left( \begin{array}{cc} 1 & 1 \\ 1 & -1 \end{array} \right) \cdot \left( \begin{array}{cc} 284 & 101 \\ 239 & 85 \end{array} \right) = \left( \begin{array}{cc} 159287 & 56648 \\ 56648 & 20146 \end{array} \right), $$
which says that
$$ 284^2 + 2 \cdot 284 \cdot 239 - 239^2 = 159287$$
and applying
$$ \left( \begin{array}{cc} 1 & 1 \\ 0 & 1 \end{array} \right) \cdot \left( \begin{array}{c} 284 \\ 239 \end{array} \right) = \left( \begin{array}{c} 523 \\ 239 \end{array} \right)
$$
tells us that $$ 523^2 - 2 \cdot 239^2 = 159287. $$