Functions $f: \mathbb{Z}^{+}\to \mathbb{R}$ satisfying $x f(y) + y f(x) = (x+y) f(x^2+y^2)$

Here is a partial result : I show below that $f(k)=f(1)$ for $1\leq k \leq 20$ (and I think it very likely that $f$ is indeed constant).

Let $a=f(1)$ and $A=\lbrace x\in{\mathbb N} | f(x)=a \rbrace$. We see from the functional equation that if any two of $x,y,x^2+y^2$ are in $A$, then the third is in $A$ also. In particular :

Rule $R_1$ : if $x,y\in A$, then $x^2+y^2\in A$. Rule $R_2$ : if $x,y,z\in A$, then $\sqrt{x^2+y^2-z^2} \in A$ (assuming this is an integer).

Since $2=1^2+1^2$, $5=1^2+2^2$, $8=2^2+2^2$, $50=5^2+5^2$, we have by rule $R_1$ that $\lbrace 1,2,5,8,50 \rbrace \subseteq A$. Since $7=\sqrt{5^2+5^2-1^2}$, $4=\sqrt{1^2+8^2-7^2}$, $17=1^2+4^2$, $20=2^2+4^2$ we deduce

$$\lbrace 1,2,4,5,7,8,17,20,50 \rbrace \subseteq A \tag{1}$$

Next, we have $f(10)=f(1^2+3^2)=\frac{3a+f(3)}{4}$ and $f(13)=\frac{3a+2f(3)}{5}$. From $2^2+11^2=5^2+10^2$, we deduce $\frac{2f(11)+11a}{13}=\frac{5f(10)+10a}{15}=\frac{f(10)+2a}{3}=\frac{f(3)+11a}{12}$. Similarly, from $7^2+11^2=1^2+13^2$, we deduce $\frac{7f(11)+11a}{18}=\frac{f(13)+13a}{14}=\frac{f(3)+34a}{35}$. We then have a Cramer $2\times 2$ system in $f(3)$ and $f(11)$ (and parameter $a$), so $f(3)=f(11)=a$. Therefore :

$$\lbrace 1,2,3,4,5,7,8,10,11,13,17,20,50 \rbrace \subseteq A \tag{2}$$

Then, from $9=\sqrt{3^2+11^2-7^2}$, $12=\sqrt{8^2+9^2-1^2}$, $18=3^2+3^2$, $6=\sqrt{2^2+9^2-7^2}$, we further deduce

$$[|1..13|]\cup\lbrace 17,18,20,50 \rbrace \subseteq A \tag{3}$$

Finally, from $14=\sqrt{10^2+10^2-2^2}$, $15=\sqrt{9^2+13^2-5^2}$, $16=\sqrt{8^2+14^2-2^2}$, $19=\sqrt{13^2+14^2-2^2}$, we obtain $[|1..20|] \subseteq A$ as announced.

We show that for any $N \ge 15$, there are $a,b,c \le N$ such that $N+1 = \sqrt{a^2+b^2-c^2}$. By Ewan Delanoy's answer, this suffices. For $N \equiv 0,1,2 \pmod{3}$, we can take $(a,b,c) = (N-1,\frac{N+9}{3},\frac{N-9}{3}),(N, \frac{N+5}{3}, \frac{N-4}{3}), (N-2, \frac{N+13}{3}, \frac{N-14}{3})$, respectively.

Moving forward with the solution idea by @Evan Delanoy, I will provide an algorithm which could lead to any number in $\mathbb{N}$ being in the set $A$. Suppose the required number is $x$.

Algorithm:

A = {1..20}  
while x is not in A:
    S = A
    for each pair (a, b) in A x A:
        S = S U {a^2 + b^2}
    for each triple (a, b, c) in A x A x A:
        if a^2 + b^2 - c^2 = d^2 where d is a positive integer:
            S = S U {d}
    A = S

We now show that this terminates.

Suppose that after the $i^\mathrm{th}$ iteration, all positive integers not exceeding $n$ are in the set $A$. We proceed to use the next few iterations to force all numbers less than or equal to $3n/2$ to lie in $A$.

The first loop gives us all positive integers less not exceeding $n^2$ which can be represented as the sum of 2 positive squares. To characterize these, let $W$ be the integers not exceeding $n^2$ in whose prime factorization, the exponent of any prime which is 3 mod 4 is even, (call such primes bad primes). Then the complement of the required set wrt $W$ is the set of numbers not exceeding $n^2$ which are perfect squares but can not be represented as a sum of 2 positive squares. For this to hold true, their square root, say $r$, is not a part of a Pythagorean triplet which is composed entirely of positive integers. We claim that a necessary and sufficient condition for such $r$ is to not have any prime factor congruent to 1 mod 4, which can be seen by applying theorem 4.4 here.

Now we try to construct numbers between $n$ and $3n/2$ in which the powers of bad primes are allowed to be odd, as well as those which are perfect squares of numbers divisible only by 2 and/or bad primes.

We give a proof using induction.

Firstly consider the case of the number not being a power of $2$.

Suppose the number is $x$, and $x = vp$, where $p$ is the product of all the bad primes which divide $x$ and have an odd exponent in the prime factorization, and $v$ is divisible by at least one prime congruent to 1 mod 4. Then $v^2p^2 + v^2 = (p^2 + 1)(a^2 + b^2) = (pb - a)^2 + (a + pb)^2$ (where $(a, b)$ satisfies $a < b$, $a^2 + b^2 = v$) implies that $vp$ is in $S$, since $v, pb - a, pa + b$ are all less than $vp$.

Now suppose the number is divisible by no prime which is 1 mod 4. This part still needs work. Note that using the proof for the two cases below, we can exhaust the cases where x is $(\mathrm{product \,of \,bad \,primes})^2 2^\alpha$ and $\alpha > 2$.

Now suppose that the number (say $x$) is a power of $2$. If the exponent of $2$ in it is odd, then we are done by using a trivial split into equal parts. So assume that $x = 2^{2k}$. Then we have $x^2 + 2^2 = 2^2(4^{2k-1} + 1)$. Now for the second factor on the right, we know that it is a sum of two squares (namely $1$ and $4^{2k-1}$), and since $5$ divides the whole thing, and $5 \equiv 1 \pmod 4$, $(4^{2k-1} + 1)/5$ should admit a sum of squares representation with 2 positive squares (else it is a power of 2 multiplied by, possibly 0, bad primes, however since -1 is a quadratic nonresidue modulo these primes, it needs to be a power of 2, which is impossible by parity reasons if the power of 2 is > 1 and size reasons if it is 1). So we have for some positive integers $a < b$, $4^{2k} + 4 = 20(a^2 + b^2) = (4b - 2a)^2 + (4a + 2b)^2$. For showing that both $j = 4b - 2a$, $k = 4a + 2b$ do not exceed $2^{2k}$, we can do the following. Suppose one of them (wlog $j$) did exceed $2^{2k}$. Then since $j \ge 2^{2k} + 1$, $j^2 \ge 4^{2k} + 2 \cdot 2^{2k} + 1 > 4^{2k} + 4$, which is a contradiction. So since by the end of the this iteration, we would have $4b - 2a$ and $4a + 2b$ both in our set $A$, the next iteration would allow us to include $2^{2k}$ in our set as well, which completes the proof.

Note: this is still not complete.