A detailed study of multiparticle azimuthal correlations is presented using pp data at s=5.02 and 13 TeV, and p+Pb data at sNN=5.02 TeV, recorded with the ATLAS detector at the CERN Large Hadron Collider. The azimuthal correlations are probed using four-particle cumulants cn{4} and flow coefficients vn{4}=(−cn{4})1/4 for n=2 and 3, with the goal of extracting long-range multiparticle azimuthal correlation signals and suppressing the short-range correlations. The values of cn{4} are obtained as a function of the average number of charged particles per event, Nch, using the recently proposed two-subevent and three-subevent cumulant methods, and compared with results obtained with the standard cumulant method. The standard method is found to be strongly biased by short-range correlations, which originate mostly from jets with a positive contribution to cn{4}. The three-subevent method, on the other hand, is found to be least sensitive to short-range correlations. The three-subevent method gives a negative c2{4}, and therefore a well-defined v2{4}, nearly independent of Nch, which implies that the long-range multiparticle azimuthal correlations persist to events with low multiplicity. Furthermore, v2{4} is found to be smaller than the v2{2} measured using the two-particle correlation method, as expected for long-range collective behavior. Finally, the measured values of v2{4} and v2{2} are used to estimate the number of sources relevant for the initial eccentricity in the collision geometry. The results based on the subevent cumulant technique provide direct evidence, in small collision systems, for a long-range collectivity involving many particles distributed across a broad rapidity interval.