The number of integer points on the curve $(7x-1)^2+(7y-1)^2=n$ The number of integral solutions to the equation
$$x^2+y^2=n$$
is defined to be $r_2(n)$ and if $n=2^ap_1^{a_1}\dots p_k^{a_k}q^2$ where $p_i\equiv 1\mod 4$ and $q$ is the product of primes which are $3\mod 4$, then
$$r_2(n)=4(a_1+1)\dots(a_k+1).$$
If we restrict $x\equiv y\equiv 1\mod 2$, that is search for solutions to
$$(2x-1)^2+(2y-1)^2=n$$
we find that if $n\equiv 2$ all solutions to the previous one are also solutions to this and otherwise none are. So the number of solutions is $r_2(n)$ or $0$ according to whether $n\equiv 2$ or not.
And if we restrict $x\equiv y\equiv 1\mod 4$, we find that each quadruple of solutions $(\pm x)^2+(\pm y)^2=n$ generates exactly one solution to the new equation. Therefore the number of solutions to $(4x-1)^2+(4y-1)^2=n$ is $\frac{r_2(n)}4$ if $n\equiv 2\mod 4$ and $0$ otherwise.
In general I want to know how many solutions there are for the equation
\begin{equation}\tag{1}
(mx-a)^2+(my-b)^2=n.
\end{equation}
The idea above solves the case $m=2,3,4,5,6$ for all values of $a,b$ but for $m=7$ it doesn't work anymore. This is because $7$ is the first integer for which something like
$$0^1+1^2\equiv 2^2+2^2\mod 7$$
happens. To be precise, $m=7$ is the least so that there are $a,b,c,d$ so that
$$a^2+b^2\equiv c^2+d^2\mod m$$
$$\{a,b\}\neq \{\pm c,\pm d\}$$
for all choices of signs. This means that the number of solutions when $m=7$ is not what you expect. For example the equation $(7x-1)^2+(7y-1)^2=9$ has no solutions, but I expected there to be$\frac{r_2(9)}4=1$. Strangely, the equation $(7x-1)^2+(7y-2)^2=n$ has exactly half as many solutions as I expected for all $n$ up to ten thousand. Since $1^2+2^2=5$, I reasoned that $5$ would behave differently from the other primes, so I excluded its exponent from the product of $r_2(n)$. Therefore I conjecture that the number of solutions to
$$(7x-1)^2+(7y-2)^2=n$$
where $n\equiv 5\mod 7$ and $n=2^a5^bp_1^{a_1}\dots p_k^{a_k}q^2$ is $\frac 12(a_1+1)\dots (a_k+1)=\frac 18r_2\left(\frac n{5^b}\right)$. This product is only odd if $n$ is a square, which is impossible so that doesn't come up.
What I want to know is how many solutions to (1) there are, and in particular the special cases above.
 A: Comment: We may use following equation as a particular case:
