IMO Shortlist, 2020

A straightforward check shows that $f(x)=x+1$ satisfies (*). We divide the proof of the converse statement into a sequence of steps.

Step 1: $f$ is injective.
Put $x=1$ in (*) and rearrange the terms to get
$$y=f(1)f(y)+1-f(1+f(y))$$
Therefore, if $f(y_1)=f(y_2)$, then $y_1=y_2$.

Step 2: $f$ is (strictly) monotone increasing.
For any fixed $y \in \mathbb{R}^{+}$, the function
$$g(x):=f(x+f(xy))=f(x)f(y)+1-y$$
is injective by Step 1. Therefore, $x_1+f(x_1y) \neq x_2+f(x_2y)$ for all $y, x_1, x_2 \in \mathbb{R}^{+}$ with $x_1 \neq x_2$. Plugging in $z_i=x_iy$, we arrive at
$$\frac{z_1-z_2}{y} \neq f(z_2)-f(z_1), \quad \text{or} \quad \frac{1}{y} \neq \frac{f(z_2)-f(z_1)}{z_1-z_2}$$
for all $y, z_1, z_2 \in \mathbb{R}^{+}$ with $z_1 \neq z_2$. This means that the right-hand side of the rightmost relation is always non-positive, i.e., $f$ is monotone non-decreasing. Since $f$ is injective, it is strictly monotone.

Step 3: There exist constants $a$ and $b$ such that $f(y)=a y+b$ for all $y \in \mathbb{R}^{+}$.
Since $f$ is monotone and bounded from below by $0$, for each $x_0 \geqslant 0$, there exists a right $\operatorname{limit} \lim_{x \searrow x_0} f(x) \geqslant 0$. Put $p=\lim_{x \searrow 0} f(x)$ and $q=\lim_{x \searrow p} f(x)$.
Fix an arbitrary $y$ and take the limit of $(*)$ as $x \searrow 0$. We have $f(xy) \searrow p$ and hence $f(x+f(xy)) \searrow q$; therefore, we obtain
$$q+y=p f(y)+1, \quad \text{or} \quad f(y)=\frac{q+y-1}{p} .$$
(Notice that $p \neq 0$, otherwise $q+y=1$ for all $y$, which is absurd.) The claim is proved.

Step 4: $f(x)=x+1$ for all $x \in \mathbb{R}^{+}$.
Based on the previous step, write $f(x)=a x+b$. Putting this relation into (*) we get
$$a(x+a x y+b)+b+y=(a x+b)(a y+b)+1$$
which can be rewritten as
$$(a-a b)x+(1-a b)y+a b+b-b^2-1=0 \quad \text{for all} \ x, y \in \mathbb{R}^{+} .$$
This identity may hold only if all the coefficients are $0$, i.e.,
$$a-a b=1-a b=a b+b-b^2-1=0 .$$
Hence, $a=b=1$.

