As one question in a somewhat bigger analysis I want to characterize the set of solutions of the generalized Pell-equation in the title: $$a^2+b^2 = 2 c^2 \tag 1$$ I'm not much fluent with the Pell-solving technology, and from some initial brute-force solutions I constructed the following scheme, which seems to me as if it could produce all solutions by changing a vector of parameters.
My question is, whether the set of solutions that I find is exhaustive or not - tests in small numbers seem to suggest that.

I introduce two matrices; $$T= \begin{bmatrix} 1&0&0\\0&3&2\\0&4&3 \end{bmatrix} \tag 2$$ $$C= \begin{bmatrix} 0&1&0\\1&0&0\\0&0&1 \end{bmatrix} \tag 3$$ The starting solution is the vector $A_0=[a,b,c]=[1,1,1]$.

With this I create further solutions by $$A_k=A_0 \cdot T^k \qquad k \in \mathbb Z \tag {4.1}$$ which gives a list with all solution having $a=1$. But from each of those solutions one can derive a new initial $A_0$ from where then again a complete list of solutions can be drawn: $$A_{k_1,k_2,k_3}=A_0 \cdot T^{k_1} \cdot C \cdot T^{k_2}\cdot C \cdot X^{k_3} \tag {4.2}$$ and so on, for arbitrarily many indexes/exponents $k_j$.
Using Pari/GP I implemented a function, which takes the vector of arbitrary length containing the indexes/exponents as argument in the following form $$ A_x = \operatorname{Pell2Sol}([k_1,k_2,...,k_n]) \tag 5$$ I don't know, whether the resulting structure can be called a tree or whether it has another name.

Q: As I said above, it seems to me that this indexing of solutions is exhaustive (assuming the number of indexes $n$ is taken from $1$ to $\infty$). Does someone know whether this is true?

Context: if this is a valid method, I'll try to extend/adapt it to some similar equations which occur in a combined set of equations in quadratics.

Here is a simple Pari/GP-implementation:
\\ Update: initial vector A0 is configurable to have a gcd-common factor
\\ vector E takes the indexes
\\ env=0 says: give only that solution, 
\\ env=1 says: printout the set of neigboured solutions towards index +- env.
\\ A0 can be configured by a common factor in the function call
{Pell2Sol(E,env=0,gcd=1)=my(A0=gcd*[1,1,1], T=[1,0,0;0,3,2;0,4,3],C=[0,1,0;1,0,0;0,0,1]);
     return(Mat(vectorv(2*env+1,r,A0 * T^(r-1-env))));}

Giving Pell2Sol([1,0],3) shows the following sequence of solutions:

  a    b    c
  7  -601  425
  7  -103   73
  7   -17   13
  7     1    5
  7    23   17
  7   137   97
  7   799  565

Here the negative signs should not be significant, since we use the squares anyway. But the signs are significant when computing the lists of solutions using the powers of the matrix T and the matrix C. I didn't supress the signs here to show that internal effect.

and Pell2Sol([1,5,1,0],3) gives

     a         b        c
  76793  -5375863  3801697
  76793   -920801   653365
  76793   -148943   118493
  76793     27143    57593
  76793    311801   227065
  76793   1843663  1304797
  76793  10750177  7601717

