We study the regularity up to the boundary of solutions to the Neumann problem for the fractional Laplacian. We prove that if $u$ is a weak solution of $(-Delta)^s u=f$ in $Omega$, $mathcal N_s u=0$ in $Omega^c$, then $u$ is $C^alpha$ up tp the boundary for some $alpha>0$. Moreover, in case $s>frac12$, we then show that $uin C^{2s-1+alpha}(overlineOmega)$. To prove these results we need, among other things, a delicate Moser iteration on the boundary with some logarithmic corrections. Our methods allow us to treat as well the Neumann problem for the regional fractional Laplacian, and we establish the same boundary regularity result. Prior to our results, the interior regularity for these Neumann problems was well understood, but near the boundary even the continuity of solutions was open.