A Class of Diophantine Equations in Three Variables This question is inspired by the problem https://projecteuler.net/problem=748
Consider the Diophantine equation
$$\frac{1}{x^2}+\frac{1}{y^2}=\frac{k}{z^2}$$
$k$ is a squarefree number. $A_k(n)$ is the number of solutions of the equation such that $1 \leq x+y+z \leq n$, $x \leq y$ and $\gcd(x,y,z)=1$. This equation has infinite solutions for $k=1$ and $k&gt;1$ that can be expressed as sum of two perfect squares.
Let $$A_k=\lim_{n \to \infty}\frac{A_k(n)}{\sqrt{n}}$$
\begin{array}{|c|c|c|c|c|}
\hline
k&amp; A_k\left(10^{12}\right)&amp; A_k\left(10^{14}\right)&amp; A_k\left(10^{16}\right)&amp; A_k\left(10^{18}\right)&amp; A_k \\ \hline
1&amp; 127803&amp; 1277995&amp; 12779996&amp; 127799963&amp; 0.12779996...\\ \hline
2&amp; 103698&amp; 1037011&amp; 10369954&amp; 103699534&amp; 0.1036995...\\ \hline
5&amp; 129104&amp; 1291096&amp; 12911049&amp; 129110713&amp; 0.129110...\\ \hline
10&amp; 90010&amp; 900113&amp; 9000661&amp; 90006202&amp; 0.0900062...\\ \hline
13&amp; 103886&amp; 1038829&amp; 10388560&amp; 103885465&amp; 0.103885...\\ \hline
17&amp; 86751&amp; 867550&amp; 8675250&amp; 86752373&amp; 0.086752...\\ \hline
\end{array}
From these data, it seems that these limits converge. I wonder if it is possible to write them in terms of known constants or some sort of infinite series. For Pythagorean triples, there are about $\frac{n}{2\pi}$ triples with hypotenuse $\leq n$.
The Python code to calculate $A_1(n)$:
def gcd(a, b):
   if b == 0:
       return a
   return gcd(b, a % b)

N = 10 ** 14

cnt = 0
for a in range(1, 22000):
   a2 = a * a
   for b in range(1, a):
       if (a + b) % 2 == 1 and gcd(a, b) == 1:
           b2 = b * b
           x = 2 * a * b * (a2 + b2)
           y = a2 * a2 - b2 * b2
           z = 2 * a * b * (a2 - b2)
           if x + y + z &gt; N:
               continue
           cnt += 1
print(cnt)

The Python code to calculate $A_2(n)$:
def gcd(a, b):
   if b == 0:
       return a
   return gcd(b, a % b)

N = 10 ** 14

cnt = 1     # (1, 1, 1)
for a in range(1, 22000):
   a2 = a * a
   a4 = a2 * a2
   for b in range(1, a):
       if gcd(a, b) == 1:
           b2 = b * b
           b4 = b2 * b2
           x = 2 * a * b * (a2 + b2) - (a4 - b4)
           y = 2 * a * b * (a2 + b2) + (a4 - b4)
           z = 6 * a2 * b2 - (a4 + b4)
           if x &gt; 0 and y &gt; 0 and z &gt; 0 and x + y + z &lt;= N and gcd(x, gcd(y, z)) == 1:
               cnt += 1
           x = (a4 - b4) - 2 * a * b * (a2 + b2)
           y = 2 * a * b * (a2 + b2) + (a4 - b4)
           z = (a4 + b4) - 6 * a2 * b2
           if x &gt; 0 and y &gt; 0 and z &gt; 0 and x + y + z &lt;= N and gcd(x, gcd(y, z)) == 1:
               cnt += 1
