How does the size of the set $$A(R) = \{(a,b) \; | \; a,b \in \mathbb{N} \times \mathbb{N}, \; \gcd(a,b) = 1, \; a^2 + b^2 \leq R^2\}$$ grow as a function of $R$?
My try: It's clear that $|A(R)| \leq \pi R^2$ (ish). Since $\mathbb{Z}_R \times \{1\} \subset A(R)$ we can also say $A(R) \geq R$ (ish). I think we can write a formula similar to $$|A(R)| = \frac{1}{2} \sum_{k \leq R} \phi(\lfloor \sqrt{k^2-R^2} \rfloor),$$ but I'm pretty sure the above is wrong. Even if it was correct, I don't have any idea how to estimate it.
It also occurs to me that if $F(n)$ is the number of ways to represent $n=a^2+b^2$ with $\gcd(a,b) = 1$, then $$A(R) = \sum_{k \leq R^2} F(k).$$ If I recall $F(n)$ is exactly known. I think $$F \left(a^2 \prod_{i=1}^m q_i \prod_{i=1}^n p_i \right) = 2^n$$ where $n/a^2$ is square free, $q_i = 3 \mod 4$ and $p_i = 1 \mod 4$. It seems this might be difficult to work with.
If we disregard the coprimality constraint for a moment, we look at the sets
$$B(R) = \{(a,b) \in\mathbb{N}^2 : a^2 + b^2 \leqslant R^2\}.$$
It is easy to see that
$$\lvert B(R)\rvert = \frac{\pi}{4} R^2 + O(R).\tag{1}$$
On the other hand, each pair $(a,b)$ with $a^2 + b^2$ can uniquely be written as $d\cdot (\alpha,\beta)$ with $\gcd(\alpha,\beta) = 1$, so we have
$$\lvert B(R)\rvert = \sum_{d=1}^R \left\lvert A\biggl( \frac{R}{d}\biggr)\right\rvert.\tag{2}$$
By generalised Möbius inversion, $(2)$ is equivalent to
$$\lvert A(R)\rvert = \sum_{d=1}^R \mu(d)\left\lvert B\biggl(\frac{R}{d}\biggr)\right\rvert.\tag{3}$$
where $\mu$ is the Möbius function. If $R$ is constrained to be an integer, we must replace $\frac{R}{d}$ with $\left\lfloor\frac{R}{d}\right\rfloor$ in $(2)$ and $(3)$, and below.
Inserting $(1)$ into $(3)$, we obtain
$$\lvert A(R)\rvert = \sum_{d=1}^R \mu(d) \left(\biggl( \frac{R}{d}\biggr)^2 + O\biggl(\frac{R}{d}\biggr)\right),\tag{4}$$
and as $\left\lvert \sum\limits_{d=1}^R\frac{\mu(d)}{d}\right\rvert \leqslant \sum\limits_{d=1}^R\frac{1}{d} = \log R + O(1)$, that becomes
$$\lvert A(R)\rvert = \frac{\pi}{4}R^2\sum_{d=1}^R \frac{\mu(d)}{d^2} + O(R\log R),\tag{5}$$
and since
$$\sum_{d=1}^\infty \frac{\mu(d)}{d^2} = \frac{1}{\zeta(2)} = \frac{6}{\pi^2}$$
and
$$\left\lvert \sum_{d=R+1}^\infty \frac{\mu(d)}{d^2}\right\rvert \leqslant \sum_{d=R+1}^\infty \frac{1}{d^2} < \frac{1}{R},$$
finally
$$\lvert A(R)\rvert = \frac{3}{2\pi}R^2 + O(R\log R).$$
If we restrict the argument of $A$ and $B$ to be integers, then in $(4)$ the dominant term is
$$\sum_{d=1}^R \mu(d)\biggl\lfloor\frac{R}{d}\biggr\rfloor^2 = \sum_{d=1}^R\mu(d)\biggl(\frac{R}{d}\biggr)^2 - \sum_{d=1}^R \mu(d)\Biggl(\biggl(\frac{R}{d}\biggr)^2 - \biggl\lfloor \frac{R}{d}\biggr\rfloor^2\Biggr),$$
and since
$$\biggl(\frac{R}{d}\biggr)^2 - \biggl\lfloor \frac{R}{d}\biggr\rfloor^2 = \underbrace{\Biggl(\frac{R}{d}-\biggl\lfloor\frac{R}{d}\biggr\rfloor\Biggr)}_{< 1}\underbrace{\Biggl(\frac{R}{d}+\biggl\lfloor\frac{R}{d}\biggr\rfloor\Biggr)}_{\leqslant 2\frac{R}{d}} < 2\frac{R}{d}$$
we see that the asymptotic behaviour doesn't change, only the constant on the second order term changes.