  We investigate quantitative properties of nonnegative solutions $u(x)\ge 0$
to the semilinear diffusion equation $\mathcal{L} u= f(u)$, posed in a bounded
domain $\Omega\subset {\mathbb R}^N$ with appropriate homogeneous Dirichlet or
outer boundary conditions. The operator $\mathcal{L}$ may belong to a quite
general class of linear operators that include the standard Laplacian, the two
most common definitions of the fractional Laplacian $(-\Delta)^s$ ($0<s<1$) in
a bounded domain with zero Dirichlet conditions, and a number of other nonlocal
versions. The nonlinearity $f$ is increasing and looks like a power function
$f(u)\sim u^p$, with $p\le 1$.
  The aim of this paper is to show sharp quantitative boundary estimates based
on a new iteration process. We also prove that, in the interior, solutions are
H\"older continuous and even classical (when the operator allows for it). In
addition, we get H\"older continuity up to the boundary.
  Particularly interesting is the behaviour of solution when the number
$\frac{2s}{1-p}$ goes below the exponent $\gamma \in(0,1]$ corresponding to the
H\"older regularity of the first eigenfunction $\mathcal{L}\Phi_1=\lambda_1
\Phi_1$. Indeed a change of boundary regularity happens in the different
regimes $\frac{2s}{1-p} \gtreqqless \gamma$, and in particular a logarithmic
correction appears in the "critical" case $\frac{2s}{1-p} = \gamma$. Indeed a
change of boundary regularity happens in the different regimes $\frac{2s}{1-p}
\gtreqqless \gamma$, and in particular a logarithmic correction appears in the
"critical" case $\frac{2s}{1-p} = \gamma$. For instance, in the case of the
spectral fractional Laplacian, this surprising boundary behaviour appears in
the range $0<s\leq \frac{1-p}{2}$.
