We calculate the triple-and quadruple-gluon inclusive distributions with arbitrary rapidity and azimuthal angle dependence in the gluon saturation regime by using glasma diagrams. Also, we predict higher dimensional ridges in triple-and quadruple-hadron correlations for p-p and p-Pb collisions at LHC, which have yet to be measured. In p-p and p-Pb collisions at the top LHC energies, gluon saturation is expected to occur since smaller Bjorken-x values are being probed. Glasma diagrams, which are enhanced at small-x, include gluon saturation effects, and they are used for calculating the long-range rapidity correlations ("ridges") and vn moments of the azimuthal distribution of detected hadrons. The glasma description reproduces the systematics of the data on both p-p and p-Pb ridges. As an alternative, relativistic hydrodynamics has also been applied to these small systems quite successfully. With the triple-and quadruple-gluon azimuthal correlations, this work aims to set the stage by going beyond the double-gluon azimuthal correlations in order to settle unambiguously the origin of "collectivity" in p-p and p-Pb collisions. We derive the triple-and quadruple-gluon azimuthal correlation functions in terms of unintegrated gluon distributions at arbitrary rapidities and azimuthal angles of produced gluons. Then, unintegrated gluon distributions from the running coupling Balitsky-Kovchegov evolution equation are used to calculate the triple-and quadruple-gluon correlations for various parameters of gluon momenta, initial scale for small-x evolution and beam energy.