Robust and efficient numerical algorithms for the discrete prolate spheroidal wave functions