Numerical algorithm for the standard pairing problem based on the Heine-Stieltjes correspondence and the polynomial approach. Nature of problem: The program calculates exact pairing energies based on the Heine-Stieltjes polynomial approach. Existing conventional exact-pairing approaches require solving systems of highly nonlinear equations, which are difficult and often impossible to solve beyond the simplest of the quantum-mechanical many-particle systems. In this study, the Heine-Stieltjes polynomial approach is employed to provide solutions for more than one or two pairs of particles residing in many energy levels. Solution method: The new Heine-Stieltjes polynomial approach transforms the pairing problem to one that involves the handling of only two matrix equations. This, combined with an efficient numerical algorithm implemented by the fast Newton-Raphson method with a Monte Carlo sampling procedure for the initial guesses, makes exact pairing solutions feasible even when more energy levels or heavy nuclei (many pairs) are considered.