Time-dependent fractional partial differential equations typically require huge amounts of memory and computational time, especially for long-time integration, which taxes computational resources heavily for high-dimensional problems. Here, we first analyze existing numerical methods of sum-of-exponentials for approximating the kernel function in constant-order fractional operators, and identify the current pitfalls of such methods. In order to overcome the pitfalls, an improved sum-of-exponentials is developed and verified. We also present several sum-of-exponentials for the approximation of the kernel function in variable-order fractional operators. Subsequently, based on the sum-of-exponentials, we propose a unified framework for fast time-stepping methods for fractional integral and derivative operators of constant and variable orders. We test the fast method based on several benchmark problems, including fractional initial value problems, the time-fractional Allen-Cahn equation in two and three spatial dimensions, and the Schrödinger equation with nonreflecting boundary conditions, demonstrating the efficiency and robustness of the proposed method. The results show that the present fast method significantly reduces the storage and computational cost especially for long-time integration problems.