The graph fractional Fourier transform (GFRFT) generalizes the graph Fourier transform (GFT) but suffers from a significant computational bottleneck: determining the optimal transform order requires expensive eigendecomposition and matrix multiplication, leading to $O(N^3)$ complexity. To address this issue, we propose a fast GFRFT (FGFRFT) algorithm for unitary GFT matrices based on Fourier series approximation and an efficient caching strategy. FGFRFT reduces the complexity of generating transform matrices to $O(2LN^2)$ while preserving differentiability, thereby enabling adaptive order learning. We validate the algorithm through theoretical analysis, approximation accuracy tests, and order learning experiments. Furthermore, we demonstrate its practical efficacy for image and point cloud denoising and present the fractional specformer, which integrates the FGFRFT into the specformer architecture. This integration enables the model to overcome the limitations of a fixed GFT basis and learn optimal fractional orders for complex data. Experimental results confirm that the proposed algorithm significantly accelerates computation and achieves superior performance compared with the GFRFT.