We have evaluated perturbation coefficients of Wilson loops up to O(g8) for the four-dimensional
twisted Eguchi-Kawai model using the numerical stochastic perturbation theory (NSPT) in [1].
In this talk we present a progress report on the higher order calculation up to O(g63),
for which we apply a fast Fourier transformation (FFT) based convolution algorithm to the multiplication
of polynomial matrices in the NSPT aiming for higher order calculation.
We compare two implementations with the CPU-only version and the GPU version of the FFT based convolution algorithm, and
find a factor 9 improvement on the computational speed of the NSPT algorithm with SU(N=225) at O(g31).
The perturbation order dependence of the computational time, we investigate it up to O(g63),
shows a mild scaling behavior on the truncation order.