It's dot products there, which are also distributive. So i⋅i + i⋅j + j⋅i + j⋅j = 1 + 0 + 0 + 1 = 2.
We got dot products from the fact that v^2 = v⋅v for any vector v (so in particular, i+j). Then dot products are linear so you can expand that out. So basically the proof compares using FOIL on geometric products to FOIL on dot products.
Note that `ij` is not a vector; it's a tensor (or "multivector" and in particular a "blade" in the GA lingo). The dot product reduction only applies to vectors from the original space. But i, j, and i+j are vectors.
For simplicity and practically this is all over real vector spaces. You can make similar definitions with other fields, but you have to be careful when working mod 2.
For simplicity I'm also using an orthonormal basis. You also need to be a bit more careful if working with e.g. special relativity where your timelike basis vector t has t⋅t=-1 (though you can see things still work).
I didn't show ij = ji = 0. I'm assuming i⋅j = j⋅i = 0 (an orthogonal basis), but i⋅j != ij and j⋅i != ji.
(i+j)^2 = (i+j)⋅(i+j) because v^2=v⋅v for any vector v. When you expand RHS, you get i⋅i + i⋅j + j⋅i + j⋅j = 1 + 0 + 0 + 1 = 2. When you expand the LHS (i+j)^2 as Clifford multiplication, you get 2 + ij + ji. Since RHS and LHS are equal, ij+ji=0.
What am I missing? I think ij = -ji is an independent axiom.