We consider rotation-based fast multipole method for the Laplace equation and its application in cases where particle interactions are governed by the Biot—Savart law. The paper presents the necessary formulas for algorithm implementation and addresses less frequently discussed topics, such as the normalization of spherical harmonics and the associated Wigner matrices normalization. The main focus of the paper is devoted to describing the details of the software implementation that significantly accelerate the performance of the code both on CPU and GPU (using CUDA technology). We provide a comprehensive explanation of proposed techniques and include code examples. We have implemented a C++ program using these methods and conducted a comparative analysis with open-source implementations of the fast multipole method, confirming the high efficiency of our approach.