If the best phase measurements are to be achieved, phase-stepping methods need algorithms that are 112 insensitive to the harmonic content of the sampled waveform and 122 insensitive to phase-shift miscalibration. A method is proposed that permits the derivation of algorithms that satisfy both requirements, up to any arbitrary order. It is based on a one-to-one correspondence between an algorithm and a polynomial. Simple rules are given to permit the generation of the polynomial that corresponds to the algorithm having the prescribed properties. These rules deal with the location and multiplicity of the roots of the polynomial. As a consequence, it can be calculated from the expansion of the products of monomials involving the roots. Novel algorithms are proposed, e.g., a six-sample one to eliminate the effects of the second harmonic and a 10-sample one to eliminate the effects of harmonics up to the fourth order. Finally, the general form of a self-calibrating algorithm that is insensitive to harmonics up to an arbitrary order is given.