The unsteady shock wave diffraction by a square cylinder in gases of arbitrary particle statistics is simulated using an accurate and direct algorithm for solving the semiclassical Boltzmann equation with relaxation time approximation in phase space. The numerical method is based on the discrete ordinate method for discretizing the velocity space of the distribution function and high-resolution method is used for evolving the solution in physical space and time. The specular reflection surface boundary condition is employed. The complete diffraction patterns including regular reflection, triple Mach reflection, slip lines, vortices and their complex nonlinear manifestations are recorded using various flow property contours. Different ranges of relaxation times corresponding to different flow regimes are considered, and the equilibrium Euler limit solution is also computed for comparison. The effects of gas particles that obey the Maxwell-Boltzmann, Bose-Einstein and Fermi-Dirac statistics are examined and depicted.