A detailed study of multiparticle azimuthal correlations is presented using pp data at root s = 5.02 and 13 TeV, and p+Pb data at root s(NN) = 5.02 TeV, recorded with the ATLAS detector at the CERN Large Hadron Collider. The azimuthal correlations are probed using four-particle cumulants c(n){4} and flow coefficients v(n){4} = (-c(n){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 c(n){4} are obtained as a function of the average number of charged particles per event, < N-ch >, 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 jetswith a positive contribution to c(n){4}. The threesubevent method, on the other hand, is found to be least sensitive to short-range correlations. The three-subevent method gives a negative c(2){4}, and therefore a well-defined v(2){4}, nearly independent of < N-ch >, which implies that the long-range multiparticle azimuthal correlations persist to events with low multiplicity. Furthermore, v(2){4} is found to be smaller than the v(2){2} measured using the two-particle correlation method, as expected for long-range collective behavior. Finally, the measured values of v(2){4} and v(2){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.