nterior solver presented in Section \[s:interior\_solver\_cylindrical\] works without modification even when $L_\phi=2\pi/P$ does not cover a full $2\pi$ domain, because the $P$-fold symmetry is automatically taken into account in Equation . This holds true also for the boundary solver (Section \[s:bc\]) as long as the periodic boundary condition with period $L_\phi$ is imposed in solving Equation . The cylindrical DGF defined as such then represents the gravitational potential from $P$ identical point masses lying along the azimuth with a uniform angular separation of $L_\phi$. When applying the boundary condition to the cylindrical DGF, therefore, one should add all the contributions from the image masses located at $\phi\in[L_\phi, 2\pi]$ to the gravitational potential of a real point mass located in the original domain with $\phi\in[0, L_\phi]$, as in Equation . We note that the inclusion of the image masses is implicit in the DGF, so that all calculations are done in the original computational domain with size $L_\phi$, enabling a factor of $P$ reduction in the computational time as well as memory compared to the cases with a full $2\pi$-periodic domain.

To verify that our method handles $P$-fold azimuthal symmetry, we consider four uniform spheres located at $\phi_0$, $\phi_0+\pi/2$, $\phi_0+\pi$, and $\phi_0+3\pi/2$ in a uniform cylindrical grid presented in Table \[tb:sphere\_test\]. We calculate the gravitational potential from the resulting mass distribution that clearly has $P=4$ symmetry over $\phi\in[0, 2\pi]$. Figure \[fig:pfold\] plots as black solid lines the mass density and the gravitational potential along the $\phi$-direction at $R=0.75$ and $z=0.0039$. We then recompute the gravitational potential of a single sphere at $\phi_0$ by reducing the computational domain to $\phi\in[0, \pi/2]$, and finally the gravitational potential of two spheres at $\phi_0+\pi$ and $\phi_0+3\pi/2$ in the domain covering $\phi\in[\pi, 2\pi]$. The resulting potentials, plotted as red and blue dashed lines, are ide