We study boundary regularity for solutions to a class of equations involving the so called regional fractional Lapacians $(-Δ)^s_Ω$, with $Ω\subset \mathbb{R}^N$. Recall that the regional fractional Laplacians are generated by symmetric stable processes which are not allowed to jump outside $Ω$. We consider weak solutions to the equation $(-Δ)^s_Ωw(x)=p.v.\int_Ω\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(Ω)$ or $H^s_0(Ω)$ 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 $Ω$ and provided $f\in L^p(Ω)$, we show that $w\in C^{2s-N/p}(\overline Ω)$ 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Ω)$. As what concerned the Dirichlet problem, we obtain ${w}/{δ^{2s-1}}\in C^{1-N/p}(\overlineΩ)$, provided $p>N$ and $s\in (1/2,1)$, where $δ(x)=\textrm{dist}(x,\partialΩ)$. To prove these results, we first classify all solutions having a certain growth at infinity when $Ω$ 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.

Minor changes has been done. To appear in the "Journal of Differential Equations"