Stochastic geometry (SG) is increasingly utilized in modeling cellular communications systems, as it yields accurate and realistic analytical results. The most common aim of such modeling is to characterize the coverage probability in both the downlink (DL) and uplink (UL) directions. To validate these results, however, researchers need to use simulation, for which efficient algorithms are by and large not available, and that is where the present article comes in. We introduce here two SG based Monte Carlo algorithms to simulate DL and UL cellular communications efficiently and accurately. The algorithms are so comprehensive that they take the parameters of the system as input and deliver its coverage probability as output, considering in the process such details as equipment deployment, signal power, interference, and channel conditions. To demonstrate their superior efficiency, the algorithms have been analyzed rigorously. To test their accuracy, they have been coded and applied to some sample systems, producing results that match those produced analytically to within 4% in DL and 5% in UL. Although they are designed for the field of communications, the algorithms can be easily adapted for other fields where SG is also applicable, such as biology, astronomy and forestry.