In this paper, we develop a quadrature framework for large-scale kernel machines via a numerical multiple integration representation. Leveraging the fact that the integration domain and measure of typical kernels, e.g., Gaussian kernels, arc-cosine kernels, are fully symmetric, we introduce a deterministic fully symmetric interpolatory rule to efficiently compute its quadrature nodes and associated weights to approximate such typical kernels. This interpolatory rule is able to reduce the number of needed nodes while retaining a high approximation accuracy. Further, we randomize the above deterministic rule such that the proposed stochastic version can generate dimension-adaptive feature mappings for kernel approximation. Our stochastic rule has the nice statistical properties of unbiasedness and variance reduction with fast convergence rate. In addition, we elucidate the relationship between our deterministic/stochastic interpolatory rules and current quadrature rules for kernel approximation, including the sparse grids quadrature and stochastic spherical-radial rule, thereby unifying these methods under our framework. Experimental results on several benchmark datasets show that our fully symmetric interpolatory rule compares favorably with other representative random features based methods.