print(cnt)
This is a partial answer. Consider the case $k=1$. The parametric solutions are given by $(x,y,z)=\left(2ab(a^2+b^2),a^4-b^4,2ab(a^2-b^2)\right)$, where $\gcd(a,b)=1$, $a&gt;b&gt;0$ and $a+b$ is odd. In the triples $(x,y,z)$ generated by this parametrization, it is not necessarily always $x\leq y$. For example, $(20,15,12)$ is a triple generated by this parametrization, and $(15,20,12)$ will not be generated. This parametrization generates only the unique solutions. So we can have two possible cases, either $0&lt;a^2-b^2\leq 2ab$ or $0&lt;2ab\leq a^2-b^2$.
The sum of the triples = $x+y+z=4a^3b+a^4-b^4 \leq n$.
Let us consider the first case. Denote the number of solutions by $A_1^1(n)$ and $A_1^1=\lim_{n \to \infty}\frac{A_1^1(n)}{\sqrt{n}}$.
$$0&lt;a^2-b^2\leq 2ab \implies \sqrt{2}-1\leq\frac{b}{a}=x&lt;1$$
Number of pairs $(a,b)$ satisfying above constraints and also $4a^3b+a^4-b^4 \leq n$ can be approximated by
\begin{align}
\begin{split}
p(n)&amp;\approx\int_{\frac{b}{a}=\sqrt{2}-1}^{1}\int_{4a^3b+a^4-b^4 \leq n}a \ da\ db\\
&amp;=\int_{x=\sqrt{2}-1}^{1}\int_{0&lt;a^4\leq \frac{n}{1+4x-x^4}}a \ da\ dx \\
&amp;=\frac{\sqrt{n}}{2}\int_{x=\sqrt{2}-1}^{1}\frac{1}{\sqrt{1+4x-x^4}}\ dx
\end{split}
\end{align}
The natural density of coprime numbers is $\frac{6}{\pi^2}$. The condition $gcd(a,b)=1$ will contribute a multiplication factor of $\frac{6}{\pi^2}$. Among all the coprime pairs $(a,b)$ such that $a&gt;b&gt;0$, fraction of which have odd value of $a+b$ is $\frac{2}{3}$. We have to multiply this factor.
\begin{align}
\begin{split}
A_1^1 &amp;=\lim_{n \to \infty}\frac{A_1^1(n)}{\sqrt{n}}\\
&amp;= \lim_{n \to \infty} \frac{2}{3} \cdot \frac{6}{\pi^2} \cdot\frac{p(n)}{\sqrt{n}}\\
&amp;= \frac{2}{\pi^2} \int_{x=\sqrt{2}-1}^{1}\frac{dx}{\sqrt{1+4x-x^4}} \\
\end{split}
\end{align}
Similarly, it can be shown that for the other case exactly the same limit exists. So, $$A_1=2A_1^1=\frac{4}{\pi^2} \int_{x=\sqrt{2}-1}^{1}\frac{dx}{\sqrt{1+4x-x^4}} \approx 0.1277999513464289...$$
Using a similar approach it can be shown that
\begin{align}
\begin{split}
A_2&amp;=\frac{2}{\pi^2} \left[ \int_{x=0}^{\sqrt{2}-1}\frac{dx}{\sqrt{3-6x^2-x^4}} +\int_{x=\sqrt{2}-1}^{1}\frac{dx}{\sqrt{4x+6x^2+4x^3-1-x^4}} \right]\\
&amp; \approx 0.1036994744684913...
\end{split}
\end{align}
\begin{align}
\begin{split}
A_{13}&amp;=\frac{52}{7\pi^2} \int_{x=(\sqrt{13}-3)/2}^{(\sqrt{26}-1)/5}\frac{dx}{\sqrt{20x+36x^2-5-7x^4}}\approx 0.1038855856479065...
\end{split}
\end{align}
To solve the Diophantine equation.
$$ \frac{ 1 }{ x^2 } +\frac{ 1 }{ y^2 } = \frac{ q }{ z^2 } $$
It is necessary to use the solution of the following equation.
$$a^2+b^2=qc^2$$
There are solutions when the coefficient can be represented as the sum of squares. $q=t^2+k^2$
$$a=-tp^2+2kps+ts^2$$
$$b=kp^2+2tps-ks^2$$
$$c=p^2+s^2$$
And the solution itself can be presented in this form.
$x=ac$
$y=bc$
$z=ab$