We provide another proof that $f(x)=x+1$ is the only function satisfying (*).
Put $a=f(1)$. Define the function $\phi: \mathbb{R}^{+} \rightarrow \mathbb{R}$ by
$$\phi(x)=f(x)-x-1$$
Then equation (*) reads as
$$\phi(x+f(xy))=f(x)f(y)-f(xy)-x-y .$$
Since the right-hand side is symmetric under swapping $x$ and $y$, we obtain
$$\phi(x+f(xy))=\phi(y+f(xy)) .$$
In particular, substituting $(x, y)=(t, 1/t)$ we get
$$\phi(a+t)=\phi\left(a+\frac{1}{t}\right), \quad t \in \mathbb{R}^{+} .$$
Notice that the function $f$ is bounded from below by a positive constant. Indeed, for each $y \in \mathbb{R}^{+}$, the relation (*) yields $f(x)f(y)>y-1$, hence
$$f(x)>\frac{y-1}{f(y)} \quad \text{for all} \ x \in \mathbb{R}^{+}$$
If $y>1$, this provides a desired positive lower bound for $f(x)$.
Now, let $M=\inf_{x \in \mathbb{R}^{+}} f(x)>0$. Then, for all $y \in \mathbb{R}^{+}$,
$$M \geqslant \frac{y-1}{f(y)}, \quad \text{or} \quad f(y) \geqslant \frac{y-1}{M}$$
Lemma 1. The function $f(x)$ (and hence $\phi(x)$) is bounded on any segment $[p, q]$, where $0<p<q<+\infty$.
Proof. $f$ is bounded from below by $M$. It remains to show that $f$ is bounded from above on $[p, q]$. Substituting $y=1$ into $(*)$, we get
$$f(x+f(x))=a f(x)$$
Take $z \in [p, q]$ and put $s=f(z)$. By the above, we have
$$f(z+s)=a s \quad \text{and} \quad f(z+s+a s)=f(z+s+f(z+s))=a^2 s .$$
Plugging in $(x, y)=\left(z, 1+\frac{s}{z}\right)$ to $(*)$ and using the previous estimate, we obtain
$$f(z+a s)=f(z+f(z+s))=s f\left(1+\frac{s}{z}\right)-\frac{s}{z} \geqslant \frac{s^2}{M z}-\frac{s}{z}$$
Now, substituting $(x, y)=\left(z+a s, \frac{z}{z+a s}\right)$ to $(*)$ and applying the above estimate and the estimate $f(y) \geqslant M$, we obtain
$$a^2 s=f(z+s+a s)=f(z+a s+f(z))=f(z+a s)f\left(\frac{z}{z+a s}\right)+1-\frac{z}{z+a s} \\ \geqslant M f(z+a s) \geqslant \frac{s^2}{z}-\frac{M s}{z} \geqslant \frac{s^2}{q}-\frac{M s}{p}$$
This yields $s \leqslant q\left(\frac{M}{p}+a^2\right)=: L$, and $f$ is bounded from above by $L$ on $[p, q]$.
Applying Lemma 1 to the segment $[a, a+1]$, we see that $\phi$ is bounded on it. By the previous symmetry, we get that $\phi$ is also bounded on $[a+1,+\infty)$, and hence on $[a,+\infty)$. Put $C=\max\{a, 3\}$.

Lemma 2. For all $x \geqslant C$, we have $\phi(x)=0$ (and hence $f(x)=x+1$).
Proof. Substituting $y=x$ to the earlier equation, we obtain
$$\phi\left(x+f\left(x^2\right)\right)=f(x)^2-f\left(x^2\right)-2x$$
hence,
$$\phi\left(x+f\left(x^2\right)\right)+\phi\left(x^2\right)=f(x)^2-(x+1)^2=\phi(x)(f(x)+x+1)$$
Since $f(x)+x+1 \geqslant C+1 \geqslant 4$, we obtain that
$$|\phi(x)| \leqslant \frac{1}{4}\left(|\phi\left(x+f\left(x^2\right)\right)|+|\phi\left(x^2\right)|\right)$$
Since $C \geqslant a$, there exists a finite supremum $S=\sup_{x \geqslant C}|\phi(x)|$. For each $x \in [C,+\infty)$, both $x+f\left(x^2\right)$ and $x^2$ are greater than $x$; hence they also lie in $[C,+\infty)$. Therefore, taking the supremum of the left-hand side over $x \in [C,+\infty)$, we obtain $S \leqslant S/2$ and hence $S=0$. Thus, $\phi(x)=0$ for all $x \geqslant C$.
It remains to show that $f(y)=y+1$ when $0<y<C$. For each $y$, choose $x>\max\{C, \frac{C}{y}\}$. Then all three numbers $x, x y$, and $x+f(x y)$ are greater than $C$, so (*) reads as
$$(x+x y+1)+1+y=(x+1)f(y)+1, \quad \text{hence} \quad f(y)=y+1$$