We develop a numerical algorithm for calculating the light-scattering properties of small particles of arbitrary shape on the basis of a method involving surface integral equations. The calculation error was estimated by performing a comparison between the proposed method and the exact Mie method with regard to the extinction efficiency factor, and the results show that the error is less than 1% when four or more nodes per wavelength are set on the surface of a spherical particle. The accuracy fluctuates in accordance with the distribution of nodal points on the particle surface with respect to the direction of propagation of the incident light. From our examinations, it is shown that the polar incidence alignment yields higher accuracy than equator incidence when a "latitude-longitude" type of mesh generation is adopted. The electric currents on the particle surface and the phase functions of all scattering directions are shown for particles shaped as spheres or hexagonal columns. It is shown that the phase function for a hexagonal column has four or eight cold spots. The phase function of a randomly oriented hexagonal column shows halolike peaks with size parameters of up to 20. This method can be applied to particles with a size parameter of up to about 20 without using the symmetry characteristic of the particle.