The Monte Carlo method has been applied to numerical modeling of an integrating sphere designed for hemispherical-directional reflectance factor measurements. It is shown that a conventional algorithm of backward ray tracing used for estimation of characteristics of the radiation field at a given point has slow convergence for small source-to-sphere-diameter ratios. A newly developed algorithm that substantially improves the convergence by calculation of direct source-induced irradiation for every point of diffuse reflection of rays traced is described. The method developed is applied to an integrating sphere reflectometer for the visible and infrared spectral ranges. Parametric studies of hemispherical radiance distributions for radiation incident onto the sample center were performed. The deviations of measured sample reflectance from the actual reflectance as a result of various factors were computed. The accuracy of the results, adequacy of the reflectance model, and other important aspects of the algorithm implementation are discussed.