Computation of singular solutions to the Helmholtz equation with high order accuracy