The classical calculation method of additional stresses in foundation soils was established based on the Flamant and Boussinesq solutions in elasticity theory. This method can only calculate the additional stresses in foundations with regular shapes and under regularly distributed base pressures. Based on the Gauss-Legendre numerical product formula, which is independent of the specific form of the product function, a suitable method for calculating the additional stresses in irregularly shaped foundation soils under non-uniform loads was established in this study. Using Simpson's formula to perform integration over the domain, Gaussian summation calculation methods for the additional stresses in foundation soil under non-uniform loads in rectangular domains, under loads in non-rectangular domains, and under irregular loads and loads in irregular integration regions were derived. The combination of the complex product algorithm and a computer program further efficiently improved the calculation accuracy. Validation examples showed that the Gaussian product calculation method of additional stresses in foundation soil was completely consistent with classical elasticity theory, and the results of the Gaussian product calculation method converged to the classical elasticity theory prediction when the complex product formula was applied. Based on validation and application examples, the method has advantages over the finite element method in terms of modeling difficulty, computational accuracy, computational volume, and computational cost when calculating additional stresses in foundation soil. The research results provide a new method and idea for the calculation of additional stresses in irregularly shaped flexible foundation soils under non-uniform loads.