High-order harmonics and phase-shift miscalibration are the two main systematic error sources that affect the accuracy of interferometry. To deal with these problems simultaneously, a multiple-beam phase shifting algorithm based on least-squares iteration is proposed. The proposed algorithm decomposes multiple-beam fringes into fundamental wave and high-order harmonics by using the Fourier series expansion. Only five frames of multiple-bean interference fringes with random phase shifts are required to accurately determine the phase distribution and phase shifts. Simulations show that the proposed algorithm reaches the error less than 0.005 (peak-valley value, PV) and 0.003 rad (root-mean-square, RMS) with 10 iterations when the reflection coefficient and the RMS of random phase shifts are less than 0.6 and 1 respectively. It also shows that the proposed algorithm exhibits higher precision than the traditional five-bucket algorithm. Experiment demonstrates that the proposed algorithm is valid and superior to two-beam phase shifting algorithm.