The stochastic Landau-Lifshitz-Gilbert-Slonczewski (s-LLGS) equation is widely used to study the temporal evolution of the macrospin subject to spin torque and thermal noise. The numerical simulation of the s-LLGS equation requires an appropriate choice of stochastic calculus and numerical integration scheme. In this paper, we comprehensively evaluate the accuracy and complexity of various numerical techniques to solve the s-LLGS equation. We focus on implicit midpoint, Heun, and Euler-Heun methods that converge to the Stratonovich solution of the s-LLGS equation. By performing numerical tests for both strong (path-wise) and weak (statistical) convergence, we quantify the accuracy of various numerical schemes used to solve the s-LLGS equation. We demonstrate a new method intended to solve Stochastic Differential Equations (SDEs) with small noise (RK4-Heun), and test its capability to handle the s-LLGS equation. We also discuss the circuit implementation of nanomagnets for large-scale SPICE-based simulations. We evaluate the efficacy of SPICE in handling the stochastic dynamics of the multiplicative noise in the s-LLGS equation. Numerical schemes such as Euler and Gear, typically used by SPICE-based circuit simulators do not yield the expected outcome when solving the Stratonovich s-LLGS equation. While the trapezoidal method in SPICE does solve for the Stratonovich solution, its accuracy is limited by the minimum time step of integration in SPICE. We implement the s-LLGS equation in both its cartesian and spherical coordinates form in SPICE and compare the stability and accuracy of the two implementations. The results in this paper will serve as guidelines for researchers to understand the tradeoffs between accuracy and complexity of various numerical methods and the choice of appropriate calculus to solve the s-LLGS equation.