Perhaps a better picture than the two above is the following, (where I also draw the motivating connection to the problem of finding magic squares-of-squares) (Screenshot):


  • 3
    $\begingroup$ Hi, Gottfried. The primitive null vectors of a ternary quadratic form can be parametrized in the same way as Pythagorean triples. Here, $ \left( \frac{a-b}{2} \right)^2 + \left( \frac{a+b}{2} \right)^2 = c^2 $ Oh, in general more than one Pythagorean type formula may be required, see Mordell, Diophantine Equations $\endgroup$
    – Will Jagy
    Dec 1, 2021 at 3:13
  • $\begingroup$ en.wikipedia.org/wiki/Tree_of_primitive_Pythagorean_triples as well as page 47, Theorem 4 in Mordell $\endgroup$
    – Will Jagy
    Dec 1, 2021 at 3:20
  • 1
    $\begingroup$ math.stackexchange.com/questions/1767109/… $\endgroup$
    – individ
    Dec 1, 2021 at 5:16
  • 1
    $\begingroup$ Hi, Will! The wikipedia entry has grown much over that version that I roughly remember - this is interesting, thank you for the mmory-refresh! Using the link to cut-the-knot in that aricle I found some more essay, which goes even a step further in details of the matrices, and this looks much familiar in the first view (Lisa Winkler in web.archive.org/web/20080516081009/http://www.geocities.com/…) And proceeds also to quadruples, which occurs as well in my eqn.system of quadratics. Something to chew on! :-) $\endgroup$ Dec 1, 2021 at 8:16
  • $\begingroup$ A minor/trivial objection would be that your method only generated positive numbers, whereas changing the sign of any number in any solution yields another solution. Also $(0,0,0)$ is not generated by this method. Perhaps you could make explicit that you are only considering positive integral solutions. $\endgroup$ Dec 2, 2021 at 22:18

2 Answers 2


I'll assume that you are only considering positive integer solutions $(a,b,c)$, as your method only generates positive integer solutions. You already note that your method preserves $\gcd(a,b,c)$, and so the seed value $(1,1,1)$ can only generated primitive triples. The question then remains whether your method generates all primitive triples.

First let's characterize all primitive positive integral solutions to your equation. If $a$, $b$ and $c$ are positive integers with $\gcd(a,b,c)=1$ such that $$a^2+b^2=2c^2,$$ then as noted in the comments, we must have $a\equiv b\pmod{2}$ and so $x=\frac{a-b}{2}$ and $y=\frac{a+b}{2}$ are integers. They are easily verified be coprime, and if $(a,b,c)\neq(1,1,1)$ then we may interchange $a$ and $b$ if necessary so that $x$ and $y$ are both positive. It is easily verified that $$x^2+y^2=c^2,$$ meaning that $(x,y,c)$ is a primitive Pythagorean triple. Then it is well known that there exist coprime positive integers $m$ and $n$, not both odd and with $m>n$, such that $$x=m^2-n^2,\qquad y=2mn,\qquad c=m^2+n^2,$$ possibly with the roles of $x$ and $y$ reversed. Conversely, any such choice of $m$ and $n$ yields a primitive Pythagorean triple, and by extension any triple $$a=m^2+2mn-n^2,\qquad b=n^2+2mn-m^2,\qquad c=m^2+n^2,$$ is a solution to your original equation, and these are all primitive solutions.

Next let's consider the solutions generated from the seed $(1,1,1)$ by your method. Because $C^2=I$, the set of all generated solutions has the structure of a binary tree, by generating the two solutions $AT$ and $ATC$ from any previous solution $A$.

From every level of the binary tree to the next, we see that $c\mapsto\ 2b+3c$. This means that for every solution $(a,b,c)$ at the $k$-th level we have $c\geq3^k$, where level $k=0$ corresponds to the seed $A=(1,1,1)$. Because these are levels of a binary tree, this shows that your method generates at most $2^{k+1}$ solutions with $c<3^k$.

To conclude that your method does not generate all primitive solutions, it now suffices to show that the number of pairs of coprime positive integers $(m,n)$ with $m>n$ that satisfy $$m^2+n^2<3^k,$$ is greater than $2^{k+1}$. For this fact I refer to this paper, that shows as an example application of a more general theorem, that the number of such solutions is asymptotically proportional to $\frac{3^k}{2\pi}$.

  • $\begingroup$ Many thanks for this, it sounds really good! It's night here and I can go through this only tomorrow. Just a short note: Using negative powers of the matrix $T$ we arrive at negative values in $(a,b,c)$ as well. The first surprise when I saw this, that the set of (only) positives $(a,b,c)$ plus the set of absolute values $| (a,b,c)| $ seemed then to give all solutions in positives . I'll look at the paper tomorrow. $\endgroup$ Dec 3, 2021 at 0:32
  • $\begingroup$ I've looked at the article you've linked to. Whew, that's too much at the moment. But I'll follow your hint about the translation to the pythagorean coordinates $(c,x,y)$. Parallelizing that coordinates as they occur from one of my lists, give a analogon to my $T$ matrix as $U=\small \begin{bmatrix} 3 & 2 & -2 \\ 2 & 2 & -1 \\ -2 & -1 & 2 \end{bmatrix}$ . I'll try whether I can find some hint that my tree of lists might be incomplete with regards to the tree of pythagorean triples... Took a long time to catch the proposal for the link to the pythagorean coordinates...(sigh) $\endgroup$ Dec 3, 2021 at 7:26
  • $\begingroup$ Adding: ah, that $C$ matrix needs a analogon too; the pythagorean version is $D=\small \begin{bmatrix} 1 & . & . \\ . & 1 & . \\ . & . & -1 \end{bmatrix}$, and we can go along with $[1,0,1]\cdot U^{k_1} \cdot D \cdot U^{k_2} \cdot D \cdots U^{k_n}$. So the question can as well be answered when proved or disproved, that the so-created tree of pythagorean coordinates is incomplete. I know the transfer-matrix-representation of the pythagorean triples; the current representation looks indeed somehow different... $\endgroup$ Dec 3, 2021 at 7:34
  • $\begingroup$ Hmm, I've looked at parallel iterations, say one list in Pell-solutions (with a fixed value in $a$-column) and the parallel list in Pythagorean solutions where I use the transfermatrices Pell2Pyth and Pyth2Pell as right multiplicators to change from one format to the other. (This keeps the $c$-column unchanged, the Pell- and Pyth.version have the same $c$-values). Then I did the opposite way: with each of the 3 iteration matrices for the pythagorean triples I constructed "pyth"-lists and from that using Pyth2Pell the parallel Pell-lists. (...) $\endgroup$ Dec 4, 2021 at 3:14
  • $\begingroup$ (...) The resulting Pell-list has a surprising pattern, but still has a constant transfermatrix from one row to the other. But aside of the reordering only familiar looking values occured. For some samples tried anyway I found the needed key for the $Pell2Sol()$ function after some trial&error - this might be refinable. But I came over no observation, which would indicate that any of the pythagorean triples, translated by Pyth2Pell into the according Pell-triples could not be decoded by the $Pell2Sol()$ key. So it would be good, if you could any counterexample (...) $\endgroup$ Dec 4, 2021 at 3:20

Partial answer conjecture: No, this method does not catch up all solutions.

  • By samples of observations: it seems, it only catches only all solutions with no common divisor in $[a,b,c]$. If we insert in the function's initialization part instead of $A_0=[1,1,1]$ a multiple like $A_0=17 \cdot [1,1,1]$ , then we get the analogue tree, but with every number triple having that same common divisor.

Remark: the existence of these solutions with common gcd are sometimes nontrivial - in the problem of "magic-squares-of-squares" ("msqosq") we have triples made from some entries with a common divisor and can only be found when $A_0=gcd \cdot [1,1,1]$ is given as starting point. For instance in the only known msqosq with $7$ quadratic entries ($2$ are not squares) there is a triple $[h,b,e]=[-289,527,425]$ which can be found by the proposed function, by changing initialization to $A_0=17 \cdot[1,1,1]$ and then perform Pell2Sol([-2,-1,1]).

I have no idea so far, whether this gcd-aspect is the only part that must be introduced into the question's indexing method.

Triggered by the doubts as user Servaes presented in his/her answer, I tried to find hints, as when the Pythagorean-triple generation rules would construct one triple, which is not/cannot be in the Pell-solution- scheme as described in my question.

The first list is using the first generating-matrix ("WP1") as shown in wikipedia, starting at $A_{py,0}=[1 ,0, 1] $ $$ \begin{array} {l} \text{Pyth WP1 } \to \quad \text{Pell} \\ \small \begin{array} {rrr|rrr} x & y & c & a & b & c \\ \hline \vdots & \vdots & \vdots & \vdots & \vdots & \vdots \\ 1 & 0 & 1 & 1 & 1 & 1 \\ 3 & 4 & 5 & 7 & -1 & 5 \\ 5 & 12 & 13 & 17 & -7 & 13 \\ 7 & 24 & 25 & 31 & -17 & 25 \\ 9 & 40 & 41 & 49 & -31 & 41 \\ \vdots & \vdots & \vdots & \vdots & \vdots & \vdots \\ \end{array} \end{array} $$

The second list is using the second generating-matrix ("WP2") as shown in wikipedia, again starting at $A_{py,0}=[1 ,0, 1] $

$$ \begin{array} {l} \text{Pyth WP2 } \to \quad \text{Pell} \\ \small \begin{array} {rrr|rrr} x & y & c & a & b & c \\ \hline \vdots & \vdots & \vdots & \vdots & \vdots & \vdots \\ 1 & 0 & 1 & 1 & 1 & 1 \\ 3 & 4 & 5 & 7 & -1 & 5 \\ 21 & 20 & 29 & 41 & 1 & 29 \\ 119 & 120 & 169 & 239 & -1 & 169 \\ 697 & 696 & 985 & 1393 & 1 & 985 \\ \vdots & \vdots & \vdots & \vdots & \vdots & \vdots \\ \end{array} \end{array} $$

The third list is using the third generating-matrix ("WP3") as shown in wikipedia, now starting at $A_{py,0}=[3 ,4, 5] $. (A starting triple $[1,0,1]$ runs into a list of zeros)

$$ \begin{array} {l} \text{Pyth WP3 } \to \quad \text{Pell} \\ \small \begin{array} {rrr|rrr} x & y & c & a & b & c \\ \hline \vdots & \vdots & \vdots & \vdots & \vdots & \vdots \\ 3 & 4 & 5 & 7 & -1 & 5 \\ 15 & 8 & 17 & 23 & 7 & 17 \\ 35 & 12 & 37 & 47 & 23 & 37 \\ 63 & 16 & 65 & 79 & 47 & 65 \\ 99 & 20 & 101 & 119 & 79 & 101 \\ \vdots & \vdots & \vdots & \vdots & \vdots & \vdots \\ \end{array} \end{array} $$

All three Pell-lists look very innocent, and trying some examples from this showed, that all generated Pell-triples could be generated by appropriate keys for the function $Pell2Sol()$


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.