This paper presents new and reviews existing algorithms for the numerical integration of multivariate functions over d-dimensional cubes using sparse grid methods. The sparse grid method, introduced by Smolyak, constructs multivariate quadrature formulas by combining tensor products of one-dimensional formulas. This approach reduces computational cost, making it almost independent of the problem's dimension if the function has bounded mixed derivatives. The authors suggest using extended Gauss (Patterson) quadrature formulas as the one-dimensional basis, which outperform previous methods like trapezoidal, Clenshaw–Curtis, and Gauss rules in numerical experiments. For path integrals, combining generalized Smolyak quadrature with the Brownian bridge construction improves results. Multivariate integrals are common in fields like statistical mechanics, financial derivatives, and partial differential equations, but traditional methods suffer from the "curse of dimension," where computational cost grows exponentially with dimension. Smolyak's construction, however, can mitigate this for functions with bounded mixed derivatives, making the number of function evaluations and accuracy independent of dimension up to logarithmic factors. It is also known as the sparse grid method and has been applied with various one-dimensional rules. Other methods for multivariate integration include Monte Carlo, Quasi-Monte Carlo, lattice rules, and neural networks. Smolyak's construction has also been applied in areas like Fourier transforms, wavelet analysis, and solving PDEs. The paper reviews Smolyak-based methods and introduces univariate extended Gauss (Patterson) formulas for higher polynomial exactness. It also discusses extensions, modifications, and a numerically stable implementation. The performance of various sparse grid quadrature formulas is compared in applications from computational physics and financial mathematics.This paper presents new and reviews existing algorithms for the numerical integration of multivariate functions over d-dimensional cubes using sparse grid methods. The sparse grid method, introduced by Smolyak, constructs multivariate quadrature formulas by combining tensor products of one-dimensional formulas. This approach reduces computational cost, making it almost independent of the problem's dimension if the function has bounded mixed derivatives. The authors suggest using extended Gauss (Patterson) quadrature formulas as the one-dimensional basis, which outperform previous methods like trapezoidal, Clenshaw–Curtis, and Gauss rules in numerical experiments. For path integrals, combining generalized Smolyak quadrature with the Brownian bridge construction improves results. Multivariate integrals are common in fields like statistical mechanics, financial derivatives, and partial differential equations, but traditional methods suffer from the "curse of dimension," where computational cost grows exponentially with dimension. Smolyak's construction, however, can mitigate this for functions with bounded mixed derivatives, making the number of function evaluations and accuracy independent of dimension up to logarithmic factors. It is also known as the sparse grid method and has been applied with various one-dimensional rules. Other methods for multivariate integration include Monte Carlo, Quasi-Monte Carlo, lattice rules, and neural networks. Smolyak's construction has also been applied in areas like Fourier transforms, wavelet analysis, and solving PDEs. The paper reviews Smolyak-based methods and introduces univariate extended Gauss (Patterson) formulas for higher polynomial exactness. It also discusses extensions, modifications, and a numerically stable implementation. The performance of various sparse grid quadrature formulas is compared in applications from computational physics and financial mathematics.