  We report an accessible and robust tool for evaluating the effects of Coulomb
collisions on a test particle in a plasma that obeys Maxwell-J\"uttner
statistics. The implementation is based on the Beliaev-Budker collision
integral which allows both the test particle and the background plasma to be
relativistic. The integration method supports adaptive time stepping, which is
shown to greatly improve the computational efficiency. The Monte Carlo method
is implemented for both the three-dimensional particle momentum space and the
five-dimensional guiding center phase space.
  Detailed description is provided for both the physics and implementation of
the operator. The focus is in adaptive integration of stochastic differential
equations, which is an overlooked aspect among existing Monte Carlo
implementations of Coulomb collision operators. We verify that our operator
converges to known analytical results and demonstrate that careless
implementation of the adaptive time step can lead to severely erroneous
results.
  The operator is provided as a self-contained Fortran 95 module and can be
included into existing orbit-following tools that trace either the full Larmor
motion or the guiding center dynamics. The adaptive time-stepping algorithm is
expected to be useful in situations where the collision frequencies vary
greatly over the course of a simulation. Examples include the slowing-down of
fusion products or other fast ions, and the Dreicer generation of runaway
electrons as well as the generation of fast ions or electrons with ion or
electron cyclotron resonance heating.
