A new extrapolation cascadic multigrid (EXCMG) method is developed to solve large sparse symmetric positive definite systems resulting from the classical 25-point finite-difference discretizations of the three-dimensional (3D) biharmonic equation on rectangular domains. We accomplish this by designing a quartic interpolation-based prolongation operator and using the symmetric successive over-relaxation (SSOR) preconditioned CG method as the multigrid smoother. For the new prolongation operator, quartic interpolations are used for the finite difference solutions on coarse and fine grids twice and once so that two approximations can be obtained on the next finer grid, and then the completed Richardson extrapolation is used for these two approximations to obtain an excellent initial guess on the next finer grid. The proposed EXCMG method with the new prolongation operator is easier to implement than the original EXCMG method. Numerical experiments demonstrate that the new EXCMG is a highly efficient solver for the 3D biharmonic equation and is considerably faster than the original EXCMG method and the aggregation-based algebraic multigrid (AGMG) method developed by Y. Notay. The proposed EXCMG method can solve discrete 3D biharmonic equations with more than 100 million unknowns in dozens of seconds.