Abstract <p>An algorithm is proposed for calculating the expansion coefficients of functions in series of classical orthogonal polynomials. The algorithm is based on simultaneous calculation of basis functions and expansion coefficients. Recurrence relations are used to calculate the orthogonal functions. Gaussian quadrature are used to calculate the expansion coefficients. The algorithm has constant complexity in terms of the memory used and computational complexity depending linearly both on the number of points in the quadrature formula and the number of expansion coefficients. An optimal vectorization scheme for the algorithm is proposed. The algorithm with a fixed vectorization depth has the maximum efficiency. With parallel calculation of the expansion coefficients on systems with shared memory, the algorithm scales linearly. The algorithm is applicable to classical orthogonal Chebyshev, Legendre, Jacobi, Sonin–Laguerre, and Hermite polynomials, as well as to trigonometric Fourier series.</p>