# What is saturated model used in glm in R?

So, I've read on stack exchange (glm summary explained) that residual deviance computed in the glm output is just the likelihood ratio chi-square stat comparing the saturated model to the reduced model. 1) But what saturated model is used? 2) How is the log likelihood evaluated at mle estimates computed for the given saturated model?

For simplicity, let's assume the glm is from the normal family and is just standard linear regression.

1)As I understand it a saturated model has as many parameters as data points-so there should be several (infinite) possible saturated models. For example, suppose you have two explanatory variables with four data points. One possible saturated model could be $$y=b_1x_1+b_2x_2+b_3x_1x_2+b_4$$ or another possible saturated model could be $$y=b_1x_1+b_2x_2+b_3x_2^2+b_4$$.

2)Regardless of the saturated model wouldn't the mle estimates for these coefficients just be found by solving the system directly using these four data points (as R naturally computes them)? Indeed, these would lead you to an infinite likelihood as residuals would be 0.