Detailed modeling of the evolution of neutral hydrogen in the intergalactic medium during the Epoch of Reionization, 5≤z≤20, is critical in interpreting the cosmological signals from current and upcoming 21-cm experiments such as the Low-Frequency Array (LOFAR) and the Square Kilometre Array (SKA). Numerical radiative transfer codes provide the most physically accurate models of the reionization process. However, they are computationally expensive as they must encompass enormous cosmological volumes while accurately capturing astrophysical processes occurring at small scales (≲Mpc). Here, we present pyC 2 Ray, an updated version of the massively parallel ray-tracing and chemistry code, C 2 -Ray, which has been extensively employed in reionization simulations. The most time-consuming part of the code is calculating the hydrogen column density along the path of the ionizing photons. Here, we present the Accelerated Short-characteristics Octahedral ray-tracing (ASORA) method, a ray-tracing algorithm specifically designed to run on graphical processing units (GPUs). We include a modern Python interface, allowing easy and customized use of the code without compromising computational efficiency. We test pyC 2 Ray on a series of standard ray-tracing tests and a complete cosmological simulation with volume size (349Mpc)3, mesh size of 2503 and approximately 106 sources. Compared to the original code, pyC 2 Ray achieves the same results with negligible fractional differences, ∼10−5, and a speedup factor of two orders of magnitude. Benchmark analysis shows that ASORA takes a few nanoseconds per source per voxel and scales linearly for an increasing number of sources and voxels within the ray-tracing radii.