Muntz theorems from classical approximation theory are used to lay the mathematical foundations for discrete spectral representations of the relaxation modulus and creep conpliance. Algorithms that can automatically generate nonnegative discrete relaxation and retardation strengths are formulated by using the nonnegative least square method.