A new iterative algorithm is established for the exact solution of the standard pairing problem, based on the Richardson-Gaudin method using the polynomial approach. It provides efficient and robust solutions for both spherical and deformed systems at a large scale. The key to its success is that the initial guess for the solutions of such a large set of the non-linear equations is provided in a physically meaningful and controllable manner. Moreover, one reduces the large-dimensional problem to a one-dimensional Monte Carlo sampling procedure, which improves the algorithm's efficiency and avoids the non-solutions and numerical instabilities that persist in most existing approaches. We calculated the ground state and low-lying excited states of equally spaced systems at different pairing strengths G. We then applied the model to study the quantum phase transitional Sm isotopes and the actinide nuclei Pu isotopes, where an excellent agreement with experimental data is obtained. Program summary Program Title: IterV1.m CPC Library link to program files: https://doi.org/10.17632/rjnbhgk2p6.1 Licensing provisions: GPLv3 Programming language: Mathematica Nature of problem: The program calculates exact pairing energies based on a new iterative algorithm. The key is the procedure of determining the initial guesses for the large-set non-linear equations involved in a controllable and physically motivated manner. It provides an efficient and robust solver for both spherical and deformed systems in super large model spaces. Solution method: The new iterative algorithm approach starts with simple systems with knucleon pairs and n = klevels, which can be solved iteratively by including one pair and one level at each step using the Newton-Raphson algorithm with a Monte Carlo sampling procedure. Then it takes the solutions of those systems as initial values and obtain the converged results for the full space by gradually adding the remaining levels. In this way, one reduces the k-dimensional Monte Carlo sampling procedure to a one-dimensional sampling, which improves the efficiency of the algorithm and avoids the non-solutions and numerical instabilities.