  We study boundary regularity for solutions to a class of equations involving
the so called regional fractional Lapacians $(-\Delta)^s_\Omega $, with
$\Omega\subset \mathbb{R}^N$. Recall that the regional fractional Laplacians
are generated by symmetric stable processes which are not allowed to jump
outside $\Omega$. We consider weak solutions to the equation
$(-\Delta)^s_\Omega w(x)=p.v.\int_{\Omega}\frac{w(x)-w(y)}{|x-y|^{N+2s}}\,
dy=f(x)$, for $s\in (0,1)$, subject to zero Neumann or Dirichlet boundary
conditions. The boundary conditions are defined by considering $w$ as well as
the test functions in the fractional Sobolev spaces $H^s(\Omega)$ or
$H^s_0(\Omega)$ respectively.
  While the interior regularity is well understood for these problems, little
is known in the boundary regularity, mainly for the Neumann problem. Under
optimal regularity assumptions on $\Omega$ and provided $f\in L^p(\Omega)$, we
show that $w\in C^{2s-N/p}(\overline \Omega)$ in the case of zero Neumann
boundary conditions. As a consequence for $2s-N/p>1$, $w\in
C^{1,2s-\frac{N}{p}-1}(\overline{\Omega})$. As what concerned the Dirichlet
problem, we obtain ${w}/{\delta^{2s-1}}\in C^{1-N/p}(\overline\Omega)$,
provided $p>N$ and $s\in (1/2,1)$, where
$\delta(x)=\textrm{dist}(x,\partial\Omega)$.
  To prove these results, we first classify all solutions having a certain
growth at infinity when $\Omega$ is a half-space and the right hand side is
zero. We then carry over a fine blow up and some compactness arguments to get
the results.

Analysis of PDEs