In this paper, we present a performance-oriented approach for calculating magnetostatic fields produced by circular coils. Those include the magnetic vector potential, magnetic flux density and magnetic field gradient, together addressed as fields. Many inductive elements used in a wide range of applications can be modeled as circular coils with a rectangular cross-section, and this is reflected in the number of papers on the subject. Instead of extensively utilising analytic methods, we focused on achieving maximum performance by efficiently implementing the Gauss-Legendre quadrature and applying it to expressions which are simple with respect to computational intensity. To properly leverage this, we implemented our approach in the C++ programming language. The CUDA framework enabled us to utilise the full power of Nvidia graphics cards, especially when multiple coils are present. Our approach can also be used from Python and MATLAB with a small performance penalty. Performance of field compute methods is primarily expressed in the number of field values calculated every second. Processor performance was tested on multiple systems and exceeded one million points per second on all of them. Graphics card performance is particularly noteworthy, with over 5 million magnetic flux density values computed every second for a system of 100 coils, placing effective performance at more than 500 million points per second on a contemporary laptop graphics card. The meticulous implementation, available on GitHub, and unique performance are the highlights of our work.