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>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& A_k\left(10^{12}\right)& A_k\left(10^{14}\right)& A_k\left(10^{16}\right)& A_k\left(10^{18}\right)& A_k \\ \hline
1& 127803& 1277995& 12779996& 127799963& 0.12779996…\\ \hline
2& 103698& 1037011& 10369954& 103699534& 0.1036995…\\ \hline
5& 129104& 1291096& 12911049& 129110713& 0.129110…\\ \hline
10& 90010& 900113& 9000661& 90006202& 0.0900062…\\ \hline
13& 103886& 1038829& 10388560& 103885465& 0.103885…\\ \hline
17& 86751& 867550& 8675250& 86752373& 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 > 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 > 0 and y > 0 and z > 0 and x + y + z <= 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 > 0 and y > 0 and z > 0 and x + y + z <= N and gcd(x, gcd(y, z)) == 1:
cnt += 1
print(cnt)
Best Answer
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>b>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<a^2-b^2\leq 2ab$ or $0<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<a^2-b^2\leq 2ab \implies \sqrt{2}-1\leq\frac{b}{a}=x<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)&\approx\int_{\frac{b}{a}=\sqrt{2}-1}^{1}\int_{4a^3b+a^4-b^4 \leq n}a \ da\ db\\ &=\int_{x=\sqrt{2}-1}^{1}\int_{0<a^4\leq \frac{n}{1+4x-x^4}}a \ da\ dx \\ &=\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>b>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 &=\lim_{n \to \infty}\frac{A_1^1(n)}{\sqrt{n}}\\ &= \lim_{n \to \infty} \frac{2}{3} \cdot \frac{6}{\pi^2} \cdot\frac{p(n)}{\sqrt{n}}\\ &= \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&=\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]\\ & \approx 0.1036994744684913... \end{split} \end{align} \begin{align} \begin{split} A_{13}&=\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}