Computing the linear response of a spatially periodic medium often requires integration with respect to a Bloch (quasi-periodicity) wavevector parameter $\mathbf k$, over the dual (reciprocal) torus for the lattice. The integrands are analytic in each coordinate of $\mathbf k$, apart from complex plane singularities of known generic type near or on the real axis. Two applications in which such integrals must be repeatedly performed are: 1) In electronic structure calculations, a Green's function $G(\mathbf k)$ must be integrated over the Brillouin zone in order to compute quantities such as density of states. 2) In time-harmonic wave scattering with a periodic medium but localized source, the Floquet-Bloch method requires integration of Helmholtz (say) solutions over their Bloch wavevector. We present high-order accurate results for 1). The integrand is the trace of the inverse of a small matrix that is analytic in $\mathbf k$. The integral is performed iteratively, dimension by dimension. For the innermost integral, we show increased efficiency using a variant of adaptive integration which explicitly handles nearby poles. For the next integral we investigate the explicit handling of the square-root type singularities, via quadratic Padé. In both cases multiple nearby singularities at unknown locations are handled robustly without contour deformations.