Is there a short proof that the set of solutions to Pell's equation form an abelian group? I'm working on a question from a past exam paper,

Show that the set of solutions $G_{d,p}$ to Pell's equation $x^2-dy^2=1$ modulo $p$ is a finite abelian group, and compute the order of this group for $p=5$ and all $d$.

Just proving associativity took 6 page-width lines of fairly meticulous working. This seems like a lot for 6/100 marks. Is there a quicker way than going through the group axioms individually?
 A: First, lets see how we can do this working over the integers instead of over $\mathbb Z/(p)$.  Consider the numbers of the form $a+b\sqrt{d}$ with $a,b\in \mathbb Z$, and define the norm $N(a+b\sqrt{d})=(a+b\sqrt{d})(a-b\sqrt{d})=a^2-db^2$.  We observe that this comes from using the automorphism of $\mathbb Z[\sqrt{d}]$ sending $\sqrt{d}$ to $-\sqrt{d}$.  Show that on numbers of this form $N(xy)=N(x)N(y)$.  Then the numbers with norm 1 are closed under multiplication, and you can construct the multiplicative inverse in a straight forward way.  Thus, the set of solutions to the Pell equation is in one to one correspondence with numbers of that form, but since numbers of that form are a subgroup of the real numbers under multiplication, they are an abelian group, and we can transport the group structure using the bijection.
Assuming that $d$ is not a perfect square mod $p$, we can use this same argument.  We would have that $\mathbb Z/(p)[\sqrt{d}]\cong \mathbb F_{p^2}$, and we have a field automorphism sending $\sqrt{d}$ to $-\sqrt{d}$.  Then the argument above works essentially unchanged.  We can compute the size of the group by noting that inside the group of units of the ring, which is of size $p^2-1$, we will have $p-1$ cosets of our group, corresponding to solutions of $x^2-dy^2=k$ where $k=1,2,\ldots, p-1$.  Thus, the cosets have size $p+1$.

If $d$ is a perfect square mod $p$, then we can still use the same formula to define the group structure. If we trace through the definitions we had before, we see that when we were working over the integers we had a binary operation defined by $(a,b)*(x,y)=(ax+byd,ay+bx)$. We can still define $N(a,b)=a^2-db^2$.  While this had an intuitive meaning over the integers, we can ignore this meaning and just remember that the multiplication is associative, and it satisfies $N((a,b)*(x,y))=N(a,b)N(x,y)$.  We can then reduce everything mod $p$ and everything is still true.  The ordered pairs with norm $1$ form a group, as do the ordered pairs with non-zero norm (mod p).  Just as before, we have that the quotient of the larger group by the subgroup is the equivalence classes of pairs with the same norm.
Since there are $p-1$ equivalence classes, we can find the size of the two groups if we can determine how many pairs satisfy $x^2-dy^2=0 \pmod p$.  If $y=0$, the only solution is $x=0$.  If $y\neq 0$, then $(x/y)^2=d \pmod p$, and so for every value of $y$ there are two choices for $x$ correseponding to the two square roots of $d$.  Therefore, there are $2p-1$ ordered pairs satisfying $x^2-dy^2=0\pmod p$, so there are $p^2-(2p-1)=(p-1)^2$ ordered pairs satisfying $x^2-dy^2\neq 0 \pmod p$, and dividing by $p-1$, we get there are $p-1$ ordered pairs satisfying $x^2-dy^2=k \pmod p$ for each $k\neq 0$.

This approach does not reveal the structure of the group as the answer of zhoraster does, but I like its simplicity.  It also explains the different answers in the size of the group based on whether $d$ is a perfect square mod $p$.  If $d$ is a square, we saw that there are $1+2(p-1)$ solutions to $x^2-dy^2=0 \pmod p$, but if $d$ is not a square, then $(0,0)$ is the only solution, so in one case we get $\frac{p^2-(2p-1)}{p-1}$ and in the other case we get $\frac{p^2-1}{p-1}$.
A: Let $A_p$ be $\mathbb Z_p$ if $d$ is a square modulo $p>2$ and the quadratic extension $\mathbb Z_p(\sqrt{d})$ otherwise; in both cases let $\alpha = \sqrt{d}$.
Observe that $(x,y)\in G_{d,p}$ iff $(x+\alpha y)(x-\alpha y) = 1$. Now consider the map from $G_{d,p}$ to $(A_p^*)^2$ given by $\psi\colon (a,y)\mapsto (x+\alpha y,x-\alpha y)$. It is easy to see that it sends the operation $*$ to multiplication in $(A_p^*)^2$. Moreover, it is easily seen to be injective. Therefore, $*$ is a group operation and $\psi$ is a monomorphism.
@Aaron's post covers the case where $d$ is a non-square. Now if $d$ is a square, I claim that any $x\in \mathbb Z_p^*$ can be uniquely represented as $a+\alpha b$ with $(a,b)\in G_{d,p}$. Indeed, to get the representation, we just need to solve the system $a+\alpha b = x$, $a-\alpha b = x^{-1}$. Therefore, the image of $G_{g,p}$ under $\psi$ is just $\{(x,x^{-1}), x\in \mathbb Z_p^*\}$ with usual multiplication, which is isomorphic to $\mathbb Z_p^*$.