$$(a-1)^2+(a+1)^2=b^2+1\space\space\space\space\space\space(1)$$
This equation can have infinitely many solution; if a and b satisfy this equation then considering following identity we can have subsequent solutions:
$$(2b+3a-1)^2+(2b+3a+1)^2=(3b+4a)^2+1\space\space\space\space\space\space(2)$$
Now we want to solve equation:
$(mx-1)^2+(my-1)^2=n\space\space\space\space\space\space(3)$
Comparing this equation with (1) we construct following system of equation:
$\begin {cases}mx-1=a-1\Rightarrow a=mx\\my-1=a+1\Rightarrow a=my-2\\ n=b^2+1\end{cases}$
For example $(a, b)=(2, 3)$ satisfies equation (1) and we have:
$(2-1)^2+(2+1)^2=3^2+1=10=2\times 5; 5\equiv 1\bmod 4$
Now we use $(a, b)=(2, 3)$ in identity (2):
$a_1=2\times 3++3\times 2=12$
$b_1=3\times 3+4\times 2=17$
Which gives:
$(12-1)^2+(12+1)^2=17^2+1=290$
Now we can rewrite :
$12-1=2\times 6-1$
$13=2\times 7-1$
Which gives $(x, y, m, n)=(6, 7, 2, 290)$
OR:
$12-1=3\times 5-4$
$13=3\times 6-5$
Which for equation $(mx-c)^2+(my-d)^2=n$ gives:
$(x, y, m, c, d, n)=(5, 6, 3, 4, 5, 290)$
We can continue this and find more solution, one more example:
$a_2=2\times 17+3\times 12=70$
$b_2=3\times 17+4\times 12=99$
$(70-1)^2+(70+1)^2=99^2+1=9802=2\cdot13^2\cdot29$
$13\equiv 1\bmod 4$ also $29\equiv 1\bmod 4$
$70-1=69=7\times 10-1$
$70+1=72=7\times 11-5$
$(x, y, m, c, d, n)=(10, 11, 7, 1, 5, 9802)$
I do not see restrictions for x, y, m and n and limit for number of solutions.
A: For a prime $p\equiv 1\mod 4$, we have a unique (up to sign) representation as a sum of squares $p=a^2+b^2$ and with the identity
$$p(x^2+y^2)=(ax-by)^2+(bx+ay)^2=(ax+by)^2+(-bx+ay)^2$$
we can define the linear transformations
$$\begin{bmatrix}
x\\
y\end{bmatrix}\mapsto\begin{bmatrix}
a & -b\\
b & a\end{bmatrix}\begin{bmatrix}
x\\
y\end{bmatrix}$$
$$\begin{bmatrix}
x\\
y\end{bmatrix}\mapsto\begin{bmatrix}
a & b\\
-b & a\end{bmatrix}\begin{bmatrix}
x\\
y\end{bmatrix}.$$
Call those matrices $M_+$ and $M_-$. Matrices of this form commute and their product is also of this form. We can classify all primes by the residues of $a,b\mod 7$. There are 9 possibilities up to sign.
\begin{array}{|c|c|}
\hline
p\mod 7& (a,b) \\ \hline
1 & (0,1),(2,2)\\ \hline
2 & (0,3),(1,1) \\ \hline
3 & (1,3) \\ \hline
4 & (0,2),(3,3) \\ \hline
5 & (1,2) \\ \hline
6 & (2,3) \\ \hline
\end{array}
So for an $n=2^ap_1^{a_1}\dots p_k^{a_k}q^2$ write each $p_i=a_i^2+b_i^2$ and define matrices $M_{i+}, M_{i-}$ as above. Each of the representations of $n$ as a sum of squares corresponds to a product
$$(M_{1+}^{t_1}M_{1-}^{a_1-t_1})(M_{2+}^{t_1}M_{2-}^{a_2-t_2})\dots (M_{k+}^{t_k}M_{k-}^{a_k-t_k})=\begin{bmatrix}
\pm x & \mp y\\
\pm y & \pm x\end{bmatrix}$$
where $0\leq t_1\leq a_1, \dots 0\leq t_k\leq a_k$. Notice that there are $4$ choices of sign and $(a_1+1)\dots(a_k+1)$ choices of $t_i$, precisely $r_2(n)$ choices in total.
This means that to find how many solutions there are to $x^2+y^2=n$ subject to $x\equiv\alpha\mod 7$, $y\equiv\beta$ we just need to classify the prime factors of $n$ into one of the nine classes and figure out how many ways there are to choose the $t_i$ such that
$$(M_{1+}^{t_1}M_{1-}^{a_1-t_1})(M_{2+}^{t_1}M_{2-}^{a_2-t_2})\dots (M_{k+}^{t_k}M_{k-}^{a_k-t_k})\equiv\begin{bmatrix}\alpha&-\beta\\ \beta&\alpha\end{bmatrix}\mod 7.$$
There are only $2\cdot 9=18$ types of matrices in this product, and they have finite order, and they commute, so it must be possible to calculate how many solutions there are as a function of the $a_i$ alone.
In other words if
$$n=2^a(p_{11}^{a_{11}}p_{12}^{a_{12}}\dots p_{1k_1}^{a_{1k_1}})(p_{21}^{a_{21}}p_{22}^{a_{22}}\dots p_{2k_1}^{a_{2k_2}})\dots (p_{91}^{a_{91}}p_{92}^{a_{92}}\dots p_{9k_1}^{a_{9k_9}})q^2$$
where each $p_{ij}$ has a representation as a sum of squares of the type $i$ modulo 7, then the number of solutions to $(7x+\alpha)^2+(7y+\beta)^2=n$ can be calculated as a function of the $a_{ij}$ alone and does not depend on the primes themselves. This is kind of the answer I wanted for the question, but I wish it was less computationally intensive.
