Accurate virial coefficients B_{N}(λ,ɛ) (where ɛ is the well depth) for the three-dimensional square-well and square-step potentials are calculated for orders N=5-9 and well widths λ=1.1-2.0 using a very fast recursive method. The efficiency of the algorithm is enhanced significantly by exploiting permutation symmetry and by storing integrands for reuse during the calculation. For N=9 the storage requirements become sufficiently large that a parallel algorithm is developed. The methodology is general and is applicable to other discrete potentials. The computed coefficients are precise even near the critical temperature, and thus open up possibilities for analysis of criticality of the system, which is currently not accessible by any other means.