Abstract: The precise measurement and error compensation technique of the main magnetic field in cyclotron accelerator are important research directions and key components in the construction of the device. Due to factors such as internal defects in the iron material and mechanical processing errors, the non-ideal magnetic field generated by the main magnet of the cyclotron accelerator usually deviates from the required isochronous magnetic field distribution and contains a certain amplitude of harmonic magnetic field. Therefore, in the process of developing a cyclotron accelerator, it is necessary to perform multiple shimming on the non-ideal magnetic field, and ultimately achieve the required distribution by correcting the actual magnetic field. Compared with traditional magnetic field compensation algorithms, this paper proposes an algorithm based on a multivariate linear regression model, incorporating the calculation of the first harmonic magnetic field. This algorithm achieves simultaneous quantitative shimming for both temporal errors in the magnetic field and first harmonic errors, eliminating the need for a separate iterative process for shimming the first harmonic magnetic field. In order to avoid increasing the finite element calculation workload, this paper uses a 1/4 magnet model in the new algorithm to obtain the shape function for the first harmonic magnetic field and the average magnetic field. Simplifying the full model finite element calculation of the main magnet to a 1/4 model can save approximately 80% of CPU time. By using the algorithm proposed in this paper, after three iterations of shimming, the isochronism error of the magnetic field in the 16 MeV cyclotron at the China Institute of Atomic Energy is reduced to the order of 10^-4, and the integral sliding phase control of particles is controlled within ±14°, while the magnitude of the first harmonic of the magnetic field is reduced to within 6 Gs. The transverse free oscillation frequency of the magnetic field was also adjusted, improving the axial focusing frequency at large radius positions, allowing the beam to pass through dangerous resonances only at low energy positions and quickly move away from resonance positions. Through the first shimming process, the magnitude of the first harmonic of the magnetic field is reduced, minimizing the coupling of the beam phase space in the transverse direction and avoiding major harmful resonances. Experimental studies have shown that this algorithm has the characteristics of low computational cost for shape function calculations, high shimming accuracy, and a small number of iterations. The algorithm in this article can be further expanded to realize the shimming calculation of any high-order harmonic magnetic field error of the cyclotron.