This is only a partial result (EDIT: full result is given below), but I'll post it anyway. Let us write $(x,y)$ if $y = f(x)$ holds. Since $f(f(x)) = x$, $(x,y)$ is equivalent to $(y,x)$ and $(x,y)\wedge (x,y')$ implies $y= y'$. What we can directly observe is that $0= f(f(0)) = f(0)$ by letting $x=y=0$. What is less trivial is $f(1) = 1$. Let $\gamma =f(1)$. Then by letting $y = 1$ and $y= \gamma$, we get
$$
(x+ \gamma + f(x),1+ f(x) +\gamma x)\wedge \;(x + 1+\gamma f(x), \gamma + f(x) +x).$$
Hence, $1+ f(x) +\gamma x = x + 1+\gamma f(x)$ for all $x$, and if $\gamma \neq 1$, then $f(x) = x$, leading to contradiction.
We next show that $$(-\frac{f(x) - f(x')}{x-x'}, -\frac{x-x'}{f(x) - f(x')}),\quad \forall x\neq x'\;\cdots(*).$$ Assume $(a,b)\wedge (x,y)\wedge (x',y').$ Then, by the FE,
$$
(a+x+by, b+y +ax)\wedge (a+x'+by', b+y' +ax').
$$ We can equate $b+y +ax$ and $ b+y' +ax'$ by letting $a= -\frac{y-y'}{x-x'}$. Then $b$ should satisfy $a+x+by = a+x'+by'$, that is, $ b= -\frac{x-x'}{y-y'}.$ Because $y = f(x), y' = f(x')$, this proves the claim.
From the previous claim, we know that $(-1,-1)$, i.e. $f(-1) = -1$. By letting $x' = 0, - 1$, we also know that $(x,y)$ implies
$$(-\frac{y}{x}, -\frac{x}{y})\wedge (-\frac{y+ 1}{x + 1}, -\frac{x+ 1}{y + 1}).
$$If we write $-\frac{y-y'}{x-x'}= k$, then $(-\frac{1/k+ 1}{k + 1}, -\frac{k+ 1}{1/k + 1}) = (-\frac{1}{k}, -k)$ also shows that $(\frac{y-y'}{x-x'}, \frac{x-x'}{y-y'})$ whenever $x\neq x'$. Next, by putting $y=-1$ in the original functional equation, we also have $(x,y)$ implies
$$(x-y-1, y-x-1).$$ Iterating this $n$ times we also have
$$(2^n(x-y)-1, 2^n(y-x)-1) \;\cdots (**).
$$
Our next claim is that if $(x,x) \vee (x, x+2^{1-N})$, then $(x+ 2^{-N}, x+2^{-N})$ ($N\geq 0$.) To show this, let us write $(x+ 2^{-N}, \alpha)$ and derive an equation about $\alpha$. If $(x,x)$, $(*)$ implies $( 2^N(x-\alpha ), \frac{1}{2^N(x-\alpha )})$. By applying $(**)$ to $(x+ 2^{-N}, \alpha)$ we also have $(2^N(x-\alpha), 2^N(\alpha-x)-2)$. Hence, $2^N(\alpha-x)-2 = \frac{1}{2^N(x-\alpha )}$ and this implies $2^N(x-\alpha )=-1$ as desired. In case $(x, x+2^{1-N})$ can be dealt with similarly by noting that $(-2^N(x-\alpha )-2, - \frac{1}{2^N(x-\alpha )+2})$.
A very similar argument can also prove that
$$(x,x) \vee (x, x-2^{1-N})\Rightarrow (x-2^{-N}, x-2^{-N}).$$
So far, all the tedious arguments show that $f(x) = x$ holds for every dyadic rationals $x=\frac{j}{2^n}$, that is, the set $F$ of fixed points of $f$ contains every dyadic rational. The following are some other facts about $F$:
(i) If $(x,y)$, then $x+y+xy$ and $x+y+1$ belongs to $F$.
(ii) If $x,y \in F$, then $x+y+xy \in F$.
(iii) If $x \in F$, then $x\pm \frac{j}{2^n}\in F$. (Above claim)
And for some $(x,f(x))$ such that $x\notin F$, this generates so many non-fixed pairs
$$(\pm \left[\frac{f(x) - q}{x-q}\right]^r,\pm \left[\frac{x- q}{f(x)-q}\right]^r), \quad (x+q+qf(x),f(x)+q+qx)_{q \neq 1},$$ for all dyadic rational $q$ and $r\in \mathbf{N}$.
I tried but failed to get further ideas about how $F$ and $\mathbf{R}\setminus F$ looks like. But I wish this helps somehow.
(Note: Actually in the above proposition, $(x, x\pm 2^{1-N})$ cannot occur since their difference cannot be in $F$.)
Following yesterday's post, I've completed proof: $f = id$ is the only solution of the equation.
I'll briefly review some facts already proved.
(i) $F + \mathbf{Q}_{dyad} = F.$
(ii) If $x,y \in F$, then $xy \in F$, or equivalently, $F\cdot F = F$. (Since $x-1, y-1 \in F$ implies $xy -1\in F$ and thus $xy \in F$.)
(iii) If $(x,y)$, then $x+y, \;xy + x+ y \in F$.
Our claim starts with: If $0\neq x \in F$, then $\frac{1}{x} \in F$. Proof is simple. Let $(\frac{1}{x},\gamma).$ Then $\gamma +\frac{1}{x}\in F$. Since $x\in F$, $\gamma x + 1 \in F$, and $\gamma x \in F$. Note that $(\gamma x, \frac{1}{\gamma x}).$(slope of $(\frac{1}{x},\gamma),(0,0).$) This shows $\gamma = \pm \frac{1}{x}.$ If it were that $(\frac{1}{x},-\frac{1}{x})$, then $\frac{1}{x^2} \in F$. This implies that $x\cdot\frac{1}{x^2}=\frac{1}{x} \in F$, as desired.
Our next claim is that $(x,y)$ implies $xy \in F$. We start from $x+y, \;x+y+xy \in F$. If $x+y =0$, it's already done. Otherwise, $\frac{x+y+xy}{x+y} = 1+ \frac{xy}{x+y} \in F$, hence $\frac{xy}{x+y}\in F$. Then this implies $$(x+y)\cdot \frac{xy}{x+y} = xy \in F.$$
Our almost final claim is that $0< y \in F$ implies $\sqrt{y} \in F$. Suppose $(\sqrt[4]{y}, \gamma)$. Then $\gamma^2 \sqrt{y} \in F$. Since $\frac{1}{y} \in F$, we have $\frac{\gamma^2}{\sqrt{y}} \in F.$ Notice that $(\frac{\gamma}{\sqrt[4]{y}}, \frac{\sqrt[4]{y}}{\gamma})$ and hence $(\frac{\gamma^2}{\sqrt{y}}, \frac{\sqrt{y}}{\gamma^2})$ Thus we have $\gamma^4 = y$, $\gamma = \pm \sqrt[4]{y}$. Suppose $\gamma = \sqrt[4]{y}$. Then, $\sqrt[4]{y} \in F$ implies $\sqrt{y} \in F$. Otherwise, $(\sqrt[4]{y}, -\sqrt[4]{y})$ implies $-\sqrt{y} \in F$. Thus $\sqrt{y} \in F.$
Finally we are ready to prove that $(x,y)$ implies $ x-y=d =0$. Assume to the contrary that $d>0$. Then, as I showed in yesterday's post, $(d-1, -d-1).$ But this implies $(-1+d)\cdot(-1-d) = 1-d^2 \in F$, and thus $d^2 \in F.$ So $d$ must be in $F$ and $d-1$ also. This leads to $d-1 = -d-1$, contradicting $d>0$. So, the only solution of the functional equation is $f(x) = x$!