An efficient and accurate 3D algorithm for dynamical simulations of many deformable drops with strong hydrodynamical interactions at zero Reynolds numbers is developed. The drop-to-medium viscosity ratio, Ξ», and the Bond number are arbitrary, and the drops are subject to gravity with stationary triply periodic boundary conditions. The algorithm, at each step, is a hybrid of boundary-integral and economical multipole techniques, with extensive use of rotational transformations and economical truncation of multipole expansions to optimize near-field interactions. A significant part of the code is the new, "best paraboloid-spline" technique for calculating the normal vectors and curvatures on drop surfaces, which greatly improves the quality of long-time simulations. Examples show the phenomenon of clustering in a concentrated sedimenting emulsion for Ξ» = 0.25 and 1, which leads to an increase in the average sedimentation velocity with time. A high efficiency of the method is demonstrated, with two orders-of-magnitude gains over the standard O(N 2 N 2 ) boundary-integral technique for N βΌ 10 2 drops in a periodic cell with N βΌ 10 3 triangular boundary elements per drop, so that typical long-time dynamical simulations can be performed in a few days or weeks on a standard workstation (as compared to the several years which would be required for the same simulations using standard boundary-integral techniques). The effects of drop triangulation and truncation of multipole expansions on dynamical simulations are assessed.