The answer is, for any subset $A\subseteq (-\infty,0)$, the function $$f(x) = \begin{cases} -x & x\in A \\ x & x\notin A. \end{cases}$$To check it, note that $f(x)=x$ for all $x\geq 0$, and $f(x)^2=x^2$ for all $x$. Thus the equation is reduced to $$2(x^2+y^2)= (x+f(y))^2 + (x-f(y))^2,$$which is indeed true. Now we restrict our attention into proving that these are all solutions. Let $P(x,y)$ denote the assertion $2f(x^2+y^2) = (x+f(y))^2 + f(x-f(y))^2$. Comparing $P(x,y)$ and $P(-x,y)$ gives $$f(x-f(y))^2 - f(-x-f(y))^2 = -4f(y)x \qquad (*)$$Since the constant function does not work, we fix two real numbers $a,b$ such that $f(a)\ne f(b)$. We use the above equation twice to obtain the equations \begin{align*} f(x+f(a)-f(b))^2 - f(-x-f(a)-f(b))^2 &= -4f(b)(x+f(a)) \\ f(x+f(b)-f(a))^2 - f(-x-f(b)-f(a))^2 &= -4f(a)(x+f(b)), \end{align*}which upon subtracting gives $$f(x+f(a)-f(b))^2 - f(x+f(b)-f(a))^2 = 4x(f(a)-f(b)).$$Setting $\delta = f(a)-f(b)$, this can be written more compactly as $$f(x+\delta)^2 - f(x-\delta)^2 = 4x\delta.$$Motivated by this, we subtract $P(x+\delta, y)$ and $P(x-\delta, y)$, which gives \begin{align*} 2f((x+\delta)^2 + y^2)&= (x+\delta + f(y))^2 + f(x+\delta - f(y))^2 \\ 2f((x-\delta)^2 + y^2) &= (x-\delta + f(y))^2 + f(x-\delta - f(y))^2 \\[4pt] 2f((x+\delta)^2 + y^2)- 2f((x-\delta)^2 + y^2) &= 4\delta(x+f(y)) + 4\delta(x-f(y)) \end{align*}Therefore if $g(x)=f(x)-x$, then $$g((x+\delta)^2+y^2) = g((x-\delta)^2+y^2)$$for all reals $x,y$. We claim that this gives $g$ eventually constant. To see why, select $u,v$ which $u,v > (1+|\delta|)^2$ and $u-v\in (0,4\delta)$. Selecting $x = \tfrac{u-v}{4\delta}$ and appropriate value of $y$, we obtain $g(u)=g(v)$. This clearly extends to any $u,v$ large enough, so $g$ is indeed eventually constant (say in $[M,\infty)$ where $M>2020$). We are nearly done. Supposethat $g(x)=c$ for all $x\geq M$. Then, for any $x,y>M$, $f(x)=x+c$, so \begin{align*} &2(x^2+y^2+c) = (x+y+c)^2 + f(x-y-c)^2 \\ \implies &f(x-y-c)^2 = x^2+y^2-2xy-2cx-2cy-c^2+2c = \underbrace{(x-y)^2 - (c^2-2c)}_{\text{fixed}}+2c(x+y) \end{align*}This forces $c=0$ and $f(t)=\pm t$ for all $t$. However, $LHS\geq 0$, hence $f(x)\geq 0$ for all $x\geq 0$. This gives the described set of solution.

Consider $a,b$ such that $f(a)\neq f(b)$ (there exists such pair since $f$ can’t be constant). Pick $s,t$ such that $s^2+a^2=t^2+b^2$ and $s-f(a)=t-f(b)$, and plug $(x,y)$ by $(s,a)$ and $(t,b)$ in the equation, we obtain $|s+f(a)|=|t+f(b)|$. Clearly, if $s+f(a)=t+f(b)$ then $f(a)=f(b)$, a contradiction. So, we have $s+f(a)+t+f(b)=0$, and thus $s=-f(b), t=-f(a)$. Hence, we deduce that if $f(a)\neq f(b)$ then $f(a)^2-a^2=f(b)^2-b^2$. Moreover, if $f(m)=f(n)$, just consider $c$ s.t. $f(c)\neq f(m)$, then $f(c)^2-c^2=f(m)^2-m^2=f(n)^2-n^2$, which stills preserves the equation $f(m)^2-m^2=f(n)^2-n^2$. Therefore, $f(x)^2-x^2=f(y)^2-y^2$ $\forall x,y$, and we conclude that $f(x)^2=x^2+c$ $\forall x$, $c$ is a constant. Note that $f(x^2)\geq 0$ $\forall x$, so $c=0$, $f(x)=x$ $\forall x\geq 0$, $f(x)=\pm x$ for each $x<0$. Just try plug back into the original equation and we can see that this answer satisfies the condition.

Oh. I've heard that Korean wintercamp test is held twice during the camp. Can I have the other one? (to any Korean users) I love those Korean questions cause they are very challenging

Really beautiful problem Here is a new way of solving the problem: We claim that the only solution is $f(x)=x \forall x\ge 0$ and $f(x)=xg(x) \forall x<0$, where $g$ is any function mapping $\mathbb{R}^-$ to $\{1,-1\}$ Verification: Note that $f(x)^2=x^2$ for all $x$. Hence $2f(x^2+y^2)=(x+f(y))^2+f(x-f(y))^2$ $\iff 2x^2+2y^2=x^2+f(y)^2+2xf(y)+(x-f(y))^2$ which is true. We now show that this is the only solution. Let $P(x,y)$ denote the assertion that $2f(x^2+y^2)=(x+f(y))^2+f(x-f(y))^2$ $P(\frac{1}{2},\frac{1}{2}): 2f(\frac{1}{2})=\frac{1}{4}+f(\frac{1}{2})^2+f(\frac{1}{2})+f(\frac{1}{2}-f(\frac{1}{2}))^2$ which is equivalent to $(\frac{1}{2}-f(\frac{1}{2}))^2+f(\frac{1}{2}-f(\frac{1}{2}))^2=0$ This is only possible when $f(\frac{1}{2})=\frac{1}{2}$ and $f(\frac{1}{2}-f(\frac{1}{2}))=0$, and thus combining the two $f(0)=0$ $P(x,0),P(0,x): 2f(x^2)=x^2+f(x)^2=f(x)^2+f(-f(x))^2$ In particular $f(x)^2=f(-x)^2 \implies x^2=f(-f(x))^2=f(f(x))^2$ Suppose $f(f(c))=-c$. Note that the original condition implies $f(\text{positive})=\text{non-negative}$ Thus if $c$ is positive then the LHS is positive while the RHS is negative, contradiction. Thus $c\le 0$. This means that $f(f(x^2))=x^2$. $P(f(y),y),P(-f(y),y): 4f(y)^2=2f(f(y)^2+y^2)=f(-2f(y))^2=f(2f(y))^2$. Take $y=f(x^2)$, $f(2x^2)^2=4x^4$ Since $f(2x^2)\ge 0$, $f(2x^2)=2x^2$ so $f(x)=x \forall x\ge 0$ Now since $f(x^2)^2=f(-x^2)^2 \implies f(-x^2)^2=x^4 \implies f(x)=\pm x \forall x<0$ Hence the general solution is $f(x)=x \text{ if } x \ge 0$ and $xg(x) \text{ if } x<0$, for any $g:\mathbb{R}^-\rightarrow \{1,-1\}$