Does these arithmetic means of Pythagorean triangles converge?
$\newcommand{\h}{\mathcal{h}}$ $\newcommand{\n}{\mathcal{n}}$ $\newcommand{\deq}{\stackrel{\text{def}}{=}}$ $\newcommand{\abs}[1]{\lvert #1 \rvert}$
Disclaimer: this answer has been edited to account for changes in the question and remove false claims.
For the first claim: parametrize Pythaogrean triples with the usual stereographic projection $$(a,b,c)=(q^2-p^2,2pq,q^2+p^2)$$ where $0<p<q$ and $p,q$ are coprime.
Then we seek the aymptotic behavior as $r\to\infty$ of the mean value of $$f(t)\deq \frac{2(t+1)}{t^2+1}$$
over all rational numbers $0<t<1$ such that
$$\h(t)<r$$
where for a rational number of the form $p/q$, $p,q$ coprime
$$\h(p/q)\deq \abs{p}\vee\abs{q}\text{.}$$
Write $\mu_{\h,r}$ for the probability measure associated over taking the mean over rationals $0<t<1$ such that $\h(t)<r$.
Now, $\mu_{\h,r}$ represents a mean with respect to rationals $t$ such that $\h(t)<r$. But since $0<t<1$, this is really a mean with respect to rationals such that the denominator $q$ is less than $r$. The sequences of such rationals are known as Farey sequences, and it is known that these are asymptotically equidistributed, so that the limiting measure is Lebesgue measure:
$$\lim_{h\to\infty}\mu_{\h,r}=\lambda\text{.}$$
Therefore the desired limiting mean value is
$$\int_0^1\frac{2(t+1)\mathrm{d}t}{t^2+1}=\frac{\pi}{2}+\log 2\text{.}$$
Edit: as @Blue pointed out in the comments, we must take into account excluding fractions for which both the numerator and denominator are odd. It is likely that these can also be shown to be equidistributed by Weyl's Criterion.
I think the limit of the average perimeter/hypotenuse values depends on the order in which the Pythagorean triples are generated.
The program in the OP generates the triples $(r^2+s^2, 2rs, r^2-s^2)$ in order of increasing $r$. The description however was (before the edit) about taking the average over those triangles with the hypotenuse below some bound $n$ (and then letting $n$ go to infinity).
This bound changes the result because if $r^2$ is close to $n$, then $s$ cannot take values almost as high as $r$ because it is bounded by $\sqrt{n−r^2}$. This leaves out some more acute triangles (with $s$ near $r$) that have a low ratio, and so increases the value of the average. If you generate the triples in order of increasing $r$, those acute triangles are structurally shifted earlier in the sequence compared to if you generated them in order of increasing hypotenuse, therefore making all the partial averages smaller.
I get a limit of about $2.2732$ instead.
Here is the straighforward C# code I used. max
is the (strict) upper bound on the hypotenuse length.
using System;
namespace test
{
/* max average
* 10^7 2.2734207124719
* 10^8 2.27329667075612
* 10^9 2.27325757481033
* 10^10 2.27324525141887
* 10^11 2.27324135532923
*/
class Msepythlimit
{
static void Main()
{
long n = 0;
double sum = 0;
double max = 10000000;
for (long r = 2; r*r <= max; r++)
{
for (long s = 1 + r % 2; s < r && s * s + r * r < max; s++)
{
if (Gcd(r, s) == 1)
{
long a = r * r + s * s;
long b = r * r - s * s;
long c = 2 * r * s;
n++;
sum += (double)(a + b + c) / a;
if (n % 100000 == 0) Console.WriteLine("{0}: {1}", n, sum / n);
}
}
}
double avg = sum / n;
Console.WriteLine(avg);
}
static long Gcd(long a, long b)
{
long x = a;
long y = b;
while (x > 0)
{
long t = y % x;
y = x;
x = t;
}
return y;
}
}
}