Average settling velocity of non-uniform hard spheres in a viscous fluid is determined by using a largescale numerical simulation that is carried out for over 103 spheres in a periodic unit cell which extends infinitely. An efficient calculation scheme is used for reducing the computation cost which steeply increases with the number of the spheres. The calculation scheme is based on a fast summation method for far-field hydrodynamic interaction among spheres. It is applied in the computation of hindered settling velocity of hard spheres with bidisperse size distribution in a viscous fluid. The simulation results are compared with the theoretical predictions by Batchelor [8] and Davis and Gecol [9]. It is found that the prediction by Davis and Gecol reasonably agrees with the numerical results.