Лема. Якщо функція $F:(0;+\infty) \rightarrow \mathbf{R}$ задовольняє на проміжку $(0;+\infty)$ функціональне рівняння Коші
$$F(u + v) = F(u) + F(v)$$
і є обмеженою на інтервалі $(0;1)$ (тобто, існує така стала $M$, що $|F(u)| \le M$ для всіх $u \in (0;1)$), то $F(u) = F(1)u$ для всіх $u \in (0;+\infty)$.

Доведення леми. Розглянемо функцію $G(u) = F(u) - F(1)u$. Легко переконатися, що $G(1) = 0$ і що функція $G$ так само, як і $F$, задовольняє на $(0;+\infty)$ функціональне рівняння Коші: для всіх $u, v \in (0;+\infty)$
$$G(u + v) = F(u + v) - (u + v)F(1) = F(u) + F(v) - uF(1) - vF(1) = G(u) + G(v).$$
Крім того, для всіх $u \in (0; +\infty)$
$$G(u + 1) = F(u + 1) - (u + 1)F(1) = F(u) + F(1) - uF(1) - F(1) = G(u).$$
Функція $G$ є обмеженою на всьому проміжку $(0;+\infty)$. Дійсно, будь-яке дійсне число $u > 0$ можна подати у вигляді $k + \alpha$, де $k$ є цілим невід'ємним числом, $0 < \alpha \le 1$, причому $|G(u)| = |G(k + \alpha)| = |G(\alpha)| \le |F(\alpha)| + \alpha|F(1)| \le M + |F(1)|$ для довільного дійсного $u > 0$. Але за індукцією легко довести, що $G(nu) = nG(u)$ для довільного дійсного $u > 0$ і довільного натурального числа $n$. Отже, нерівність $|G(u)| = \frac{1}{n} |G(nu)| \le \frac{1}{n} (M + |F(1)|)$ виконується для довільного натурального $n$ і довільного дійсного $u > 0$. З цього випливає, що $G(u) = 0$ на $(0; +\infty)$.

Перейдемо до розв'язання задачі.
З умови виводимо, що
$$\begin{align*} f(x+y+z) - f(x) - f(y) - f(z) &= \\ &= f(x+y+z) - f(x+y) - f(z) + f(x+y) - f(x) - f(y) \\ &= f(xz + yz + 1) - f(xz + yz) - f(1) + f(xy + 1) - f(xy) - f(1) \end{align*}$$
для будь-яких $x, y, z > 0$.
Звідки
$$\begin{align*} f(xz + xy + 1) - f(xz + xy) - f(1) + f(yz + 1) - f(yz) - f(1) &= \\ &= f(xz + yz + 1) - f(xz + yz) - f(1) + f(xy + 1) - f(xy) - f(1). \end{align*}$$
Отже, для довільних $a,b,c > 0$ — якщо покласти $x = \sqrt{ca/b}$, $y = \sqrt{bc/a}$, $z = \sqrt{ab/c}$ —
одержуємо рівність
$$f(a+b+1)-f(a+b)-f(b+1)+f(b)=f(a+c+1)-f(a+c)-f(c+1)+f(c).$$
Розглянемо першу різницю функції $f$ як функцію $\Delta f$, що визначається на $[0;+\infty)$ співвідношенням $\Delta f(u) = f(u+1) - f(u)$. Тоді, як нескладно отримати, для всіх $a,b,c > 0$ виконується рівність
$$\Delta f(b+a) - \Delta f(b) = \Delta f(c+a) - \Delta f(c).$$
Права частина цієї тотожності не залежить від $b$, тому й ліва її частина від $a$ не залежить. Отже, $\Delta f(b+a) - \Delta f(b) = \phi(a)$ для певної дійсної функції $\phi$, визначеної на проміжку $(0;+\infty)$. Помінявши місцями $a$ та $b$, отримаємо, що $\Delta f(a+b) - \Delta f(a) = \phi(b)$, звідки випливає, що $\Delta f(a) - \phi(a) = \Delta f(b) - \phi(b)$ для всіх $a, b > 0$. Знову ж таки, права частина останньої рівності не залежить від $a$, отже, ліва також не залежить, тому є сталою: $\Delta f(a) - \phi(a) = C_1$ для якогось дійсного $C_1$. Таким чином, $\Delta f(a+b) = \Delta f(b) + \Delta f(a) - C_1$ для всіх $a, b > 0$. Легко переконатися, що функція $g(u) = \Delta f(u) - C_1$ задовольняє умови леми (обмеженість $g$ на $(0;1]$ випливає з монотонності $f$), тому $g(u) = g(1)u$, і функція $\Delta f$ є лінійною на $(0;+\infty)$: $\Delta f(u) = C_2u + C_1$, де $C_2 = g(1)$. Тепер з вихідної рівності для довільних $x, y > 0$ виводимо, що $f(x+y)-f(x)-f(y)=C_2xy+C_3$, де $C_3 = C_1 - f(1)$. Функція $h(u) = f(u) - \frac{1}{2}C_2u^2 + C_3$, як нескладно довести, задовольняє умови леми. Тому $h(u) = h(1)u$, і на проміжку $(0; +\infty)$ функція $f$ має вигляд $f(u) = Au^2 + Bu + C$, де $A = \frac{1}{2}C_2$, $B = h(1)$, $C = -C_4$. Згідно з умовою задачі маємо також, що
$$f(0) = f(3) - 3f(2) + 3f(1) = (9A + 3B + C) - 3(4A + 2B + C) + 3(A + B + C) = C.$$
Отже, $f(u) = Au^2 + Bu + C$ для всіх $u \ge 0$. Підстановка такого тричлена в вихідну рівність дає тотожність для будь-яких дійсних $A, B$ і $C$. Залишається відкинути випадок $AB < 0$, коли $f$ не є монотонною функцією на $[0; +\infty)$.