\documentclass[12pt,reqno]{article}

\usepackage[usenames]{color}
\usepackage{amssymb}
\usepackage{graphicx}
\usepackage{amscd}

\usepackage[colorlinks=true,
linkcolor=webgreen,
filecolor=webbrown,
citecolor=webgreen]{hyperref}

\definecolor{webgreen}{rgb}{0,.5,0}
\definecolor{webbrown}{rgb}{.6,0,0}

\usepackage{color}
\usepackage{fullpage}
\usepackage{float}

\usepackage{psfig}
\usepackage{graphics,amsmath,amssymb}
\usepackage{amsthm}
\usepackage{amsfonts}
\usepackage{latexsym}
\usepackage{epsf}

\setlength{\textwidth}{6.5in}
\setlength{\oddsidemargin}{.1in}
\setlength{\evensidemargin}{.1in}
\setlength{\topmargin}{-.1in}
\setlength{\textheight}{8.4in}

\newcommand{\seqnum}[1]{\href{http://oeis.org/#1}{\underline{#1}}}

\begin{document}

\begin{center}
\epsfxsize=4in
\leavevmode\epsffile{logo129.eps}
\end{center}

\theoremstyle{plain}
\newtheorem{theorem}{Theorem}
\newtheorem{corollary}[theorem]{Corollary}
\newtheorem{lemma}[theorem]{Lemma}
\newtheorem{proposition}[theorem]{Proposition}

\theoremstyle{definition}
\newtheorem{definition}[theorem]{Definition}
\newtheorem{example}[theorem]{Example}
\newtheorem{conjecture}[theorem]{Conjecture}

\theoremstyle{remark}
\newtheorem{remark}[theorem]{Remark}

\begin{center}
\vskip 1cm{\LARGE\bf 
Lagrange's Algorithm Revisited: \\
\vskip .02in
Solving $at^2+btu+cu^2=n$ in the  \\
\vskip .1in
Case of Negative Discriminant
}
\vskip 1cm
\large
Keith R. Matthews\\
School of Mathematics and Physics\\
University of Queensland\\
Brisbane, QLD 4072 \\
Australia \\
and\\
Centre for Mathematics and its Applications\\
Australian National University\\
Canberra, ACT 0200 \\
Australia\\
\href{mailto:keithmatt@gmail.com}{\tt keithmatt@gmail.com} \\
\end{center}

\vskip .2 in

\begin{abstract}
We make more accessible a neglected continued fraction algorithm of
Lagrange for solving the equation $at^2+btu+cu^2=n$ in relatively prime integers
$t,u$, where $a>0$, $\gcd(a,n)=1$, and $D=b^2-4ac<0$. The cases $D=-4$ and
$D=-3$ present a consecutive convergents phenomenon which aids the search
for solutions.
\end{abstract}

\bigskip

\section{Introduction}

\newcommand{\BZ}{\mbox{$\mathbb Z$}}
At the end of a memoir in 1770, Lagrange \cite[pp.\ 717--726]{lagrange} gave a method for finding the solutions of a positive definite binary form equation 
\begin{equation}
at^2+btu+cu^2=n,\label{eq:-1}
\end{equation}
where $\gcd(t,u)=1$, $\gcd(a,n)=1$, $b^2-4ac<0$, $a>0$, and $n>0$.  For such a solution, $\gcd(u,n)=1$ and hence the congruence $\theta{}u \equiv t$ (mod $n$) has a unique solution $\theta$ in the range $-n/2<\theta\leq n/2$. Then 
\begin{align}
at^2+btu+cu^2&\equiv 0\pmod{n}\notag\\
a(\theta{}u)^2+b(\theta{}u)u+cu^2&\equiv 0\pmod{n}\notag\\
a\theta^2+b\theta+c&\equiv 0\pmod{n}.\label{eq:congruence}
\end{align}
Lagrange \cite[p.\ 700]{lagrange} used the transformation 
\begin{equation}
t=\theta{}u-ny\label{eq:0}
\end{equation}
 to convert equation \eqref{eq:-1} to
\begin{equation}
Pu^2+Quy+Ry^2=1,\label{eq:01}
\end{equation}
where $P=(a\theta^2+b\theta+c)/n, Q=-(2a\theta+b), R=na$.

(We remark that, conversely, if $(u,y)$ is a solution of \eqref{eq:01}, then $(t,u)=(\theta{}u-ny,u)$  is a solution of \eqref{eq:-1} with $\gcd(t,u)=1$.)

We note that $D=b^2-4ac=Q^2-4PR$. Also if $(u,y)$ is a solution of \eqref{eq:01}, so is $(-u,-y)$.

Lagrange \cite{lagrange} proved in  Sections $20, 35$ and $39$ of his paper that $u/y, y>0$ is a convergent to $-Q/2P$. His proof was long and hard to follow. The aim of this paper is to give a short proof that Lagrange's assertion holds, apart from certain exceptional cases. If $D\neq -3$, this is done in Section \ref{section:Dneq3}, where we use the following standard test due to Lagrange
\cite[Satz $2.11$, p.\ 39]{perron}:
\begin{lemma}\label{lemma:lagrange}
If a rational $x/y$ with $\gcd(x,y)=1$ and $y\geq1$ has the property that $|\omega-x/y|<1/2y^2$, then $x/y$ is a convergent of the continued fraction expansion of $\omega$.
\end{lemma}
If $D=-3$, more care is needed. In Section \ref{section:Dequalsminus3}, we use the following criterion from \cite[Theorem $172$, p.\ 140]{hardy}:
\begin{lemma}\label{lemma:unimodular} If $\omega=\frac{P\zeta+R}{Q\zeta +S}$, 
where $\zeta>1$ and $P,Q,R,S$ are integers such that $Q>S>0$ and $PS-QR=\pm 1$, then $R/S=A_{n-1}/B_{n-1}$ and $P/Q=A_n/B_n$ are consecutive convergents of $\omega$.  Also $\zeta$ is the $(n+1)$th complete convergent of $\omega$.
\end{lemma}
(The author \cite{matthews} used this approach successfully in an earlier paper \cite{matthews} on Lagrange's work, when $D>0$.)

Lagrange gave solution bounds 
\begin{equation}
u\leq\sqrt{4R/(4PR-Q^2)}, \quad y\leq\sqrt{4P/(4PR-Q^2)},\label{lagrange_bounds}
\end{equation}
which are easy to derive by completion of the square in equation \eqref{eq:01}.

We note that in a series of papers, K. S. Williams \cite{williams} also considered congruence \eqref{eq:congruence}, but did not consider equation \eqref{eq:01} and instead examined the continued fraction of $\theta/n$, thereby generalizing a method of Hermite and Cornacchia.  The algorithm presented in Section \ref{section:algorithm} of the present paper is quite different and is also easy to program.

We use the continued fraction notation $[a_0,\ldots,a_n]$.

It is well-known (see \cite[Theorem $59$, p.\ 75]{dickson}) that equation \eqref{eq:01}, when soluble, has two solutions if $D<-4$, four solutions if $D=-4$ and six solutions if $D=-3$. 

It is easy to show that if $D=-4,Q=2N$ and $(u,y)$ is a solution of \eqref{eq:01}, then
$(-(Nu+Ry),Pu+Ny)$ is also a solution. Whereas if $D=-3, Q=2N+1$ and $(u,y)$ is a solution of \eqref{eq:01}, then $(-(Nu+Ry),Pu+(N+1)y)$ and $(-(u(N+1)+yR), Pu+Ny)$ are also solutions.  In sections \ref{section:Dequalsminus4} and \ref{section:Dequalsminus3}, if $D=-4$ or $-3$, it is shown that apart from certain exceptional cases, these solutions of \eqref{eq:01}, apart from sign, arise from consecutive convergents to $-Q/2P$. This important fact is used in the algorithm of Section \ref{section:algorithm} and was not mentioned by Lagrange.

The algorithm is available for online experimentation at \cite{bcmath}.
\begin{remark} The assumption that $\gcd(a,n)>1$ involves no loss of generality. For we can assume 
that $\gcd(a,b,c)=1$.  Then as pointed out by Gauss \cite[p.\ 221]{gauss} (also see \cite[Lemma 2, pp.\ 311--312]{hua}), there exist relatively prime integers 
$\alpha, \gamma$ such that $a\alpha^2+b\alpha\gamma+c\gamma^2=A$, where $\gcd(A,n)=1$. The construction uses the factorization of $n$. Then if $\alpha\delta-\beta\gamma=1$, the transformation 
$t=\alpha T+\beta U, u=\gamma T+\delta U$ converts $at^2+btu+cu^2$ to $AT^2+BTU+CU^2$, with the two forms representing the same integers.
\end{remark}
\section{Exceptional cases}\label{section:exceptional}
We first list some exceptional cases where the solutions $(u,y)$ of \eqref{eq:01} are easily found.
\begin{enumerate}
\item[(a)] $D<-4$ and $P=1$. Then the solutions are $(u,y)=\pm(1,0)$.
 \item[(b)] $D=-4$. Then $Q=2N$.
  \begin{enumerate}
      \item[(i)] If $P=1$, then $R=N^2+1$ and the solutions $(u,y)$ are $\pm(1,0)$ and $\pm(-N,1)$. Here $(-N,1)=(A_0,B_0)$.
      \item[(ii)] If $P=2$, then $R=(N^2+1)/2$, where $N$ is odd and the solutions $(u,y)$ are $\pm(\frac{(-N+1)}{2},1)$ and $\pm(\frac{-(N+1)}{2},1)$. Here $(-(N+1)/2,1)=(A_1,B_1)$.
   \end{enumerate}
 \item[(c)] $D=-3$. Then $Q=2N+1$.
  \begin{enumerate}
      \item[(i)] If $P=1$, then $R=N^2+N+1$ and the solutions $(u,y)$ are $\pm(1,0), \pm(-N,1)$ and $\pm(-(N+1),1)$. Here $(-(N+1),1)=(A_0,B_0)$.
\item[(ii)] If $P=3$, then $R=(N^2+N+1)/3$, where $N \equiv 1$ (mod ${3}$) and solutions $(u,y)$ are $\pm(\frac{(-N+1)}{3},1), \pm(\frac{-(2N+1)}{3},2)$ and $\pm(\frac{-(N+2)}{3},1)$.  Here $(\frac{-(2N+1)}{3},2)=(A_1,B_1)$ and $(\frac{-(N+2)}{3},1)=(A_0,B_0)$.
  \end{enumerate}
\end{enumerate}
From now on, we exclude these cases.
\section{The case \texorpdfstring{$D\neq -3$}{D not equal to minus 3}}\label{section:Dneq3}
\begin{theorem}\label{theorem:lagrange_lemma} Let $u$ and $y>0$ be integers satisfying \eqref{eq:01},
where $D=Q^2-4PR<0$ and $P, Q, R$ are integers, $P>0, D\neq-3$ and $P\neq 2$ if $D=-4$.  Then $u/y$ is a convergent to $\omega=-Q/2P$. 
\end{theorem}
\begin{proof} We derive the inequality
\begin{equation}
\left|\omega-\frac{u}{y}\right|<\frac{1}{2y^2}.\label{eq:2}
\end{equation}
Then Lemma~\ref{lemma:lagrange} shows that $u/y$ is a convergent to $\omega=-Q/2P$.

\bigskip

(a) Let $Q=2N$.  Then $\omega=-N/P$ and 
\begin{align}
\left|\omega-\frac{u}{y}\right|&=\left|-\frac{N}{P}-\frac{u}{y}\right|<\frac{1}{2y^2}\notag\\
                    &\iff |Pu+Ny|<\frac{P}{2y}.\label{eq:3}
\end{align}
From \eqref{eq:01}, with $\Delta=-D/4=PR-N^2$, we have
\begin{equation}
u=\frac{-Ny\pm\sqrt{P-\Delta{}y^2}}{P}\label{eq:4}
\end{equation}
and hence
\begin{equation}
Pu+Ny=\pm\sqrt{P-\Delta{}y^2}.
\end{equation}
Then \eqref{eq:3} becomes
\begin{equation}
\sqrt{P-\Delta{}y^2}<\frac{P}{2y},
\end{equation}
which reduces, on cross-multiplying, to
\begin{equation}
(P-2y^2)^2+4(\Delta-1)y^4>0.\label{eq:5}
\end{equation}
However \eqref{eq:5} holds if $\Delta>1$ or if $\Delta=1$ and $P\neq 2y^2$. But if $\Delta=1$ and $P=2y^2$, then $y=1$, as $\gcd(P,y)=1$. Hence $P=2$ and this case was excluded.

\bigskip

(b) Let $Q=2N+1$ and $\Delta=-D=4PR-(2N+1)^2>0$. Then
\begin{align}
\left|\omega-\frac{u}{y}\right|&=\left|-\frac{2N+1}{2P}-\frac{u}{y}\right|<\frac{1}{2y^2}\notag\\
                    &\iff |2Pu+(2N+1)y|<\frac{P}{y}.\label{eq:6}
\end{align}
From \eqref{eq:01}, we have
\begin{equation}
u=\frac{-(2N+1)y\pm\sqrt{4P-\Delta{}y^2}}{2P}\label{eq:7}
\end{equation}
and hence
\begin{equation}
2Pu+(2N+1)y=\pm\sqrt{4P-\Delta{}y^2}.\label{eq:70}
\end{equation}
Then \eqref{eq:6} becomes
\begin{equation}
\sqrt{4P-\Delta{}y^2}<\frac{P}{y},
\end{equation}
which reduces, on cross-multiplying, to
\begin{equation}
(P-2y^2)^2+(\Delta-4)y^4>0.\label{eq:8}
\end{equation}
This inequality holds, as $\Delta \equiv -1$ (mod $4$) and so $\Delta>4$ if $\Delta\neq 3$.
\end{proof}
\section{The case \texorpdfstring{$D=-4$}{D equal to -4}: Finer detail}\label{section:Dequalsminus4}
The next result has the useful computational aspect that once we find a convergent that gives a solution of \eqref{eq:01}, we know that the next convergent will also give a solution and complete the search for that value of $\theta$.
\begin{theorem}\label{theorem;1}  Let $D=-4, P\neq 1,2$ and $(u,y), y>0$ be a solution of \eqref{eq:01}. Then $Q=2N$ and $\frac{u}{y}$ and $\frac{-t(Nu+Ry)}{t(Pu+Ny)}$ are consecutive convergents to $-Q/2P$, where $t=sgn(Pu+Ny)$.
\end{theorem}
\begin{proof} 
We have the identity
\begin{equation}
\frac{-Q}{2P}=\frac{-N}{P}=\frac{u\xi-Nu-Ry}{y\xi+Pu+Ny}, \label{eq:9}
\end{equation}
where $\xi=\frac{y}{Pu+Ny}$. From equation \eqref{eq:4} with $\Delta=1$, we have
\begin{equation}
Pu+Ny=\pm\sqrt{P-y^2},\label{eq:10}
\end{equation}
where $P>y^2$. (We note that $P=y^2$ would imply $y=1=P$, as $\gcd(P,y)=1$ and this is excluded.) 

\bigskip

\noindent
{\bf Case} 1. Assume $Pu+Ny=\sqrt{P-y^2}$. Then $\xi=y/\sqrt{P-y^2}>0$.  Then $\xi\neq1$, as otherwise $\sqrt{P-y^2}=y, P=2y^2$ and $y=1, P=2$, as $\gcd(P,y)=1$; however this case was excluded.
   \begin{enumerate}
      \item[(i)] Assume $2y^2>P$.  Then $\xi>1$.  For
\begin{align*}
\xi>1&\iff y>\sqrt{P-y^2}\\
     &\iff 2y^2>P.
\end{align*} Also
\[
\frac{-N}{P}=\frac{u\xi+(-Nu-Ry)}{y\xi+(Pu+Ny)}.
\]
Then as $y>Pu+Ny>0$,  Lemma \ref{lemma:unimodular} implies that $\frac{u}{y}=\frac{A_m}{B_m}$ and $\frac{-Nu-Ry}{Pu+Ny}=\frac{A_{m-1}}{B_{m-1}}$ are consecutive convergents to $-N/P$.
      \item[(ii)] Assume $2y^2<P$. Then $0<\xi<1$.  Also
\[
\frac{-N}{P}=\frac{u-(Nu-Ry)(\xi^{-1})}{y+(Pu+Ny)(\xi^{-1})}.
\]
Then as $y<Pu+Ny$, Lemma \ref{lemma:unimodular} implies that $\frac{u}{y}=\frac{A_{m-1}}{B_{m-1}}$ and $\frac{-Nu-Ry}{Pu+Ny}=\frac{A_m}{B_m}$ are consecutive convergents to $-N/P$.
   \end{enumerate}

\bigskip
\noindent
{\bf Case} 2. Assume $Pu+Ny=-\sqrt{P-y^2}$. Then $\xi=y/(-\sqrt{P-y^2})<0$. We cannot have $\xi=-1$ as otherwise $P=2$.
   \begin{enumerate}
      \item[(i)] Assume $2y^2>P$.  Then $|\xi|>1$ and hence $-\xi>1$. Also
\[
\frac{-N}{P}=\frac{u(-\xi)+Nu+Ry}{y(-\xi)-(Pu+Ny)},
\]
and $y>-(Pu+Ny)>0$.  Hence $\frac{u}{y}=\frac{A_m}{B_m}$ and $\frac{Nu+Ry}{-(Pu+Ny)}=\frac{A_{m-1}}{B_{m-1}}$ are consecutive convergents to $-N/P$.
      \item[(ii)] Assume $2y^2<P$. Then $|\xi|<1$ and hence $-\xi^{-1}>1$. Also
\[
\frac{-N}{P}=\frac{u+(Nu+Ry)(-\xi^{-1})}{y-(Pu+Ny)(-\xi^{-1})},
\]
where $y<-(Pu+Ny)$.  Hence $\frac{u}{y}=\frac{A_{m-1}}{B_{m-1}}$ and $\frac{Nu+Ry}{-(Pu+Ny)}=\frac{A_m}{B_m}$ are consecutive convergents to $-N/P$.
   \end{enumerate}
\end{proof}

\section{The case \texorpdfstring{$D=-3$}{D equal to -3}}\label{section:Dequalsminus3}
The case $D=-3$ was excluded from Theorem \ref{theorem:lagrange_lemma} and
we discuss it now.

The next result has the useful computational aspect that once we find a convergent that gives a solution of \eqref{eq:01}, we know that the next two convergents will give two further solutions and complete the search for that value of $\theta$.
\begin{theorem}\label{theorem:1}  Let $D=-3, P\neq 1,3$ and $(u,y), y>0$ be a solution of\eqref{eq:01}.
Then $Q=2N+1$ and the rational numbers
\[
\frac{u}{y}, \quad\frac{-s(Nu+Ry)}{s(Pu+(N+1)y)}, \quad\frac{t(u(N+1)+yR)}{t(Pu+Ny)}
\]
 are consecutive convergents in some order to $-Q/2P$, where $s=sgn(Pu+(N+1)y$ and $t=sgn(Pu+Ny)$.
\end{theorem}
\begin{proof} We have the identity
\begin{equation}
\frac{-Q}{2P}=\frac{-(2N+1)}{2P}=\frac{u\xi-(N+1)u-Ry}{y\xi+Pu+Ny}, \label{eq:11}
\end{equation}
where $\xi=\frac{Pu+(N+2)y}{2Pu+(2N+1)y}$. From equation \eqref{eq:70} with $\Delta=3$, we have
\begin{equation}
2Pu+(2N+1)y=\pm\sqrt{4P-3y^2},\label{eq:12}
\end{equation}
where $4P>3y^2$.  (We have $4P-3y^2\neq0$, as otherwise $4P=3y^2$ and hence $y=2, P=3$, which was excluded.)

\bigskip

\noindent
{\bf Case} 1. Assume $2Pu+(2N+1)y=\sqrt{4P-3y^2}$.  Then 
\[
Pu+(N+2)y=\frac{\sqrt{4P-3y^2}+3y}{2}>0
\]
and hence $\xi>0$.  Also $\xi\neq1$.   For
\begin{align*}
\xi=1&\implies\frac{\sqrt{4P-3y^2}+3y}{2}=\sqrt{4P-3y^2}\\
     &\implies 3y=\sqrt{4P-3y^2}\\
     &\implies 3y^2=P\\
     &\implies y=1, P=3,
\end{align*}
which was excluded.

We note that $Pu+Ny>0\iff\sqrt{4P-3y^2}>y\iff P>y^2$.
    \begin{enumerate}
        \item[(i)] Assume $3y^2<P$.  Then $0<\xi<1$. For
\begin{align*}
\xi<1&\iff Pu+(N+2)y<2Pu+(2N+1)y\\
     &\iff y<Pu+Ny\\
     &\iff y<\frac{\sqrt{4P-3y^2}-y}{2}\\
     &\iff 3y<\sqrt{4P-3y^2}\\
     &\iff 3y^2<P.
\end{align*}
 Then  \eqref{eq:11} gives
\begin{equation}
\frac{-Q}{2P}=\frac{-(2N+1)}{2P}=\frac{u-((N+1)u+Ry)\xi^{-1}}{y+(Pu+Ny)\xi^{-1}}. \label{eq:13}
\end{equation}
Also we have $y<Pu+Ny$ and $\xi^{-1}>1$.  Then Lemma \ref{lemma:unimodular} applied to \eqref{eq:13} implies that $\frac{u}{y}=\frac{A_{m-1}}{B_{m-1}}$ and $\frac{-(N+1)u-Ry}{Pu+Ny}=\frac{A_m}{B_m}$ are consecutive convergents of $-Q/2P$.
    \item[(ii)] Assume $3y^2>P>y^2$.  Then we have $\xi>1, y>Pu+Ny>0$ and by Lemma \ref{lemma:unimodular} applied to equation \eqref{eq:13}, it follows that $\frac{u}{y}=\frac{A_r}{B_{r}}$ and $\frac{-(N+1)u-Ry}{Pu+Ny}=\frac{A_{r-1}}{B_{r-1}}$ are consecutive convergents to $-Q/2P$.
    \item[(iii)] Assume $P<y^2$.  Then we have $\xi>2$. For
\begin{align*}
\xi>2&\iff Pu+(N+2)y>2(2Pu+(2N+1)y)\\
     &\iff 0>Pu+Ny\\
     &\iff P<y^2.
\end{align*}
We rewrite equation \eqref{eq:11} as
\begin{equation}
\frac{-(2N+1)}{2P}=\frac{u(\xi-1)-(Nu+Ry)}{y(\xi-1)+Pu+(N+1)y}. \label{eq:14}
\end{equation}
Then we have $\xi-1>1, y>Pu+(N+1)y>0$ and by Lemma \ref{lemma:unimodular} applied to equation \eqref{eq:14}, it follows that $\frac{u}{y}=\frac{A_s}{B_{s}}$ and $\frac{-Nu-Ry}{Pu+(N+1)y}=\frac{A_{s-1}}{B_{s-1}}$ are consecutive convergents to $-Q/2P$.
    \end{enumerate}

We now link up each pair of solution convergents found in Cases 1(i)--(iii) with a third solution convergent.  We start by employing the equations
\begin{align}
\xi_1&=\frac{-Pu-(N-1)y}{2Pu+(2N+1)y},\label{eq:15}\\
\frac{-(2N+1)}{2P}&=\frac{u\xi_1-Nu-Ry}{y\xi_1+Pu+(N+1)y}.\label{eq:16}
\end{align}
We find a pair of convergents which we list corresponding to Cases 1(i) and (ii):
\begin{align*}
{\rm(i)}\quad\frac{u}{y}=\frac{A_{m-1}}{B_{m-1}},\ \frac{-(N+1)u-Ry}{Pu+Ny}=\frac{A_m}{B_m} &\ \frac{-Nu-Ry}{Pu+(N+1)y}=\frac{A_{m+1}}{B_{m+1}},\\
{\rm(ii)}\quad\frac{-(N+1)u-Ry}{Pu+Ny}=\frac{A_{r-1}}{B_{r-1}},\ \frac{u}{y}=\frac{A_r}{B_r},&\ \frac{-Nu-Ry}{Pu+(N+1)y}=\frac{A_{r+1}}{B_{r+1}}.
\end{align*}
For Case 1(iii), we employ the equations
\begin{align}
\xi_2&=\frac{-Pu-(N-1)y}{Pu+(N+2)y},\label{eq:17}\\
\frac{-(2N+1)}{2P}&=\frac{((N+1)u+Ry)\xi_2-Nu-Ry}{-(Pu+Ny)\xi_2+Pu+(N+1)y}.\label{eq:18}
\end{align}
We then find a pair of convergents that is listed with the pair found in Case 1(iii):
\[
{\rm(iii)}\quad\frac{(N+1)u+Ry}{-(Pu+Ny)}=\frac{A_{s-1}}{B_{s-1}} \mbox { and }\frac{-Nu-Ry}{Pu+(N+1)y}=\frac{A_s}{B_s},\ \frac{u}{y}=\frac{A_{s+1}}{B_{s+1}}.
\]
This finishes Case 1.

\bigskip
\noindent
{\bf Case} 2. Assume $2(Pu+Ny)+y=-\sqrt{4P-3y^2}$.  Summarising, we find after tedious calculation, the following three results:
\begin{enumerate}
\item[(i)] $P>3y^2$: 
\[
\frac{u}{y}=\frac{A_{m-1}}{B_{m-1}},\ \frac{Nu+Ry}{-Pu-(N+1)y}=\frac{A_m}{B_m},\ \frac{(N+1)u+Ry}{-Pu-Ny}=\frac{A_m}{B_{m+1}};
\]
\item[(ii)] $3y^2>P>y^2$: 
\[
\frac{Nu+Ry}{-Pu-(N+1)y}=\frac{A_{m-1}}{B_{m-1}},\ 
\frac{u}{y}=\frac{A_m}{B_m},\ 
\frac{(N+1)u+Ry}{-Pu-Ny}=\frac{A_{m+1}}{B_{m+1}};
\]
\item[(iii)] $y^2>P$: 
\[
\frac{-Nu-Ry}{Pu+(N+1)y}=\frac{A_{m-1}}{B_{m-1}},\ 
\frac{(N+1)u+Ry}{-Pu-Ny}=\frac{A_m}{B_m},\ 
\frac{u}{y}=\frac{A_{m+1}}{B_{m+1}}. 
\]
\end{enumerate}
\end{proof}
\section{The continued fraction based algorithm}\label{section:algorithm}

The following algorithm finds all solutions $(t,u)$ with $\gcd(t,u)=1$,  of equation \eqref{eq:-1}, when $\gcd(a,n)=1$. 
The most time-consuming part involves solving the quadratic congruence \eqref{eq:congruence}, as this depends on finding the prime factorization of $n$. This is also a feature of Williams' algorithm, as well as Gauss' algorithm in \cite[p.\ 75]{dickson}.  This dependence means that for practial purposes, $n$ is restricted to less than $200$ digits. The present algorithm, like that of Williams, 
also uses Euclid's algorithm, whereas Gauss' algorithm uses reduction of positive definite forms, 
and these have similar running times. 

\bigskip

\noindent {\bf Input}: Integers $a,b,c,n, b^2-4ac<0, n>0, \gcd(a,n)=1$.

\noindent {\bf Output}: All solutions, if any, of $at^2+btu+cu^2=n, \gcd(t,u)=1$.

\noindent Solve $a\theta^2+b\theta+c \equiv 0$ (mod $n$), $-n/2 < \theta \leq n/2$.

\noindent {\bf If} there are no solutions, {\bf exit}.

\noindent Let $\theta_0,\ldots,\theta_{s-1}$ be the congruence solutions in the range $(-n/2,n/2]$.

\noindent $D :=b^2-4ac$.

\noindent {\bf for} $k=0,\ldots,s-1$ {\bf do}

\noindent \hspace{.5cm} $P: =(a\theta_k^2+b\theta_k+c)/n, Q :=2a\theta_k+b$;

\noindent \hspace{.5cm} {\bf if} $D<-4$ {\bf and} $P=1$ {\bf then} $(u,y) :=\pm(1,0)$;

\noindent \hspace{.5cm} {\bf if} $D = -4$, {\bf then} $N :=Q/2$; 

\noindent \hspace{1cm}    {\bf if} $P=1$, {\bf then} $(u,y) :=\pm(1,0), \pm(-N,1)$;

\noindent \hspace{1cm}    {\bf if} $P=2$, {\bf then} $(u,y) :=\pm(-N+1)/2,1), \pm(-(N+1)/2,1)$.

\noindent \hspace{.5cm} {\bf if} $D=-3$, {\bf then} $N :=(Q-1)/2$.

\noindent \hspace{1cm}    {\bf if} $P=1$, {\bf then} $(u,y) :=\pm(1,0), \pm(-N,1), \pm(-(N+1),1)$.

\noindent \hspace{1cm}    {\bf if} $P=3$, {\bf then} 

\noindent \hspace{1.5cm}     $(u,y) :=\pm((-N+1)/3,1), \pm(-(2N+1)/3,2), \pm(-(N+2)/3,1)$.

\noindent \hspace{.5cm} {\bf print} exceptional solutions $(t,u) :=(\theta_ku-ny,u)$;

\noindent \hspace{.5cm} {\bf continue} to next $k$;

\noindent \hspace{.5cm} $i :=0$;

\noindent \hspace{.5cm} bound := $\sqrt{4P/(-D)}$;

\noindent \hspace{.5cm} calculate convergent $A_0/B_0$ of $-Q/2P$;

\noindent \hspace{.5cm} {\bf while} $(B_i\leq {\rm bound})$ {\bf do}

\noindent \hspace{1.5cm}    {\bf if} $aA_i^2+bA_iB_i+cB_i^2=1$ {\bf then}

\noindent \hspace{2cm}      {\bf print} solutions $(t,u) :=\pm(\theta_kA_i-nB_i,A_i)$;
%If no $A_i/B_i$ satisfies $f_i=bA_i^2+cA_iB_i+dB_i^2=1$, {\bf continue} to next $k$.

\noindent \hspace{2cm}       {\bf if} $D<-4$, {\bf continue} to next $k$;

\noindent \hspace{2cm}       {\bf if} $D = -4$ {\bf or} $-3$ {\bf then}

\noindent \hspace{2.5cm}           calculate convergent $A_{i+1}/B_{i+1}$ of $-Q/2P$;

\noindent \hspace{2.5cm}           {\bf print} solutions $(t,u) :=\pm(\theta_kA_{i+1}-nB_{i+1},A_{i+1})$; 

\noindent \hspace{2.5cm}           {\bf if} $D = -4$, {\bf continue} to next $k$;

\noindent \hspace{2.5cm}           {\bf if} $D = -3$, {\bf then} 

\noindent \hspace{3cm}                calculate convergent $A_{i+2}/B_{i+2}$ of $-Q/2P$;

\noindent \hspace{3cm}                {\bf print} solutions $(t,u) :=\pm(\theta_kA_{i+2}-nB_{i+2},A_{i+2})$;

\noindent \hspace{3cm}                {\bf continue} to next $k$;
 
\noindent \hspace{1.5cm}  $i :=i+1$;

\noindent \hspace{1.5cm}  calculate convergent $A_i/B_i$ of $-Q/2P$;

\noindent \hspace{.5cm} {\bf end while} loop;

\noindent {\bf end for} loop.


\begin{example}
Find all solutions of $7t^2-9tu+3u^2=19, \gcd(t,u)=1$.   Here $D=-3$.

The congruence $7\theta^2-9\theta+3 \equiv 0$ (mod $19$) has solutions
$\theta = -1$ and $5$.

\bigskip

\noindent $\theta=-1$: The transformation $t=-u-19y$ converts $7t^2-9tu+3u^2=19$ to $u^2+23uy+133y^2=1$.  
This is one of the exceptional cases from Section \ref{section:exceptional}, with solutions $(u,y)=(\pm1,0), \pm(-11,1), \pm(-12,1)$. 

These produce primitive solutions $(t,u) = \pm(1,-1), \pm(8,-11), \pm(7,-12)$.

\bigskip

\noindent $\theta=5$:  The transformation $t=5u-19y$ converts $7t^2-9tu+3u^2=19$ to $7u^2-61uy+133y^2=1$.

 Then $-Q/2P=61/14=[4,2,1,4]$ and we find that convergents $A_0/B_0=4/1, A_1/B_1=9/2, A_2/B_2=13/3$ give solutions $(u,y)$.  These in turn give $(t,u)=(1,4), (7,9), (8,13)$. Hence we get primitive solutions $\pm(1,4), \pm(7,9), \pm(8,13)$ of $7t^2-9tu+3u^2=19$.

Hence the equation $7t^2-9tu+3u^2=19$ has $12$ primitive solutions.
\end{example}

\section{Acknowledgment}
I am grateful to the referee whose comments improved the presentation of the manuscript.


\begin{thebibliography}{9}

\bibitem{dickson}{L. E. Dickson},
\textit{Introduction to the Theory of Numbers}, 
Dover, 1957.

\bibitem{gauss} C. F. Gauss,
\textit{Disquisitiones Arithmeticae},
Yale University Press, 1966.
 
\bibitem{hardy}{G. H. Hardy and E. M. Wright},
\textit{An Introduction to the Theory of Numbers}, 
Oxford University Press, 1960.

\bibitem{hua} L. K. Hua,
\textit{Introduction to Number Theory},
Springer, 1982.

\bibitem{matthews}{K. R. Matthews},
The Diophantine equation $ax^2+bxy+cy^2=N, D=b^2-4ac>0$,
\textit{J. Th\'{e}or. Nombres Bordeaux} {\bf14} (2002) 257--270.

\bibitem{bcmath}{K. R. Matthews},
Finding primitive solutions of the diophantine equation 
$ax^2+bxy+cy^2=n$, where $b^2-4ac < 0$, $\gcd(a,b,c)=1$, and $a > 0$, available
at \url{http://www.numbertheory.org/php/posrep2.html}.

\bibitem{perron}{O. Perron},
\textit{Die Lehre von den Kettenbr\"{u}chen},
Band 1, B. G. Teubner, 1954.

\bibitem{lagrange} {J. A. Serret ed.},
\textit{Oeuvres de Lagrange},
Gauthiers-Villars, 1877.

\bibitem{williams}{K. S. Williams},
On finding the solutions of $n=au^2+buv+cv^2$ in integers $u$ and $v$, 
\textit{Util. Math}. {\bf 46} (1994), 3--19.
\end{thebibliography}

\bigskip
\hrule
\bigskip

\noindent 2010 {\it Mathematics Subject Classification}:
Primary 11A55; Secondary 11J70, 11D09.

\noindent \emph{Keywords: } diophantine equation, positive definite binary quadratic form, continued fraction, convergent.

\bigskip
\hrule
\bigskip

\vspace*{+.1in}
\noindent
Received August 26 2014;
revised version received September 22 2014.  
Published in {\it Journal of Integer Sequences}, November 5 2014.

\bigskip
\hrule
\bigskip

\noindent
Return to
\htmladdnormallink{Journal of Integer Sequences home page}{http://www.cs.uwaterloo.ca/journals/JIS/}.
\vskip .1in


\end{document}

                                                                                


.
