Suppose I wanted to solve (for x and y) an equation of the form $x^2-dp^2y^2=1$ where d is squarefree and p is a prime. Of course I could simply generate the solutions to the Pell equation $x^2-dy^2=1$ and check if the value of y was divisible by $p^2,$ but that would be slow. Any better ideas?
It would be useful to be able to distinguish cases where the equation is solvable from cases where it is unsolvable, even without finding an explicit solution.
On
A different and important way to view this equation is still as a norm equation, only in a non-maximal order of the quadratic field - namely, we are looking at the ring $\mathbb{Z}[p\sqrt{d}]$.
Since we are still in the realm of algebra, we can prove results algebraically. For example:
Proposition. Let $p$ be a prime that doesn't divide $d$, and $e=x+y\sqrt{d}$ be a unit with integer $x,y$. Then $e^{p-(d/p)} \in \mathbb{Z}[p\sqrt{d}]$, where $(d/p)$ is the Legendre symbol.
Proof. If the Legendre symbol is $-1$, this means that $(p)$ doesn't split in $O$ (which will denote the maximal order from here on), so that $O/(p)$ is isomorphic to the finite field of size $p^2$. Now every number in $\mathbb{F}_{p^2}$ has order dividing $p^2-1$, so $e^{p+1} \pmod{p}$ has order diving $p-1$. But the numbers of order diving $p-1$ are exactly those in the subfield $\mathbb{F}_p$, and that means that $e^{p+1} = z+pw\sqrt{d}$, proving the claim in this case.
If $(p)$ splits into $\pi_1 \pi_2$, then $O/(p)$ is isomorphic to $O/\pi_1 \pi_2$ which is in turn, by chinese remainder theorem, isomorphic to $\mathbb{F}_p \oplus \mathbb{F}_p$. By inspection we see that the galois action permutes the two coordinates. So, since under the isomorphism $e^{p-1}\pmod{\pi_1}=e^{p-1}\pmod{\pi_2}=1$, we see that $e^{p-1}\pmod{p}$ is invariant. Hence it is of the form $z+pw\sqrt{d}$, proving the claim in this case as well. $\Box$
The above proposition shows how to get solutions to the non-squarefree Pell equation from solutions of the squarefree one - simply power appropriately.
On
For what it is worth, my lecture notes on the Pell Equation are in the context of $d$ any positive, nonsquare integer. Also see problem 9.2 here where I ask the students to exploit the fact that $d$ does not need to be squarefree to show that one can find solutions $(x,y)$ to the Pell equation satisfying an additional congruence $y \equiv 0 \pmod M$ for any $M \in \mathbb{Z}^+$.
From the perspective of algebraic number theory, this comes down to a version of the Dirichlet Unit Theorem for nonmaximal orders in a number field $K$. Recall that a $\mathbb{Z}$-order $R$ in $K$ is a subring of $\mathbb{Z}_K$ such that the additive group $(R,+)$ is finitely generated as a $\mathbb{Z}$-module. The last condition is equivalent to the finiteness of the index $[\mathbb{Z}_K:R]$. The usual statement of the Dirichlet Unit Theorem is that the unit group $\mathbb{Z}_K^{\times}$ is a finitely generated abelian group with rank equal to $r_1 + r_2 - 1$, where if $K \cong \mathbb{Q}[t]/(P(t))$, the polynomial $P$ has $r_1$ real roots and $r_2$ pairs of complex conjugate non-real roots.
But if I am not mistaken (and please let me know if I am!), the standard proof of the Dirichlet Unit Theorem works to show that exactly the same is true for the unit group $R^{\times}$ of any nonmaximal order. (Certainly $R^{\times}$ is finitely generated, being a subgroup of the finitely generated abelian group $\mathbb{Z}_K^{\times}$; the claim is that its rank is no less than that of $\mathbb{Z}_K^{\times}$.)
Using the structure theory of finitely generated abelian groups, one easily deduces the following relative version of the Dirichlet Unit Theorem: for any order $R$ in $\mathbb{Z}_K^{\times}$, the quotient group $\mathbb{Z}_K^{\times}/R^{\times}$ is finite.
William LeVeque's Topics in Number Theory, Volume I treats Pell's equation with nonsquare $d$ in Chapter 8, without the assumption that $d$ is squarefree.
In Section 8.2 (page 139 in the first half of the Dover edition), he considers $x^2-dy^2 = 1$ as follows:
If $d$ is negative, then when $d=-1$ the only solutions are $(\pm 1,0)$ and $(0,\pm 1)$. If $d\lt -1$, then the only solutions are $(\pm 1,0)$.
If $d$ is a square, $d=(d')^2$, then we can rewrite the equation as $x^2 - (d'y)^2 = 1$; the only squares that differ by $1$ are $0$ and $1$, so again the only solutions are $(\pm 1,0)$.
For positive nonsquare $d$, we start with a lemma on approximation and then proceed.
Lemma. If $\xi$ is a real number and $t$ is a positive integer, then there are integers $x$ and $y$ such that $$\left|\xi - \frac{x}{y}\right| \leq \frac{1}{y(t+1)},\qquad 1\leq y\leq t.$$
Proof. The $t+1$ numbers $$0\cdot \xi - \lfloor 0\cdot \xi\rfloor,\quad 1\cdot\xi - \lfloor 1\cdot \xi\rfloor, \quad\cdots\quad t\xi - \lfloor t\xi\rfloor$$ are all in the interval $0\leq u \lt 1$. In increasing order of magnitude, call them $a_0$, $a_1,\ldots,a_t$. Mark the numbers on a circle of perimeter $1$. Then the $t+1$ differences $$a_1-a_0,\quad a_2-a_1,\quad\ldots,\quad a_t-a_{t-1},\quad a_0-a_t+1,$$ are the lengths of the arcs between successive $a$s, so they are nonnegative, and $$(a_1-0) + (a_2-a_1) + \cdots + (1-a_t) = 1.$$ So at least one of these $t+1$ differences is no more than $\frac{1}{t+1}$. But the each difference is of the form $g_1\xi - g_2\xi - m$, with $m$ an integer, so take $y=|g_1-g_2|$, $x=\pm m$. $\Box$
Theorem. For any irrational $\xi$, the inequality $$|x - \xi y| \lt \frac{1}{y}$$ has infinitely many integer solutions.
Proof. For each positive integer $t$, $0\lt |x-\xi y|\lt \frac{1}{t}$, $1\leq y\leq t$ has solutions. Taking $t=1$ gives a solution $(x_1,y_1)$ to the original equation. For sufficiently large $t_1$ we have $|x_1 - \xi y_1|\gt \frac{1}{t_1}$, so taking $t=t_1$ gives a new solution $(x_2,y_2)$ to the original equation. Lather, rinse, repeat. $\Box$
The irrationality of $\xi$ is used to ensure that you have strict inequalities when you apply the lemma, so that you can set up the recursion and get infinitely many solutions.
Theorem. If $d$ is positive and not a square, then there are infinitely many integer solutions to the equation $$x^2 - dy^2 = k$$ for positive integers $x,y$ for some $k$ with $|k|\lt 1+2\sqrt{d}$.
Proof. Pick a solution to $|x - \sqrt{d}y|\lt \frac{1}{y}$. Then \begin{align*} |x+y\sqrt{d}| &= |x-y\sqrt{d} + 2y\sqrt{d}|\\ &\lt \frac{1}{y}+2y\sqrt{d}\\ &\leq (1+2\sqrt{d})y \end{align*} and so $|x^2-dy^2| \lt \frac{1}{y}(1+2\sqrt{d})y = 1+2\sqrt{d}$.
Since there are infinitely many pairs $(x,y)$ you can use, but only finitely many integers that are positive and smaller than $1+2\sqrt{d}$, infinitely many of the values of $x^2-dy^2$ must coincide, which is the theorem. $\Box$
Notice that the only requirement here is that $\sqrt{d}$ be real and irrational, i.e., that $d$ is positive and not a perfect square.
Theorem. If $d\gt 0$ is not a square, then the equation $$x^2 - dy^2 = 1$$ has at least one solution with $y\neq 0$.
Proof. Take the infinitely many solutions to $x^2-dy^2 = k$ (for some $k$ with $|k|\lt 1+2\sqrt{d}$) and divide them into $k^2$ equivalence classes, where $(x_1,y_1)\sim (x_2,y_2)$ if and only if $x_1\equiv x_2\pmod{k}$ and $y_1\equiv y_2\pmod{k}$. Some class contains more than one solution, say $(x_1,y_1)$ and $(x_2,y_2)$ with $x_1x_2\gt 0$. Let $$x = \frac{x_1x_2 - dy_1y_2}{k},\qquad y=\frac{x_1y_2 - x_2y_1}{k};$$ then $x$ and $y$ are integers with $y\neq 0$ and $x^2-dy^2 = 1$.
Indeed, $x_1y_2\equiv x_2y_1\pmod{k}$, so $y$ is an integer. Also, $$x_1x_2 - dy_1y_2 \equiv x_1^2 - dy_1^2 = k \equiv 0 \pmod{k}$$ so $x$ is an integer. Also \begin{align*} x^2 - dy^2 &= \frac{1}{k^2}\left( (x_1x_2 - dy_1y_2)^2 - d(x_1y_2-x_2y_1)^2\right)\\ &= \frac{1}{k^2}\left(x_1^2x_2^2 - dx_1^2y_2^2 + d^2y_1^2y_2^2 - dx_2^2y_1^2\right)\\ &= \frac{1}{k^2}\left(x_1^2-dy_1^2\right)\left(x_2^2 - dy_2^2\right) = 1. \end{align*} Finally, if $y=0$, then $x_1y_2=x_2y_1$, so $x_1=ax_2$ and $y_1=ay_2$ for some $a$; but plugging into $x^2- dy^2 = k$ we get $a=1$, contradicting that $(x_1,y_1)$ and $(x_2,y_2)$ are distinct solutions. $\Box$
Theorem. If $(x_1,y_1)$ and $(x_2,y_2)$ are solutions to $x^2-dy^2 = 1$, then so are the integers $x$ and $y$ defined by $$(x_1 + y_1\sqrt{d})(x_2+y_2\sqrt{d}) = x+y\sqrt{d}.$$